text:
from fuse.server.immunespace.dispatcher import GetObject
import json
import os
import pytest
import numpy as np

# this takes about 20s to return
# go get a session id and group objectIdfrom immunespace for  user for this to work:
# https://www.immunespace.org/security/externalToolsView.view?returnUrl=%2Fproject%2FStudies%2Fbegin.view%3F
#g_debug = True
g_debug = False
def test_GetObject():

    if os.getenv('TEST_LIBRARY') == "0":
        pytest.skip("Only testing docker container")

    objectId = os.getenv('GROUP')
    username = os.getenv('USERNAME')
    sess = os.getenv('APIKEY')

    obj = {
        "id": objectId,
        "resourceType": "eset",
        "resource": GetObject(objectId,sess,username)
    }

    with open('tests/expected/test_1.json', 'r', encoding='utf-8') as f:
        expected = json.load(f)

    #make smaller chunks for easier debugging
    if(g_debug):
        max_subjs=3
        max_pheno=4
        max_genes=5
        obj["resource"]["exprs"] = np.array(obj["resource"]["exprs"])[0:max_genes,0:max_subjs].tolist() # 3 genes, 2 subjects
        obj["resource"]["featureNames"] = np.array(obj["resource"]["featureNames"])[0:max_genes].tolist()
        obj["resource"]["pData"] = np.array(obj["resource"]["pData"])[0:max_pheno,0:max_subjs].tolist() # 4 phenoetypes, 2 subjects

        expected["resource"]["exprs"] = np.array(expected["resource"]["exprs"])[0:max_genes,0:max_subjs].tolist() # 3 genes, 2 subjects
        expected["resource"]["featureNames"] = np.array(expected["resource"]["featureNames"])[0:max_genes].tolist()
        expected["resource"]["pData"] = np.array(expected["resource"]["pData"])[0:max_pheno,0:max_subjs].tolist() # 4 phenoetypes, 2 subjects

    # Uncomment this to capture output:
    #with open('tests/test_1.out.json', 'w', encoding='utf-8') as f:
    #     json.dump(obj, f, ensure_ascii=False, indent=4, sort_keys=True)

    objs = json.dumps(obj, ensure_ascii=False, indent=4, sort_keys=True)
    expecteds = json.dumps(expected, ensure_ascii=False, indent=4, sort_keys=True)

    if(g_debug):
        print("obj:")
        print(obj["resource"]["exprs"])
        #print("expected:")
        #print(expected["resource"]["exprs"])

    # xxx sort the keyes, then copy this to test_func.py
    assert objs == expecteds

ext: py

max_stars_repo_name: RENCI/fuse-server-immunespace

max_issues_repo_name: RENCI/fuse-server-immunespace

max_forks_repo_name: RENCI/fuse-server-immunespace

reason: import numpy


text:
import numpy as np

a = [1, 6, 4, 2, 6, 4, 1]
b = [1, 5, 9]
c = [6, 4, 1]

# a = 3,3,3,2,2,2,1,1,1
# b = 3,3,2,2
# c = [4]

data = [a, b, c]

def u1b1(data):
    n = len(data)
    index_n = n
    values = [d[0] + np.sqrt(2 * np.log(n)) for d in data]
    compare_values = [[values[i], i] for i in range(n)]
    already_calculate = [[data[i][0]] for i in range(n)]
    lengths = [len(d) for d in data]
    maxn = sum(lengths)
    data = [d[1:] for d in data]

    while n < maxn:
        n += 1
        i = max(compare_values)[1]
        for k in range(index_n):
            if i == k:
                already_calculate[k].append(data[k][0])
                data[k] = data[k][1:]
                if not data[k]:
                    compare_values[k][0] = -float('inf')
            values[k] = np.mean(already_calculate[k]) + np.sqrt(2*np.log(n) / len(already_calculate[k]))
            if compare_values[k][0] != -float('inf'):
                compare_values[k][0] = values[k]

    print(values)

u1b1(data)


text:
import h5py
import cv2
import numpy as np
import copy
import os
import imageio
import io
from multiprocessing import Pool, Manager
from visual_mpc.utils.sync import ManagedSyncCounter
import random
import functools
from tqdm import tqdm

MANDATORY_KEYS = ['camera_configuration', 'policy_desc', 'bin_type', 'bin_insert', 'contains_annotation',
                                'robot', 'gripper', 'background', 'action_space', 'object_classes', 'primitives', 'camera_type']

def serialize_image(img):
    assert img.dtype == np.uint8, "Must be uint8!"
    return cv2.imencode('.jpg', img)[1]

def serialize_video(imgs, temp_name_append):
    mp4_name = './temp{}.mp4'.format(temp_name_append)
    try:
        assert imgs.dtype == np.uint8, "Must be uint8 array!"
        assert not os.path.exists(mp4_name), "file {} exists!".format(mp4_name)
        # this is a hack to ensure imageio succesfully saves as a mp4 (instead of getting encoding confused)
        writer = imageio.get_writer(mp4_name)
        [writer.append_data(i[:, :, ::-1]) for i in imgs]
        writer.close()

        f = open(mp4_name, 'rb')
        buf = f.read()
        f.close()
    finally:
        if os.path.exists(mp4_name):
            os.remove(mp4_name)

    return np.frombuffer(buf, dtype=np.uint8)

def save_dict(data_container, dict_group, video_encoding, t_index):
    for k, d in data_container.items():
                    if 'images' == k:
                        T, n_cams = d.shape[:2]
                        dict_group.attrs['n_cams'] = n_cams

                        for n in range(n_cams):
                            cam_group = dict_group.create_group("cam{}_video".format(n))
                            if video_encoding == 'mp4':
                                data = cam_group.create_dataset("frames", data=serialize_video(d[:, n], t_index))
                                data.attrs['shape'] = d[0, n].shape
                                data.attrs['T'] = d.shape[0]
                                data.attrs['image_format'] = 'RGB'
                            elif video_encoding == 'jpeg':
                                for t in range(T):
                                    data = cam_group.create_dataset("frame{}".format(t), data=serialize_image(d[t, n]))
                                    data.attrs['shape'] = d[t, n].shape
                            else:
                                raise ValueError
                    elif 'image' in k:
                        data = dict_group.create_dataset(k, data=serialize_image(d))
                        data.attrs['shape'] = d.shape
                    else:
                        dict_group.create_dataset(k, data=d)

def save_hdf5(filename, env_obs, policy_out, agent_data, meta_data, video_encoding='mp4', t_index=None):
    if t_index is None:
        t_index = random.randint(0, 9999999)
    # meta-data includes calibration "number", policy "type" descriptor, environment bounds
    with h5py.File(filename, 'w') as f:
        [save_dict(data_container, f.create_group(name), video_encoding, t_index) for data_container, name in zip([env_obs, agent_data], ['env', 'misc'])]

        policy_dict = {}
        first_keys = list(policy_out[0].keys())
        for k in first_keys:
            assert all([k in p for p in policy_out[1:]]), "hdf5 format requires keys must be uniform across time!"
            policy_dict[k] = np.concatenate([p[k][None] for p in policy_out], axis=0)
        save_dict(policy_dict, f.create_group('policy'), video_encoding, t_index)

        meta_data_group = f.create_group('metadata')
        for mandatory_key in MANDATORY_KEYS:
            meta_data_group.attrs[mandatory_key] = meta_data.pop(mandatory_key)

        for k in meta_data.keys():
            meta_data_group.attrs[k] = meta_data[k]

def save_worker(traj_data, cntr, group_name=''):
    t_index = random.randint(0, 9999999)
    t, meta_data = traj_data

    try:
        env_obs = pkl.load(open('{}/obs_dict.pkl'.format(t), 'rb'), encoding='latin1')
        if meta_data['contains_annotation']:
            env_obs['bbox_annotations'] = pkl.load(open('{}/annotation_array.pkl'.format(t), 'rb'), encoding='latin1')
        n_cams = len(glob.glob('{}/images*'.format(t)))
        if n_cams:
            T = min([len(glob.glob('{}/images{}/*.jpg'.format(t, i))) for i in range(n_cams)])
            height, width = cv2.imread('{}/images0/im_0.jpg'.format(t)).shape[:2]
            env_obs['images'] = np.zeros((T, n_cams, height, width, 3), dtype=np.uint8)

            for n in range(n_cams):
                for time in range(T):
                    env_obs['images'][time, n] = cv2.imread('{}/images{}/im_{}.jpg'.format(t, n, time))

        policy_out = pkl.load(open('{}/policy_out.pkl'.format(t), 'rb'), encoding='latin1')
        agent_data = pkl.load(open('{}/agent_data.pkl'.format(t), 'rb'), encoding='latin1')

        def store_in_metadata_if_exists(key):
            if key in agent_data:
                meta_data[key] = agent_data.pop(key)
        [store_in_metadata_if_exists(k) for k in ['goal_reached', 'term_t']]

        c = cntr.ret_increment
        save_hdf5('{}/{}traj{}.hdf5'.format(args.output_folder, group_name, c), env_obs, policy_out, agent_data, meta_data, video_encoding, t_index)
        return True
    except (FileNotFoundError, NotADirectoryError):
        return False

if __name__ == '__main__':
    import argparse
    import glob
    import json
    import random
    import sys
    import os
    import shutil
    import math
    if sys.version_info[0] == 2:
        import cPickle as pkl
        input_fn = raw_input
    else:
        import pickle as pkl
        input_fn = input

    parser = argparse.ArgumentParser(description="converts dataset from pkl format to hdf5")
    parser.add_argument('input_folder', type=str, help='where raw files are stored')
    parser.add_argument('output_folder', type=str, help='where to save')
    parser.add_argument('--output_group_name', type=str, default='', help='name to prepend in front of trajs')
    parser.add_argument('--video_jpeg_encoding', action='store_true', default=False, help='uses jpeg encoding for video frames instead of mp4')
    parser.add_argument('--counter', type=int, help='where to start counter', default=0)
    parser.add_argument('--n_workers', type=int, help='number of multi-threaded workers', default=1)
    args = parser.parse_args()

    assert args.n_workers >= 1, "can't have less than 1 worker thread!"
    args.input_folder, args.output_folder = [os.path.expanduser(x) for x in (args.input_folder, args.output_folder)]
    if not os.path.exists(args.output_folder):
        os.makedirs(args.output_folder)
    elif input_fn('path {} exists, should folder be deleted? (y/n): '.format(args.output_folder)) == 'y':
        shutil.rmtree(args.output_folder)
        os.makedirs(args.output_folder)

    if args.video_jpeg_encoding:
        video_encoding = 'jpeg'
    else:
        video_encoding = 'mp4'
        if len(glob.glob('./temp*.mp4')) != 0:
            print("Please delete all temp*.mp4 files! (needed for saving)")
            raise EnvironmentError

    traj_groups = glob.glob(args.input_folder + "*/*/*/*")
    print('found {} traj groups!'.format(len(traj_groups)))

    trajs, annotations_loaded = [], 0
    for group in traj_groups:
        meta_data_dict = json.load(open('{}/hparams.json'.format(group), 'r'))
        group_trajs = glob.glob('{}/*'.format(group))
        for t in group_trajs:
            traj_meta_data = copy.deepcopy(meta_data_dict)
            traj_meta_data['object_batch'] = group
            if os.path.exists('{}/annotation_array.pkl'.format(t)):
                traj_meta_data['contains_annotation'] = True
                annotations_loaded += 1
            else:
                traj_meta_data['contains_annotation'] = False

            if isinstance(traj_meta_data['object_classes'], str):
                traj_meta_data['object_classes'] = traj_meta_data['object_classes'].split("+")

            assert all([k in traj_meta_data for k in MANDATORY_KEYS]), 'metadata for {} is missing keys!'.format(t)
            assert isinstance(traj_meta_data['object_classes'], list), "did not split object classes!"
            assert all([isinstance(x, str) for x in traj_meta_data['object_classes']]), 'object classes is not a string!'

            trajs.append((t, traj_meta_data))
    random.shuffle(trajs)

    print('Loaded {} trajectories with {} annotations!'.format(len(trajs), annotations_loaded))

    cntr = ManagedSyncCounter(Manager(), args.counter)
    if args.n_workers == 1:
        saved = 0
        for t in tqdm(trajs):
            saved += save_worker(t, cntr, args.output_group_name)

        print('saved {} total trajs'.format(saved))
    else:
        map_fn = functools.partial(save_worker, cntr=cntr, group_name=args.output_group_name)
        p = Pool(args.n_workers)
        print('saved {} total trajs'.format(sum(tqdm(p.imap_unordered(map_fn, trajs), total=len(trajs)))))


text:
import matplotlib.pyplot as plt
import numpy as np
import os

def plot_metrics(
    comp_dist, gt_dist, braking_action, p_values, fontsize=10,
    figsize=(8, 8), save=False, save_path=None, caption='fig', mode='in'
):
    title = 'In-distribution episode' if mode == 'in' else 'Out-of-distribution episode'
    fontdict = {
        'fontsize': fontsize
    }
    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=figsize)
    fig.suptitle(title, fontsize=12)
    fig.tight_layout(pad=3.5)

    # Plot distances
    axes[0, 0].plot(comp_dist, label='Predicted Distance')
    axes[0, 0].plot(gt_dist, label='GroundTruth Distance')
    axes[0, 0].set_xlabel('Simulation steps')
    axes[0, 0].set_ylabel('Distance(m)')
    axes[0, 0].set_title('Distance Comparison (Computed vs GT)', fontdict=fontdict)
    axes[0, 0].legend(loc='upper right')

    # Plot the absolute error between the distances and Ground truth
    axes[0, 1].set_yticks(np.linspace(0, 50, num=5))
    axes[0, 1].set_autoscaley_on(False)
    axes[0, 1].plot(np.abs(comp_dist - gt_dist))
    axes[0, 1].set_xlabel('Simulation steps')
    axes[0, 1].set_ylabel('Distance(m)')
    axes[0, 1].set_title('Absolute distance Error', fontdict=fontdict)

    # Plot the braking action
    axes[1, 0].plot(braking_action)
    axes[1, 0].set_xlabel('Simulation steps')
    axes[1, 0].set_ylabel('Brake intensity')
    axes[1, 0].set_title('Braking action policy', fontdict=fontdict)

    # Plot the p-values
    axes[1, 1].set_yticks(np.linspace(0, 1, num=5))
    axes[1, 1].set_autoscaley_on(False)
    axes[1, 1].scatter(np.arange(p_values.shape[-1]), p_values)
    axes[1, 1].set_xlabel('Simulation steps')
    axes[1, 1].set_ylabel('P-values')
    axes[1, 1].set_title('P-values', fontdict=fontdict)

    plt.show()

    if save:
        if not save_path:
            raise ValueError('Save path must be provided when saving the figure')
        os.makedirs(save_path, exist_ok=True)
        plt.savefig(os.path.join(save_path, f'{caption}.png'), bbox_inches='tight')


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 27 08:04:07 2019

@author: abibeka
# Get Travel time on LPGA and I-95 (After Validation)
"""

#0.0 Housekeeping. Clear variable space
from IPython import get_ipython  #run magic commands
ipython = get_ipython()
ipython.magic("reset -f")
ipython = get_ipython()

import os
import pandas as pd
import numpy as np
import subprocess
import seaborn as sns
import matplotlib.pyplot as plt

def TTSegName(x):
    TTSeg= {1: 'EB',
            2: 'WB',
            3: 'EB LPGA (Tomoka Rd to I-95 SB Ramp)',
            4: 'EB LPGA (I-95 SB Ramp to I-95 NB Ramp)',
            5: 'EB LPGA (I-95 NB Ramp to Technology Blvd)',
            6: 'EB LPGA (Technology Blvd to Willamson Blvd)',
            7: 'EB LPGA (Willamson Blvd to Clyde-Morris Blvd)',
            8: 'WB LPGA (Clyde-Morris Blvd to Willamson Blvd)',
            9: 'WB LPGA (Willamson Blvd to Technology Blvd)',
            10:'WB LPGA (Technology Blvd to I-95 NB Ramp)',
            11:'WB LPGA (I-95 NB Ramp to I-95 SB Ramp)',
            12:'WB LPGA (I-95 SB Ramp to Tomoka Rd)',
            13:'SB I-95',
            14:'NB I-95',
            15:'SB I-95 (SR-40 to SB Off-Ramp)',
            16:'SB I-95 (SB Off-Ramp to SB Loop-Ramp)',
            17:'SB I-95 (SB Loop-Ramp to SB On-Ramp)',
            18:'SB I-95 (SB On-Ramp to US-92)',
            19:'NB I-95 (US-92 to NB Off-Ramp)',
            20:'NB I-95 (NB Off-Ramp to NB Loop-Ramp)',
            21:'NB I-95 ( NB Loop-Ramp to NB On-Ramp)',
            22:'NB I-95 (NB On-Ramp to SR-40)'}
    if x < 23:
        Nm = TTSeg[x]
    else:
        Nm = None
    return Nm

def PreProcessVissimTT(file, SegIN):
    '''
    file : VISSIM Results file
    SegIN : Specify segments for which results are needed
    Summarize Vissim Travel time results
    '''
    ExistingAMDat=pd.read_csv(file,sep =';',skiprows=17)
    ExistingAMDat.columns
    ExistingAMDat.rename(columns={'TRAVTM(ALL)':'VissimTT','VEHS(ALL)':'Veh','DISTTRAV(ALL)':'Len'},inplace=True)
    mask=ExistingAMDat["$VEHICLETRAVELTIMEMEASUREMENTEVALUATION:SIMRUN"]=="AVG"
    ExistingAMDat = ExistingAMDat[mask]
    ExistingAMDat["TTSegNm"]=ExistingAMDat['VEHICLETRAVELTIMEMEASUREMENT'].apply(TTSegName)
    ExistingAMDat.TIMEINT = pd.Categorical(ExistingAMDat.TIMEINT,['900-1800','1800-2700','2700-3600','3600-4500',
                                        '4500-5400','5400-6300','6300-7200','7200-8100',
                                        '8100-9000','9000-9900','9900-10800','10800-11700'])

    mask = ExistingAMDat.TTSegNm.isin(SegIN)
    ExistingAMDat = ExistingAMDat[mask]
    ExistingAMDat.loc[:,"VissimSMS"] = (ExistingAMDat['Len']/ExistingAMDat['VissimTT']/1.47).round(1)
    mask1 = ~ExistingAMDat['VEHICLETRAVELTIMEMEASUREMENT'].isin([13,14])
    TTData = ExistingAMDat[mask1]
    TTData.loc[:,"VissimSMS"] = TTData.groupby(['TTSegNm'])['VissimSMS'].transform(lambda v: v.ffill())

    return(TTData)

# Get the Weighted Density by Segments
#Debug:
#file = NoBuildAMfi
#SegKeyVal = TTSegLaneDat
#SegIN = SegIn1
def PreProcessVissimDensity(file, SegKeyVal):
    '''
    file : VISSIM Results file
    SegKeyVal : Key value pair for segment # and TT segment name
    Summarize Vissim Density results
    '''
    ExistingAMDat=pd.read_csv(file,sep =';',skiprows=17)
    ExistingAMDat.columns
    ExistingAMDat.rename(columns={'TRAVTM(ALL)':'VissimTT','VEHS(ALL)':'Veh','DISTTRAV(ALL)':'Len'},inplace=True)
    mask=ExistingAMDat["$VEHICLETRAVELTIMEMEASUREMENTEVALUATION:SIMRUN"]=="AVG"
    ExistingAMDat = ExistingAMDat[mask]
    ExistingAMDat = SegKeyVal.merge(ExistingAMDat,left_on=['SegNO'],right_on = ['VEHICLETRAVELTIMEMEASUREMENT'],how= 'left')
    ExistingAMDat.TIMEINT = pd.Categorical(ExistingAMDat.TIMEINT,['900-1800','1800-2700','2700-3600','3600-4500',
                                        '4500-5400','5400-6300','6300-7200','7200-8100',
                                        '8100-9000','9000-9900','9900-10800','10800-11700'])
    ExistingAMDat.loc[:,"VissimSMS"] = (ExistingAMDat['Len']/ExistingAMDat['VissimTT']/1.47).round(1)
    #Get flow rate and density
    ExistingAMDat.loc[:,'FlowRate'] = ExistingAMDat.Veh *4
    ExistingAMDat.loc[:,'FlowRate_PerLane'] = ExistingAMDat.Veh *4/ ExistingAMDat.NumLanes
    ExistingAMDat.loc[:,"VissimSMS"] = (ExistingAMDat['Len']/ExistingAMDat['VissimTT']/1.47).round(1)
    ExistingAMDat.loc[:,'Density'] = (ExistingAMDat.FlowRate/ ExistingAMDat.VissimSMS).round(1)
    ExistingAMDat.loc[:,'Density_Len_Lanes'] = ExistingAMDat.Density *ExistingAMDat.Len
    ExistingAMDat.loc[:,'Len_Lanes'] = ExistingAMDat.Len * ExistingAMDat.NumLanes
    ExistingAMDat.columns
    DensityData = ExistingAMDat.groupby(['TIMEINT','SegName']).aggregate({'Len_Lanes':'sum','Density_Len_Lanes':'sum','VissimSMS':'mean','FlowRate_PerLane':'mean'}).reset_index()
    #DensityData = ExistingAMDat.groupby(['TIMEINT','SegName'])['Len_Lanes','Density_Len_Lanes'].sum().reset_index()
    DensityData.loc[:,'WeightedDensity'] = (DensityData.Density_Len_Lanes/ DensityData.Len_Lanes).round(1)
    DensityData.loc[:,"WeightedDensity"] = DensityData.groupby(['SegName'])['WeightedDensity'].transform(lambda v: v.ffill())
    DensityData.loc[:,"VissimSMS"] = DensityData.groupby(['SegName'])['VissimSMS'].transform(lambda v: v.ffill())
    DensityData = DensityData[['TIMEINT','SegName','WeightedDensity']]
    return(DensityData)

I95_Segs = [     'NB I-95 (US-92 to NB Off-Ramp)',
                'NB I-95 (NB Off-Ramp to NB Loop-Ramp)',
                'NB I-95 ( NB Loop-Ramp to NB On-Ramp)',
                'NB I-95 (NB On-Ramp to SR-40)',
                'SB I-95 (SR-40 to SB Off-Ramp)',
                'SB I-95 (SB Off-Ramp to SB Loop-Ramp)',
                'SB I-95 (SB Loop-Ramp to SB On-Ramp)',
                'SB I-95 (SB On-Ramp to US-92)'
               ]

LPGASeg=['EB LPGA (Tomoka Rd to I-95 SB Ramp)',
    'EB LPGA (I-95 SB Ramp to I-95 NB Ramp)',
    'EB LPGA (I-95 NB Ramp to Technology Blvd)',
    'EB LPGA (Technology Blvd to Willamson Blvd)',
    'EB LPGA (Willamson Blvd to Clyde-Morris Blvd)',
    'WB LPGA (Clyde-Morris Blvd to Willamson Blvd)',
    'WB LPGA (Willamson Blvd to Technology Blvd)',
    'WB LPGA (Technology Blvd to I-95 NB Ramp)',
    'WB LPGA (I-95 NB Ramp to I-95 SB Ramp)',
    'WB LPGA (I-95 SB Ramp to Tomoka Rd)']

# VISSIM File
#*********************************************************************************
PathToFile= r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\VISSIM - V2'
ExistingPMfi = r'Existing\20834_Existing_PM--C1C2C3C4C5C6C7C8_Vehicle Travel Time Results.att'
ExistingPMfi = os.path.join(PathToFile,ExistingPMfi)
ExistingAMfi =r'Existing\20834_Existing_AM--C1C2aC3C4C5C6C7C8C9_Vehicle Travel Time Results.att'
ExistingAMfi = os.path.join(PathToFile,ExistingAMfi)

NoBuildAMfi = r'No Build\20834_NoBuild_AM--C1C2C3C4C5C6C7C8C9_Vehicle Travel Time Results.att'
NoBuildPMfi = r'No Build\20834_NoBuild_PM---C1C2C3C4C5C6C7C8_Vehicle Travel Time Results.att'
BuildDCDI_AMfi = r'Build\DCDI\20834_Build1_6LnMod_AM---C1C2C3C4C5C6C7C8C9C10_Vehicle Travel Time Results.att'
BuildDCDI_PMfi = r'Build\DCDI\20834_Build1_6LnMod_PM---C1C2C3C4C5C6C7C8C9_Vehicle Travel Time Results.att'
BuildDCMI_AMfi = r'Build\DCMI\20834_DCMI_AM_20180808--C1C2C3C4C5C6C7C8_Vehicle Travel Time Results.att'
BuildDCMI_PMfi = r'Build\DCMI\20834_DCMI_PM_20180808--C1C2C3C4C5C6C7C8_Vehicle Travel Time Results.att'

NoBuildAMfi = os.path.join(PathToFile,NoBuildAMfi)
NoBuildPMfi = os.path.join(PathToFile,NoBuildPMfi)
BuildDCDI_AMfi = os.path.join(PathToFile,BuildDCDI_AMfi)
BuildDCDI_PMfi = os.path.join(PathToFile,BuildDCDI_PMfi)
BuildDCMI_AMfi = os.path.join(PathToFile,BuildDCMI_AMfi)
BuildDCMI_PMfi = os.path.join(PathToFile,BuildDCMI_PMfi)

#*******************************************************************************
#*******************************************************************************
# Which File are we analysising. Only make changes here and done:
FileAM = ExistingAMfi
FileAM  = NoBuildAMfi
FileAM  = BuildDCDI_AMfi
FileAM  = BuildDCMI_AMfi
#File Being Processed
#*******************************************************************************
FileAM = BuildDCMI_AMfi
#*******************************************************************************

def switch_PMFile(argument):
    switcher = {
        ExistingAMfi: ExistingPMfi,
        NoBuildAMfi: NoBuildPMfi,
        BuildDCDI_AMfi: BuildDCDI_PMfi,
        BuildDCMI_AMfi: BuildDCMI_PMfi
    }
    return switcher.get(argument, "Invalid Name")

FilePM = switch_PMFile(FileAM)

def switch_OutputTag(argument):
    switcher = {
        ExistingAMfi: "Existing",
        NoBuildAMfi: "No-Build",
        BuildDCDI_AMfi: "DCDI",
        BuildDCMI_AMfi: "DCMI"
    }
    return switcher.get(argument, "Invalid Name")

OutputFolderName = switch_OutputTag(FileAM)

if((FileAM == BuildDCDI_AMfi) | (FileAM == BuildDCMI_AMfi)):
    AreResultsForBuild = True
else: AreResultsForBuild = False
# Get the time keys to convert Vissim intervals to Hours of the day
PathToVis = r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files'
TimeConvFi  = 'TravelTimeKeyValuePairs.xlsx'
TimeConvFi = pd.ExcelFile(os.path.join(PathToVis,TimeConvFi))
TimeConvFi.sheet_names
TimeKeys = TimeConvFi.parse('TimeVissimKey')

if(AreResultsForBuild):
    TTSegLaneDat = TimeConvFi.parse('Build')
    TTSegLaneDat = TTSegLaneDat.dropna()
else:
    # Get The Travel time (TT) segment and # of lane data
    TTSegLaneDat = TimeConvFi.parse('Existing')
    TTSegLaneDat = TTSegLaneDat.dropna()

#Summarize the Travel time for AM and PM. Concatenate the 2 datasets.
SegIn1 = I95_Segs+ LPGASeg
TTAM = PreProcessVissimTT(FileAM,SegIn1)
TTPM = PreProcessVissimTT(FilePM,SegIn1)
DenAM = PreProcessVissimDensity(FileAM,TTSegLaneDat)
DenPM = PreProcessVissimDensity(FilePM,TTSegLaneDat)
TTAM = TTAM.merge(DenAM,left_on=['TIMEINT','TTSegNm'],right_on=['TIMEINT','SegName'],how='left')
TTPM = TTPM.merge(DenPM,left_on=['TIMEINT','TTSegNm'],right_on=['TIMEINT','SegName'],how='left')
TTAM.loc[:,'AnalysisPeriod'] = 'AM'
TTPM.loc[:,'AnalysisPeriod'] = 'PM'
Vissim_TT = pd.concat([TTAM,TTPM],sort=False)
Vissim_TT.columns
Vissim_TT = Vissim_TT[['AnalysisPeriod','TIMEINT','TTSegNm','VissimSMS','WeightedDensity','Veh','Len']]

#Correction-Oversaturated-Conditions
# For 10 to 20 mph if the density > 45, just use it else replace with 160
mask = (Vissim_TT.VissimSMS < 20) & (Vissim_TT.WeightedDensity < 45)
Vissim_TT.loc[mask,"WeightedDensity"] = 160 # 33 ft. spacing
# For < 10 mph replace all values with 160
mask = (Vissim_TT.VissimSMS < 10) & (Vissim_TT.WeightedDensity < 200)
Vissim_TT.loc[mask,"WeightedDensity"] = 160 # 33 ft. spacing
Debug = Vissim_TT.loc[mask,"WeightedDensity"]

# Convert Vissim intervals to Hours of the day
lkey_ = ['AnalysisPeriod','TIMEINT']
Vissim_TT = Vissim_TT.merge(TimeKeys,left_on = lkey_, right_on =lkey_,how = 'inner')

#Sort the Travel Segments in Correct Order
Vissim_TT.TTSegNm = pd.Categorical(Vissim_TT.TTSegNm,SegIn1)
Vissim_TT.Time = pd.Categorical(Vissim_TT.Time,TimeKeys.Time)
Vissim_TT.TIMEINT = pd.Categorical(Vissim_TT.TIMEINT,['900-1800','1800-2700','2700-3600','3600-4500',
                                        '4500-5400','5400-6300','6300-7200','7200-8100',
                                        '8100-9000','9000-9900','9900-10800','10800-11700'])

Vissim_TT.loc[:,'Dir'] = Vissim_TT.TTSegNm.apply(lambda x: x.split(' ')[0])
Vissim_TT.rename(columns={'WeightedDensity':'Density'},inplace=True)

Vissim_TT.sort_values(['AnalysisPeriod','TIMEINT','Dir','TTSegNm'],inplace=True)
Vissim_TT.set_index(['TIMEINT','Dir','TTSegNm','AnalysisPeriod'],inplace=True)
Vissim_TT = Vissim_TT.unstack()
Vissim_TT = Vissim_TT.swaplevel(0,1,axis=1)

mux = pd.MultiIndex.from_product([['AM','PM'],['Time','VissimSMS','Density','Veh','Len','FlowRate'],
        ], names=Vissim_TT.columns.names)
Vissim_TT = Vissim_TT.reindex(mux,axis=1)
idx = pd.IndexSlice
Vissim_I95 = Vissim_TT.loc[idx[:,:,I95_Segs], idx[:,['VissimSMS','Density']]]
Vissim_LPGA =Vissim_TT.loc[idx[:,:,LPGASeg], idx[:,['VissimSMS']]]

TimePer = ['AM','PM']
Direction = ['NB','SB','EB','WB']
PerfMeasure = ['VissimSMS','Density']
t = 'AM'
d = 'SB'
PerfMes ='Density'
(not((d=='EB' or d=='WB') & (PerfMes=='Density')))
if not os.path.exists(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName)):
    os.makedirs(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName))
os.chdir(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName))
os.getcwd()
DictFigs = {}
idx = pd.IndexSlice
#Color Palettes
#https://matplotlib.org/users/colormaps.html
for t in TimePer:
    for d in Direction:
        for PerfMes in PerfMeasure:
            if(not((d=='EB' or d=='WB') & (PerfMes=='Density'))):
                HeatMapDat = Vissim_TT.loc[idx[:,d,:], idx[t,['Time',PerfMes]]]
                HeatMapDat = HeatMapDat.droplevel(0,axis=1).reset_index()
                HeatMapDat = HeatMapDat.pivot('Time','TTSegNm',PerfMes)
                HeatMapDat = HeatMapDat.astype('int32')
                plt.figure()
                if PerfMes=='Density':
                    title_ = "Density (veh/mile/lane)"
                    colorBar_ = 'viridis_r'
                    vmax_ = 160; vmin_ = 0
                    g = sns.heatmap(HeatMapDat,cmap = colorBar_, linewidths=.5,vmin=vmin_, vmax=vmax_, annot=True,fmt="d")
                else:
                    title_ = "Space Mean Speed (mph)"
                    colorBar_ = 'viridis'
                    if(d=='SB') | (d=='NB'):
                        vmax_ = 70; vmin_ = 0
                    else:
                        vmax_ = 45; vmin_ = 0
                    g = sns.heatmap(HeatMapDat,cmap = colorBar_, linewidths=.5, vmin=vmin_, vmax=vmax_, annot=True,fmt="d")

                # figure size in inches
                #sns.set(rc={'figure.figsize':(6,10)})

                g.set_xticklabels(rotation=30,labels = g.get_xticklabels(),ha='right')
                g.set_title('{}—{}'.format(t,title_))
                g.set_xlabel("")
                g.set_ylabel("Time Interval")
                fig = g.get_figure()
                fig.savefig("".join([t,d,PerfMes,'.jpg']),bbox_inches="tight")

if not os.path.exists(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName)):
    os.makedirs(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName))
os.chdir(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName))
os.getcwd()
t = 'AM'
d = 'SB'
PerfMes ='VissimSMS'
for t in TimePer:
    for d in Direction:
        for PerfMes in PerfMeasure:
            os.chdir(r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files\Figs\{}'.format(OutputFolderName))
            if not os.path.exists('{}/{}/{}'.format(t,d,PerfMes)):
                os.makedirs('{}/{}/{}'.format(t,d,PerfMes))
            os.chdir('{}/{}/{}'.format(t,d,PerfMes))
            if(not((d=='EB' or d=='WB') & (PerfMes=='Density'))):
                LinePlotDat = Vissim_TT.loc[idx[:,d,:], idx[t,['Time',PerfMes]]]
                LinePlotDat = LinePlotDat.droplevel(0,axis=1).reset_index()
                for Tseg in LinePlotDat.TTSegNm.unique():
                    TempDat = LinePlotDat[LinePlotDat.TTSegNm==Tseg]
                    sns.set(rc={'figure.figsize':(6,4)})
                    sns.set_style("ticks")
                    fig,ax1 = plt.subplots()
                    sns.lineplot('Time',PerfMes,data = TempDat,ax=ax1)
                    if(PerfMes == 'VissimSMS'):
                        ax1.set_ylim(0,70)
                        ylab = "Space Mean Speed (mph)"
                    else:
                        ax1.set_ylim(0,170)
                        ylab = "Density (veh/mile/lane)"

                    ax1.set_title('{}—{}'.format(t,Tseg), y =1.05)
                    ax1.set_ylabel(ylab)
                    ax1.set_xlabel("Time Interval")
                    # Must draw the canvas to position the ticks
                    fig.canvas.draw()
                    ax1.set_xticklabels(rotation=30,labels = ax1.get_xticklabels(),ha='right')
                    fig.savefig("_".join([Tseg,d,PerfMes,'.jpg']),bbox_inches="tight")
                    plt.close()
#Write Results to Output File
PathToOutFi = r'C:\Users\abibeka\OneDrive - Kittelson & Associates, Inc\Documents\LPGA\VISSIM-Files'
OutFi = "AfterValidationRes.xlsx"
OutFi = os.path.join(PathToOutFi,OutFi)
if(OutputFolderName =='Existing'):
    writer=pd.ExcelWriter(OutFi)
else:
    writer=pd.ExcelWriter(OutFi, mode ='a')
Vissim_LPGA.to_excel(writer, '{}-LPGA-SMS-Res'.format(OutputFolderName),na_rep='NaN')
Vissim_I95.to_excel(writer,'{}-I95-SMS&Density-Res'.format(OutputFolderName),na_rep='NaN')
writer.save()

if(OutputFolderName == 'DCMI'):
    subprocess.Popen([OutFi],shell=True)


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Letter-color Consistency test
O.Colizoli 2020
Each letter of the alphabet in random order x 2
Color wheel opens at a randomized color on each trial (but does not turn)
Python 2..7
"""
# data saved in ~/LogFiles/sub-XXX

# Import necessary modules
import random
import numpy as np
import pandas as pd
import os, time  # for paths and data
from IPython import embed as shell
try:
    import Tkinter as tk # py27
    from tkColorChooser import askcolor
except:
    import tkinter as tk
    from tkinter.colorchooser import askcolor

# Get subject number via tkinter (command line doesn't work in PsychoPy)
subject_ID = []
session = []
## INPUT WINDOW
class GetInput():
    def __init__(self):
        self.root2 = tk.Tk()
        self.root2.title("Subject and Session")
        # always put in same location
        w = 400 # width for the Tk root
        h = 200 # height for the Tk root
        # get screen width and height
        ws = self.root2.winfo_screenwidth() # width of the screen
        hs = self.root2.winfo_screenheight() # height of the screen
        # calculate x and y coordinates for the Tk root window
        x = (ws/6) - (w/6)
        y = (hs/6) - (h/6)
        self.root2.geometry('%dx%d+%d+%d' % (w, h, x, y))
        # Subject
        self.e = tk.Entry(self.root2)
        self.e.insert(0, 'Subject Number')
        self.e.pack()
        self.e.focus_set()
        # Session
        self.e2 = tk.Entry(self.root2)
        self.e2.insert(0, 'Session')
        self.e2.pack()
        self.e2.focus_set()

        txt='If each letter of the alphabet\
        \nwere to have a unique color,\
        \nwhat color would it have?\
        \n\nThere are no right or wrong answers.'
        # instructions
        self.instr = tk.Label(self.root2, bg='white', text=txt, font=("Helvetica", 14))
        self.instr.pack()

        b = tk.Button(self.root2,text='OK',command=self.get_input)
        b.pack(side='bottom')

        self.root2.mainloop()

    def get_input(self):
        subj_str = self.e.get()
        sess_str = self.e2.get()
        subject_ID.append(subj_str)
        session.append(sess_str)
        self.root2.destroy()

## ASK INPUT
app = GetInput()   # subject and session
subject_ID = int(subject_ID[0])
session = int(session[0])

## Create LogFile folder cwd/LogFiles
cwd = os.getcwd()
logfile_dir = os.path.join(cwd,'LogFiles','sub-{}'.format(subject_ID),'sess-{}'.format(session),'behav')
if not os.path.isdir(logfile_dir):
    os.makedirs(logfile_dir)
timestr = time.strftime("%Y%m%d-%H%M%S")
output_alphabet = os.path.join(logfile_dir,'sub-{}_sess-{}_task-consistency_events_{}.tsv'.format(subject_ID,session,timestr))

### CONSISTENCY TASK ###
alphabet = ['a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z']
#alphabet = ['a','b','c']

REPS = 2 # number of times to repeat whole alphabet

RGBS = [] # save output
L = '2'  # place holder

class Test():
    def __init__(self):
        self.counter = 1
        self.root = tk.Tk()
        self.root.title("Subject {} Session {}".format(subject_ID, session))
        # always put in same location
        # get screen width and height
        ws = self.root.winfo_screenwidth() # width of the screen
        hs = self.root.winfo_screenheight() # height of the screen
        # open in full screen
        self.root.geometry('%dx%d+%d+%d' % (ws, hs, 0, 0))
        self.open1 = tk.Button(self.root, text='Pick a color:', command=self.pick_a_color, font=('Helvetica', '36'),padx=5, pady=5)
        self.open1.pack(fill=tk.X, expand=False)
        self.letter = tk.Label(self.root, bg='white', text=L, font=("Helvetica", 90))
        self.letter.pack()
        self.root.mainloop()

    def quit(self):
        RGBS.append( [L ,self.RGB, self.HEX, abc] )
        self.root.destroy()

    def pick_a_color(self,):
        # GET COLOR CHOOSER NOT OPEN ON TOP OF ROOT
        self.RGB,self.HEX = askcolor((random.randint(0,255), random.randint(0,255), random.randint(0,255)), parent=None, title='Pick a color: {}'.format(L) )
        self.letter.configure(fg = self.HEX)
        if self.counter:
            exit_button = tk.Button(self.root, text='FINISHED', command=self.quit, font=('Helvetica', '28'))
            exit_button.pack()
            self.counter = 0
        self.root.mainloop()

# MAIN LOOP
abc = 1 # round
for R in np.arange(REPS):
    random.shuffle(alphabet)
    # Open a new GUI per letter
    for L in alphabet:
        app = Test()
        # save colors on each trial to prevent losing data

        DFS = pd.DataFrame(RGBS)
        print(RGBS)

        try:
            DFS.columns = ["letter","rgb","hex","choice"]
            DFS['subject'] = np.repeat(subject_ID,len(DFS))
            DFS['r']            = [c[0] for c in DFS['rgb']]
            DFS['g']            = [c[1] for c in DFS['rgb']]
            DFS['b']            = [c[2] for c in DFS['rgb']]
        except:
            # clicked window away
            pass
        DFS.to_csv(output_alphabet, sep='\t') # save all alphabet/preferences for both groups (also in case it goes wrong)
    abc+=1

####################################
## SAVE OUTPUT & determine conditions
print(RGBS)
print('consistency test - success!')

##### OUTPUT FIGURE WITH COLORS #####
# Sort and show letters x 2 side by side
del tk  # py27
del askcolor
import matplotlib.pyplot as plt # doesn't work together with tkinter
import seaborn as sns
fig = plt.figure(figsize=(10,5))

# Sort so the same letters go side by side for each choice
try:
    DFS.sort_values(by=['choice', 'letter'],inplace=True)
except:
    DFS = DFS.sort(['choice', 'letter'])

DFS.reset_index(inplace=True)
for i,A in enumerate(alphabet):
    ax = fig.add_subplot(6,5,i+1)
    ax.text(0.5, 0.5, DFS['letter'][i], color=DFS['hex'][i],fontsize=18)
    ax.text(0.25, 0.5, DFS['letter'][i+len(alphabet)], color=DFS['hex'][i+len(alphabet)],fontsize=18)
    ax.set_axis_off()

sns.despine(offset=10, trim=True)
plt.tight_layout()
fig.savefig(os.path.join(cwd,'LogFiles','sub-{}'.format(subject_ID),'sess-{}'.format(session),'behav','sub-{}_sess-{}_colors.pdf'.format(subject_ID,session)))
print('success: sub-{}_sess-{}_colors.pdf'.format(subject_ID,session))


text:
#!/usr/bin/env python3

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import os
import argparse

import numpy as np
import torch as th
import torchvision
from tqdm import tqdm
from tensorboardX import SummaryWriter

import utils

def main(args):
    writer = SummaryWriter(os.path.join("logs", args.dir))
    trainloader, testloader = get_loaders(args.batch_size)
    epsilon = 0.5
    beta = 1.0
    alpha1 = 0.1
    alpha2 = 0.05

    a = np.sqrt(2.0 / (784 + 500))
    W1 = th.randn(784, 500, device=device) * a
    b1 = th.randn(500, device=device) * a

    a = np.sqrt(2.0 / (500 + 10))
    W2 = th.randn(500, 10, device=device) * a
    b2 = th.randn(10, device=device) * a

    states = [(th.rand(args.batch_size, 500, device=device), \
            th.rand(args.batch_size, 10, device=device)) for _ in range(len(trainloader))]
    states_test = [(th.rand(args.batch_size, 500, device=device), \
            th.rand(args.batch_size, 10, device=device)) for _ in range(len(testloader))]

    for epoch in range(args.epochs):
        running_loss = running_energy = running_true_positive = 0.

        for i, (x, labels) in enumerate(tqdm(trainloader)):
            x, labels = x.view(-1, 784).to(device), labels.to(device)
            h, y = states[i]

            # Free phase
            for j in range(20):
                dh = d_rho(h) * (x @ W1 + y @ W2.t() + b1) - h
                dy = d_rho(y) * (h @ W2 + b2) - y

                h = rho(h + epsilon * dh)
                y = rho(y + epsilon * dy)
                '''
                energy = (h.pow(2).sum() + y.pow(2).sum() \
                    - (W1 * (x.t() @ h)).sum() - (W2 * (h.t() @ y)).sum() \
                    - 2 * (h @ b1).sum() - 2 * (y @ b2).sum()).item() / 2
                print(energy, dh.norm().item())
                '''

            h_free, y_free = h, y
            states[i] = h_free, y_free

            t = th.zeros(x.shape[0], 10, device=device)
            t.scatter_(1, labels.unsqueeze(1), 1)

            # Weakly clamped
            for j in range(4):
                dy = d_rho(y) * (h @ W2 + b2) - y + beta * (t - y)
                dh = d_rho(h) * (x @ W1 + y @ W2.t() + b1) - h

                h = rho(h + epsilon * dh)
                y = rho(y + epsilon * dy)
                '''
                energy = (h.pow(2).sum() + y.pow(2).sum() \
                    - (W1 * (x.t() @ h)).sum() - (W2 * (h.t() @ y)).sum() \
                    - 2 * (h @ b1).sum() - 2 * (y @ b2).sum()).item() / 2
                print(energy, dh.norm().item())
                '''

            h_clamped = h
            y_clamped = y
            W1 += alpha1 / beta * (rho(x.t()) @ (rho(h_clamped) - rho(h_free))) / args.batch_size
            W2 += alpha2 / beta * (rho(h_clamped.t()) @ rho(y_clamped) - rho(h_free.t()) @ rho(y_free)) / args.batch_size
            b1 += alpha1 / beta * (rho(h_clamped) - rho(h_free)).mean(0)
            b2 += alpha2 / beta * (rho(y_clamped) - rho(y_free)).mean(0)

            running_energy += (h_free.pow(2).sum() / 2 + y_free.pow(2).sum() \
                - ((W1 * (x.t() @ h_free)) / 2).sum() - ((W2 * (h_free.t() @ y_free)) / 2).sum() \
                - (h_free @ b1).sum() - (y_free @ b2).sum()).item()
            running_loss += (t - y_free).pow(2).sum().item()
            running_true_positive += (y_free.argmax(1) == labels).sum().item()

        energy_train = running_energy / (len(trainloader) * args.batch_size)
        accuracy_train = running_true_positive / (len(trainloader) * args.batch_size)
        loss_train = running_loss / (len(trainloader) * args.batch_size)
        print(f"Energy: {energy_train}, Accuracy: {accuracy_train}, Loss: {loss_train}")
        writer.add_scalar(f"loss", loss_train, epoch)
        writer.add_scalar(f"energy", energy_train, epoch)
        writer.add_scalar(f"accuracy", accuracy_train, epoch)

        running_true_positive = 0.
        for i, (x, labels) in enumerate(tqdm(testloader)):
            x, labels = x.view(-1, 784).to(device), labels.to(device)
            h, y = states_test[i]

            # Free phase
            for j in range(20):
                dh = d_rho(h) * (x @ W1 + y @ W2.t() + b1) - h
                dy = d_rho(y) * (h @ W2 + b2) - y

                h = rho(h + epsilon * dh)
                y = rho(y + epsilon * dy)

            h_free, y_free = h, y
            states_test[i] = h_free, y_free

            running_true_positive += (y_free.argmax(1) == labels).sum().item()
        accuracy_test = running_true_positive / (len(testloader) * args.batch_size)
        writer.add_scalar(f"accuracy_test", accuracy_test, epoch)
        print(f"Energy: {energy_train}, Accuracy: {accuracy_train}, Test Accuracy: {accuracy_test}, Loss: {loss_train}")

def rho(x):
    return x.clamp(0., 1.)

def d_rho(x):
    return ((x >= 0.) * (x <= 1.)).float()

def get_loaders(batch_size, fashion=False):
    mnist = torchvision.datasets.MNIST
    if fashion:
        mnist = torchvision.datasets.FashionMNIST

    transform = torchvision.transforms.Compose(
        [torchvision.transforms.ToTensor(),])

    trainloader = th.utils.data.DataLoader(
        mnist(root="./data", train=True, download=True, transform=transform),
        batch_size=batch_size,
        shuffle=True,
        num_workers=2)
    testloader = th.utils.data.DataLoader(
        mnist(root="./data", train=False, download=True, transform=transform),
        batch_size=batch_size,
        shuffle=False,
        num_workers=2)
    return trainloader, testloader

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("--epochs", type=int, default=25)
    parser.add_argument("--batch-size", type=int, default=20)

    parser.add_argument("--fashion", action="store_true", default=False,
        help="use fashion mnist")
    parser.add_argument("--dir", default=utils.timestamp(),
        help="name of output log directory")
    parser.add_argument("--no-cuda", action="store_true", default=False)
    args = parser.parse_args()
    device = th.device("cpu" if (not th.cuda.is_available() or args.no_cuda) else "cuda")
    device = th.device("cpu")
    main(args)

max_stars_repo_name: jiangdaniel/dl-papers

max_issues_repo_name: jiangdaniel/ml-implementations

max_forks_repo_name: jiangdaniel/ml-implementations


text:
#!/usr/bin/env python3
"""
Generate GSE64913
"""

__author__ = "Enrico Maiorino"
__version__ = "0.1.0"
__license__ = "MIT"

import logging
import GEOparse
import argparse
import pandas as pd
from funcs import utils
from os.path import join
import numpy as np

#def append_postfix(filename,postfix):
#    return "{0}_{2}.{1}".format(*filename.rsplit('.', 1) + postfix)

def main(args):
    logging.basicConfig(level=logging.INFO,
                        format='%(module)s:%(levelname)s:%(asctime)s:%(message)s',
                            handlers=[logging.FileHandler("../logs/report.log")])
    logging.info(args)

    utils.create_dir_if_not_exist(args.out_expr_dir)
    utils.create_dir_if_not_exist(join(args.out_expr_dir,'raw'))
    utils.create_dir_if_not_exist(join(args.out_expr_dir,'processed'))

    gse = GEOparse.get_GEO(geo='GSE64913', destdir=join(args.out_expr_dir,'raw'))
    annotated = gse.pivot_and_annotate('VALUE', gse.gpls['GPL570'], 'ENTREZ_GENE_ID')
    annotated2 = annotated[~pd.isnull(annotated.ENTREZ_GENE_ID)]
    annotated2 = annotated2.loc[~annotated2.isnull().values.all(axis=1)]
    annotated2['ENTREZ_GENE_ID'] = annotated2.ENTREZ_GENE_ID.str.split('///').str[0].astype(int)
    annotated2 = annotated2.set_index('ENTREZ_GENE_ID')

    classes = {}
    classes['healthy_cae'] = ['diagnosis: Healthy','cell type: Central airway epithelium']
    classes['healthy_pae'] = ['diagnosis: Healthy', 'cell type: Peripheral airway epithelium']
    classes['asthma_cae'] = ['diagnosis: Severe Asthmatic', 'cell type: Central airway epithelium']
    classes['asthma_pae'] = ['diagnosis: Severe Asthmatic', 'cell type: Peripheral airway epithelium']

    logging.info(classes)
    gsms = {cls: [gsm for gsm in gse.gsms if gse.gsms[gsm].metadata['characteristics_ch1'][1] == classes[cls][0] and
                  gse.gsms[gsm].metadata['characteristics_ch1'][5] == classes[cls][1]] for cls in classes}

    logging.info(' '.join(['{} GSM:{}'.format(cls, len(gsms[cls])) for cls in classes]))
    utils.create_dir_if_not_exist(args.out_expr_dir)
    utils.write_expr(join(args.out_expr_dir, 'processed', 'expr.tsv'), annotated2)
    for cls in classes:
        utils.write_text(join(args.out_expr_dir, 'processed', '{}_gsms.txt'.format(cls)), gsms[cls])

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='Process GSE64913')
    parser.add_argument('out_expr_dir', type=str, help='Output directory for expression data file and GSM lists')
    args = parser.parse_args()
    main(args)

max_stars_repo_name: reemagit/flowcentrality

max_issues_repo_name: reemagit/flowcentrality

max_forks_repo_name: reemagit/flowcentrality


text:
"""
@brief      test log(time=2s)
"""
import unittest
import warnings
from logging import getLogger
import numpy
from pyquickhelper.pycode import ExtTestCase
from sklearn.neighbors import RadiusNeighborsRegressor
from sklearn.datasets import make_regression
from skl2onnx.algebra.onnx_ops import (  # pylint: disable=E0611
    OnnxMul, OnnxAdd)
from mlprodict.onnxrt import OnnxInference
from mlprodict.onnx_conv import to_onnx
from mlprodict.tools.asv_options_helper import (
    get_ir_version_from_onnx, get_opset_number_from_onnx)

class TestOnnxrtOnnxRuntimeRuntime(ExtTestCase):

    def setUp(self):
        logger = getLogger('skl2onnx')
        logger.disabled = True

    def test_onnxt_runtime_add(self):
        idi = numpy.identity(2, dtype=numpy.float32)
        onx = OnnxAdd('X', idi, output_names=['Y'],
                      op_version=get_opset_number_from_onnx())
        model_def = onx.to_onnx({'X': idi.astype(numpy.float32)})
        X = numpy.array([[1, 2], [3, 4]], dtype=numpy.float32)

        model_def.ir_version = get_ir_version_from_onnx()
        oinf = OnnxInference(model_def, runtime='onnxruntime1')
        got = oinf.run({'X': X})
        self.assertEqual(list(sorted(got)), ['Y'])
        self.assertEqualArray(idi + X, got['Y'], decimal=6)

        oinf = OnnxInference(model_def, runtime='onnxruntime2')
        got = oinf.run({'X': X})
        self.assertEqual(list(sorted(got)), ['Y'])
        self.assertEqualArray(idi + X, got['Y'], decimal=6)

        oinf = OnnxInference(model_def, runtime='onnxruntime1')
        got = oinf.run({'X': X}, intermediate=True)
        self.assertEqual(list(sorted(got)), ['Ad_Addcst', 'X', 'Y'])
        self.assertEqualArray(idi + X, got['Y'], decimal=6)

    def test_onnxt_runtime_add_raise(self):
        idi = numpy.identity(2)
        onx = OnnxAdd('X', idi, output_names=['Y'],
                      op_version=get_opset_number_from_onnx())
        model_def = onx.to_onnx({'X': idi.astype(numpy.float32)})
        self.assertRaise(lambda: OnnxInference(model_def, runtime='onnxruntime-1'),
                         ValueError)

    def test_onnxt_runtime_add1(self):
        idi = numpy.identity(2, dtype=numpy.float32)
        onx = OnnxAdd('X', idi, output_names=['Y'],
                      op_version=get_opset_number_from_onnx())
        model_def = onx.to_onnx({'X': idi.astype(numpy.float32)})
        X = numpy.array([[1, 2], [3, 4]], dtype=numpy.float32)
        model_def.ir_version = get_ir_version_from_onnx()
        oinf = OnnxInference(model_def, runtime='onnxruntime1')
        got = oinf.run({'X': X})
        self.assertEqual(list(sorted(got)), ['Y'])
        self.assertEqualArray(idi + X, got['Y'], decimal=6)

    def test_onnxruntime_bug(self):
        rnd = numpy.random.randn(2, 20, 20).astype(numpy.float32)
        bni = (numpy.random.random((20, 20)).astype(  # pylint: disable=E1101
            numpy.float32) >= 0.7).astype(numpy.float32)
        mul = rnd * bni
        isn = any(numpy.isnan(mul.ravel()))
        self.assertFalse(isn)

        node = OnnxMul('X', bni, output_names=['Y'],
                       op_version=get_opset_number_from_onnx())
        onx = node.to_onnx({'X': rnd})
        for rt in ['python', 'onnxruntime1']:
            with self.subTest(runtime=rt):
                oinf = OnnxInference(onx, runtime=rt)
                y = oinf.run({'X': rnd})['Y']
                self.assertEqualArray(mul, y)

    def test_onnxruntime_knn_radius(self):
        def _get_reg_data(self, n, n_features, n_targets, n_informative=10):
            X, y = make_regression(  # pylint: disable=W0632
                n, n_features=n_features, random_state=0,
                n_targets=n_targets, n_informative=n_informative)
            return X, y

        def _fit_model(model, n_targets=1, label_int=False,
                       n_informative=10):
            X, y = _get_reg_data(20, 4, n_targets, n_informative)
            if label_int:
                y = y.astype(numpy.int64)
            model.fit(X, y)
            return model, X

        model, X = _fit_model(RadiusNeighborsRegressor())
        model_onnx = to_onnx(
            model, X[:1].astype(numpy.float32),
            target_opset=get_opset_number_from_onnx(),
            options={id(model): {'optim': 'cdist'}})
        oinf = OnnxInference(model_onnx, runtime='onnxruntime1')
        X = X[:7]
        got = oinf.run({'X': X.astype(numpy.float32)})['variable']
        exp = model.predict(X.astype(numpy.float32))
        if any(numpy.isnan(got.ravel())):
            # The model is unexpectedly producing nan values
            # sometimes.
            res = oinf.run({'X': X.astype(numpy.float32)}, intermediate=True)
            rows = ['--EXP--', str(exp), '--GOT--', str(got),
                    '--EVERY-OUTPUT--']
            for k, v in res.items():
                rows.append('-%s-' % k)
                rows.append(str(v))
            if any(map(numpy.isnan, res["variable"].ravel())):
                raise AssertionError('\n'.join(rows))
            # onnxruntime and mlprodict do not return the same
            # output
            warnings.warn('\n'.join(rows))
            return
        self.assertEqualArray(exp, got, decimal=4)

if __name__ == "__main__":
    unittest.main()

max_stars_repo_name: xadupre/mlprodict

max_issues_repo_name: xadupre/mlprodict

max_forks_repo_name: xadupre/mlprodict


text:
import logging
import xarray
import numpy as np
from typing import Tuple
from pymap3d.haversine import anglesep
import pymap3d as pm
from pymap3d.vincenty import vdist
import histutils.findnearest as fnd

try:
    import scipy.ndimage as ndi
except ImportError:
    ndi = None

def getimgind(imgs: xarray.Dataset, lla: np.ndarray, az: np.ndarray, el: np.ndarray) -> np.ndarray:
    """find pixels in images according to lat,lon,alt or az,el spec"""
    if lla is not None:
        lla = np.atleast_2d(lla)
        assert lla.ndim == 2
        assert lla.shape[1] == 3
        N = lla.shape[0]
    elif az is not None and el is not None:
        az = np.atleast_1d(az)
        el = np.atleast_1d(el)
        assert az.ndim == 1 and el.ndim == 1
        N = az.size
    else:
        raise ValueError("must specify LLA or az/el")
    # %% work
    ind = np.empty((N, 2), dtype=int)

    if lla is not None:
        az, el, _ = pm.geodetic2aer(
            lla[:, 0], lla[:, 1], lla[:, 2] * 1000, imgs.lat, imgs.lon, imgs.alt_m
        )

    for i, (a, e) in enumerate(zip(az, el)):
        ind[i, :] = fnd.findClosestAzel(imgs.az, imgs.el, a, e)

    ind = np.unique(ind, axis=0)

    assert ind.ndim == 2
    assert ind.shape[1] == 2  # row, column

    return ind

def projected_coord(imgs: xarray.Dataset, ind: np.ndarray, lla: Tuple[float, float, float]):
    az = imgs.az[ind[:, 0], ind[:, 1]].values.squeeze()
    el = imgs.el[ind[:, 0], ind[:, 1]].values.squeeze()

    alt_m = lla[2] * 1000 if lla is not None else 100e3

    plat, plon, palt_m = pm.aer2geodetic(
        az, el, alt_m / np.sin(np.radians(el)), imgs.lat, imgs.lon, imgs.alt_m
    )

    return az, el, plat, plon, palt_m

def line2plane(cam: xarray.Dataset) -> xarray.Dataset:
    """
    previously, you find the row, col pixels along the line from the other camera to magnetic zenith.
    In this function, you find the indices of the plane passing through that line and both cameras.

    Inputs:
    -------
    row: row indices of line
    col: column indices of line
    nPlane: number of samples (pixels) to take for plane.
      If the plane is horizontal in the images, it would be approximately the number of x-pixels.
      If diagonal, it would be about sqrt(2)*Nx.
      If less, the image is sampled more sparsely, which is in effect recducing the resolution of the image.
      This parameter is not critical.
    """
    # %% linear regression
    polycoeff = np.polyfit(cam["cols"], cam["rows"], deg=3, full=False)
    # %% columns (x)  to cut from picture (have to pick either rows or cols, arbitrarily I use cols)
    # NOT range, NEED dtype= for arange api

    # %% rows (y) to cut from picture
    cam["cutrow"] = np.rint(np.polyval(polycoeff, cam["cutcol"])).astype(int)
    assert (cam["cutrow"] >= 0).all() and (
        cam["cutrow"] < cam.y.size
    ).all(), "impossible least squares fit for 1-D cut\n is your video orientation correct? are you outside the FOV?"

    # DONT DO THIS: cutrow.clip(0,self.supery,cutrow)
    # %% angle from magnetic zenith corresponding to those pixels
    anglesep_deg = anglesep(
        cam.Bel,
        cam.Baz,
        # NEED .values or you'll get 2-D instead of 1-D!
        cam.el.values[cam.cutrow, cam.cutcol],
        cam.az.values[cam.cutrow, cam.cutcol],
    )

    if cam.verbose:
        from maatplotlib.pyplot import figure

        ax = figure().gca()
        ax.plot(cam.cutcol, anglesep_deg, label="angle_sep from MZ")
        ax.plot(cam.cutcol, cam.el.values[cam.cutrow, cam.cutcol], label="elevation angle [deg]")
        ax.legend()

    assert anglesep_deg.ndim == 1
    if not np.isfinite(anglesep_deg).all():
        logging.error(f"did you pick areas outside the {cam.filename} FOV?")
    # %% assemble angular distances from magnetic zenith--these are used in the tomography algorithm
    cam["angle_deg"], cam["angleMagzenind"] = sky2beam(anglesep_deg, cam.cutcol.size)

    return cam

def sky2beam(anglesep_deg, Ncol: int):
    angle_deg = np.empty(Ncol, float)
    # whether minimum angle distance from MZ is slightly positive or slightly negative, this should be OK
    MagZenInd = anglesep_deg.argmin()

    angle_deg[MagZenInd:] = 90.0 + anglesep_deg[MagZenInd:]
    angle_deg[:MagZenInd] = 90.0 - anglesep_deg[:MagZenInd]
    # %% LSQ
    col = np.arange(Ncol, dtype=int)
    polycoeff = np.polyfit(col, angle_deg, deg=1, full=False)

    return np.polyval(polycoeff, col), MagZenInd

def mergefov(w0: xarray.Dataset, w1: xarray.Dataset, projalt: float = 110e3, method: str = None):
    """
    inputs:
    -------
    w0: wide FOV data, particularly az/el
    w1: other camera FOV data contained in w0
    projalt: projection altitude METERS
    ofn: plot filename
    fovrow: to vastly speedup intiial overlap exploration, pick row(s) of pixels to examing--it could take half an hour + otherwise.

    find the ECEF x,y,z, at 110km altitude for narrow camera outer pixel
    boundary, then find the closest pixels in the wide FOV to those points.
    Remember, it can be (much) faster to brute force this calculation than to use
    k-d tree.

    """
    if projalt < 1e3:
        logging.warning(
            f"this function expects meters, you picked projection altitude {projalt/1e3} km"
        )

    # %% print distance from wide camera to narrow camera (just for information)
    print(
        f"intercamera distance with {w0.site}:  {vdist(w0.lat,w0.lon, w1.lat,w1.lon)[0]/1e3:.1f} kilometers"
    )
    # %% ENU projection from cam0 to cam1
    e1, n1, u1 = pm.geodetic2enu(w1.lat, w1.lon, w1.alt_m, w0.lat, w0.lon, w0.alt_m)
    # %% find the ENU of narrow FOV pixels at 110km from narrow FOV
    w1 = pixelmask(w1, method)
    if method is not None and method.lower() == "mzslice":
        w0 = pixelmask(w0, method)
        azSlice0, elSlice0, rSlice0 = pm.ecef2aer(
            w1.x2mz, w1.y2mz, w1.z2mz, w0.lat, w0.lon, w0.alt_m
        )
        azSlice1, elSlice1, rSlice1 = pm.ecef2aer(
            w0.x2mz, w0.y2mz, w0.z2mz, w1.lat, w1.lon, w1.alt_m
        )
        # find image indices (mask) corresponding to slice az,el
        w0["rows"], w0["cols"] = fnd.findClosestAzel(
            w0["az"].where(w0["fovmask"]), w0["el"].where(w0["fovmask"]), azSlice0, elSlice0
        )
        w0.attrs["Brow"], w0.attrs["Bcol"] = fnd.findClosestAzel(
            w0["az"], w0["el"], w0.Baz, w0.Bel
        )

        w1["rows"], w1["cols"] = fnd.findClosestAzel(
            w1["az"].where(w1["fovmask"]), w1["el"].where(w1["fovmask"]), azSlice1, elSlice1
        )
        w1.attrs["Brow"], w1.attrs["Bcol"] = fnd.findClosestAzel(
            w1["az"], w1["el"], w1.Baz, w1.Bel
        )
    else:
        # csc(x) = 1/sin(x)
        slantrange = projalt / np.sin(
            np.radians(np.ma.masked_invalid(w1["el"].where(w1["fovmask"])))
        )
        assert (slantrange >= projalt).all(), "slantrange must be >= projection altitude"

        e0, n0, u0 = pm.aer2enu(w1["az"], w1["el"], slantrange)
        # %% find az,el to narrow FOV from ASI FOV
        az0, el0, _ = pm.enu2aer(e0 - e1, n0 - n1, u0 - u1)
        assert (el0 >= 0).all(), "FOVs may not overlap, negative elevation from cam0 to cam1"
        # %% nearest neighbor brute force
        w0["rows"], w0["cols"] = fnd.findClosestAzel(w0["az"], w0["el"], az0, el0)

    return w0, w1

def pixelmask(data: xarray.Dataset, method: str = None) -> xarray.Dataset:
    """
    Use list because image may not be square

    returns mask of chosen pixels
    """
    if method is None or method == "rect":  # outermost edge of image (trivial)
        mask = np.zeros(data["az"].shape, dtype=bool)
        mask[:, 0] = True
        mask[0, :] = True
        mask[:, -1] = True
        mask[-1, :] = True
    elif method == "perimeter":  # perimeter for arbitrary shapes  e.g all-sky cameras
        if ndi is None:
            raise ImportError("pip install scipy")

        mask = ndi.distance_transform_cdt(~np.isnan(data["az"]), "taxicab") == 1
    elif method.lower() == "mzslice":
        """
        Assuming the imagers are all-sky, we arbitrarily discard pixels of low elevation as distortion is high low to horizon.
        """
        MIN_EL = 5  # degrees, arbitrary
        if data.srpts is None:
            raise ValueError("must include slant range points")

        mask = np.zeros(data["az"].shape, dtype=bool)
        mask[data["el"] >= MIN_EL] = True

        data["fovmask"] = (("y", "x"), mask)

        data.attrs["x2mz"], data.attrs["y2mz"], data.attrs["z2mz"] = pm.aer2ecef(
            data.attrs["Baz"],
            data.attrs["Bel"],
            data.srpts,
            data.attrs["lat"],
            data.attrs["lon"],
            data.attrs["alt_m"],
        )
        return data
    else:
        raise ValueError(f"unknown mask {method}")
    # %% sanity check
    if ~mask.any():
        raise ValueError("no FOV overlap found")

    data["fovmask"] = (("y", "x"), mask)

    return data

max_stars_repo_name: scienceopen/themis-asi-reader

max_issues_repo_name: scienceopen/themis-asi-reader

max_forks_repo_name: scivision/themisasi

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
"""
Cheryl Ngo

11/24/19
"""
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Code Flow
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
'''
The true/false variables in the Variables Used section are to be manipulated
to determine what type of model should be trained. For each model, the data
loaded and split into training and test data. The image shape is saved for
use as parameters in the model. According to the true/false variables set,
the parameters number of epochs for which the model is trained is varied.
The data is reshaped to be represented as sparse, binary vectors.  After this,
the data may be augmented.

Depending on the true/false parameters specified, the convolutional model is
built with one or many convolutional and max pooling layers and possible
dropout followed by two dense layers with relu activation and a final dense
layer with 10 units and softmax activation.

The model is then trained for the specified number of epochs and the
results in terms of the loss function and accuracy are plotted on a graph
for both the training and test data.  The start and stop time for the model
training are recorded and the duration is printed along with the final test
loss and test accuracy.
'''
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Import Libraries Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
import matplotlib.pyplot as plt
from keras import models
from keras import layers
from keras.layers import Dropout
import datetime
from keras.utils import to_categorical

from keras.preprocessing.image import ImageDataGenerator
import numpy as np

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Variables Used
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Set as True if the model should include data augmentation
dataAugmentation=True
#Set as True if the model has at least three convolutional layers (no dropout)
threeConvLayer=True
#Set as True if the model has four convolutional layers (no dropout)
fourConvLayer=False
#Set as True if the model has dropout (three convolutional layers)
dropOut=False

'''
num_nodes (first and second) are used to set the number of hidden units in the
     two convolutional layers
num_epochs is used to determine how many epochs the model should be trained for
size_batch is used to set the number of samples used to train the network
     at a time
verbose_setting is used to generate output that allows the user to visualize
     the progress of the training
num_classes is used to signify that there are 10 classes in the data

train/test_images are used to store the normalized input vector
train/test_labels are used to store the model factors
train/test_number_images, train/test_x_image_size, and
     train/test_y_image_size are used in the convolutional layers to specify
     the shape of the inputs
train_datagen is used to perform data augmentation

acc, val_acc, loss, and val_loss are used to store the results of the model
     training for use in a plot to compare loss and metrics for the training and
     test data
'''
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Objective
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
'''The goal of this program is to create a ConvNet that accurately classifies
data trained from the MNIST fashion dataset.  The code is split into
sections and subtasks are commented out and classified in the applicable
section
'''
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Load Data Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Set a seed to allow for replication
np.random.seed(7)

from keras.datasets import fashion_mnist
#Split into test and training data
(train_images, train_labels), (test_images, test_labels) = fashion_mnist.load_data()

#Capture the time to allow us to calculate training duration
start_time = datetime.datetime.now()
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Parameters Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Set a variable to represent the first dimension--number of images
train_number_images = train_images.shape[0]
test_number_images = test_images.shape[0]

#Set a variable to represent the second dimension--pixel value length
train_x_image_size = train_images.shape[1]
train_y_image_size = train_images.shape[2]

#Set a variabelt o represent the third dimension--pixel value width
test_x_image_size = test_images.shape[1]
test_y_image_size = test_images.shape[2]

#Set the number of hidden units in the first dense layer
num_nodes_first = 448
#Set the number of hidden units in the second dense layer
num_nodes_second=112 #128

#Train the dropout model through 16 epochs
if dropOut:
    num_epochs=16
#Train the data augmentation model through 8 epochs
elif dataAugmentation:
    num_epochs=100
#Train the four convolutional layer model through 6 epochs
elif fourConvLayer:
    num_epochs=6
#Train the three convolutional layer model through 6 epochs
elif threeConvLayer:
    num_epochs=6 #Four conv layer trained through 15 epochs
#Train the one convolutional layer model through 4 epochs
else:
    num_epochs=4

#Set the number of samples processed before the model is updated as 30
size_batch = 30

#Return a visual progress indicator
verbose_setting = 1

#Set the number of classes as 10 since there are 10 types of clothing
num_classes = 10
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Pretreat Data Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Reshape the training images
train_images = train_images.reshape((train_number_images, train_x_image_size , train_y_image_size, 1))
#Convert training images to a float
train_images = (train_images.astype('float32'))

#Reshape the test images
test_images = test_images.reshape((test_number_images, test_x_image_size , test_y_image_size, 1))
#Convert test images to a float
test_images = (test_images.astype('float32'))

#Use one hot encoding to transform the integer into a catagorical variable (binary vector)
train_labels = to_categorical(train_labels)
test_labels = to_categorical(test_labels)

#Perform if the model is data augmentation
if dataAugmentation:
    #Create a data generator for the training data
    train_datagen = ImageDataGenerator(
        #Rescale the images
        rescale=1./255,
        #Include copies of the images rotated by up to 40 degrees
        rotation_range=25,
        #Shift the image horizontally by up to 25% of width
        width_shift_range=0.15,
        #Shift the image vertically by up to 25% of width
        height_shift_range=0.15,
        #Set the shear angle intensity
        shear_range=0.1,
        #Set the amount of zoom
        zoom_range=(.9,1.1),
        horizontal_flip=True,)
    # The validation data does not get augmented
    test_datagen = ImageDataGenerator(rescale=1./255)

    #Augment the training data
    train_generator = train_datagen.flow(
            train_images, train_labels,
            batch_size=32)
    #Do not augment the validation data
    validation_generator = test_datagen.flow(
            test_images, test_labels,
            batch_size=32)

else:
    #Rescale test and training images to [0,1]
    train_images = (train_images / 255)-.5
    test_images = (test_images/ 255)-.5

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Define Model Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Start the sequential model
network = models.Sequential()

#Add a convolutional layer of 32 hidden units with relu activation
network.add(layers.Conv2D(32, (3, 3), activation='relu', input_shape=(train_x_image_size , train_y_image_size, 1)))
#Add a max pooling layer to reduce dimensionality
network.add(layers.MaxPooling2D((2, 2)))
#Add dropout to reduce overfitting
if dropOut:
    network.add(Dropout(0.2))

if threeConvLayer:
    #Add a convolutional layer of 64 hidden units with relu activation
    network.add(layers.Conv2D(64, (3, 3), activation='relu', input_shape=(train_x_image_size , train_y_image_size, 1),padding='same'))
    #Add a max pooling layer to reduce dimensionality
    network.add(layers.MaxPooling2D((2, 2)))
#Add dropout to reduce overfitting
    if dropOut:
        network.add(Dropout(0.2))

    #Add a convolutional layer of 128 hidden units with relu activation
    network.add(layers.Conv2D(128, (3, 3), activation='relu', input_shape=(train_x_image_size , train_y_image_size, 1),padding='same'))
    #Add a max pooling layer to reduce dimensionality
    network.add(layers.MaxPooling2D((2, 2)))
#Add dropout to reduce overfitting
    if dropOut:
        network.add(Dropout(0.2))

if fourConvLayer:
    #Add a convolutional layer of 128 hidden units with relu activation
    network.add(layers.Conv2D(128, (3, 3), activation='relu', input_shape=(train_x_image_size , train_y_image_size, 1),padding='same'))
    #Add a max pooling layer to reduce dimensionality
    network.add(layers.MaxPooling2D((2, 2)))

#Reshape the tensor to one dimension
network.add(layers.Flatten())

#Add a dense layer of num_nodes_first hidden units with relu activation
network.add(layers.Dense(num_nodes_first, activation='relu'))
#Add dropout to reduce overfitting
if dropOut:
    network.add(Dropout(0.2))

#Add a dense layer of num_nodes_second hidden units with relu activation
network.add(layers.Dense(num_nodes_second, activation='relu'))
if dropOut:
    network.add(Dropout(0.2))

#Add the final hidden layer used to classify the images in one of num_classes classes
network.add(layers.Dense(num_classes, activation='softmax'))

#Set the compiler to have the adam optimizer with categorical crossentropy loss and accuracy as the metric
network.compile(optimizer='adam',
                loss='categorical_crossentropy',
                metrics=['accuracy'])

#View a detail summary of the model structure
network.summary()

"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Fit Model Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Use fit_generator to train the model if data augmentation was used
if dataAugmentation:
    history=network.fit_generator(train_generator, steps_per_epoch=len(train_images)/32, epochs=num_epochs, validation_data=validation_generator, validation_steps=len(test_images)/32)
#Use .fit to train the model
else:
    history = network.fit(train_images, train_labels, epochs=num_epochs, batch_size=size_batch, verbose = verbose_setting, validation_split=0.2)
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
Show output Section
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Save test loss and accuracy from the non-data augmentation results
if not dataAugmentation:
    test_loss, test_acc = network.evaluate(test_images, test_labels)
    #Print the test loss
    print('test_loss:', test_loss)
    #Print the test accuracy
    print('test_acc:', test_acc)

#Save the stop time to calculate the run time
stop_time = datetime.datetime.now()
#Print the run time
print ("Time required for training:",stop_time - start_time)

#Save the test accuracy of the model per epoch to plot
acc = history.history['acc']
#Save the validation accuracy of the model per epoch to plot
val_acc = history.history['val_acc']
#Save the test loss of the model per epoch to plot
loss = history.history['loss']
#Save the validation loss of the model per epoch to plot
val_loss = history.history['val_loss']
#Save the number of epochs for use in our metric plot
epochs = range(len(acc))

#Create a plot of training accuracy
plt.plot(epochs, acc, 'bo', label='Training acc')
#Add to the plot validation accuracy
plt.plot(epochs, val_acc, 'b', label='Validation acc')
#Set a title for the plot
plt.title('Training and validation accuracy')
#Set a y-axis label for the plot
plt.ylabel('accuracy')
#Set an x-axis label for the plot
plt.xlabel('epoch')
#Put a legend on the plot
plt.legend()

#Create figure
plt.figure()
#Create a plot of training loss
plt.plot(epochs, loss, 'bo', label='Training loss')
#Add to the plot validation loss
plt.plot(epochs, val_loss, 'b', label='Validation loss')
#Set a title for the plot
plt.title('Training and validation loss')
#Set a y-axis label for the plot
plt.ylabel('loss')
#Set an x-axis label for the plot
plt.xlabel('epoch')
#Put a legend on the plot
plt.legend()

#Show the plots
plt.show()


text:
# Purpose: get experiment ready

import os
import glob
import numpy as np
from subprocess import call
import time
import nilearn
from scipy import stats
import scipy.io as sio
import pickle
import nibabel as nib
import argparse
import random
import sys
from datetime import datetime
from dateutil import parser

# WHEN TESTING - COMMENT OUT
currPath = os.path.dirname(os.path.realpath(__file__))
rootPath = os.path.dirname(os.path.dirname(currPath))
sys.path.append(rootPath)
#WHEN TESTING: UNCOMMENT TO ADD PATH TO RT-CLOUD
#sys.path.append('/Users/amennen/github/rt-cloud/')
#sys.path.append('/jukebox/norman/amennen/github/brainiak/rt-cloud/')

import rtCommon.utils as utils
from rtCommon.fileClient import FileInterface
import rtCommon.projectUtils as projUtils
from rtCommon.structDict import StructDict
#from rtCommon.dicomNiftiHandler import getTransform
from rtCommon.imageHandling import getTransform
from projects.amygActivation.initialize import initialize
from projects.amygActivation.amygActivation import getRunFilename

# obtain the full path for the configuration toml file
defaultConfig = os.path.join(currPath, 'conf/amygActivation.toml')

def finalize(cfg, args):
	# first find the number of runs completed
	run_path = os.path.join(cfg.local.subject_full_day_path, 'run*')
	if args.filesremote:
		run_path = os.path.join(cfg.server.subject_full_day_path, 'run*')
	nRuns_completed = len(glob.glob(run_path))
	return nRuns_completed

def main(argv=None):
	"""
	This is the main function that is called when you run 'finialize.py'.

	Here, you will load the configuration settings specified in the toml configuration
	file, initiate the class fileInterface, and set up some directories and other
	important things through 'finalize()'
	"""

	# define the parameters that will be recognized later on to set up fileIterface
	argParser = argparse.ArgumentParser()
	argParser.add_argument('--config', '-c', default=defaultConfig, type=str,
						   help='experiment config file (.json or .toml)')
	# This parameter is used for projectInterface
	argParser.add_argument('--commpipe', '-q', default=None, type=str,
						   help='Named pipe to communicate with projectInterface')
	argParser.add_argument('--filesremote', '-x', default=False, action='store_true',
						   help='retrieve files from the remote server')
	argParser.add_argument('--addr', '-a', default='localhost', type=str,
			   help='server ip address')
	argParser.add_argument('--runs', '-r', default='', type=str,
					   help='Comma separated list of run numbers')
	argParser.add_argument('--scans', '-s', default='', type=str,
					   help='Comma separated list of scan number')
	args = argParser.parse_args(argv)

	# load the experiment configuration file
	cfg = utils.loadConfigFile(args.config)
	cfg = initialize(cfg, args)
	print(args.config)
	nRunsCompleted = finalize(cfg,args)
	# copy subject folders from server to local
	# subject-specific folder
	# everything in temp/convertedNiftis
	if args.filesremote:

		# open up the communication pipe using 'projectInterface'
		projectComm = projUtils.initProjectComm(args.commpipe, args.filesremote)

		# initiate the 'fileInterface' class, which will allow you to read and write
		#   files and many other things using functions found in 'fileClient.py'
		#   INPUT:
		#       [1] args.filesremote (to retrieve dicom files from the remote server)
		#       [2] projectComm (communication pipe that is set up above)
		fileInterface = FileInterface(filesremote=args.filesremote, commPipes=projectComm)

		# we don't need the tmp/convertedNiftis so first remove those
		tempNiftiDir = os.path.join(cfg.server.dataDir, 'tmp/convertedNiftis/')
		if os.path.exists(tempNiftiDir):
			projUtils.deleteFolder(tempNiftiDir)
			print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
			print('deleting temporary convertedNifti folder: ', tempNiftiDir)
			print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
		# next, go through each run and put each run data into local run folder
		for r in np.arange(nRunsCompleted):
			runNum = r + 1 # run numbers start at 1
			runId = 'run-{0:02d}'.format(runNum)
			runFolder = os.path.join(cfg.server.subject_full_day_path, runId, '*')
			listOfFiles = glob.glob(runFolder)
			runFolder_local = os.path.join(cfg.local.subject_full_day_path, runId)
			projUtils.downloadFilesFromList(fileInterface, listOfFiles, runFolder_local)
			print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
			print('downloading data to local computer: ', runFolder)
		# next delete the entire subject folder on the cloud
		# MAKE SURE THIS IS CORRECT FOR YOUR EXPERIMENT BEFORE YOU RUN
		subject_dir = os.path.join(cfg.server.dataDir, cfg.bids_id)
		print('FOLDER TO DELETE ON CLOUD SERVER: ', subject_dir)
		print('IF THIS IS CORRECT, GO BACK TO THE CONFIG FILE USED ON THE WEB SERBER COMPUTER AND CHANGE THE FLAG FROM false --> true IN [server] deleteAfter')
		if cfg.server.deleteAfter:
			print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
			print('DELETING SUBJECT FOLDER ON CLOUD SERVER: ', subject_dir)
			print('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
			if os.path.exists(subject_dir):
				projUtils.deleteFolder(subject_dir)

	return 0

if __name__ == "__main__":
	"""
	If 'finalize.py' is invoked as a program, then actually go through all of the
	portions of this script. This statement is not satisfied if functions are called
	from another script using "from finalize.py import FUNCTION"
	"""
	main()
	sys.exit(0)

max_stars_repo_name: Chibee/rt-cloud

max_issues_repo_name: Chibee/rt-cloud

max_forks_repo_name: Chibee/rt-cloud

reason: import numpy,import scipy,from scipy


text:
"""
Evaluation code for multimodal-ranking
"""
import numpy

from app.sandi.quotes.visual_semantic_embedding.datasets import load_dataset
from app.sandi.quotes.visual_semantic_embedding.tools import encode_sentences, encode_images

def evalrank(model, data, split='dev'):
    """
    Evaluate a trained model on either dev or test
    data options: f8k, f30k, coco
    """
    print 'Loading dataset'
    if split == 'dev':
        X = load_dataset(data, load_train=False)[1]
    else:
        X = load_dataset(data, load_train=False)[2]

    print 'Computing results...'
    ls = encode_sentences(model, X[0])
    lim = encode_images(model, X[1])

    (r1, r5, r10, medr) = i2t(lim, ls)
    print "Image to text: %.1f, %.1f, %.1f, %.1f" % (r1, r5, r10, medr)
    (r1i, r5i, r10i, medri) = t2i(lim, ls)
    print "Text to image: %.1f, %.1f, %.1f, %.1f" % (r1i, r5i, r10i, medri)

def i2t(images, captions, npts=None):
    """
    Images->Text (Image Annotation)
    Images: (5N, K) matrix of images
    Captions: (5N, K) matrix of captions
    """
    if npts == None:
        npts = images.shape[0] / 5
    index_list = []

    ranks = numpy.zeros(npts)
    for index in range(npts):

        # Get query image
        im = images[5 * index].reshape(1, images.shape[1])

        # Compute scores
        d = numpy.dot(im, captions.T).flatten()
        inds = numpy.argsort(d)[::-1]
        index_list.append(inds[0])

        # Score
        rank = 1e20
        for i in range(5*index, 5*index + 5, 1):
            tmp = numpy.where(inds == i)[0][0]
            if tmp < rank:
                rank = tmp
        ranks[index] = rank

    # Compute metrics
    r1 = 100.0 * len(numpy.where(ranks < 1)[0]) / len(ranks)
    r5 = 100.0 * len(numpy.where(ranks < 5)[0]) / len(ranks)
    r10 = 100.0 * len(numpy.where(ranks < 10)[0]) / len(ranks)
    medr = numpy.floor(numpy.median(ranks)) + 1
    return (r1, r5, r10, medr)

def t2i(images, captions, npts=None):
    """
    Text->Images (Image Search)
    Images: (5N, K) matrix of images
    Captions: (5N, K) matrix of captions
    """
    if npts == None:
        npts = images.shape[0] / 5
    ims = numpy.array([images[i] for i in range(0, len(images), 5)])

    ranks = numpy.zeros(5 * npts)
    for index in range(npts):

        # Get query captions
        queries = captions[5*index : 5*index + 5]

        # Compute scores
        d = numpy.dot(queries, ims.T)
        inds = numpy.zeros(d.shape)
        for i in range(len(inds)):
            inds[i] = numpy.argsort(d[i])[::-1]
            ranks[5 * index + i] = numpy.where(inds[i] == index)[0][0]

    # Compute metrics
    r1 = 100.0 * len(numpy.where(ranks < 1)[0]) / len(ranks)
    r5 = 100.0 * len(numpy.where(ranks < 5)[0]) / len(ranks)
    r10 = 100.0 * len(numpy.where(ranks < 10)[0]) / len(ranks)
    medr = numpy.floor(numpy.median(ranks)) + 1
    return (r1, r5, r10, medr)


text:
import torch
import numpy as np
import random
import torch.backends.cudnn as cudnn
import sys
import os

class Logger(object):
    def __init__(self,filename="result.txt"):
        self.terminal=sys.stdout
        self.log=open(filename,"a")
    def write(self,message):
        self.terminal.write(message)
        self.log.write(message)
    def flush(self):
        pass

class counter():
    def __init__(self):
        self.value = 0
        self.ave = 0
        self.sum = 0
        self.count = 0

    def update(self, val, n=1):
        self.value = val
        self.sum += val * n
        self.count += n
        self.ave = self.sum / self.count

def topk_accuracy(output, labels, k):
    batch_size = labels.size(0)
    _, pred = output.topk(k, 1, True, True)
    pred = pred.t()
    judge = pred.eq(labels.view(1, -1).expand_as(pred))
    acc_top_k = judge[:k].view(-1).float().sum(0)*100.0/batch_size
    return acc_top_k

def init_seeds(seed=0):
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    if seed == 0:
        cudnn.deterministic = True
        cudnn.benchmark = False

max_stars_repo_name: cititude/Media-and-Cognition-Homework

max_issues_repo_name: cititude/Media-and-Cognition-Homework

max_forks_repo_name: cititude/Media-and-Cognition-Homework


text:
#!/usr/local/sci/bin/python

#***************************************
# 22 March 2014 KMW - v1
# Plots station locations for each variable including removed sub/super sats if desired
#
#************************************************************************
#                                 START
#************************************************************************
# USE python2.7
# python2.7 PlotStationLocs_JUL2015.py
#
# REQUIRES
#
#************************************************************************
# Set up python imports
import matplotlib.pyplot as plt
import numpy as np
import sys, os
import scipy.stats
import struct
from mpl_toolkits.basemap import Basemap
import datetime as dt
from matplotlib.dates import date2num,num2date
from netCDF4 import Dataset

# RESTART VALUE
Restarter='------'				#'------'		#'681040'

Plottype='SATSSUBS' #'JUSTGOODS','SATSSUBS'

# Set up initial run choices
param='t'	# tw, q, e, rh, t, td, dpd
param2='T'	# Tw, q, e, RH, T, Td, DPD
homogtype='IDPHA'	# 'IDPHA','PHA','PHADPD'
nowmon='JUL'
nowyear='2015'
thenmon='JAN'
thenyear='2015'
version='2.0.1.2014p'
version2='2.0.1.2014p_SELECTmax1_most3'

# Set up directories and files
INDIR='/data/local/hadkw/HADCRUH2/UPDATE2014/LISTS_DOCS/'
OUTDIR='/data/local/hadkw/HADCRUH2/UPDATE2014/IMAGES/MAPS/'

#INLIST='Posthomog'+homogtype+param+'_goodsHadISDH.'+version+'_'+thenmon+thenyear+'.txt'
INLIST='HadISDH.land'+param2+'.2.0.1.2014p_SELECTmax1_most3_JUL2015.txt'
INSATS='Posthomog'+homogtype+param+'_satsHadISDH.'+version+'_'+thenmon+thenyear+'.txt'
INSUBS='Posthomog'+homogtype+param+'_subzerosHadISDH.'+version+'_'+thenmon+thenyear+'.txt'
OUTPLOT='StationLocations.land'+param2+'.'+version2+'_'+nowmon+nowyear

# Set up variables
ngoods=0	#set once file read in
nsats=0		#set once file read in
nsubs=0		#set once file read in

GoodWMOs=[]
GoodLats=[]
GoodLons=[]
SatWMOs=[]
SatLats=[]
SatLons=[]
SubWMOs=[]
SubLats=[]
SubLons=[]

#************************************************************************
# Subroutines
#************************************************************************
# READDATA
def ReadData(FileName,typee,delimee):
    ''' Use numpy genfromtxt reading to read in all rows from a complex array '''
    ''' Need to specify format as it is complex '''
    ''' outputs an array of tuples that in turn need to be subscripted by their names defaults f0...f8 '''
    return np.genfromtxt(FileName, dtype=typee,delimiter=delimee) # ReadData

#************************************************************************
# PlotNiceDotsMap
def PlotNiceDotsMap(TheFile,TheGLats,TheGLons,TheSTLats,TheSTLons,TheSBLats,TheSBLons,
                    TheGCount,TheSTCount,TheSBCount,TheLetter,TheNamee,ThePlot):
    ''' Plot all good lats and lons TheGLats, TheGLons '''
    ''' If there are sats (TheSTCount > 0) plot all sat lats TheSTLats, TheSTLons '''
    ''' If there are subzeros (TheSBCount > 0) plot all sat lats TheSBLats, TheSBLons '''
    ''' Label plot with totals '''
    ''' Save as png and eps '''

    # set up plot

    plt.clf()
    plt.figure(1,figsize=(8,3))
    plt.axes([0.05,0.05,0.9,0.9])

    # plot map without continents and coastlines
    m = Basemap(projection='kav7',lon_0=0)
    # draw map boundary, transparent
    m.drawmapboundary()
    m.drawcoastlines()
    # draw paralells and medians, no labels
    m.drawparallels(np.arange(-90,90.,30.))
    m.drawmeridians(np.arange(-180,180.,60.))

    # plot black dots for the goods
    x,y = m(TheGLons,TheGLats)
    m.scatter(x,y,s=7,marker='o',color="CornflowerBlue")

    # plot blue dots for the supersats
    if TheSTCount > 0:
        x,y = m(TheSTLons,TheSTLats)
        m.scatter(x,y,s=7,marker='o',color="Violet")

    # plot red dots for the subs
    if TheSBCount > 0:
        x,y = m(TheSBLons,TheSBLats)
        m.scatter(x,y,s=7,marker='o',color="Gold")

    # add labals and watermark
#    watermarkstring="/".join(os.getcwd().split('/')[4:])+'/'+os.path.basename( __file__ )+"   "+dt.datetime.strftime(dt.datetime.now(), "%d-%b-%Y %H:%M")
#    plt.figtext(0.01,0.01,watermarkstring,size=6)
    plt.figtext(0.05,0.9,TheLetter,size=16)
    plt.figtext(0.5,0.9,TheNamee,size=16,ha='center')

    print(ngoods,nsats,nsubs)
    totalgoods=TheGCount-(TheSTCount+TheSBCount)
    plt.annotate(str(totalgoods)+" good stations",xy=(0.52,0.265),xytext=None, xycoords='axes fraction',color="black")
    x,y = m(5,-40)
    m.scatter(x,y,s=15,marker='o',color='CornflowerBlue',)

    if (ThePlot == 'SATSSUBS'):
        plt.annotate(str(TheSTCount)+" supersaturated error stations",xy=(0.52,0.22),xytext=None, xycoords='axes fraction',color="black")
        x,y = m(5,-48)
        m.scatter(x,y,s=15,marker='o',color='Violet',)

        plt.annotate(str(TheSBCount)+" subzero error stations",xy=(0.52,0.175),xytext=None, xycoords='axes fraction',color="black")
        x,y = m(5,-56)
        m.scatter(x,y,s=15,marker='o',color='Gold',)

#    plt.show()
#    stop()

    plt.savefig(TheFile+".eps")
    plt.savefig(TheFile+".png")

    return #PlotNiceDotsMap

#************************************************************************
# MAIN PROGRAM
#************************************************************************
# read in station list
#MyTypes=("|S6","|S5","float","float","float","|S4","|S30","|S7")
MyTypes=("|S6","|S5","float","float","float","|S4","|S30")
#MyDelimiters=[6,5,9,10,7,4,30,7]
MyDelimiters=[6,5,9,10,7,4,30]
MyFile=INDIR+INLIST
RawData=ReadData(MyFile,MyTypes,MyDelimiters)
GoodWMOs=np.array(RawData['f0'])
GoodLats=np.array(RawData['f2'])
GoodLons=np.array(RawData['f3'])
ngoods=len(GoodWMOs)

# If Plottype='SATSSUBS' then look at sats and subs too
nsats=0
nsubs=0

if (Plottype == 'SATSSUBS'):
    # get longs and lats of goods, supersats and subzeros if they exist
    if param != 't':
        print("Looking for sats...")
        MyTypes=("|S6","|S5","|S21")
        MyDelimiters=[6,5,21]
        MyFile=INDIR+INSATS
        RawData=ReadData(MyFile,MyTypes,MyDelimiters)
        SatWMOs=np.array(RawData['f0'])
        totsats=len(SatWMOs)
        SatLats=[]
        SatLons=[]
        for id in range(totsats):
            for findee in range(ngoods):
	        if SatWMOs[id] == GoodWMOs[findee]:
	            SatLats.append(GoodLats[findee])
		    SatLons.append(GoodLons[findee])
		    nsats=nsats+1
		    break

    if param == 'q' or param == 'e':
        print("Looking for subs...")
        MyTypes=("|S6","|S5","|S21")
        MyDelimiters=[6,5,21]
        MyFile=INDIR+INSUBS
        RawData=ReadData(MyFile,MyTypes,MyDelimiters)
        SubWMOs=np.array(RawData['f0'])
        totsubs=len(SubWMOs)
        SubLats=[]
        SubLons=[]
        for id in range(totsubs):
            for findee in range(ngoods):
	        if SubWMOs[id] == GoodWMOs[findee]:
	            SubLats.append(GoodLats[findee])
		    SubLons.append(GoodLons[findee])
		    nsubs=nsubs+1
		    break

# pass to plotter

if param == 't':
    Letty='a)'
    Namey='temperature'
if param == 'tw':
    Letty='a)'
    Namey='wet bulb temperature'
if param == 'td':
    Letty='b)'
    Namey='dew point temperature'
if param == 'q':
    Letty='c)'
    Namey='specific humidity'
if param == 'e':
    Letty='b)'
    Namey='vapour pressure'
if param == 'rh':
    Letty='d)'
    Namey='relative humidity'
if param == 'dpd':
    Letty='c)'
    Namey='dew point depression'

MyFile=OUTDIR+OUTPLOT
PlotNiceDotsMap(MyFile,GoodLats,GoodLons,
                        SatLats,SatLons,
			SubLats,SubLons,
			ngoods,nsats,nsubs,
			Letty,Namey,Plottype)

#    stop()

print("And, we are done!")


text:
# from wordcloud import WordCloud
# import PIL
# import re
# import matplotlib.pyplot as plt
# import plotly.graph_objects as go
# import plotly.io as pio
# import numpy as np
# from wrapplotly.util import as_list
# wc = WordCloud(font_path = '/Library/Fonts/Microsoft/Arial.ttf',
#                prefer_horizontal = 0.2,
#                random_state=89860).generate_from_frequencies({'hi':5, 'lo': 1, 'medium':3})
# plt.imshow(wc)
# plt.show()
# dir(plt.imshow(wc))
# [w for w in dir(plt.imshow(wc).axes) if 'lim' in w]
# plt.imshow(wc).axes.get_xlim()
# plt.imshow(wc).axes.get_ylim()

# for item in wc.layout_:
#     print(item)

# unnested = [dict(text = t, freq = f, fontsize = fs, x = x, y = y, orientation = o, col = c)
#             for (t, f), fs, (y, x), o, c in wc.layout_]
# # note the swapping of x and y

# font_family = as_list(PIL.ImageFont.truetype(wc.font_path).getname())[0]
# #f = PIL.ImageFont.truetype(wc.font_path)
# #f.getname() #('Droid Sans Mono', 'Regular')

# fig = go.Figure()
# for word in unnested[:]:
#     if word['orientation'] is None:
#         textangle = 0
#     elif word['orientation'] == 2:
#         textangle = 270

#     #text = word['text']
#     #xanchor = 'left'
#     #yanchor = 'top' # might change if we flip the y-axis
#     #showarrow = False

#     font = dict(size = word['fontsize'] * 1,
#                 color = word['col'],
#                 family = font_family+', sans-serif')
#     hovertext = ('Word: {}
Relative Frequency: {}'
#                  .format(word['text'], word['freq']))
#     hoverlabel = dict(bgcolor = word['col'],
#                       bordercolor = word['col'],
#                       font_family = font_family+', sans-serif',
#                       font_color = wc.background_color
#                       )
#     fig = fig.add_annotation(text = word['text'],
#                        x = word['x'],
#                        y = word['y'],
#                        borderpad = wc.margin,
#                        textangle = textangle,
#                        font = font,
#                        hovertext = hovertext,
#                        hoverlabel = hoverlabel,
#                        xanchor = 'left',
#                        yanchor = 'top',
#                        showarrow = False)

# x_max = wc.width
# y_max = wc.height
# fig = fig.update_layout(yaxis_showgrid = False, xaxis_showgrid = False,
#                         yaxis_zeroline = False, xaxis_zeroline = False,
#                         yaxis_range = [wc.height,0], yaxis_constrain = 'domain',
#                         yaxis_scaleanchor = 'x', yaxis_scaleratio = 1,
#                         xaxis_range = [0,wc.width], xaxis_constrain = 'domain',
#                         plot_bgcolor = wc.background_color,
#                         #paper_bgcolor = wc.background_color
#                         margin={"l": 0, "r": 0, "t": 0, "b": 0},
#                         width = wc.width, height = wc.height,
#                         xaxis_showticklabels = False,
#                         yaxis_showticklabels = False
#                         )
# #fig.show()
# fig.write_image('eg.png', scale = 1)

# img_obj = PIL.Image.open('eg.png')
# img_obj.show()
# fig = go.Figure().update_layout(yaxis_range = [wc.height,0], yaxis_constrain = 'domain',
#                         yaxis_scaleanchor = 'x', yaxis_scaleratio = 1,
#                         xaxis_range = [0,wc.width], xaxis_constrain = 'domain',)
# fig = fig.add_layout_image(
#         x=0,
#         sizex=img_obj.width,
#         y=0,
#         sizey=img_obj.height,
#         xref="x",
#         yref="y",
#         opacity=1.0,
#         layer="below",
#         source=wc.to_image()
# )

# fig.show()

# go.Figure(go.Image(source = img_obj)).show()

# [attr for attr in dir(wc) if re.search('_$', attr) is not None]
# wc.layout_
# wc.words_

max_stars_repo_name: jill-augustine/wrapplotly

max_issues_repo_name: jill-augustine/wrapplotly

max_forks_repo_name: jill-augustine/wrapplotly