text:
import pandas as pd
import numpy as np
import random
import csv
import io
import os
import pickle as pkl

DATA = "data/sentences.csv"
DATA_TRAIN = "data/wili-2018/x_train.txt"
DATA_TEST = "data/wili-2018/x_test.txt"
DATA_PATH = 'data'
STOP_CHAR = '\u0003' # U+0003 \x03

def read_data_taboeta(datafile):
    sent_df = pd.read_csv(datafile, sep="\t", encoding="utf-8", names=['num', 'lang', 'sent'])
    return trim_taboeta_sentences(sent_df)

def trim_taboeta_sentences(data_df):
    languages = data_df['lang'].unique()
    char_set = set()
    sum = 0.0
    sent = 0.0

    num_less_than = 0

    for i, lang in enumerate(languages):
        sentence_df = data_df.loc[data_df['lang'] == lang]
        data = sentence_df['sent'].values
        random.shuffle(data)

        for line in data:
            sum += len(line)
            sent += 1

            if len(line) < 30:
                num_less_than += 1

            char_set.update(c for c in line)

        if (i + 1) % 10 == 0:
            print("Completed: {0} out of {1}".format(i + 1, len(languages)))
    print("Num unique chars in Taboeta: {}".format(len(char_set)))

    print("AVG = ", sum / sent)
    print("NUM less = ", num_less_than)
    return char_set

def read_data_wili(datafile):
    data = io.open(datafile, 'r', encoding='utf-8')
    char_set = set()
    all_data = list(data)
    print("LEN OF FILE: ", len(all_data))
    random.shuffle(all_data)

    for i, sentence in enumerate(all_data):
        char_set.update(c for c in sentence.strip())
        if (i + 1) % 20000 == 0:
            print("Completed: {0} out of {1}".format(i + 1, len(all_data)))
    print("Num unique chars in wili: {}".format(len(char_set)))
    return char_set

def make_mapping(char_set):
    char_to_num_map = dict()
    num_to_char_map = dict()
    for i, c in enumerate(char_set):
        num_to_char_map[i] = c
        char_to_num_map[c] = i

    return char_to_num_map, num_to_char_map

def save_map(map, filepath):
    with open(filepath, 'wb') as file:
        pkl.dump(map, file)

def main():
    universal_char_set = set()
    universal_char_set = universal_char_set | read_data_taboeta(DATA)
    universal_char_set = universal_char_set | read_data_wili(DATA_TRAIN)
    universal_char_set = universal_char_set | read_data_wili(DATA_TEST)

    universal_char_set.add(STOP_CHAR)

    print("Num unique chars: {}".format(len(universal_char_set)))

    char_to_num_map, num_to_char_map = make_mapping(universal_char_set)
    save_map(char_to_num_map, os.path.join(DATA_PATH, 'char_to_num.pkl'))
    save_map(num_to_char_map, os.path.join(DATA_PATH, 'num_to_char.pkl'))

if __name__ == '__main__':
    main()


text:
import torch
import random
import numpy as np

def get_device(device_num=0):
    # If there's a GPU available...
    if torch.cuda.is_available():
        # Tell PyTorch to use the GPU.
        device = torch.device(f"cuda:{device_num}")
        print("There are %d GPU(s) available." % torch.cuda.device_count())
        print("We will use the GPU:", torch.cuda.get_device_name(device_num))
    # If not...
    else:
        print("No GPU available, using the CPU instead.")
        device = torch.device("cpu")
    return device

def set_seed(seed_val=888):
    random.seed(seed_val)
    np.random.seed(seed_val)
    torch.manual_seed(seed_val)
    torch.cuda.manual_seed_all(seed_val)

ext: py

max_stars_repo_name: KRR-Oxford/BERTMap

max_issues_repo_name: KRR-Oxford/BERTMap

max_forks_repo_name: KRR-Oxford/BERTMap

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding:utf-8 -*-

# Copyright (c) 2020 Carnegie Mellon University, Wenshan Wang <wenshanw@andrew.cmu.edu>
# For License information please see the LICENSE file in the root directory.
import math
import numpy as np
from evaluator_base import ATEEvaluator, RPEEvaluator, KittiEvaluator, transform_trajs, quats2SEs
from os.path import isdir, isfile
from data_prepare import writeTimeStamp2PoseGt,get_sr,get_match_list,get_match_pose,get_gt_est_file,get_gt_timeLength,get_est_file_timelength

# from trajectory_transform import timestamp_associate

class TartanAirEvaluator:
    def __init__(self, scale = False, round=1):
        self.ate_eval = ATEEvaluator()
        self.rpe_eval = RPEEvaluator()
        self.kitti_eval = KittiEvaluator()

    def evaluate_one_trajectory(self, gt_traj_name, est_traj_name, scale=False):
        """
        scale = True: calculate a global scale
        """
        # load trajectories
        gt_traj = np.loadtxt(gt_traj_name)
        est_traj = np.loadtxt(est_traj_name)

        if gt_traj.shape[0] != est_traj.shape[0]:
            raise Exception("POSEFILE_LENGTH_ILLEGAL")
        if gt_traj.shape[1] != 7 or est_traj.shape[1] != 7:
            raise Exception("POSEFILE_FORMAT_ILLEGAL")

        # transform and scale
        gt_traj_trans, est_traj_trans, s = transform_trajs(gt_traj, est_traj, scale)
        gt_SEs, est_SEs = quats2SEs(gt_traj_trans, est_traj_trans)

        ate_score, gt_ate_aligned, est_ate_aligned = self.ate_eval.evaluate(gt_traj, est_traj, scale)
        rpe_score = self.rpe_eval.evaluate(gt_SEs, est_SEs)#(rot_error_mean, trans_error_mean)
        kitti_score = self.kitti_eval.evaluate(gt_SEs, est_SEs)#np.mean(rot), np.mean(tra)

        return {'ate_score': ate_score,
                'rpe_score': rpe_score,
                'kitti_score': kitti_score}

if __name__ == "__main__":

    # scale = True for monocular track, scale = False for stereo track
    aicrowd_evaluator = TartanAirEvaluator()
    groundtruthFile,estimateFileList= get_gt_est_file()#estimateFileList是一个数组，有5个文件
    #groundtruthFile='pose_left.txt'
    #输入的是tartanair(有或)没有时间戳的pose_left.txt文件，和tum格式(time tx ty tz qx qy qz qw)的estimate_traj文件
    writeTimeStamp2PoseGt(groundtruthFile)#为groundtruth写时间戳
    #SR=get_sr(groundtruthFile,estimateFile)#利用两个的时间戳计算sr
    gt_time_length= get_gt_timeLength(groundtruthFile)
    result_dict_list=[]
    for estimateFile in estimateFileList:
        est_time_length=get_gt_timeLength(estimateFile)
        if abs(est_time_length-0.0)<0.01:#完全没有跟上的
            continue
        gt_pose_file,est_pose_file=get_match_pose(groundtruthFile,estimateFile)#将gt和est文件生成pose_gt.txt,pose_est.txt，并且里面的位姿一一对应
        result = aicrowd_evaluator.evaluate_one_trajectory(gt_pose_file, est_pose_file, scale=True)
        result['sr']=est_time_length/gt_time_length
        with open('result.txt',"a+") as f:
            f.write(estimateFile)
            f.write('\n')
            for key in result:
                f.write(str(key)+": "+str(result[key])+'\n')#如果一个score有两个数据，则前者表示旋转平均误差，后者表示平移平均误差
            f.write('\n\n\n')
        #print(result)
        result_dict_list.append(result)
    kitti_score_rot=0.0
    kitti_score_trans=0.0
    rpe_score_rot=0.0
    rpe_score_trans=0.0
    ate_score=0.0
    sr=0.0
    nan_time=0#无效数据的个数
    with open('result.txt',"a+") as f:#将所有的值都计算一个平均值
        for result in result_dict_list:
            for key in result:
                if str(key)=="kitti_score":
                    if math.isnan(float(list(result[key])[0])):
                        nan_time+=1
                        print("HAVE AN NAN DATA")
                        break
                    if math.isnan(float(list(result[key])[1])):
                        nan_time+=1
                        print("HAVE AN NAN DATA")
                        break
                if str(key)=="rpe_score":
                    if math.isnan(float(list(result[key])[0])):
                        nan_time+=1
                        print("HAVE AN NAN DATA")
                        break
                    if math.isnan(float(list(result[key])[1])):
                        nan_time+=1
                        print("HAVE AN NAN DATA")
                        break

                if str(key)=="ate_score":
                    if math.isnan(float(result[key])):
                        nan_time+=1
                        print("HAVE AN NAN DATA")
                        break
                if str(key)=="sr":
                    if math.isnan(float(result[key])):
                        nan_time+=1
                        print("HAVE AN NAN DATA")
                        break

                if str(key)=="kitti_score":
                    kitti_score_rot=kitti_score_rot+float(list(result[key])[0])
                    kitti_score_trans=kitti_score_trans+float(list(result[key])[1])
                if str(key)=="rpe_score":
                    rpe_score_rot=rpe_score_rot+float(list(result[key])[0])
                    rpe_score_trans=rpe_score_trans+float(list(result[key])[1])
                if str(key)=="ate_score":
                    ate_score=ate_score+float(result[key])
                if str(key)=="sr":
                    sr=sr+float(result[key])

        kitti_score_rot=kitti_score_rot/(len(result_dict_list)-nan_time)
        kitti_score_trans=kitti_score_trans/(len(result_dict_list)-nan_time)
        rpe_score_rot=rpe_score_rot/(len(result_dict_list)-nan_time)
        rpe_score_trans=rpe_score_trans/(len(result_dict_list)-nan_time)
        ate_score=ate_score/(len(result_dict_list)-nan_time)
        sr=sr/(len(result_dict_list)-nan_time)

        av_kitti_score=[kitti_score_rot,kitti_score_trans]
        av_rpe_score=[rpe_score_rot,rpe_score_trans]
        av_ate_score=ate_score

        new_result={}
        new_result["average_kitti_score"]=av_kitti_score
        new_result["average_rpe_score"]=av_rpe_score
        new_result["average_ate_score"]=av_ate_score
        new_result["average_sr"]=sr

        with open('result.txt',"a+") as f:
            f.write('\n\n')
            f.write("average score:\n")
            for key in new_result:
                f.write(str(key)+": "+str(new_result[key])+'\n')#如果一个score有两个数据，则前者表示旋转平均误差，后者表示平移平均误差
        print(new_result)


text:
import numpy as np
from queue import PriorityQueue
with open("day15_test.txt") as f:
    caves = np.array([list(map(int, line.rstrip())) for line in f.readlines()])

def traverse_caves(caves):
    xdim, ydim = caves.shape
    exit = (xdim-1, ydim-1)
    points = PriorityQueue()
    points.put(((0, 0), (0,0)))
    visited = {(0,0): 0}
    while not points.empty():
        (c, _), p = points.get()
        if p == exit: return c
        x, y = p
        neighs = [
            (max(0, x-1), y),
            (min(xdim-1, x+1), y),
            (x, max(0, y-1)),
            (x, min(ydim-1, y+1))]

        for n in neighs:
            n_cost = c+caves[n]
            if n in visited and n_cost >= visited[n]:
                continue
            visited[n] = n_cost
            dist = (exit[0]-n[0])**2+(exit[1]-n[1])**2
            points.put(((n_cost, dist), n))
print(traverse_caves(caves))

max_stars_repo_name: boostjanbjorge/adventofcode

max_issues_repo_name: boostjanbjorge/adventofcode

max_forks_repo_name: boostjanbjorge/adventofcode


text:
#!/usr/bin/env python

import pickle
import yaml

import numpy as np
import rospy
import sklearn
from sklearn.preprocessing import LabelEncoder
import tf

from geometry_msgs.msg import Pose
from std_msgs.msg import Float64
from std_msgs.msg import Int32
from std_msgs.msg import String
from visualization_msgs.msg import Marker

from pr2_robot.srv import *
from rospy_message_converter import message_converter
from sensor_stick.srv import GetNormals
from sensor_stick.features import compute_color_histograms
from sensor_stick.features import compute_normal_histograms
from sensor_stick.marker_tools import *
from sensor_stick.msg import DetectedObjectsArray
from sensor_stick.msg import DetectedObject
from sensor_stick.pcl_helper import *

def get_normals(cloud):
    get_normals_prox = rospy.ServiceProxy('/feature_extractor/get_normals',
                                          GetNormals)
    return get_normals_prox(cloud).cluster

def make_yaml_dict(test_scene_num, arm_name, object_name, pick_pose, place_pose):
    yaml_dict = {}
    yaml_dict["test_scene_num"] = test_scene_num.data
    yaml_dict["arm_name"]  = arm_name.data
    yaml_dict["object_name"] = object_name.data
    yaml_dict["pick_pose"] = \
            message_converter.convert_ros_message_to_dictionary(pick_pose)
    yaml_dict["place_pose"] = \
            message_converter.convert_ros_message_to_dictionary(place_pose)
    return yaml_dict

def send_to_yaml(yaml_filename, dict_list):
    data_dict = {"object_list": dict_list}
    with open(yaml_filename, 'w') as outfile:
        yaml.dump(data_dict, outfile, default_flow_style=False)

def pcl_callback(pcl_msg):
    # Convert ROS msg to PCL data (XYZRGB)
    cloud = ros_to_pcl(pcl_msg)
    pcl.save(cloud, './misc/pipeline_0_raw_cloud.pcd')

    # Statistical outlier filtering
    outlier_filter = cloud.make_statistical_outlier_filter()
    # Set the number of neighboring points to analyze for any given point
    outlier_filter.set_mean_k(20)
    # Any point with a mean distance larger than global will be considered out
    outlier_filter.set_std_dev_mul_thresh(0.1)
    cloud_filtered = outlier_filter.filter()
    pcl.save(cloud_filtered, './misc/pipeline_1_outlier_removal_filter.pcd')

    # Voxel Grid Downsampling
    vox = cloud_filtered.make_voxel_grid_filter()
    LEAF_SIZE = 0.01
    vox.set_leaf_size(LEAF_SIZE, LEAF_SIZE, LEAF_SIZE)
    cloud_filtered = vox.filter()
    pcl.save(cloud_filtered, './misc/pipeline_2_voxel_grid_filter.pcd')

    # PassThrough Filter to remove the areas on the side of the table
    passthrough_y = cloud_filtered.make_passthrough_filter()
    passthrough_y.set_filter_field_name('y')
    y_axis_min = -0.4
    y_axis_max = 0.4
    passthrough_y.set_filter_limits(y_axis_min, y_axis_max)
    cloud_filtered = passthrough_y.filter()

    # PassThrough Filter to isolate only the objects on the table surface
    passthrough_z = cloud_filtered.make_passthrough_filter()
    passthrough_z.set_filter_field_name('z')
    z_axis_min = 0.61
    z_axis_max = 0.9
    passthrough_z.set_filter_limits(z_axis_min, z_axis_max)
    cloud_filtered = passthrough_z.filter()
    pcl.save(cloud_filtered, './misc/pipeline_3_passthrough_filter.pcd')

    # RANSAC Plane Segmentation to identify the table from the objects
    seg = cloud_filtered.make_segmenter()
    seg.set_model_type(pcl.SACMODEL_PLANE)
    seg.set_method_type(pcl.SAC_RANSAC)
    max_distance = 0.01
    seg.set_distance_threshold(max_distance)
    inliers, coefficients = seg.segment()

    # Extract inliers (table surface) and outliers (objects)
    cloud_table = cloud_filtered.extract(inliers, negative=False)
    cloud_objects = cloud_filtered.extract(inliers, negative=True)
    pcl.save(cloud_table, './misc/pipeline_4_extracted_inliers.pcd')
    pcl.save(cloud_objects, './misc/pipeline_5_extracted_outliers.pcd')

    # Convert the xyzrgb cloud to only xyz since k-d tree only needs xyz
    white_cloud = XYZRGB_to_XYZ(cloud_objects)
    # k-d tree decreases the computation of searching for neighboring points.
    #   PCL's euclidian clustering algo only supports k-d trees
    tree = white_cloud.make_kdtree()

    # Euclidean clustering used to build individual point clouds for each
    #   object. The Cluster-Mask point cloud allows each cluster to be
    #   visualized separately.
    ec = white_cloud.make_EuclideanClusterExtraction()
    ec.set_ClusterTolerance(0.05)
    ec.set_MinClusterSize(100)
    ec.set_MaxClusterSize(3000)
    ec.set_SearchMethod(tree)       # Search the k-d tree for clusters
    # Extract indices for each found cluster. This is a list of indices for
    #   each cluster (list of lists)
    cluster_indices = ec.Extract()

    # Assign a random color to each isolated object in the scene
    cluster_color = get_color_list(len(cluster_indices))

    # Store the detected objects and labels in these lists
    detected_objects_labels = []
    detected_objects = []
    color_cluster_point_list = []

    # Iterate through each detected object cluster for object recognition
    for index, pts_list in enumerate(cluster_indices):

        # Store the object's cloud in this list
        object_cluster = []

        # Create an individual cluster just for the object being processed
        for i, pts in enumerate(pts_list):
            # Retrieve cloud values for the x, y, z, rgb object
            object_cluster.append([cloud_objects[pts][0],
                                   cloud_objects[pts][1],
                                   cloud_objects[pts][2],
                                   cloud_objects[pts][3]])

            # Retrieve cloud values for the x, y, z object, assigning a
            #   preidentified color to all cloud values
            color_cluster_point_list.append([white_cloud[pts][0],
                                             white_cloud[pts][1],
                                             white_cloud[pts][2],
                                             rgb_to_float(cluster_color[index])])

        # Convert list of point cloud features (x,y,z,rgb) into a point cloud
        pcl_cluster = pcl.PointCloud_PointXYZRGB()
        pcl_cluster.from_list(object_cluster)

        # Convert the cluster from pcl to ROS using helper function
        ros_cloud = pcl_to_ros(pcl_cluster)
        #pcl_objects_pub.publish(ros_cloud)

        # Extract histogram features (similar to capture_features.py)
        histogram_bins = 128
        chists = compute_color_histograms(ros_cloud,
                                          nbins=histogram_bins,
                                          using_hsv=True)
        normals = get_normals(ros_cloud)
        nhists = compute_normal_histograms(normals,
                                           nbins=histogram_bins)
        feature = np.concatenate((chists, nhists))

        # Make the prediction, retrieve the label for the result and add it
        #   to detected_objects_labels list
        prediction = clf.predict(scaler.transform(feature.reshape(1, -1)))
        label = encoder.inverse_transform(prediction)[0]
        detected_objects_labels.append(label)

        # Publish a label into RViz
        label_pos = list(white_cloud[pts_list[0]])
        label_pos[2] += .4
        object_markers_pub.publish(make_label(label, label_pos, index))

        # Add the detected object to the list of detected objects.
        do = DetectedObject()
        do.label = label
        do.cloud = ros_cloud
        detected_objects.append(do)

    rospy.loginfo('Detected {} objects: {}'.format(len(detected_objects_labels), detected_objects_labels))

    # Create new cloud containing all clusters, each with a unique color
    cluster_cloud = pcl.PointCloud_PointXYZRGB()
    cluster_cloud.from_list(color_cluster_point_list)

    # Convert PCL data to ROS messages
    ros_cloud_object_cluster = pcl_to_ros(cluster_cloud)
    ros_cloud_objects = pcl_to_ros(cloud_objects)
    ros_cloud_table = pcl_to_ros(cloud_table)

    # Publish ROS messages of the point clouds and detected objects
    pcl_objects_cloud_pub.publish(ros_cloud_object_cluster) # solid color objects
    pcl_objects_pub.publish(ros_cloud_objects)      # original color objects
    pcl_table_pub.publish(ros_cloud_table)          # table cloud
    detected_objects_pub.publish(detected_objects)  # detected object labels

    try:
        # Add some logic to determine whether or not your object detections
        #   are robust enough before calling pr2_mover()
        pr2_mover(detected_objects)
    except rospy.ROSInterruptException:
        pass

def pr2_mover(object_list):
    # Initialize output list that'll store multiple object yaml dictionaries
    output_list = []

    # Load the parameters from the YAML files located in /pr2_robot/config/
    object_list_param = rospy.get_param('/object_list')
    dropbox_list_param = rospy.get_param('/dropbox')

    # TODO: Get this PR2 rotation working
    ## Rotate PR2 in place to capture side tables for the collision map. The
    ##   world_joint_controller controls the revolute joing world_joint between
    ##   the robot's base_footprint and world coordinate frames
    #pub_j1 = rospy.Publisher('/pr2/world_joint_controller/command',
    #                         Float64, queue_size=10)
    #pub_j1.publish(np.pi/2)     # Rotate the PR2 counter-clockwise 90 deg
    #pub_j1.publish(-np.pi)      # Rotate the PR2 clockwise 180 deg
    #pub_j1.publish(np.pi/2)     # Rotate the PR2 back to neutral position

    # Iterate through all objects that should be moved
    for object_params in object_list_param:
        object_name = object_params['name']
        object_group = object_params['group']

        # Check if the object to be moved was found in the perception analysis,
        #   populating the pick_pose message if it was
        for object_i, object_val in enumerate(object_list):
            if object_name != object_val.label:
                # Skip this object b/c it doesn't match the object to be moved
                continue

            # Assign the scene number
            ros_scene_num = Int32()
            # TODO: Figure out what parameter holds the scene data
            #ros_scene_num.data = rospy.get_param('/world_name')
            test_num = 3
            ros_scene_num.data = test_num

            # Assign the object name
            ros_object_name = String()
            ros_object_name.data = object_name

            # Assign the arm that'll be used to pickup the object
            ros_arm_to_use = String()
            if object_group == 'green':
                # The green bin is on the robot's right
                ros_arm_to_use.data = 'right'
            else:
                # The red bin is on the robot's left
                ros_arm_to_use.data = 'left'

            # Get the PointCloud for the object and obtain it's centroid
            #   (the average position of all points in the object cloud).
            #   Convert the cloud to an array, then calculate the average
            #   of the array.
            points_array = ros_to_pcl(object_val.cloud).to_array()
            centroid_numpy = np.mean(points_array, axis=0)[:3]
            # Convert the numpy float64 to native python floats
            centroid = [np.asscalar(x) for x in centroid_numpy]

            # Assign the object pick pose
            ros_pick_pose = Pose()
            ros_pick_pose.position.x = centroid[0]
            ros_pick_pose.position.y = centroid[1]
            ros_pick_pose.position.z = centroid[2]

            # Find the correct dropbox's position
            box_pos = [0, 0, 0]     # Set a default box position
            for box_params in dropbox_list_param:
                if box_params['group'] == object_group:
                    box_pos = box_params['position']
                    break
            # TODO: Add a random offset to the dropbox's position

            # Assign the dropbox pose
            ros_place_pos = Pose()
            ros_place_pos.position.x = box_pos[0]
            ros_place_pos.position.y = box_pos[1]
            ros_place_pos.position.z = box_pos[2]

            # Add the object's yaml dict to the output_list
            obj_yaml_dict = make_yaml_dict(ros_scene_num, ros_arm_to_use,
                                           ros_object_name, ros_pick_pose,
                                           ros_place_pos)
            output_list.append(obj_yaml_dict)
            print('processed %s' % ros_object_name.data)

            ## Wait for 'pick_place_routine' service to come up
            #rospy.wait_for_service('pick_place_routine')
            ## Uses /srv/PickPlace.srv file
            #pick_place_routine = rospy.ServiceProxy('pick_place_routine',
            #                                        PickPlace)
            #try:
            #    # Insert the message variables to be sent as a service request
            #    resp = pick_place_routine(ros_scene_num, ros_object_name,
            #                              ros_arm_to_use, ros_pick_pose,
            #                              ros_place_pos)

            #    print("Response: ", resp.success)
            #except rospy.ServiceException, e:
            #    print("Service call failed: %s" % e)

            # Remove the object from object_list to indicate it was picked up
            del object_list[object_i]

            # Stop looking through the other identified objects
            break

    # Output your request parameters into an output yaml file
    send_to_yaml('output_%i.yaml' % test_num, output_list)

if __name__ == '__main__':

    # ROS node initialization
    rospy.init_node('recognition', anonymous=True)

    # Create Subscriber to receive the published data coming from the
    #   pcl_callback() function that will be processing the point clouds
    pcl_sub = rospy.Subscriber('/pr2/world/points', pc2.PointCloud2,
                               pcl_callback, queue_size=1)

    # Create Publishers
    object_markers_pub = rospy.Publisher('/object_markers', Marker,
                                         queue_size=1)
    detected_objects_pub = rospy.Publisher('/detected_objects',
                                           DetectedObjectsArray,
                                           queue_size=1)
    # Isolated object point cloud with the object's original colors
    pcl_objects_pub = rospy.Publisher('/pcl_objects', PointCloud2, queue_size=1)
    # Isolated object point cloud with random colors
    pcl_objects_cloud_pub = rospy.Publisher('/pcl_objects_cloud', PointCloud2,
                                            queue_size=1)
    # Table point cloud without the objects
    pcl_table_pub = rospy.Publisher('/pcl_table', PointCloud2, queue_size=1)

    # Load model from disk
    model = pickle.load(open('model.sav', 'rb'))
    clf = model['classifier']
    encoder = LabelEncoder()
    encoder.classes_ = model['classes']
    scaler = model['scaler']

    # Initialize color_list
    get_color_list.color_list = []

    # Spin while node is not shutdown
    while not rospy.is_shutdown():
        rospy.spin()

max_stars_repo_name: daemonslayer/robond

max_issues_repo_name: daemonslayer/robond

max_forks_repo_name: daemonslayer/robond


text:
import pandas as pd
import numpy as np

# concatenating
df1 = pd.DataFrame(np.ones((3, 4)) * 0, columns=['a', 'b', 'c', 'd'])
df2 = pd.DataFrame(np.ones((3, 4)) * 1, columns=['a', 'b', 'c', 'd'])
df3 = pd.DataFrame(np.ones((3, 4)) * 2, columns=['a', 'b', 'c', 'd'])
# print(df1, '\n', df2, '\n', df3)

# ignore_index=True 忽略 index 的影响
res = pd.concat([df1, df2, df3], axis=0, ignore_index=True)
print(res)

##############

# join, ['inner', 'outer']
df1 = pd.DataFrame(np.ones((3, 4)) * 0, columns=['a', 'b', 'c', 'd'], index=[1, 2, 3])
df2 = pd.DataFrame(np.ones((3, 4)) * 1, columns=['b', 'c', 'd', 'e'], index=[2, 3, 4])
# print(df1, '\n', df2)

# print(pd.concat([df1, df2])) # 默认的 join 为 outer，没有的数据用 NAN 填充
print(pd.concat([df1, df2]))  # 用 inner 会把没有数据的列裁剪掉，只保留共有列

##############

# join, ['inner', 'outer']
df1 = pd.DataFrame(np.ones((3, 4)) * 0, columns=['a', 'b', 'c', 'd'], index=[1, 2, 3])
df2 = pd.DataFrame(np.ones((3, 4)) * 1, columns=['b', 'c', 'd', 'e'], index=[2, 3, 4])
# print(df1, '\n', df2)
# 左右合并，即行的方向
# res = pd.concat([df1, df2], axis=1)
# 指定 join_axes 指定按照哪个 df 的下标为准
res = pd.concat([df1, df2], axis=1, join_axes=[df1.index])
print(res)

##############

# append

df1 = pd.DataFrame(np.ones((3, 4)) * 0, columns=['a', 'b', 'c', 'd'])
df2 = pd.DataFrame(np.ones((3, 4)) * 1, columns=['a', 'b', 'c', 'd'])
df3 = pd.DataFrame(np.ones((3, 4)) * 1, columns=['b', 'c', 'd', 'e'], index=[2, 3, 4])
# res = df1.append(df2, ignore_index=True)
# res = df1.append([df2, df3])

s1 = pd.Series([1, 2, 3, 4], index=['a', 'b', 'c', 'd'])
# 追加一行数据
res = df1.append(s1, ignore_index=True)
print(res)

max_stars_repo_name: InnoFang/misc-code

max_issues_repo_name: InnoFang/playground

max_forks_repo_name: InnoFang/playground


text:
import sklearn.decomposition
import numpy as np

class PCA:
    def __init__(self):
        self.m = None
        self.V = None

    def _number_of_components(self, contribution_rate, eigenvalues):
        total = np.sum(eigenvalues)
        increment = 0.0
        for m, value in enumerate(eigenvalues):
            increment += value
            if increment / total > contribution_rate:
                # Starts at 0
                return m + 1

    def replace_principal_individuals(self, X, lb, ub):
        self._custom_pca(X, return_pca=False)
        V = self.V
        for i in range(self.m):
            X[i] = np.clip(np.dot(V[i], X[i]), lb, ub)
        return X

    def _custom_pca(self, matrix, return_pca=True):
        if len(matrix.shape) >= 3:
            Exception("Dimensions >=3 not supported in custom_pca.")
        # calculate the mean of each column
        M = np.mean(matrix.T, axis=1)
        # center columns by subtracting column means
        C = matrix - M
        # calculate covariance matrix of centered matrix
        self.V = np.cov(C.T)
        # eigendecomposition of covariance matrix
        eigenvalues, eigenvectors = np.linalg.eig(self.V)
        self.m = self._number_of_components(0.85, eigenvalues)
        if return_pca:
            # project data
            P = eigenvectors.T.dot(C.T)
            return P.T[:self.m]
        return None

    def _industry_pca(self, matrix):
        if len(matrix.shape > 3):
            Exception("Input dimensions >=3 not supported in industry_pca.")
        if len(matrix.shape) == 3:
            samples, nx, ny = matrix.shape
            matrix = matrix.reshape(samples, nx*ny)
        pca = sklearn.decomposition.PCA()
        pca = pca.fit(matrix)
        return pca.transform(matrix)


text:
import pandas as pd
import numpy as np
import glob

def score_calibration(y_pred):
    thr = 0.932107 #Calculated by validation set for minimize ACER
    a = 0.5/thr
    b = 1
    y_pred = np.array([a*x if x<thr else b*x for x in y_pred])
    return y_pred

def main():
	test_df = pd.read_csv('data/lists/phase2_list.txt')
	exp_groups = [['exp1_2stage_seed1'], ['exp2_seed1', 'exp2_seed2'], ['exp3b_seed1', 'exp3b_seed2'], ['exp3c_seed1', 'exp3c_seed2']]
	tta_columns = ['crop0', 'crop0_hflip']
	test_res = np.empty((0,len(test_df)))
	for exp_seeds in exp_groups:
		seed_res = np.empty((0,len(test_df)))
		for exp_name in exp_seeds:
			exp_res = np.empty((0,len(test_df)))
			files = glob.glob(f'data/opts/{exp_name}/fold*/tta_epoch31_test_score.txt')
			for fname in files:
				curr_df = pd.read_csv(fname)
				curr_res = curr_df[tta_columns].mean(1).values
				exp_res = np.vstack((exp_res, curr_res))
			exp_res = exp_res.mean(0)
			seed_res = np.vstack((seed_res, exp_res))

		seed_res = seed_res.mean(0)
		test_res = np.vstack((test_res, seed_res))

	test_res = test_res.mean(0)
	test_res = score_calibration(test_res)
	test_df['score'] = np.around(test_res, decimals=6)
	test_df.to_csv('data/predict_sample.txt', header=False, index=False, sep=' ')

if __name__ == '__main__':
	main()


text:
# coding: utf-8
import numpy as np
from kerasy.Bio.string import ALPHABETS
from kerasy.Bio.suffix import SAIS, BWT_create, reverseBWT

len_string = 100000

def get_test_data():
    rnd = np.random.RandomState(1)
    string = "".join(rnd.choice(ALPHABETS, len_string))
    return string

def test_SAIS():
    input = "abaabababbabbb"
    answer = np.asarray([14,  2,  0,  3,  5,  7, 10, 13,  1,  4,  6,  9, 12,  8, 11])
    SA = SAIS(input)
    assert np.all(SA == answer)

def test_burrows_wheeler_transform():
    string = get_test_data()
    SA = SAIS(string)
    bwt = BWT_create(string, SA)
    reversed_bwt = reverseBWT(bwt)
    assert reversed_bwt == string

max_stars_repo_name: iwasakishuto/Keras-Imitation

max_issues_repo_name: iwasakishuto/Keras-Imitation

max_forks_repo_name: iwasakishuto/Keras-Imitation


text:
import numpy

class X:
    counter = 0

    def __repr__(self):
        string = str(X.counter)
        X.counter = X.counter + 1
        return string

    # __str__() isn't invoked by numpy
    __str__ = __repr__

shape = numpy.asarray([50, 50])
# print X()

flat = numpy.zeros(shape.prod(), dtype=numpy.object)
for i in xrange(shape.prod()):
    flat[i] = X()
shaped = flat.reshape(shape)

print shaped

max_stars_repo_path: studies/2015-08-28.1502 Printing/e01 Copy of b02.py

max_stars_repo_name: friedrichromstedt/upy

max_issues_repo_name: friedrichromstedt/upy

max_forks_repo_name: friedrichromstedt/upy


text:
# Copyright The PyTorch Lightning team.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import collections
import inspect
import os
import time
from typing import Optional

import numpy
import pytorch_lightning as pl
import torch
from pytorch_lightning import Callback
from pytorch_lightning.utilities import rank_zero_info
from pytorch_lightning.utilities.exceptions import MisconfigurationException
from torch import Tensor

from lightning_transformers.utilities.imports import _BOLTS_AVAILABLE

if _BOLTS_AVAILABLE:
    from pl_bolts.callbacks import SparseMLCallback

    class TransformerSparseMLCallback(SparseMLCallback):
        def __init__(self, output_dir, recipe_path):
            self.output_dir = output_dir
            super().__init__(recipe_path=recipe_path)

        def on_fit_start(self, trainer: "pl.Trainer", pl_module: "pl.LightningModule") -> None:
            optimizer = trainer.optimizers

            if len(optimizer) > 1:
                raise MisconfigurationException("SparseML only supports training with one optimizer.")
            optimizer = optimizer[0]

            loggers = trainer.logger

            if not isinstance(loggers, list):
                loggers = [loggers]

            self.manager.initialize(pl_module, epoch=0.0, logger=loggers)
            self.manager.initialize_loggers(loggers)

            optimizer = self.manager.modify(
                pl_module, optimizer, steps_per_epoch=self._num_training_steps_per_epoch(trainer), epoch=0
            )

            trainer.optimizers = [optimizer]

        @staticmethod
        def export_to_sparse_onnx(
            model: "pl.LightningModule", output_dir: str, sample_batch: Optional[Tensor] = None, **kwargs
        ) -> None:
            """Exports the model to ONNX format."""
            import onnxruntime
            from sparseml.pytorch.utils import ModuleExporter

            with model._prevent_trainer_and_dataloaders_deepcopy():
                exporter = ModuleExporter(model.model, output_dir=output_dir)
                sample_batch = sample_batch if sample_batch is not None else model.example_input_array
                if sample_batch is None:
                    raise MisconfigurationException(
                        "To export the model, a sample batch must be passed via "
                        "``SparseMLCallback.export_to_sparse_onnx(model, output_dir, sample_batch=sample_batch)`` "
                        "or an ``example_input_array`` property within the LightningModule"
                    )

                # the following is adapted from @natuan and @spacemanidol
                sample_inputs = os.path.join(output_dir, "sample-inputs")
                sample_outputs = os.path.join(output_dir, "sample-outputs")
                os.makedirs(sample_inputs, exist_ok=True)
                os.makedirs(sample_outputs, exist_ok=True)

                forward_args_spec = inspect.getfullargspec(exporter._module.__class__.forward)
                try:
                    # assume sample_batch is a callable dictionary-type object
                    one_sample_input = collections.OrderedDict(
                        [
                            (f, sample_batch[f][0].long().reshape(1, -1))
                            for f in forward_args_spec.args
                            if f in sample_batch
                        ]
                    )
                except RuntimeError:
                    # assume sample_batch is a tensor
                    one_sample_input = sample_batch

                try:
                    exporter.export_onnx(sample_batch=one_sample_input, convert_qat=True, **kwargs)
                    exporter.export_onnx(
                        sample_batch=one_sample_input,
                        name="small_model.onnx",
                        export_params=False,
                        **kwargs,
                    )
                    onnx_file = os.path.join(output_dir, "model.onnx")

                except RuntimeError:
                    raise RuntimeError("Error exporting ONNX models and/or inputs/outputs")

                sess = onnxruntime.InferenceSession(onnx_file)

                num_samples = 0
                # add additional files for testing since this feature is very new
                if isinstance(one_sample_input, (collections.OrderedDict, dict)):
                    input_names = list(sample_batch.keys())
                    output_names = [o.name for o in sess.get_outputs()]
                    for input_vals in zip(*sample_batch.values()):
                        input_feed = {k: v.long().numpy() for k, v in zip(input_names, input_vals)}
                        output_vals = sess.run(output_names, {k: input_feed[k].reshape(1, -1) for k in input_feed})
                        output_dict = {name: numpy.squeeze(val) for name, val in zip(output_names, output_vals)}
                        file_idx = f"{num_samples}".zfill(4)
                        numpy.savez(f"{sample_inputs}/inp-{file_idx}.npz", **input_feed)
                        numpy.savez(f"{sample_outputs}/out-{file_idx}.npz", **output_dict)
                        num_samples += 1
                else:
                    output_names = [o.name for o in sess.get_outputs()]
                    input_feed = {"input": sample_batch.numpy()}
                    output_vals = sess.run(output_names, input_feed)
                    output_dict = {name: numpy.squeeze(val) for name, val in zip(output_names, output_vals)}
                    file_idx = f"{num_samples}".zfill(4)
                    numpy.savez(f"{sample_inputs}/inp-{file_idx}.npz", **input_feed)
                    numpy.savez(f"{sample_outputs}/out-{file_idx}.npz", **output_dict)

        def teardown(self, trainer: "pl.Trainer", pl_module: "pl.LightningModule", stage: Optional[str] = None) -> None:
            sample_batch = next(iter(trainer.train_dataloader))
            # if asked for output names, bert's ModelOutput gives two names
            # but when run, this the model only gives one output
            # workaround is just to force onnx to realize there is only one output
            output_names = ["logits"]
            self.export_to_sparse_onnx(
                output_dir=self.output_dir, model=pl_module, sample_batch=sample_batch, output_names=output_names
            )

class CUDACallback(Callback):
    def on_train_epoch_start(self, trainer, pl_module):
        # Reset the memory use counter
        torch.cuda.reset_peak_memory_stats(trainer.root_gpu)
        torch.cuda.synchronize(trainer.root_gpu)
        self.start_time = time.time()

    def on_train_epoch_end(self, trainer, pl_module, outputs):
        torch.cuda.synchronize(trainer.root_gpu)
        max_memory = torch.cuda.max_memory_allocated(trainer.root_gpu) / 2 ** 20
        epoch_time = time.time() - self.start_time

        max_memory = trainer.training_type_plugin.reduce(max_memory)
        epoch_time = trainer.training_type_plugin.reduce(epoch_time)

        rank_zero_info(f"Average Epoch time: {epoch_time:.2f} seconds")
        rank_zero_info(f"Average Peak memory {max_memory:.2f}MiB")

max_stars_repo_name: mariomeissner/lightning-transformers

max_issues_repo_name: mathemusician/lightning-transformers

max_forks_repo_name: mathemusician/lightning-transformers


text:
# These are all the modules we'll be using later. Make sure you can import them
# before proceeding further.
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
import tarfile
import tensorflow as tf
import cv2
from rotationInvariantNetwork import RICNN

import imageHelper as ih
from tensorflow.examples.tutorials.mnist import input_data

batch_size = 50
raw_image_size = 28
image_size = 46
num_labels = 10

mnist = input_data.read_data_sets('MNIST_data', one_hot=True)

# I don't think we need this, unless we're calling this code from some other python code
# that might have defined tf nodes? I guess this is just here to be safe? Because I think
# all this does is clear out previously constructed tf nodes, and I don't see any in the very few lines
# defined above this call
tf.reset_default_graph()

tf_train_dataset = tf.placeholder(
    tf.float32, shape=[None, image_size * image_size])
tf_train_labels = tf.placeholder(tf.float32, shape=[None, num_labels])
keep_prob = tf.placeholder(tf.float32)

#pad our valid dataset with zeros
valid_dataset = ih.padImages(mnist.validation.images, raw_image_size, image_size)
test_dataset = ih.padImages(mnist.test.images, raw_image_size, image_size)
valid_labels = mnist.validation.labels
test_labels = mnist.test.labels

numOutputClasses = 10
fullyConnectedLayerWidth = 1024
nn = RICNN([10, 100], fullyConnectedLayerWidth, numOutputClasses, 7)

logits = nn.setupNodes(tf_train_dataset, keep_prob, image_size, iterativelyMaxPool=False)

loss = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=logits, labels=tf_train_labels))
correct_prediction = tf.equal(tf.argmax(logits,1), tf.argmax(tf_train_labels,1))
accuracy = tf.reduce_sum(tf.cast(correct_prediction, tf.float32))

optimizer = tf.train.AdamOptimizer(0.0005).minimize(loss)
init_op = tf.global_variables_initializer()
saver = tf.train.Saver()

cols=image_size
rows=image_size
rawimgs = np.reshape(valid_dataset,(-1,image_size,image_size))
valid_rotimg = np.zeros(rawimgs.shape)
for ind, rawimg in enumerate(rawimgs):
  ang = np.random.randint(1, 360)
  M = cv2.getRotationMatrix2D((cols/2 - 0.5,rows/2 - 0.5),ang,1)
  valid_rotimg[ind] = cv2.warpAffine(rawimg,M,(rows,cols))
valid_rotimg = np.reshape(valid_rotimg, (-1,image_size**2))

cols=image_size
rows=image_size
rawimgs = np.reshape(test_dataset,(-1,image_size,image_size))
test_rotimg = np.zeros(rawimgs.shape)
for ind, rawimg in enumerate(rawimgs):
  ang = np.random.randint(1, 360)
  M = cv2.getRotationMatrix2D((cols/2 - 0.5,rows/2 - 0.5),ang,1)
  test_rotimg[ind] = cv2.warpAffine(rawimg,M,(rows,cols))
test_rotimg = np.reshape(test_rotimg, (-1,image_size**2))

with tf.Session() as sess:
  sess.run(init_op)

  for i in range(20000):
      # batch = [np.reshape(b, [28, 28]) for b in mnist.train.next_batch(batch_size=batch_size)[0]]
      batch_in, batch_out = mnist.train.next_batch(batch_size=batch_size)
      #pad our training dataset with zeros
      batch_in = ih.padImages(batch_in, raw_image_size, image_size)
      sess.run(optimizer, feed_dict = {tf_train_dataset: batch_in, tf_train_labels: batch_out, keep_prob: 0.8})

      if not i % 5000:
        print("checking confusion")
        confus = sess.run(tf.contrib.metrics.confusion_matrix(tf.argmax(logits,1),tf.argmax(valid_labels,1)), feed_dict={tf_train_dataset: valid_rotimg, tf_train_labels: valid_labels, keep_prob: 1.0})
        print(confus)
      if not i % 5000:
          print("checking accuracy...")
          numValid = len(valid_labels)
          accsum = sess.run(accuracy, feed_dict={tf_train_dataset: valid_dataset, tf_train_labels: valid_labels, keep_prob: 1.0})
          rotaccsum = sess.run(accuracy, feed_dict={tf_train_dataset: valid_rotimg, tf_train_labels: valid_labels, keep_prob: 1.0})

          accbat = accsum / numValid
          rotaccbat = rotaccsum / numValid
          print(i, accbat, rotaccbat)
          save_path = saver.save(sess, "/tmp/ricnn.ckpt")
          print("Model saved in file: %s" % save_path)
          """ls, d, i_ls, mu, sigm = sess.run([loss, dec, img_loss, mn, sd], feed_dict = {X_in: batch.images, Y: batch.labels, keep_prob: 1.0})
          plt.imshow(np.reshape(batch[0], [28, 28]), cmap='gray')
          plt.show()
          plt.imshow(d[0], cmap='gray')
          plt.show()
          print(i, ls, np.mean(i_ls))"""
  accsum = 0.0
  rotaccsum = 0.0
  numValid = 10000
  v_batch_size = 10
  for batchNum in range(numValid/v_batch_size):
    test_batch = range(batchNum * v_batch_size, min((batchNum + 1) * v_batch_size, numValid))
    accsum += sess.run(accuracy, feed_dict={tf_train_dataset: test_dataset[test_batch,:], tf_train_labels: test_labels[test_batch], keep_prob: 1.0})
    rotaccsum += sess.run(accuracy, feed_dict={tf_train_dataset: test_rotimg[test_batch,:], tf_train_labels: test_labels[test_batch], keep_prob: 1.0})

  accbat = accsum / numValid
  rotaccbat = rotaccsum / numValid
  print("For real homiez this is it! Our test data performance!")
  print(i, accbat, rotaccbat)

max_stars_repo_name: travers-rhodes/ricnn

max_issues_repo_name: travers-rhodes/ricnn

max_forks_repo_name: travers-rhodes/ricnn


text:
# Copyright (c) 2021 Jacob F. Valdez. Released under the MIT license.
"""Agent implementations for reinforcement learning."""

from __future__ import annotations

import random
from typing import Tuple, List, Mapping, Union, Callable

import numpy as np
from numpy.core.defchararray import _use_unicode

import jnumpy.core as jnp
import jnumpy.nn as jnn
from jnumpy.rl.types import BatchStep, Traj, reshape_traj_batch_size

class Agent:
    """Base class for agents.
    See docstrings on `forward`, `reward`, `train` for more information."""

    def __init__(self, policy: Callable, name: str = "Agent"):
        self.policy = policy
        self.name = jnp.NameScope.get_name("AGENT_NAMES", name)

    def forward(self, step: BatchStep) -> np.ndarray:
        """Generates an action for a given observation using `self.policy`.
        Override if you want to give your policy more information such as
        recurrent state or previous reward.

        Args:
            step (Step): last step output by the environment. This means the agent
                should feed `step.next_obs`, not `step.obs` to its policy. If the
                environment is multi-agent, then the `reward` attribute has already
                updated to reflect the reward for this agent by the driver.

        Returns:
            np.ndarray: The action to take
        """
        with jnp.NameScope(f"{self.name}.forward()"):
            action = self.policy(step.next_obs)
        return action

    def reward(self, traj: Traj) -> float:
        """Evaluates the mean reward of the agent on a trajectory.

        If your agent uses intrinsic rewards, be sure to add them in here.
        Do not introduce Q-values or predicted rewards here.

        Args:
            traj (Traj): timestep trajectory.

        Returns:
            float: Cumulative (sum) reward over the entire sequence.
        """
        return np.mean(sum(step.reward for step in traj) / len(traj))

    def train(self, traj: Traj):
        """Train your agent on a sequence of Timestep's.

        Args:
            traj (Traj): timestep trajectory.
        """
        with jnp.NameScope(f"{self.name}.train()"):
            self._train(traj)

    def _train(self, traj: Traj):
        """Private method implemented by subclasses to train on a sequence of Timestep's.

        Args:
            traj (Traj): timestep trajectory.
        """
        pass

    def __repr__(self) -> str:
        return self.__str__()

    def __str__(self) -> str:
        return self.name

class RandomAgent(Agent):
    """Takes a random action on each timestep.

    Example:
        >>> agent = RandomAgent(num_actions=5)
        >>> obs = jnp.Var(np.array([[1, 2], [3, 4]]))
        >>> step = Step(obs=None, next_obs=obs, action=None, reward=None, done=None, info=None)
        >>> agent.forward(step)

    Args:
        num_actions (int): The number of actions the agent can take.
    """

    def __init__(self, num_actions: int, name: str = "RandomAgent"):
        super(RandomAgent, self).__init__(policy=self._policy, name=name)
        self.num_actions = num_actions

    def _policy(self, obs: np.ndarray) -> np.ndarray:
        choices = np.random.randint(0, self.num_actions, (obs.shape[0],))
        onehots = np.eye(self.num_actions)[choices]
        return onehots

class RealDQN(Agent):
    """'Classic' Deep Q-learning agent.
    Implements the approach in https://arxiv.org/pdf/1312.5602.pdf.

    NOTE: This agent expects its encoder to output a per-column vector.
        I.E.: [B, H, W, C] --encoder--> [B, W, d_enc]

    This agent uses the following hyperparameters:
        - activation: activation function to use for the hidden layers.
        - hidden_size: hidden layer size.
        - discount: discount factor.
        - optimizer: optimizer to use.
        - epsilon_start: initial epsilon value.
        - min_epsilon: minimum epsilon value.
        - epsilon_decay: epsilon decay rate.
        - epoch: current epoch.
        - train_batch_size: batch size for training.

    Example:
        >>> test_batch_size = 5
        >>> test_num_actions = 7
        >>> test_encoder = jnn.Sequential([
                jnn.Conv2D(32, 3, 2, 'same', jnp.Relu),
                jnn.Conv2D(64, 3, 2, 'same', jnp.Relu),
                jnn.GlobalMaxPooling(1),
            ])  # [B, H, W, 2] -> [B, W, d_enc]
        >>> test_step = Step(
                obs=np.random.rand(test_batch_size, test_num_actions, test_num_actions, 2),
                next_obs=np.random.rand(test_batch_size, test_num_actions, test_num_actions, 2),
                action=np.random.rand(test_batch_size, test_num_actions),
                reward=np.random.rand(test_batch_size),
                done=False,
                info={})
        >>> test_hparams = {
                'hidden_size': 256,
                'activation': jnp.Relu,
                'optimizer': jnp.SGD(1e-3),
                'epsilon_start': 1.0,
                'epsilon_decay': 0.9,
                'min_epsilon': 0.1,
                'discount': 0.99,
                'epoch': 0,
                'train_batch_size': 32,}
        >>> test_agent = RealDQN(num_actions=test_num_actions, encoder=test_encoder, hparams=test_hparams)
        >>> test_agent.forward(test_step)

    Args:
        num_actions (int): The number of actions the agent can take.
        encoder (jnn.Layer): An image encoder [B, H, W, C] -> [B, W, d_enc]
        hparams (dict): Hyperparameters.
    """

    def __init__(
        self,
        num_actions: int,
        encoder: jnn.Layer,
        hparams: dict,
        name: str = "RealDQNAgent",
    ):

        super(RealDQN, self).__init__(policy=self._policy, name=name)

        self.num_actions = num_actions
        self.encoder = encoder  # [B, H, W, C] -> [B, W, d_enc]
        self.neck = jnn.Flatten(name="neck")
        self.head = jnn.Sequential(
            [
                jnn.Dense(
                    hparams["hidden_size"],
                    hparams["activation"],
                    name="hidden",
                ),
                jnn.Dense(1, jnp.Linear, name="output"),
            ],
            name="head",
        )  # [B, L+|A|] -> [B, 1]
        self.hparams = hparams

    def _policy(self, obs: np.ndarray, use_epsilon: bool = True) -> np.ndarray:

        B = obs.shape[0]

        # Maybe take greedy step
        epsilon = (
            self.hparams["epsilon_start"]
            * self.hparams["epsilon_decay"] ** self.hparams["epoch"]
        )
        epsilon = max(epsilon, self.hparams["min_epsilon"])
        if use_epsilon and random.random() < epsilon:
            indeces = np.random.randint(0, self.num_actions, (B,))
            return np.eye(self.num_actions)[indeces]

        # Otherwise take the action with the highest Q-value
        q_vals = np.zeros((B, self.num_actions))  # [B, self.num_actions]

        # prepare inputs
        obs_T = jnp.Var(obs)  # [B, H, W, 2]
        enc_T = self.neck(self.encoder(obs_T))  # [B, W*d_enc]

        for action_index in range(self.num_actions):

            # prepare action
            action_enc = np.eye(self.num_actions)[action_index]  # [num_actions]
            action_T = jnp.Var(
                np.repeat(action_enc[None, :], repeats=B, axis=0),
            )  # [B, 1]

            # run the network
            cat_T = jnp.Concat([enc_T, action_T], axis=1)  # [B, W*d_enc+num_actions]
            q_T = self.head(cat_T)  # [B, 1]

            # store q-values
            q_vals[:, action_index] = q_T.val[:, 0]

        # select the action with the highest Q-value
        action_indeces = np.argmax(q_vals, axis=1)  # [B]
        onehots = np.eye(self.num_actions)[action_indeces]  # [B, num_actions]
        return onehots

    def _train(self, traj: Traj):
        """Train your agent on a sequence of Timestep's.

        Args:
            traj (Traj): batched timestep trajectory.
        """

        optimizer = self.hparams["optimizer"]
        discount_T = jnp.Var(self.hparams["discount"], name="discount")  # []
        for step in reshape_traj_batch_size(traj, self.hparams["train_batch_size"]):

            obs_T = jnp.Var(step.obs, name="obs")  # [B, H, W, 2]
            action_T = jnp.Var(step.action, name="action")  # [B, A]
            next_obs_T = jnp.Var(step.next_obs, name="next_obs")  # [B, H, W, 2]
            next_action_T = jnp.Var(
                self.policy(step.next_obs, use_epsilon=False),
                name="next_action",
            )  # [B, A]
            r_T = jnp.Var(step.reward, name="reward")  # [B]

            # compute previous Q value using the actual (not necesarily optimal) action selected
            enc_T = self.neck(self.encoder(obs_T))  # [B, W*d_enc]
            cat_T = jnp.Concat([enc_T, action_T], axis=1)  # [B, d_enc+A]
            Qnow_T = self.head(cat_T)[:, 0]  # [B]
            reg_loss_now_T = self.encoder.loss + self.head.loss  # []

            # compute the maximum possible next step Q-value
            next_enc_T = self.neck(self.encoder(next_obs_T))  # [B, W*d_enc]
            next_cat_T = jnp.Concat([next_enc_T, next_action_T], axis=1)  # [B, d_enc+A]
            Qnext_T = self.head(next_cat_T)[:, 0]  # [B]
            reg_loss_next_T = self.encoder.loss + self.head.loss  # []

            # train the current policy to estimate new Q-value
            # i.e.: approx Qnew_T <= (1-lr)*Qnow_T + lr*(r_T+discount_T*Qnext_T)  # [B, 1]
            # using gradient descent (let lr=1 in the above eq; small updates are
            # ensured by small SGD lr instead)
            loss_T = (
                jnp.ReduceSum(
                    ((r_T + discount_T * jnp.StopGrad(Qnext_T)) - Qnow_T) ** 2, axis=0
                )
                + reg_loss_now_T
                + reg_loss_next_T
            )  # []
            optimizer.minimize(loss_T)

    def q_eval(self, obs: np.ndarray, action: np.ndarray) -> np.ndarray:
        """Evaluate the Q-value of a given state-action pair.

        Args:
            obs (np.ndarray): observation.
            action (np.ndarray): action.

        Returns:
            q_val (np.ndarray): Q-value of the given state-action pair.
        """

        # prepare inputs
        obs_T = jnp.Var(obs)  # [B, H, W, 2]
        action_T = jnp.Var(action)  # [B, A]
        enc_T = self.neck(self.encoder(obs_T))  # [B, W*d_enc]
        cat_T = jnp.Concat([enc_T, action_T], axis=1)  # [B, W*d_enc+A]
        Qnow_T = self.head(cat_T)[:, 0]  # [B]

        return Qnow_T.val

class CategoricalDQN(Agent):
    """Categorical deep Q-network agent.
    I never read the paper for this architecture, so my implementation may
    be different from the origonal researchers.

    NOTE: This agent expects its encoder to output a per-column vector.
        I.E.: [B, H, W, C] --encoder--> [B, W, d_enc]

    This agent uses the following hyperparameters:
        - activation: activation function to use for the hidden layers.
        - categorical_hidden_size: hidden layer size.
        - discount: discount factor.
        - optimizer: optimizer to use.
        - epsilon_start: initial epsilon value.
        - min_epsilon: minimum epsilon value.
        - epsilon_decay: epsilon decay rate.
        - epoch: current epoch.
        - train_batch_size: batch size for training.

    Example:
        >>> test_batch_size = 5
        >>> test_num_actions = 7
        >>> test_encoder = jnn.Sequential([
                jnn.Conv2D(32, 3, 2, 'same', jnp.Relu),
                jnn.Conv2D(64, 3, 2, 'same', jnp.Relu),
                jnn.GlobalMaxPooling(1),
            ])  # [B, H, W, 2] -> [B, W, d_enc]
        >>> test_step = Step(
                obs=np.random.rand(test_batch_size, test_num_actions, test_num_actions, 2),
                next_obs=np.random.rand(test_batch_size, test_num_actions, test_num_actions, 2),
                action=np.random.rand(test_batch_size, test_num_actions),
                reward=np.random.rand(test_batch_size),
                done=False,
                info={})
        >>> test_hparams = {
                'categorical_hidden_size': 32,
                'activation': jnp.Relu,
                'optimizer': jnp.SGD(1e-3),
                'epsilon_start': 1.0,
                'epsilon_decay': 0.9,
                'min_epsilon': 0.1,
                'discount': 0.99,
                'epoch': 0,
                'train_batch_size': 32,}
        >>> test_agent = CategoricalDQN(num_actions=test_num_actions, encoder=test_encoder, hparams=test_hparams)
        >>> test_agent.forward(test_step)

    Args:
        num_actions (int): The number of actions the agent can take.
        encoder (jnn.Layer): An image encoder [B, H, W, C] -> [B, W, d_enc]
        hparams (dict): Hyperparameters.
    """

    def __init__(
        self,
        num_actions: int,
        encoder: jnn.Layer,
        hparams: dict,
        name: str = "CategoricalDQN",
    ):

        self.num_actions = num_actions
        self.encoder = encoder  # [B, H, W, C] -> [B, W, d_enc]
        self.neck = jnn.Flatten(name="neck")
        self.head = jnn.Sequential(
            [
                jnn.Dense(hparams["categorical_hidden_size"], hparams["activation"]),
                jnn.Dense(1, jnp.Linear),
            ],
            name="head",
        )  # [B, W, d_enc] -> [B, A, 1]
        self.hparams = hparams

        super(CategoricalDQN, self).__init__(policy=self._policy, name=name)

    def _policy(self, obs: np.ndarray) -> np.ndarray:

        B = obs.shape[0]

        # Maybe take greedy step
        epsilon = (
            self.hparams["epsilon_start"]
            * self.hparams["epsilon_decay"] ** self.hparams["epoch"]
        )
        epsilon = max(epsilon, self.hparams["min_epsilon"])
        if random.random() < epsilon:
            indeces = np.random.randint(0, self.num_actions, (B,))
            return np.eye(self.num_actions)[indeces]

        # Otherwise estimate Q-values for all actions
        obs_T = jnp.Var(obs)  # [B, H, W, C]
        enc_T = self.encoder(obs_T)  # [B, W, d_enc]
        qvals_T = self.head(enc_T)  # [B, W, 1]
        return qvals_T[..., 0].val  # [B, W]

    def _train(self, traj: Traj):
        """Train your agent on a sequence of Timestep's.

        Args:
            traj (Traj): batched timestep trajectory.
        """
        optimizer = self.hparams["optimizer"]
        discount_T = jnp.Var(self.hparams["discount"], trainable=False)  # []
        for step in reshape_traj_batch_size(traj, self.hparams["train_batch_size"]):

            obs_T = jnp.Var(step.obs, name="obs")  # [B, H, W, 2]
            action_indeces = np.argmax(step.action, axis=1)  # [B]
            action_mask = np.eye(self.num_actions)[action_indeces]  # [B, A]
            action_mask_T = jnp.Var(action_mask, name="action_mask")  # [B, A]
            next_obs_T = jnp.Var(step.next_obs, name="next_obs")  # [B, H, W, 2]
            r_T = jnp.Var(step.reward, name="reward")  # [B]

            # compute previous Q value using the actual (not necesarily optimal) action selected
            enc_T = self.encoder(obs_T)  # [B, W, d_enc]
            qvals_T = self.head(enc_T)[..., 0]  # [B, W]
            Q_now_T = jnp.ReduceSum(qvals_T * action_mask_T, axis=1)  # [B]
            reg_loss_now_T = self.encoder.loss + self.head.loss  # []

            # compute the maximum possible next step Q-value
            enc_next_T = self.encoder(next_obs_T)  # [B, W, d_enc]
            qvals_T = self.head(enc_next_T)[..., 0]  # [B, W]
            Q_next_T = jnp.ReduceMax(qvals_T, axis=1)  # [B]
            # equivalent to: Q_next_T = ReduceMax(qvals_T, axis=1)
            reg_loss_next_T = self.encoder.loss + self.head.loss  # []

            # train the current policy to estimate new Q-value
            # i.e.: approx Qnew_T <= (1-lr)*Qnow_T + lr*(r_T+discount_T*Qnext_T)  # [B, 1]
            # using gradient descent (let lr=1 in the above; small updates are handled in the SGD step)
            # but only update the targets that were actually selected for action at `step_now`.
            loss_T = (
                jnp.ReduceSum(
                    ((r_T + discount_T * jnp.StopGrad(Q_next_T)) - Q_now_T) ** 2, axis=0
                )
                + reg_loss_now_T
                + reg_loss_next_T
            )  # []
            optimizer.minimize(loss_T)

    def q_eval(self, obs: np.ndarray, action: np.ndarray) -> np.ndarray:
        """Evaluate the Q-value of a given state-action pair.

        Args:
            obs (np.ndarray): observation.
            action (np.ndarray): action.

        Returns:
            q_val (np.ndarray): Q-value of the given state-action pair.
        """

        # prepare inputs
        obs_T = jnp.Var(obs)  # [B, H, W, 2]
        action_indeces = np.argmax(action, axis=1)  # [B]
        action_mask = np.eye(self.num_actions)[action_indeces]  # [B, A]
        action_mask_T = jnp.Var(action_mask, name="action_mask")  # [B, A]

        enc_T = self.encoder(obs_T)  # [B, W, d_enc]
        qvals_T = self.head(enc_T)[..., 0]  # [B, W]
        Q_T = jnp.ReduceSum(qvals_T * action_mask_T, axis=1)  # [B]

        return Q_T.val  # [B]

max_stars_repo_name: JacobFV/jumpy

max_issues_repo_name: JacobFV/jumpy

max_forks_repo_name: JacobFV/jumpy

reason: import numpy,from numpy


text:
# coding: utf-8

# ___
#
#
# ___

# # Tensorflow Project Exercise
# Let's wrap up this Deep Learning by taking a a quick look at the effectiveness of Neural Nets!
#
# We'll use the [Bank Authentication Data Set](https://archive.ics.uci.edu/ml/datasets/banknote+authentication) from the UCI repository.
#
# The data consists of 5 columns:
#
# * variance of Wavelet Transformed image (continuous)
# * skewness of Wavelet Transformed image (continuous)
# * curtosis of Wavelet Transformed image (continuous)
# * entropy of image (continuous)
# * class (integer)
#
# Where class indicates whether or not a Bank Note was authentic.
#
# This sort of task is perfectly suited for Neural Networks and Deep Learning! Just follow the instructions below to get started!

# ## Get the Data
#
# ** Use pandas to read in the bank_note_data.csv file **

# In[2]:

import pandas as pd
import numpy as np
df = pd.read_csv('bank_note_data.csv')

# ** Check the head of the Data **

# In[3]:

df.head()

# ## EDA
#
# We'll just do a few quick plots of the data.
#
# ** Import seaborn and set matplolib inline for viewing **

# In[4]:

import matplotlib.pyplot as plt
import seaborn as sns
get_ipython().run_line_magic('matplotlib', 'inline')

# ** Create a Countplot of the Classes (Authentic 1 vs Fake 0) **

# In[6]:

sns.countplot(x='Class',data=df)

# ** Create a PairPlot of the Data with Seaborn, set Hue to Class **

# In[8]:

sns.pairplot(df,hue='Class')

# ## Data Preparation
#
# When using Neural Network and Deep Learning based systems, it is usually a good idea to Standardize your data, this step isn't actually necessary for our particular data set, but let's run through it for practice!
#
# ### Standard Scaling
#
#

# In[9]:

from sklearn.preprocessing import StandardScaler

# **Create a StandardScaler() object called scaler.**

# In[10]:

scaler = StandardScaler()

# **Fit scaler to the features.**

# In[12]:

scaler.fit(df.drop('Class',axis=1))

# **Use the .transform() method to transform the features to a scaled version.**

# In[14]:

scaled_features = scaler.fit_transform(df.drop('Class',axis=1))

# **Convert the scaled features to a dataframe and check the head of this dataframe to make sure the scaling worked.**

# In[22]:

df_feat = pd.DataFrame(scaled_features,columns=df.columns[:-1])
df_feat.head()

# ## Train Test Split
#
# ** Create two objects X and y which are the scaled feature values and labels respectively.**

# In[24]:

X = df_feat

# In[25]:

y = df['Class']

# ** Use SciKit Learn to create training and testing sets of the data as we've done in previous lectures:**

# In[27]:

from sklearn.model_selection import train_test_split

# In[28]:

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

# # Tensorflow

# In[29]:

import tensorflow as tf

# ** Create a list of feature column objects using tf.feature.numeric_column() as we did in the lecture**

# In[30]:

df_feat.columns

# In[31]:

image_var = tf.feature_column.numeric_column("Image.Var")
image_skew = tf.feature_column.numeric_column('Image.Skew')
image_curt = tf.feature_column.numeric_column('Image.Curt')
entropy =tf.feature_column.numeric_column('Entropy')

# In[32]:

feat_cols = [image_var,image_skew,image_curt,entropy]

# ** Create an object called classifier which is a DNNClassifier from learn. Set it to have 2 classes and a [10,20,10] hidden unit layer structure:**

# In[33]:

classifier = tf.estimator.DNNClassifier(hidden_units=[10, 20, 10], n_classes=2,feature_columns=feat_cols)

# ** Now create a tf.estimator.pandas_input_fn that takes in your X_train, y_train, batch_size and set shuffle=True. You can play around with the batch_size parameter if you want, but let's start by setting it to 20 since our data isn't very big. **

# In[35]:

input_func = tf.estimator.inputs.pandas_input_fn(x=X_train,y=y_train,batch_size=40,shuffle=True)

# ** Now train classifier to the input function. Use steps=500. You can play around with these values if you want!**
#
# *Note: Ignore any warnings you get, they won't effect your output*

# In[36]:

classifier.train(input_fn=input_func,steps=700)

# ## Model Evaluation

# ** Create another pandas_input_fn that takes in the X_test data for x. Remember this one won't need any y_test info since we will be using this for the network to create its own predictions. Set shuffle=False since we don't need to shuffle for predictions.**

# In[42]:

pred_fn = tf.estimator.inputs.pandas_input_fn(x=X_test,batch_size=len(X_test),shuffle=False)

# ** Use the predict method from the classifier model to create predictions from X_test **

# In[43]:

note_predictions = list(classifier.predict(input_fn=pred_fn))

# In[44]:

note_predictions[0]

# In[45]:

final_preds  = []
for pred in note_predictions:
    final_preds.append(pred['class_ids'][0])

# ** Now create a classification report and a Confusion Matrix. Does anything stand out to you?**

# In[46]:

from sklearn.metrics import classification_report,confusion_matrix

# In[47]:

print(confusion_matrix(y_test,final_preds))

# In[48]:

print(classification_report(y_test,final_preds))

# ## Optional Comparison
#
# ** You should have noticed extremely accurate results from the DNN model. Let's compare this to a Random Forest Classifier for a reality check!**
#
# **Use SciKit Learn to Create a Random Forest Classifier and compare the confusion matrix and classification report to the DNN model**

# In[49]:

from sklearn.ensemble import RandomForestClassifier

# In[50]:

rfc = RandomForestClassifier(n_estimators=200)

# In[51]:

rfc.fit(X_train,y_train)

# In[52]:

rfc_preds = rfc.predict(X_test)

# In[53]:

print(classification_report(y_test,rfc_preds))

# In[54]:

print(confusion_matrix(y_test,rfc_preds))

# ** It should have also done very well, possibly perfect! Hopefully you have seen the power of DNN! **

# # Great Job!

max_stars_repo_path: 04-Tensorflow Project Exercise.py

max_issues_repo_path: 04-Tensorflow Project Exercise.py

max_forks_repo_path: 04-Tensorflow Project Exercise.py


text:
"""
Scripts looks at resampling methods to understanding the number of ensembles
for statistically significant responses

Notes
-----
    Author : Zachary Labe
    Date   : 16 August 2019
"""

### Import modules
import numpy as np
import matplotlib.pyplot as plt
import datetime
import read_MonthlyData as MO
import string

### Define time
now = datetime.datetime.now()
currentmn = str(now.month)
currentdy = str(now.day)
currentyr = str(now.year)
currenttime = currentmn + '_' + currentdy + '_' + currentyr
titletime = currentmn + '/' + currentdy + '/' + currentyr
print('\n' '----Plotting Monthly Map Comparison- %s----' % titletime)

### Alott time series (300 ensemble members)
year1 = 1701
year2 = 2000
years = np.arange(year1,year2+1,1)

###############################################################################
###############################################################################
###############################################################################
### Call arguments
directoryfigure = '/home/zlabe/Desktop/STRATOVARI/'
varnames = ['U10','Z50','U200','Z500','SLP','THICK','T2M','RNET']
experi = np.array([r'\textbf{$\bf{\Delta}$Pi}',r'\textbf{$\bf{\Delta}$Cu}'])
letters = list(string.ascii_lowercase)

### Functions
def readPolarCapData(varnames,level,levq,sliceq,period):
    """
    Read in all data for polar cap (>65N)
    """
    if varnames == 'U10': # data needed to calculate polar vortex definition
        lat,lon,lev,varf = MO.readExperiAll(varnames,'Future','surface')
        lat,lon,lev,varcu = MO.readExperiAll(varnames,'Current','surface')
        lat,lon,lev,varpi = MO.readExperiAll(varnames,'Past','surface')
    else: # averaged over the polar cap
        lat,lon,lev,varf = MO.readExperiAveVar(varnames,'Future','polar',level)
        lat,lon,lev,varcu = MO.readExperiAveVar(varnames,'Current','polar',level)
        lat,lon,lev,varpi = MO.readExperiAveVar(varnames,'Past','polar',level)

    if sliceq == True: # data from the polar cap
        print('\nSlicing grid at level=%s!' % levq)
        levqq = np.where(lev == levq)[0]
        varf = varf[:,:,levqq].squeeze()
        varcu = varcu[:,:,levqq].squeeze()
        varpi = varpi[:,:,levqq].squeeze()
    elif varnames == 'U10':
        latq = np.where((lat >= 59.5) & (lat <= 60.5))[0]
        varf = np.nanmean(varf[:,:,latq,:].squeeze(),axis=2)
        varcu = np.nanmean(varcu[:,:,latq,:].squeeze(),axis=2)
        varpi = np.nanmean(varpi[:,:,latq,:].squeeze(),axis=2)
    else:
        print('\nSurface variable!')

    ### Calculate time mean [returns 1D array]
    print('-------Period of time: %s!-------\n' % period)
    if period == 'NDJFM':
        varwf = np.nanmean(np.append(varf[:,-2:],varf[:,:3],axis=1),axis=1)
        varwcu = np.nanmean(np.append(varcu[:,-2:],varcu[:,:3],axis=1),axis=1)
        varwpi = np.nanmean(np.append(varpi[:,-2:],varpi[:,:3],axis=1),axis=1)
    elif period == 'DJF':
        varwf = np.nanmean(np.append(varf[:,-1:],varf[:,:2],axis=1),axis=1)
        varwcu = np.nanmean(np.append(varcu[:,-1:],varcu[:,:2],axis=1),axis=1)
        varwpi = np.nanmean(np.append(varpi[:,-1:],varpi[:,:2],axis=1),axis=1)
    elif period == 'JFM':
        varwf = np.nanmean(varf[:,0:3],axis=1)
        varwcu = np.nanmean(varcu[:,0:3],axis=1)
        varwpi = np.nanmean(varpi[:,0:3],axis=1)
    elif period == 'JF':
        varwf = np.nanmean(varf[:,0:2],axis=1)
        varwcu = np.nanmean(varcu[:,0:2],axis=1)
        varwpi = np.nanmean(varpi[:,0:2],axis=1)
    elif period == 'OND':
        varwf = np.nanmean(varf[:,-3:],axis=1)
        varwcu = np.nanmean(varcu[:,-3:],axis=1)
        varwpi = np.nanmean(varpi[:,-3:],axis=1)
    elif period == 'ND':
        varwf = np.nanmean(varf[:,-2:],axis=1)
        varwcu = np.nanmean(varcu[:,-2:],axis=1)
        varwpi = np.nanmean(varpi[:,-2:],axis=1)
    elif period == 'MAM':
        varwf = np.nanmean(varf[:,2:5],axis=1)
        varwcu = np.nanmean(varcu[:,2:5],axis=1)
        varwpi = np.nanmean(varpi[:,2:5],axis=1)
    elif period == 'Annual':
        varwf = np.nanmean(varf[:,:],axis=1)
        varwcu = np.nanmean(varcu[:,:],axis=1)
        varwpi = np.nanmean(varpi[:,:],axis=1)
    elif period == 'Dec':
        varwf = varf[:,-1]
        varwcu = varcu[:,-1]
        varwpi = varpi[:,-1]
    elif period == 'Jan':
        varwf = varf[:,0]
        varwcu = varcu[:,0]
        varwpi = varpi[:,0]
    elif period == 'Feb':
        varwf = varf[:,1]
        varwcu = varcu[:,1]
        varwpi = varpi[:,1]
    elif period == 'Mar':
        varwf = varf[:,2]
        varwcu = varcu[:,2]
        varwpi = varpi[:,2]

    ### Check for missing data! this occurs between 100,200,300 ensembles
    varwf[np.where(varwf < -1e20)] = np.nan
    varwcu[np.where(varwcu < -1e20)] = np.nan
    varwpi[np.where(varwpi < -1e20)] = np.nan

    ### Calculate anomalies [1D arrays]
    anomcu = varwf - varwcu
    anompi = varwf - varwpi

    return lat,lon,anomcu,anompi

def calcSubEns(size,var):
    """
    Calculate random subsamples from large ensemble
    """
    subens = np.empty((var.shape[0]))
    subens[:] = np.nan
    for i in range(2,var.shape[0]): # min ensemble size of 2
        superens = np.random.choice(var,size=(size,i),replace=True)
        supermean = np.nanmean(superens,axis=1)
        subens[i] = abs(np.nanmax(supermean) - np.nanmin(supermean))
    print('Completed: calcSubEns function!')
    return subens

def readSubEns(variable,level,levq,levslice,period):
    """
    Process all data
    """
    lat,lon,anomcu,anompi = readPolarCapData(variable,level,levq,
                                             levslice,period)

    ### Perform combinations (100,000 random samples)
    subens_cu = calcSubEns(100000,anomcu)
    subens_pi = calcSubEns(100000,anompi)

    print('\n<<<<Completed: finished readings %s data combos!>>>>' % variable)
    return subens_cu,subens_pi

###############################################################################
###############################################################################
###############################################################################
### Read all functions
periodm = 'Annual'
pv_cu,pv_pi = readSubEns('U10','surface',None,False,periodm)
u10_cu,u10_pi = readSubEns('U','profile',10,True,periodm)
z50_cu,z50_pi = readSubEns('GEOP','profile',50,True,periodm)

slp_cu,slp_pi = readSubEns('SLP','surface',None,False,periodm)
z500_cu,z500_pi = readSubEns('Z500','surface',None,False,periodm)

t2_cu,t2_pi = readSubEns('T2M','surface',None,False,periodm)
p_cu,p_pi = readSubEns('P','surface',None,False,periodm)

### Create subplot files and arguement
suball = np.array([[u10_cu,u10_pi],[z50_cu,z50_pi],[slp_cu,slp_pi],
                   [z500_cu,z500_pi],[t2_cu,t2_pi],[p_cu,p_pi]])
xaxis = np.array([np.arange(0,301,50),np.arange(0,301,50),np.arange(0,301,50),
                  np.arange(0,301,50),np.arange(0,301,50),np.arange(0,301,50)])
yaxis = np.array([np.arange(0,51,10),np.arange(0,501,100),
                  np.arange(0,11,2),np.arange(0,201,25),
                  np.arange(0,8,1),np.round(np.arange(0,0.51,0.1),2)])
variablelist = np.array(['U10','Z50','SLP','Z500','T2M','P'])
ylabels = np.array([r'm/s',r'm',r'hPa',r'm',r'$\bf{^{\circ}}$C',r'mm/day'])
explist = np.array([[r'{\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pi}'],
                    [r'{\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pi}'],
                    [r'{\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pi}'],
                    [r'{\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pi}'],
                    [r'{\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pi}'],
                    [r'{\textbf{$\bf{\Delta}$Cu}',r'\textbf{$\bf{\Delta}$Pi}']])

###############################################################################
###############################################################################
###############################################################################
### Create subplots of ensembles
plt.rc('text',usetex=True)
plt.rc('font',**{'family':'sans-serif','sans-serif':['Avant Garde']})

def adjust_spines(ax, spines):
    for loc, spine in ax.spines.items():
        if loc in spines:
            spine.set_position(('outward', 5))
        else:
            spine.set_color('none')
    if 'left' in spines:
        ax.yaxis.set_ticks_position('left')
    else:
        ax.yaxis.set_ticks([])

    if 'bottom' in spines:
        ax.xaxis.set_ticks_position('bottom')
    else:
        ax.xaxis.set_ticks([])

fig = plt.figure(figsize=(8,5))
for i in range(len(suball)):
    ax = plt.subplot(2,3,i+1)

    if i < 3:
        adjust_spines(ax, ['left', 'bottom'])
        ax.spines['top'].set_color('none')
        ax.spines['right'].set_color('none')
        ax.spines['left'].set_color('dimgrey')
        ax.spines['bottom'].set_color('w')
        ax.spines['left'].set_linewidth(2)
        ax.spines['bottom'].set_linewidth(0)
        ax.tick_params('y',length=4,width=2,which='major',color='dimgrey',
                       pad=0.3)
        ax.tick_params('x',length=0,width=0,which='major',color='w',
                       pad=0,labelbottom=False)
    elif i >= 3:
        adjust_spines(ax, ['left', 'bottom'])
        ax.spines['top'].set_color('none')
        ax.spines['right'].set_color('none')
        ax.spines['left'].set_color('dimgrey')
        ax.spines['bottom'].set_color('dimgrey')
        ax.spines['left'].set_linewidth(2)
        ax.spines['bottom'].set_linewidth(2)
        ax.tick_params('both',length=4,width=2,which='major',color='dimgrey',
                       pad=0.3)

    plt.axvline(100,color='dimgrey',linestyle='--',dashes=(1,0.3),
                linewidth=2)
    plt.hlines(suball[i][0][100],xmin=100,xmax=300,color='dimgrey',
              linestyle='--',linewidth=0.8)
    plt.plot(suball[i][0],linestyle='-',linewidth=0.8,color='rebeccapurple')
    plt.plot(suball[i][1],linestyle='-',linewidth=0.8,color='teal')

    if i == 0:
       plt.plot(pv_pi,linestyle='-',linewidth=0.8,color='darkorange')
       plt.plot(pv_cu,linestyle='-',linewidth=0.8,color='dodgerblue')

    plt.xticks(xaxis[i],list(map(str,xaxis[i])),fontsize=6)
    plt.xlim([0,xaxis[i].max()])
    plt.yticks(yaxis[i],list(map(str,yaxis[i])),fontsize=6)
    plt.ylim([0,yaxis[i].max()])

    if i == 4:
        plt.xlabel(r'\textbf{Number of Ensembles}',color='k',fontsize=10)
    ax.annotate(r'\textbf{[%s]}' % variablelist[i],xy=(0,0),
            xytext=(0.75,0.9),xycoords='axes fraction',
            color='dimgrey',fontsize=17)
    plt.ylabel(r'\textbf{Uncertainty [%s]}' % (ylabels[i]),
                         color='k',fontsize=10)

plt.tight_layout()

plt.savefig(directoryfigure + 'Ensemble_Subsampling_%s.png' % periodm,dpi=900)

max_stars_repo_name: zmlabe/StratoVari

max_issues_repo_name: zmlabe/StratoVari

max_forks_repo_name: zmlabe/StratoVari