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

# In[1]:

# Essentials
import os, sys, glob
import pandas as pd
import numpy as np
import nibabel as nib
import scipy.io as sio
from tqdm import tqdm

# Stats
import scipy as sp
from scipy import stats
import statsmodels.api as sm
import pingouin as pg

# Plotting
import seaborn as sns
import matplotlib.pyplot as plt
plt.rcParams['svg.fonttype'] = 'none'

# In[2]:

from networkx import from_numpy_matrix, degree_centrality, closeness_centrality, betweenness_centrality, subgraph_centrality
# from sklearn.kernel_ridge import KernelRidge

# In[3]:

sys.path.append('/Users/lindenmp/Google-Drive-Penn/work/research_projects/neurodev_cs_predictive/1_code/')
from func import set_proj_env, my_get_cmap, get_fdr_p, get_fdr_p_df, node_strength, ave_control

# In[4]:

parc_str = 'schaefer' # 'schaefer' 'lausanne' 'glasser'
parc_scale = 200 # 200/400 | 125/250 | 360
edge_weight = 'streamlineCount' # 'streamlineCount' 'volNormStreamline'
parcel_names, parcel_loc, drop_parcels, num_parcels = set_proj_env(parc_str = parc_str, parc_scale = parc_scale, edge_weight = edge_weight)

# In[5]:

# output file prefix
run_hemi = ''
if run_hemi == 'ipsi':
    outfile_prefix = parc_str+'_'+str(parc_scale)+'_'+edge_weight+'_ipsi_'
elif run_hemi == 'contra':
    outfile_prefix = parc_str+'_'+str(parc_scale)+'_'+edge_weight+'_contra_'
else:
    outfile_prefix = parc_str+'_'+str(parc_scale)+'_'+edge_weight+'_'
outfile_prefix

# In[6]:

# we want to calculate conn features including subcortex
# drop brainstem but retain subcortex.
# here because conn features have already been calculate, we actually drop the brainstem from parcel_loc
if parc_str == 'lausanne':
    parcel_loc = parcel_loc[parcel_loc != 2]
    num_parcels = len(parcel_loc)
print(num_parcels)

# ### Setup directory variables

# In[7]:

print(os.environ['PIPELINEDIR'])
if not os.path.exists(os.environ['PIPELINEDIR']): os.makedirs(os.environ['PIPELINEDIR'])

# In[8]:

figdir = os.path.join(os.environ['OUTPUTDIR'], 'figs')
print(figdir)
if not os.path.exists(figdir): os.makedirs(figdir)

# In[9]:

phenos = ['Overall_Psychopathology','Psychosis_Positive','Psychosis_NegativeDisorg']
phenos_label = ['Overall Psychopathology','Psychosis (Positive)','Psychosis (Negative)']
phenos_short = ['Ov. Psy.','Psy. (pos)','Psy. (neg)']

print(phenos)

# metrics = ['str', 'ac']
metrics = ['str', 'ac', 'bc', 'cc', 'sgc']
algs = ['rr', 'krr_rbf']
scores = ['corr', 'rmse']

# ## Setup plots

# In[10]:

if not os.path.exists(figdir): os.makedirs(figdir)
os.chdir(figdir)
sns.set(style='white', context = 'paper', font_scale = 1)
sns.set_style({'font.family':'sans-serif', 'font.sans-serif':['Public Sans']})
cmap = my_get_cmap('pair')

# ## Load data

# In[11]:

df = pd.read_csv(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'store', outfile_prefix+'df.csv'))
df.set_index(['bblid', 'scanid'], inplace = True)

df_node = pd.read_csv(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'out', outfile_prefix+'X.csv'))
df_node.set_index(['bblid', 'scanid'], inplace = True)

df_node_ac_overc = pd.read_csv(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'out', outfile_prefix+'X_ac_c.csv'))
df_node_ac_overc.set_index(['bblid', 'scanid'], inplace = True)

df_node_ac_i2 = pd.read_csv(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'out', outfile_prefix+'X_ac_i2.csv'))
df_node_ac_i2.set_index(['bblid', 'scanid'], inplace = True)

df_pheno = pd.read_csv(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'out', outfile_prefix+'y.csv'))
df_pheno.set_index(['bblid', 'scanid'], inplace = True)

c = pd.read_csv(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'out', outfile_prefix+'c.csv'))
c.set_index(['bblid', 'scanid'], inplace = True); print(c.shape)
c.columns

# In[12]:

print(np.all(df_node.filter(regex = 'ac') == df_node_ac_i2))
print(np.sum(df_node.filter(regex = 'ac').corrwith(df_node_ac_i2, method='spearman') < 0.99))

# In[13]:

f, ax = plt.subplots()
f.set_figwidth(5)
f.set_figheight(5)

df_x = df_node.filter(regex = 'str'); df_x.columns = np.arange(df_x.shape[1])
df_y = df_node.filter(regex = 'ac'); df_y.columns = np.arange(df_y.shape[1])
xy_corr = df_x.corrwith(df_y, method = 'pearson')
print(np.round(xy_corr.mean(),2))
sns.distplot(xy_corr, ax=ax, label='str_ac')

df_x = df_node.filter(regex = 'str'); df_x.columns = np.arange(df_x.shape[1])
df_y = df_node.filter(regex = 'mc'); df_y.columns = np.arange(df_y.shape[1])
xy_corr = df_x.corrwith(df_y, method = 'pearson')
print(np.round(xy_corr.mean(),2))
sns.distplot(xy_corr, ax=ax, label='str_mc')

df_x = df_node.filter(regex = 'ac'); df_x.columns = np.arange(df_x.shape[1])
df_y = df_node.filter(regex = 'mc'); df_y.columns = np.arange(df_y.shape[1])
xy_corr = df_x.corrwith(df_y, method = 'pearson')
print(np.round(xy_corr.mean(),2))
sns.distplot(xy_corr, ax=ax, label='ac_mc')

f.legend()

# In[14]:

# # replace ac with initial condition = 2
# ac_labels = ['ac_' + str(i) for i in range(num_parcels)]
# df_node.loc[:,ac_labels]

# ### Nuisance regression

# In[15]:

# Note, this nuisance regression made no difference.
# Leaving off for simplicity of reporting.

# nuis_reg = KernelRidge(kernel='rbf'); nuis_reg.fit(c, df_node)
# X_pred = nuis_reg.predict(c); df_node = df_node - X_pred

# nuis_reg = KernelRidge(kernel='rbf'); nuis_reg.fit(c, df_node_ac_overc)
# X_pred = nuis_reg.predict(c); df_node_ac_overc = df_node_ac_overc - X_pred

# In[16]:

g = 0
gradient = np.loadtxt(os.path.join(os.environ['PIPELINEDIR'], '2_compute_gradient', 'out', parc_str+'_'+str(parc_scale)+'_'+'pnc_grads_template.txt'))[:,g]

# In[17]:

if parc_str == 'lausanne' and len(gradient) == num_parcels:
    print('cutting subcortex')
    gradient = gradient[parcel_loc==1]

# In[18]:

# load spin test outputs
if parc_scale != 250:
    spin_indices = pd.read_csv(os.path.join(os.environ['PROJDIR'],'rotated_ind_'+parc_str+str(parc_scale)+'.csv'),sep=',' , header=None).values
    num_perms = spin_indices.shape[1]

# In[19]:

# if run_ipsi == True and parc_str == 'schaefer' and parc_scale == 200: # reting ipsilateral hemisphere (optional)
#     num_parcels = int(num_parcels/2)
#     gradient = gradient[:num_parcels]
#     spin_indices = spin_indices[:num_parcels,:]

# In[20]:

A = np.load(os.path.join(os.environ['PIPELINEDIR'], '1_compute_node_features', 'store', outfile_prefix+'A.npy'))
A = np.mean(A, axis = 2)
np.any(np.isnan(A))

# In[21]:

# density
np.count_nonzero(np.triu(A))/((A.shape[0]**2-A.shape[0])/2)

# In[22]:

# threshold to mean density of cohort
print(df.loc[:,'network_density'].mean())
thresh = np.percentile(A,100-(df.loc[:,'network_density'].mean()*100))
A[A < thresh] = 0

# In[23]:

# new density
np.count_nonzero(np.triu(A))/((A.shape[0]**2-A.shape[0])/2)

# In[24]:

r_vals = pd.DataFrame(index = metrics, columns = phenos)
p_vals = pd.DataFrame(index = metrics, columns = phenos)

for metric in metrics:
    for pheno in phenos:
        X = df_node.filter(regex = metric).copy()
        if parc_str == 'lausanne':
            X = X.iloc[:,parcel_loc==1]
        y = df_pheno.loc[:,pheno].copy()

        X_y_corr = X.corrwith(y, method = 'spearman').values

        my_corr = sp.stats.spearmanr(gradient, X_y_corr)
        r_vals.loc[metric,pheno] = my_corr[0]
        p_vals.loc[metric,pheno] = my_corr[1]

p_vals = get_fdr_p_df(p_vals, rows = True)

# In[25]:

r_vals

# In[26]:

r_vals[p_vals < .05]

# In[27]:

my_ac = ave_control(A)
my_str = node_strength(A)
G = from_numpy_matrix(A)
my_bc = np.array(list(betweenness_centrality(G, normalized=False).values()))
my_cc = np.array(list(closeness_centrality(G).values()))
my_sgc = np.array(list(subgraph_centrality(G).values()))

my_ac = sp.stats.boxcox(my_ac)[0]
my_str = sp.stats.boxcox(my_str)[0]

# In[28]:

f = sns.jointplot(x = my_str, y = my_ac)
f.fig.set_figwidth(2.5)
f.fig.set_figheight(2.5)
# f.annotate(sp.stats.spearmanr)
f.plot_joint(sns.kdeplot, zorder = 0, n_levels = 20, alpha = 0.1)
f.plot_joint(plt.scatter, c = "k", s = 10, linewidth = 2, marker = ".", alpha = 0.5)
f.ax_joint.collections[0].set_alpha(0)
f.ax_joint.set_xlabel('Regional Strength')
f.ax_joint.set_ylabel('Regional Average Controllability')
f.ax_joint.tick_params(pad = -2)
f.savefig(outfile_prefix+'regional_str_vs_ac.png', dpi = 300, bbox_inches = 'tight')

# In[29]:

sp.stats.pearsonr(my_ac,my_str)

# In[30]:

sp.stats.spearmanr(my_ac,my_str)

# In[31]:

c_params = np.array([10, 100, 1000, 10000])
c_params

# In[32]:

ac_orig_r = np.zeros(c_params.shape,)
ac_str_r = np.zeros(c_params.shape,)

for i, c_param in enumerate(c_params):
    ac_tmp = ave_control(A, c = c_param)
    ac_tmp = sp.stats.boxcox(ac_tmp)[0]

    ac_orig_r[i] = sp.stats.spearmanr(my_ac,ac_tmp)[0]
    ac_str_r[i] = sp.stats.spearmanr(my_str,ac_tmp)[0]

# In[33]:

ac_orig_r[-1]

# In[34]:

ac_str_r[-1]

# In[35]:

f = sns.jointplot(x = c_params, y = ac_orig_r)
f.fig.set_figwidth(2.5)
f.fig.set_figheight(2.5)
f.plot_joint(plt.scatter, c = "k", s = 5, linewidth = 2, marker = ".", alpha = 0.3)
f.ax_joint.collections[0].set_alpha(0)
f.ax_joint.set_xlabel('C params')
f.ax_joint.set_ylabel('corr(ac_orig,ac)')
f.ax_joint.tick_params(pad = -2)
f.savefig(outfile_prefix+'regional_c_vs_corr(ac_orig,ac).png', dpi = 300, bbox_inches = 'tight')

# In[36]:

f = sns.jointplot(x = c_params, y = ac_str_r)
f.fig.set_figwidth(2.5)
f.fig.set_figheight(2.5)
f.plot_joint(plt.scatter, c = "k", s = 5, linewidth = 2, marker = ".", alpha = 0.3)
f.ax_joint.collections[0].set_alpha(0)
f.ax_joint.set_xlabel('C params')
f.ax_joint.set_ylabel('corr(str,ac)')
f.ax_joint.tick_params(pad = -2)
f.savefig(outfile_prefix+'regional_c_vs_corr(str,ac).png', dpi = 300, bbox_inches = 'tight')

# ## Correlation between metric gradient value

# In[37]:

# plot_me = my_cc.copy()

# f, ax = plt.subplots()
# f.set_figwidth(2)
# f.set_figheight(2)
# sns.regplot(x = gradient, y = plot_me, color = 'gray', scatter = False)
# plt.scatter(x = gradient, y = plot_me, marker = 'o', s = 10, c = gradient, cmap = 'viridis', alpha=0.75, linewidth = 1)
# ax.set_xlabel('Regional gradient value')
# # ax.set_ylabel(metric)
# # ax.set_ylim([0.25,0.9])
# # ax.set_xlim([-0.13,0.115])
# ax.tick_params(pad = -2)

# textstr = 'r = {:.2f}'.format(sp.stats.pearsonr(gradient,plot_me)[0])
# ax.text(0.95, 0.95, textstr, transform=ax.transAxes,
#         verticalalignment='top', horizontalalignment='right')

# # f.savefig(outfile_prefix+'regional_gradient_'+str(g)+'_vs_corr(str,ac).svg', dpi = 300, bbox_inches = 'tight')

# ## Correlation between (x,y) cross-subject similarity and gradient value

# In[38]:

# num_perms = spin_indices.shape[1]
num_perms = 10000

# In[39]:

if parc_str == 'schaefer' and parc_scale == 200 and not run_hemi:
    metric_pairs = ['str_ac', 'str_bc', 'str_cc', 'str_sgc',
                    'ac_bc', 'ac_cc', 'ac_sgc',
                    'bc_cc', 'bc_sgc',
                    'cc_sgc',
                    'str_mc', 'ac_mc']
else:
#     metric_pairs = ['str_ac', 'bc_cc']
    metric_pairs = ['str_ac',]
metric_pairs

# In[40]:

# # age split

# # age_split = c['ageAtScan1'] < c['ageAtScan1'].mean()
# age_split = c['ageAtScan1'] > c['ageAtScan1'].mean()

# c = c.loc[age_split,:]
# c.drop(labels=['ageAtScan1'], axis=1, inplace=True)

# df_node = df_node.loc[age_split,:]
# df_pheno = df_pheno.loc[age_split,:]

# In[41]:

emp_corr = np.zeros(len(metric_pairs))
p_val = np.zeros(len(metric_pairs))

for m, metric_pair in enumerate(metric_pairs):
    print(metric_pair)
    metric_x = metric_pair.split('_')[0]
    metric_y = metric_pair.split('_')[1]

    df_x = df_node.filter(regex = metric_x)
    df_y = df_node.filter(regex = metric_y)
    if parc_str == 'lausanne':
        df_x = df_x.iloc[:,parcel_loc==1]
        df_y = df_y.iloc[:,parcel_loc==1]

    df_x.columns = np.arange(df_x.shape[1])
    df_y.columns = np.arange(df_y.shape[1])
    xy_corr = df_x.corrwith(df_y, method = 'pearson')

    emp_corr[m] = sp.stats.pearsonr(gradient,xy_corr)[0]

    if parc_scale != 250:
        spun_corr = np.zeros(num_perms)
        for i in tqdm(np.arange(num_perms)):
            spun_corr[i] = sp.stats.pearsonr(gradient,xy_corr[spin_indices[:,i]-1])[0]

        if emp_corr[m] < 0:
            p_val[m] = np.sum(spun_corr < emp_corr[m])/num_perms
        else:
            p_val[m] = np.sum(spun_corr > emp_corr[m])/num_perms
    else:
        p_val[m] = sp.stats.pearsonr(gradient,xy_corr)[1]

# p_val = get_fdr_p(p_val)

# In[42]:

sns.set(style='white', context = 'paper', font_scale = 0.8)
sns.set_style({'font.family':'sans-serif', 'font.sans-serif':['Public Sans']})

# In[43]:

for m, metric_pair in enumerate(metric_pairs):
    metric_x = metric_pair.split('_')[0]
    metric_y = metric_pair.split('_')[1]

    df_x = df_node.filter(regex = metric_x)
    df_y = df_node.filter(regex = metric_y)
    if parc_str == 'lausanne':
        df_x = df_x.iloc[:,parcel_loc==1]
        df_y = df_y.iloc[:,parcel_loc==1]

    df_x.columns = np.arange(df_x.shape[1])
    df_y.columns = np.arange(df_y.shape[1])
    xy_corr = df_x.corrwith(df_y, method = 'pearson')

    f, ax = plt.subplots()
#     f.set_figwidth(1.25)
#     f.set_figheight(1.25)
    f.set_figwidth(1)
    f.set_figheight(1)
    sns.regplot(x = gradient, y = xy_corr, color = 'gray', scatter = False)
    plt.scatter(x = gradient, y = xy_corr, marker = 'o', s = 10, c = gradient, cmap = 'viridis', alpha=0.75, linewidth = 1)
    ax.set_xlabel('Regional gradient value')
    ax.set_ylabel('Pearson('+metric_x+','+metric_y+')')
    # ax.set_ylim([0.25,0.9])
#     ax.set_xlim([-0.13,0.115])
    ax.set_xlim([gradient.min()+gradient.min()*.25,
                 gradient.max()+gradient.max()*.25])
    ax.tick_params(pad = -2)

    if p_val[m] < 0.05:
        textstr = 'r = {:.2f}\npFDR < 0.05'.format(emp_corr[m])
    else:
        textstr = 'r = {:.2f}\npFDR = {:.2f}'.format(emp_corr[m], p_val[m])
    ax.text(0.95, 0.95, textstr, transform=ax.transAxes,
            verticalalignment='top', horizontalalignment='right')

    f.savefig(outfile_prefix+'regional_gradient_'+str(g)+'_vs_corr('+metric_x+','+metric_y+').svg', dpi = 300, bbox_inches = 'tight')

# ## Correlation between (ac,str) cross-subject similarity and gradient value

# In[44]:

m=0
metric_pair = metric_pairs[m]
metric_x = metric_pair.split('_')[0]
metric_y = metric_pair.split('_')[1]

df_x = df_node.filter(regex = metric_x)
df_y = df_node.filter(regex = metric_y)
if parc_str == 'lausanne':
    df_x = df_x.iloc[:,parcel_loc==1]
    df_y = df_y.iloc[:,parcel_loc==1]

df_x.columns = np.arange(df_x.shape[1])
df_y.columns = np.arange(df_y.shape[1])
xy_corr = df_x.corrwith(df_y, method = 'pearson')

f, ax = plt.subplots()
f.set_figwidth(2)
f.set_figheight(2)
sns.regplot(x = gradient, y = xy_corr, color = 'gray', scatter = False)
plt.scatter(x = gradient, y = xy_corr, marker = 'o', s = 10, c = gradient, cmap = 'viridis', alpha=0.75, linewidth = 1)
ax.set_xlabel('Regional gradient value')
ax.set_ylabel('Pearson(str,ac)')
ax.set_ylim([0.25,0.9])
ax.set_xlim([-0.13,0.115])
ax.tick_params(pad = -2)

if p_val[m] < 0.05:
    textstr = 'r = {:.2f}\npFDR < 0.05'.format(emp_corr[m])
else:
    textstr = 'r = {:.2f}\npFDR = {:.2f}'.format(emp_corr[m], p_val[m])
ax.text(0.95, 0.95, textstr, transform=ax.transAxes,
        verticalalignment='top', horizontalalignment='right')

f.savefig(outfile_prefix+'figure3_regional_gradient_'+str(g)+'_vs_corr('+metric_x+','+metric_y+').svg', dpi = 300, bbox_inches = 'tight')

# ## Correlation between (ac_over_c,str) cross-subject similarity and gradient value

# In[45]:

c_params = np.array([10, 100, 1000, 10000])
c_params

# In[46]:

emp_corr = np.zeros(len(metric_pairs))
p_val = np.zeros(len(metric_pairs))

for m, c_param in enumerate(c_params):
    print(c_param)
    df_x = df_node.filter(regex = 'str')
    df_y = df_node_ac_overc.filter(regex = 'ac_c'+str(c_param)+'_')
    if parc_str == 'lausanne':
        df_x = df_x.iloc[:,parcel_loc==1]
        df_y = df_y.iloc[:,parcel_loc==1]

    df_x.columns = np.arange(df_x.shape[1])
    df_y.columns = np.arange(df_y.shape[1])
    xy_corr = df_x.corrwith(df_y, method = 'pearson')

    emp_corr[m] = sp.stats.pearsonr(gradient,xy_corr)[0]

    spun_corr = np.zeros(num_perms)
    for i in tqdm(np.arange(num_perms)):
        spun_corr[i] = sp.stats.pearsonr(gradient,xy_corr[spin_indices[:,i]-1])[0]

    if emp_corr[m] < 0:
        p_val[m] = np.sum(spun_corr < emp_corr[m])/num_perms
    else:
        p_val[m] = np.sum(spun_corr > emp_corr[m])/num_perms

p_val = get_fdr_p(p_val)

# In[47]:

for m, c_param in enumerate(c_params):
    print(c_param)
    df_x = df_node.filter(regex = 'str')
    df_y = df_node_ac_overc.filter(regex = 'ac_c'+str(c_param)+'_')
    if parc_str == 'lausanne':
        df_x = df_x.iloc[:,parcel_loc==1]
        df_y = df_y.iloc[:,parcel_loc==1]

    df_x.columns = np.arange(df_x.shape[1])
    df_y.columns = np.arange(df_y.shape[1])
    xy_corr = df_x.corrwith(df_y, method = 'pearson')

    f, ax = plt.subplots()
#     f.set_figwidth(1.25)
#     f.set_figheight(1.25)
    f.set_figwidth(1)
    f.set_figheight(1)
    sns.regplot(x = gradient, y = xy_corr, color = 'gray', scatter = False)
    plt.scatter(x = gradient, y = xy_corr, marker = 'o', s = 10, c = gradient, cmap = 'viridis', alpha=0.75, linewidth = 1)
    ax.set_xlabel('Regional gradient value')
#     ax.set_xlabel('')
    ax.set_ylabel('Pearson(str,ac)')
    ax.set_ylim([0.25,1])
    ax.set_xlim([-0.12,0.12])
    ax.set_xlim([gradient.min()+gradient.min()*.25,
                 gradient.max()+gradient.max()*.25])
    ax.tick_params(pad = -2)

    ax.set_title('c = ' + str(c_param))

    if p_val[m] < 0.05:
        textstr = 'r = {:.2f}\npFDR < 0.05'.format(emp_corr[m])
    else:
        textstr = 'r = {:.2f}\npFDR = {:.2f}'.format(emp_corr[m], p_val[m])
    ax.text(0.95, 0.95, textstr, transform=ax.transAxes,
            verticalalignment='top', horizontalalignment='right')

    f.savefig(outfile_prefix+'regional_gradient_'+str(g)+'_vs_corr(str,ac)_c_'+str(c_param)+'.svg', dpi = 300, bbox_inches = 'tight')

reason: import numpy,import scipy,from scipy,import statsmodels,from networkx


text:
# ******************************************************************************
#
#   File:   metal/views.py
#   Rev:    a-2
#   Date:   07/28/2017
#
#   Developed for the U.S. Government under contract(s):
#           HR001117C0099
#
# ******************************************************************************
#
#   View set definitions for models
#
# ******************************************************************************
#
#   Modification Log:
#
#   a-1:    07/14/2017  pcharasala
#           : Initial version
#   a-2:    07/28/2017  cstarkey
#           : Added risk areas, colors and map data
#
# ******************************************************************************

from django.http import HttpResponse
from django.contrib.auth.models import User
from django.http import Http404
from rest_framework.decorators import detail_route
from corsheaders import defaults

from .serializers import ScenarioSerializer
from .serializers import SiteSerializer
from .serializers import ResourceSerializer
from .serializers import AssetSerializer
from .serializers import AssetResourceSerializer
from .serializers import RouteSerializer
from .serializers import RouteSegmentSerializer
from .serializers import AssetRouteAssignmentSerializer
from .serializers import RiskTypeSerializer
from .serializers import RiskAreaSerializer
from .serializers import RiskAreaVertexSerializer
from .serializers import MapDataSerializer
from .serializers import TimeToFailureDistributionSerializer

from rest_framework import status
from rest_framework import viewsets
from rest_framework.response import Response
from rest_framework.views import APIView

from .models import Asset
from .models import AssetResource
from .models import AssetRouteAssignment
from .models import Perspective
from .models import Resource
from .models import Route
from .models import RouteSegment
from .models import Scenario
from .models import Site
from .models import RiskType
from .models import RiskArea
from .models import RiskAreaVertex
from .models import TimeToFailureDistribution
from .serializers import PerspectiveSerializer
from .serializers import UserSerializer

def index(request):

    return build_scenarios()

def build_scenarios():

    import os
    import json

    scenario_dir = os.path.expanduser("~/metal_scenarios")
    # Get the list of scenarios in the database
    scenario_list = Scenario.objects.all()

    # Return list of directories in the scenario directory
    scenario_files = os.listdir(scenario_dir)
    for file in scenario_files:
        if file.endswith(".json"):
            # This file is a potential scenario file, see if it has been previously loaded
            fullfilename = os.path.join(scenario_dir, file)

            file_modified_time = os.path.getmtime(fullfilename)
            found_in_database = False
            for scenario in scenario_list:
                if scenario.file_name == fullfilename:
                    found_in_database = True
                    if file_modified_time > scenario.date_modified:
                        # The file is more recent than the database entry
                        scenario.date_modified = file_modified_time
                        json_string = open(fullfilename, 'r').read()
                        scenario.json_file = json_string
                        data = json.loads(json_string)
                        scenario.name = data["name"]
                        scenario.file_name = fullfilename
                        scenario.save()
                        break

            if not found_in_database:
                json_string = open(fullfilename, 'r').read()
                data = json.loads(json_string)
                scenario = Scenario(name=data["name"],
                                    json_file=json_string,
                                    file_name=fullfilename,
                                    date_modified=file_modified_time)
                scenario.save()

    scenario_list = Scenario.objects.all()
    outstring = ""
    for scenario in scenario_list:
        outstring += scenario.name + ", " + scenario.file_name + ', ' + str(scenario.date_modified) + '\n'

    return HttpResponse(outstring)

def run_model(scenario_key):

    import json
    import json_to_model_system as gda_json
    import analyze_lp as gda_lp

    scenario = Scenario.objects.filter(name=scenario_key).first()
    # Delete all records associated with this scenario
    resources_to_delete = Resource.objects.filter(scenario=scenario)
    assets_to_delete = Asset.objects.filter(scenario=scenario)
    for resource in resources_to_delete:
        for asset in assets_to_delete:
            asset_resources_to_delete = AssetResource.objects.filter(asset=asset, resource=resource)
            for asset_resource in asset_resources_to_delete:
                asset_resource.delete()
        resource.delete()

    for asset in assets_to_delete:
        asset.delete()

    site_to_delete = Site.objects.filter(scenario=scenario)
    for site in site_to_delete:
        site.delete()

    routes_to_delete = Route.objects.filter(scenario=scenario)
    for route in routes_to_delete:
        route_segments_to_delete = RouteSegment.objects.filter(route=route)
        for route_segment in route_segments_to_delete:
            route_segment.delete()
        route.delete()

    asset_route_assignments_to_delete = AssetRouteAssignment.objects.filter(scenario=scenario)
    for asset_route_assignment in asset_route_assignments_to_delete:
        asset_route_assignment.delete()

    risktypes_to_delete = RiskType.objects.filter(scenario=scenario)
    for risktype in risktypes_to_delete:
        risktype.delete()

    riskareas_to_delete = RiskArea.objects.filter(scenario=scenario)
    for riskarea in riskareas_to_delete:
        riskareavertex_to_delete = RiskAreaVertex.objects.filter(riskarea=riskarea)
        for riskareavertex in riskareavertex_to_delete:
            riskareavertex.delete()
        riskarea.delete()

    ttfdist_to_delete = TimeToFailureDistribution.objects.filter(scenario=scenario)
    for ttfdist in ttfdist_to_delete:
        ttfdist.delete()

    # Break up the input file into it's components to be written for display
    scenario_pnode = json.loads(scenario.json_file)

    # Write scenario details to the database -- start by determining the resources that are played
    resources_found = dict()
    assets_raw = scenario_pnode['assets']
    for asset_dict in assets_raw:
        asset = Asset(name=asset_dict['name'],
                      speed=(asset_dict['attributes'])['max_speed'],
                      htmlcolor=(asset_dict['attributes']).get('htmlcolor','black'),
                      scenario=scenario)
        asset.save()
        resource_dict = (asset_dict['attributes'])['resources']
        for resource_name, resource_def in resource_dict.items():
            if resource_name not in resources_found:
                resource = Resource(name=resource_name,
                                    scenario=scenario)
                resource.save()
                resources_found[resource_name] = resource
            else:
                resource = resources_found[resource_name]

            transport_capacity = 0.0
            if isinstance(resource_def['capacity'], dict):
                consumption_capacity = (resource_def['capacity'])['bunker']
                transport_capacity = (resource_def['capacity'])['deliverable']
            else:
                consumption_capacity = resource_def['capacity']

            contested_consumption = 0.0
            uncontested_consumption = 0.0
            if 'consumption' in resource_def.keys():
                contested_consumption = (resource_def['consumption'])['contested']
                uncontested_consumption = (resource_def['consumption'])['uncontested']

            asset_resource = AssetResource(asset=asset,
                                           resource=resource,
                                           transport_capacity=transport_capacity,
                                           consumption_capacity=consumption_capacity,
                                           contested_consumption=contested_consumption,
                                           uncontested_consumption=uncontested_consumption)
            asset_resource.save()

    sites_raw = scenario_pnode['sites']
    for site_dict in sites_raw:
        # FIXME - Note the input file permits multiple assets per site but now we are only taking one
        asset = Asset.objects.filter(name=(site_dict['needed_asset_types'])[0], scenario=scenario)
        site = Site(name=site_dict['name'],
                    scenario=scenario,
                    asset=asset[0],
                    latitude=(site_dict['location'])[0],
                    longitude=(site_dict['location'])[1])
        site.save()

    risktype_dict = dict()
    risk_raw = scenario_pnode.get('risktypes', [])
    for risk_dict in risk_raw:
        risktype_name = risk_dict['name']
        risktype = RiskType(name=risktype_name,
                      htmlcolor=risk_dict.get('htmlcolor','white'),
                      scenario=scenario)
        risktype.save()
        risktype_dict[risktype_name] = risktype

    riskarea_raw = scenario_pnode.get('riskareas', [])
    for riskarea_dict in riskarea_raw:
        risktype_name = riskarea_dict['risktype']
        risktype = risktype_dict[risktype_name]
        riskarea = RiskArea(name=riskarea_dict['name'],
                      risktype=risktype,
                      scenario=scenario)
        riskarea.save()
        for vertex_array in riskarea_dict['region']:
            riskareavertex = RiskAreaVertex(riskarea=riskarea,
			    	latitude=vertex_array[0],
			    	longitude=vertex_array[1])
            riskareavertex.save()

    # budget_raw = scenario_pnode['budget']
    # composites_raw = scenario_pnode['composite_relations']
    # services_raw = scenario_pnode['service_relations']

    # Call GDA function to create model system
    model_system, model_route = gda_json.json_to_model_system(scenario.json_file)

    # Create the model route and route segment tables
    for route_dict in model_route.routes:
        route = Route(scenario=scenario,
                      name=route_dict['name'],
                      distance=route_dict['total_dist'])
        route.save()

        for segment_dict in route_dict['segments']:
            start_site = Site.objects.filter(name=segment_dict['points'][0])
            end_site = Site.objects.filter(name=segment_dict['points'][1])
            segment = RouteSegment(route=route,
                                   start_latitude=start_site[0].latitude,
                                   start_longitude=start_site[0].longitude,
                                   end_latitude=end_site[0].latitude,
                                   end_longitude=end_site[0].longitude,
                                   distance=segment_dict['dist'])
            segment.save()

    print("Beginning Run of LP for ", scenario.name)
    # Construct LP from model system and model routes
    gda_lp.assemble_problem(model_system, model_route)

    # Solve LP
    result = gda_lp.solve_LP_per_route(model_system, model_route)
    print("Completed Run of LP for ", scenario.name)

    composite_asset_list = model_system.assets
    comp_asset_based_results = result.lp_res_counts_df
    for df_index, row in comp_asset_based_results.iterrows():
        if row['count'] > 0:
            composite_asset_name = df_index[0]
            route_name = df_index[1]
            for comp_asset in composite_asset_list:
                if comp_asset.symbol_composite == composite_asset_name:
                    for asset_name in comp_asset.contains:
                        asset_list = Asset.objects.filter(scenario=scenario, name=asset_name)
                        route_list = Route.objects.filter(scenario=scenario, name=route_name)

                        asset_route_assign = AssetRouteAssignment(scenario=scenario,
                                                                  asset=asset_list[0],
                                                                  route=route_list[0],
                                                                  count=row['count'],
                                                                  utilization=row['utilization'])
                        asset_route_assign.save()
                    break

    print("Completed Asset Route Assignment writing for ", scenario.name)
    # FIXME - add a fake time to failure distribution until the model is updated
    data = "[{"
    import numpy as np
    s = np.random.normal(6, 1.5, 1000)
    counts, edges = np.histogram(s, bins=10)
    for day in range(1,11):
        data += "\"label\": \" " + str(day) + "\", \"value\":\" " + str((counts[day - 1]/1000))
        if day == 10:
            data += "\"}"
        else:
            data += "\"},{"
    data += "]"

    ttfdist = TimeToFailureDistribution(scenario=scenario,
                                        key=(scenario.name + "- Time to Failure Distribution"),
                                        x_axis_label='Days to Failure',
                                        y_axis_label='Probability',
                                        data=data)
    ttfdist.save()
    print("Completed building Time to failure distribution for ", scenario.name)
    print("Data = ", data)
    # print_database()

    #return HttpResponse(result)

def print_database():

    # Print Scenario Database
    scenario_list = Scenario.objects.all()
    for scenario in scenario_list:
        print("Scenario: ", scenario.name)
        print("\tFilename: ", scenario.file_name)
        print("\tModified Date: ", scenario.date_modified)

        resource_list = Resource.objects.filter(scenario=scenario)
        print("\tResources:")
        for resource in resource_list:
            print("\t\t", resource.name)

        site_list = Site.objects.filter(scenario=scenario)
        for site in site_list:
            print("\tSite: ", site.name)
            print("\t\tLocation: ", site.latitude, ", ", site.longitude)
            print("\t\tAsset:", site.asset.name)

        asset_list = Asset.objects.filter(scenario=scenario)
        for asset in asset_list:
            print("\tAsset: ", asset.name)
            print("\t\tSpeed: ", asset.speed)
            for resource in resource_list:
                asset_resource_list = AssetResource.objects.filter(asset=asset, resource=resource)
                print("\t\tResource:", resource.name)
                for asset_resource in asset_resource_list:
                    print("\t\t\tConsumption Capacity:", asset_resource.consumption_capacity)
                    print("\t\t\tTransport Capacity:", asset_resource.transport_capacity)
                    print("\t\t\tContested Consumption:", asset_resource.consumption_capacity)
                    print("\t\t\tUncontested Consumption:", asset_resource.uncontested_consumption)

        print("\tRoutes:")
        route_list = Route.objects.filter(scenario=scenario)
        for route in route_list:
            print("\t\tName:", route.name)
            print("\t\tDistance:", route.distance)
            print("\t\tSegments:")
            segment_list = RouteSegment.objects.filter(route=route)
            for segment in segment_list:
                print("\t\t\tStart Latitude:", segment.start_latitude)
                print("\t\t\tStart Longitude:", segment.start_longitude)
                print("\t\t\tEnd Latitude:", segment.end_latitude)
                print("\t\t\tEnd Longitude:", segment.end_longitude)

        asset_route_assign_list = AssetRouteAssignment.objects.filter(scenario=scenario)
        print("\tAsset Route Assignment")
        for asset_route_assign in asset_route_assign_list:
            print("\t\tAsset:", asset_route_assign.asset.name)
            print("\t\tRoute:", asset_route_assign.route.name)
            print("\t\tCount:", asset_route_assign.count)

        print("\tRiskAreas:")
        riskarea_list = RiskArea.objects.filter(scenario=scenario)
        for riskarea in riskarea_list:
            print("\t\tName:", riskarea.name)
            print("\t\tRisk Type:", riskarea.risktype.name)
            print("\t\tVerteces:")
            vertex_list = RiskAreaVertex.objects.filter(riskarea=riskarea)
            for vertex in vertex_list:
                print("\t\t\tLatitude:", vertex.latitude)
                print("\t\t\tLongitude:", vertex.longitude)

        dist_list = TimeToFailureDistribution.objects.filter(scenario=scenario)
        for dist in dist_list:
            print("\tDistribution:")
            print("\t\tKey:", dist.key)
            print("\t\tX-Axis Label:", dist.x_axis_label)
            print("\t\tY-Axis Label:", dist.y_axis_label)
            print("\t\tData:", dist.data)

class UserViewSet(viewsets.ViewSet):
    """
    A simple ViewSet for listing or retrieving users.
    """
    def list(self, request):
        queryset = User.objects.all()
        serializer = UserSerializer(queryset, many=True)
        return Response(serializer.data)

    def retrieve(self, request, pk=None):
        queryset = User.objects.all()
        user = get_object_or_404(queryset, pk=pk)
        serializer = UserSerializer(user)
        return Response(serializer.data)

class UserList(APIView):
    """
        List all users, or create a new user.
    """
    def get(self, request, format=None):
        users = User.objects.all()
        serializer = UserSerializer(users, many=True)
        return Response(serializer.data)

    def post(self, request, format=None):
        serializer = UserSerializer(data=request.DATA)
        if serializer.is_valid():
            serializer.save()
            return Response(serializer.data, status=status.HTTP_201_CREATED)
        return Response(serializer.errors, status=status.HTTP_400_BAD_REQUEST)

    def delete(self, request, pk, format=None):
        user = self.get_object(pk)
        user.delete()
        return Response(status=status.HTTP_204_NO_CONTENT)

class UserDetail(APIView):
    """
    Retrieve, update or delete a user instance.
    """
    def get_object(self, pk):
        try:
            return User.objects.get(pk=pk)
        except User.DoesNotExist:
            raise Http404

    def get(self, request, pk, format=None):
        user = self.get_object(pk)
        user = UserSerializer(user)
        return Response(user.data)

    def put(self, request, pk, format=None):
        user = self.get_object(pk)
        serializer = UserSerializer(user, data=request.DATA)
        if serializer.is_valid():
            serializer.save()
            return Response(serializer.data)
        return Response(serializer.errors, status=status.HTTP_400_BAD_REQUEST)

    def delete(self, request, pk, format=None):
        user = self.get_object(pk)
        user.delete()
        return Response(status=status.HTTP_204_NO_CONTENT)

class PerspectivesList(APIView):
    """
        List all perspectives, or create a new perspective.
    """
    def get(self, request, format=None):
        perspectives = Perspective.objects.all()
        serializer = PerspectiveSerializer(perspectives, many=True)
        return Response(serializer.data)

    def post(self, request, format=None):
        serializer = PerspectiveSerializer(data=request.DATA)
        if serializer.is_valid():
            serializer.save()
            return Response(serializer.data, status=status.HTTP_201_CREATED)
        return Response(serializer.errors, status=status.HTTP_400_BAD_REQUEST)

    def delete(self, request, pk, format=None):
        perspective = self.get_object(pk)
        perspective.delete()
        return Response(status=status.HTTP_204_NO_CONTENT)

class PerspectiveViewSet(viewsets.ModelViewSet):
    serializer_class = PerspectiveSerializer
    queryset = Perspective.objects.all()

class ScenarioViewSet(viewsets.ModelViewSet):
    build_scenarios()
    serializer_class = ScenarioSerializer
    queryset = Scenario.objects.all()

    @detail_route(methods=['get'], url_path='run_model')
    def run_model(self, request, pk=None):
        scenario = Scenario.objects.filter(name=pk).first()
        if(scenario != None):
            return run_model(scenario)
        response = HttpResponse({})
        response['Access-Control-Allow-Origin'] = 'http://localhost:8091'
        response['Access-Control-Allow-Methods'] = 'GET, POST, OPTIONS'
        return response

    @detail_route(methods=['get'], url_path='routes')
    def get_routes(self, request, pk=None):
        routes = Route.objects.filter(scenario=pk)
        serializer = RouteSerializer(routes, many=True)
        response = Response(serializer.data, status=status.HTTP_200_OK)
        response['Access-Control-Allow-Origin'] = 'http://localhost:8091'
        response['Access-Control-Allow-Methods'] = 'GET, POST, OPTIONS'
        return response

class ScenarioVewSet2(viewsets.ViewSet):
    """
    Viewset that calls index view too
    """
    def list(selfself, request):
        index(request)
        queryset = Scenario.objects.all()
        serializer = ScenarioSerializer(queryset, many=True)
        return Response(serializer.data)

class ResourceViewSet(viewsets.ModelViewSet):
    serializer_class = ResourceSerializer
    queryset = Resource.objects.all()

class SiteViewSet(viewsets.ModelViewSet):
    serializer_class = SiteSerializer
    queryset = Site.objects.all()

class AssetViewSet(viewsets.ModelViewSet):
    serializer_class = AssetSerializer
    queryset = Asset.objects.all()

class AssetResourceViewSet(viewsets.ModelViewSet):
    serializer_class = AssetResourceSerializer
    queryset = AssetResource.objects.all()

class RouteViewSet(viewsets.ModelViewSet):
    serializer_class = RouteSerializer
    queryset = Route.objects.all()

class RouteSegmentViewSet(viewsets.ModelViewSet):
    serializer_class = RouteSegmentSerializer
    queryset = RouteSegment.objects.all()

class AssetRouteAssignmentViewSet(viewsets.ModelViewSet):
    serializer_class = AssetRouteAssignmentSerializer
    queryset = AssetRouteAssignment.objects.all()

class RiskTypeViewSet(viewsets.ModelViewSet):
    serializer_class = RiskTypeSerializer
    queryset = RiskType.objects.all()

class RiskAreaViewSet(viewsets.ModelViewSet):
    serializer_class = RiskAreaSerializer
    queryset = RiskArea.objects.all()

class RiskAreaVertexViewSet(viewsets.ModelViewSet):
    serializer_class = RiskAreaVertexSerializer
    queryset = RiskAreaVertex.objects.all()

class MapViewSet(viewsets.ModelViewSet):
    serializer_class = MapDataSerializer
    queryset = Scenario.objects.all()

class TimeToFailureDistributionViewSet(viewsets.ModelViewSet):
    serializer_class = TimeToFailureDistributionSerializer
    queryset = TimeToFailureDistribution.objects.all()

    def get_object(self):
        if self.action == 'retrieve':
            scenario = self.kwargs['pk']
            run_model(scenario)
        return super().get_object()

ext: py

max_stars_repo_name: stephenst/base

max_issues_repo_name: stephenst/base

max_forks_repo_name: stephenst/base

reason: import numpy


text:
"""
Export colormaps from Python / matplotlib to JavaScript.
"""

# -----------------------------------------------------------------------------
# IMPORTS
# -----------------------------------------------------------------------------

import json

from matplotlib.colors import Colormap

import matplotlib.cm as cm
import matplotlib.colors as colors
import numpy as np

# -----------------------------------------------------------------------------
# MAIN CODE
# -----------------------------------------------------------------------------

if __name__ == "__main__":

    # Loop over all matplotlib colormaps and store them in a dictionary. This
    # dictionary contains the colors of the colormap as list of lists (= RGB
    # tuples), and whether or not the colormap should be interpolated (= false
    # for qualitative colormaps).
    colormaps = {}
    for name in dir(cm):

        # Skip reversed colormaps (they don't contain additional information)
        if name.endswith("_r"):
            continue

        # If `cmap` is a colormap, we can store the association information
        if isinstance(cmap := getattr(cm, name), Colormap):

            # Evaluate colormap on the grid to get colors, drop alpha channel
            # information, and round to a reasonable precision
            colors = np.around(cmap(np.linspace(0, 1, cmap.N))[:, 0:3], 4)

            # Store relevant colormap information
            colormaps[cmap.name] = {
                "interpolate": cmap.N >= 256,
                "colors": colors.tolist(),
            }

    # Save colormap data and shortcuts to data.js. The contents of this file
    # need to be copied manually to js-colormaps.js.
    with open("data.js", "w") as json_file:

        # Write the data dictionary to data.js
        json_file.write(f"const data = {json.dumps(colormaps)};")
        json_file.write('\n\n')

        # Write partial function applications to data.js so that we can use
        # a colormap by its name --- e.g., call viridis(0.5) to evaluate the
        # viridis colormap at a value of 0.5.
        for name in colormaps.keys():
            json_file.write(f"const {name} = partial('{name}');\n")
            json_file.write(f"const {name}_r = partial('{name}_r');\n")

    # Final words
    print("\nExported data to data.js, please copy to js-colormaps.js!\n")

max_stars_repo_name: bensnell/js-colormaps

max_issues_repo_name: bensnell/js-colormaps

max_forks_repo_name: bensnell/js-colormaps


text:
from tensorforce.environments import Environment
import numpy as np
import math
import signal
import time
import logging

eps = 1e-12

# Load UnityEnvironment and my wrapper
from mlagents.envs import UnityEnvironment

class UnityEnvWrapper(Environment):
    def __init__(self, game_name = None, no_graphics = True, seed = None, worker_id=0, size_global = 8, size_two = 5,
                 with_local = True, size_three = 3, with_stats=True, size_stats = 1,
                 with_previous=True, manual_input = False, config = None, curriculum = None, verbose = False,
                 agent_separate = False, agent_stats = 6, action_size = 19,
                 _max_episode_timesteps = 100,
                 with_class=False, with_hp = False, size_class = 3,
                 # Double agent
                 use_double_agent=False, double_agent = None, play_with_random=0.33,
                 # Reward model
                 reward_model = None):

        self.probabilities = []

        self.size_global = size_global
        self.size_two = size_two
        self.with_local = with_local
        self.size_three = size_three
        self.with_stats = with_stats
        self.size_stats = size_stats
        self.action_size = action_size

        self.manual_input = manual_input
        self.with_previous = with_previous
        self.config = config
        self.curriculum = curriculum
        self.verbose = verbose

        self.agent_separate = agent_separate
        self.agent_stats = agent_stats
        self.with_class = with_class
        self.with_hp = with_hp
        self.size_class = size_class
        self.use_double_agent = use_double_agent

        self.game_name = game_name

        self.no_graphics = no_graphics
        self.seed = seed

        self.worker_id = worker_id
        self.unity_env = self.open_unity_environment(game_name, no_graphics, seed, worker_id)
        self.default_brain = self.unity_env.brain_names[0]

        self.input_channels = 6
        self.one_hot = True
        self.reward_model = reward_model
        self._max_episode_timesteps = _max_episode_timesteps

        # Adversarial play
        if(self.use_double_agent):
            self.double_brain = self.unity_env.brain_names[1]
            self.double_agent = double_agent
            self.play_with_random = play_with_random

        self.global_timesteps = 0

        self.with_transformer = False

    count = 0

    def to_one_hot(self, a, channels):
        return (np.arange(channels) == a[..., None]).astype(float)

    def get_input_observation(self, env_info, action=None, transformer=False):
        size = self.size_global * self.size_global * self.input_channels

        global_in = env_info.vector_observations[0][:size]
        global_in = np.reshape(global_in, (self.size_global, self.size_global, self.input_channels))
        if self.one_hot:
            global_in_one_hot = self.to_one_hot(global_in[:,:,0], 7)
            for i in range(1, self.input_channels):
                global_in_one_hot = np.append(global_in_one_hot, self.to_one_hot(global_in[:,:,i], 9), axis = 2)

        if self.with_local:
            size_local = self.size_two * self.size_two * self.input_channels
            local_in = env_info.vector_observations[0][size:(size + size_local)]
            local_in = np.reshape(local_in, (self.size_two, self.size_two, self.input_channels))
            if self.one_hot:
                local_in_one_hot = self.to_one_hot(local_in[:, :, 0], 7)
                for i in range(1, self.input_channels):
                    local_in_one_hot = np.append(local_in_one_hot, self.to_one_hot(local_in[:, :, i], 9), axis=2)

            size_local_two = self.size_three * self.size_three * self.input_channels
            local_in_two = env_info.vector_observations[0][(size + size_local):(
                    size + size_local + size_local_two)]
            local_in_two = np.reshape(local_in_two, (self.size_three, self.size_three, self.input_channels))
            if self.one_hot:
                local_in_two_one_hot = self.to_one_hot(local_in_two[:, :, 0], 7)
                for i in range(1, self.input_channels):
                    local_in_two_one_hot = np.append(local_in_two_one_hot, self.to_one_hot(local_in_two[:, :, i], 9), axis=2)

        if self.with_local and self.with_stats:
            stats = env_info.vector_observations[0][
                    (size + size_local + size_local_two):
                    (size + size_local + size_local_two + self.size_stats)]

        if self.with_local and self.with_stats and self.with_hp:
            hp = env_info.vector_observations[0][
                 (size + (self.size_two * self.size_two) + (self.size_three * self.size_three)):
                 (size + (self.size_two * self.size_two) + (self.size_three * self.size_three) + 4)]
            stats = env_info.vector_observations[0][
                    (size + (self.size_two * self.size_two) + (self.size_three * self.size_three) + 4):
                    (size + (self.size_two * self.size_two) + (
                                self.size_three * self.size_three) + self.size_stats + 4)]
            self.size_stats += 4

        if self.with_local and self.with_stats and self.with_class:
            agent_class = env_info.vector_observations[0][
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats):
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats + self.size_class)]

            enemy_class = env_info.vector_observations[0][
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats + self.size_class):
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats + self.size_class * 2)]

        observation = {
            'global_in': global_in,
            'local_in': local_in
        }

        if self.with_local:
            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two
            }
        if self.with_local and self.with_stats:
            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two,
                'stats': stats
            }

        if self.with_local and self.with_stats and self.with_previous:
            action_vector = np.zeros(19)
            if action != None:
                action_vector[action] = 1

            observation = {
                # Global View
                'global_in': global_in_one_hot,
                # 'attr_global_1': global_in[:, :, 1],
                # 'attr_global_2': global_in[:, :, 2],
                # 'attr_global_3': global_in[:, :, 3],
                # 'attr_global_4': global_in[:, :, 4],
                # 'attr_global_5': global_in[:, :, 5],
                #'global_in_attributes': global_in[:, :, 1:]/5.,

                # Local View
                'local_in': local_in_one_hot,
                # 'attr_local_1': local_in[:, :, 1],
                # 'attr_local_2': local_in[:, :, 2],
                # 'attr_local_3': local_in[:, :, 3],
                # 'attr_local_4': local_in[:, :, 4],
                # 'attr_local_5': local_in[:, :, 5],
                #'local_in_attributes': local_in[:, :, 1:]/5.,

                # Local Two View
                'local_in_two': local_in_two_one_hot,
                # 'attr_local_two_1': local_in_two[:, :, 1],
                # 'attr_local_two_2': local_in_two[:, :, 2],
                # 'attr_local_two_3': local_in_two[:, :, 3],
                # 'attr_local_two_4': local_in_two[:, :, 4],
                # 'attr_local_two_5': local_in_two[:, :, 5],
                #'local_in_two_attributes': local_in_two[:, :, 1:],

                # Stats
                'agent_stats': stats[:16],
                'target_stats': stats[16:],
                'prev_action': action_vector
            }

            if transformer:
                x_map = np.asarray([[0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
                                     [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]])

                y_map = np.asarray([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                                     [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                                     [2, 2, 2, 2, 2, 2, 2, 2, 2, 2],
                                     [3, 3, 3, 3, 3, 3, 3, 3, 3, 3],
                                     [4, 4, 4, 4, 4, 4, 4, 4, 4, 4],
                                     [5, 5, 5, 5, 5, 5, 5, 5, 5, 5],
                                     [6, 6, 6, 6, 6, 6, 6, 6, 6, 6],
                                     [7, 7, 7, 7, 7, 7, 7, 7, 7, 7],
                                     [8, 8, 8, 8, 8, 8, 8, 8, 8, 8],
                                     [9, 9, 9, 9, 9, 9, 9, 9, 9, 9]])

                x_map = np.reshape(x_map, [10, 10, 1])
                y_map = np.reshape(y_map, [10, 10, 1])

                x_map = 2 * (x_map - 0) / (9 - 0) - 1
                y_map = 2 * (y_map - 0) / (9 - 0) - 1

                x_local_two_map = np.asarray([[0, 1, 2],
                                              [0, 1, 2],
                                              [0, 1, 2]])

                y_local_two_map = np.asarray([[0, 0, 0],
                                              [1, 1, 1],
                                              [2, 2, 2]])

                x_local_two_map = np.reshape(x_local_two_map, [3, 3, 1])
                y_local_two_map = np.reshape(y_local_two_map, [3, 3, 1])

                x_local_two_map = 2 * (x_local_two_map - 0) / (3 - 0) - 1
                y_local_two_map = 2 * (y_local_two_map - 0) / (3 - 0) - 1

                observation['global_positions'] = np.concatenate([x_map, y_map], axis=2)
                observation['local_two_positions'] = np.concatenate([x_local_two_map, y_local_two_map], axis=2)

        if self.with_local and self.with_stats and self.with_previous and self.with_class:
            action_vector = np.zeros(17)
            if action != None:
                action_vector[action] = 1

            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two,
                'stats': stats,
                'agent_class': agent_class,
                'enemy_class': enemy_class,
                'action': action_vector
            }

        if self.with_local and self.with_stats and self.with_previous and self.with_hp:
            action_vector = np.zeros(17)
            if action != None:
                action_vector[action] = 1

            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two,
                'hp': hp,
                'stats': stats,
                'prev_action': action_vector
            }

        return observation

    def get_input_observation_adapter(self, env_info, action = None):
        size = self.size_global * self.size_global * self.input_channels

        global_in = env_info.vector_observations[0][:size]
        global_in = np.reshape(global_in, (self.size_global, self.size_global, self.input_channels))
        if self.one_hot:
            global_in_one_hot = self.to_one_hot(global_in[:,:,0], 7)
            for i in range(1, self.input_channels):
                global_in_one_hot = np.append(global_in_one_hot, self.to_one_hot(global_in[:,:,i], 9), axis = 2)

        if self.with_local:
            size_local = self.size_two * self.size_two * self.input_channels
            local_in = env_info.vector_observations[0][size:(size + size_local)]
            local_in = np.reshape(local_in, (self.size_two, self.size_two, self.input_channels))
            if self.one_hot:
                local_in_one_hot = self.to_one_hot(local_in[:, :, 0], 7)
                for i in range(1, self.input_channels):
                    local_in_one_hot = np.append(local_in_one_hot, self.to_one_hot(local_in[:, :, i], 9), axis=2)

            size_local_two = self.size_three * self.size_three * self.input_channels
            local_in_two = env_info.vector_observations[0][(size + size_local):(
                    size + size_local + size_local_two)]
            local_in_two = np.reshape(local_in_two, (self.size_three, self.size_three, self.input_channels))
            if self.one_hot:
                local_in_two_one_hot = self.to_one_hot(local_in_two[:, :, 0], 7)
                for i in range(1, self.input_channels):
                    local_in_two_one_hot = np.append(local_in_two_one_hot, self.to_one_hot(local_in_two[:, :, i], 9), axis=2)

        if self.with_local and self.with_stats:
            stats = env_info.vector_observations[0][
                    (size + size_local + size_local_two):
                    (size + size_local + size_local_two + self.size_stats)]

        if self.with_local and self.with_stats and self.with_hp:
            hp = env_info.vector_observations[0][
                 (size + (self.size_two * self.size_two) + (self.size_three * self.size_three)):
                 (size + (self.size_two * self.size_two) + (self.size_three * self.size_three) + 4)]
            stats = env_info.vector_observations[0][
                    (size + (self.size_two * self.size_two) + (self.size_three * self.size_three) + 4):
                    (size + (self.size_two * self.size_two) + (
                                self.size_three * self.size_three) + self.size_stats + 4)]
            self.size_stats += 4

        if self.with_local and self.with_stats and self.with_class:
            agent_class = env_info.vector_observations[0][
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats):
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats + self.size_class)]

            enemy_class = env_info.vector_observations[0][
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats + self.size_class):
                          (size + (self.size_two * self.size_two) + (
                                      self.size_three * self.size_three) + self.size_stats + self.size_class * 2)]

        observation = {
            'global_in': global_in,
            'local_in': local_in
        }

        if self.with_local:
            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two
            }
        if self.with_local and self.with_stats:
            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two,
                'stats': stats
            }

        if self.with_local and self.with_stats and self.with_previous:
            action_vector = np.zeros(19)
            if action != None:
                action_vector[action] = 1

            observation = {
                # Global View
                'global_in': global_in_one_hot,
                # 'attr_global_1': global_in[:, :, 1],
                # 'attr_global_2': global_in[:, :, 2],
                # 'attr_global_3': global_in[:, :, 3],
                # 'attr_global_4': global_in[:, :, 4],
                # 'attr_global_5': global_in[:, :, 5],
                #'global_in_attributes': global_in[:, :, 1:]/5.,

                # Local View
                'local_in': local_in_one_hot,
                # 'attr_local_1': local_in[:, :, 1],
                # 'attr_local_2': local_in[:, :, 2],
                # 'attr_local_3': local_in[:, :, 3],
                # 'attr_local_4': local_in[:, :, 4],
                # 'attr_local_5': local_in[:, :, 5],
                #'local_in_attributes': local_in[:, :, 1:]/5.,

                # Local Two View
                'local_in_two': local_in_two_one_hot,
                # 'attr_local_two_1': local_in_two[:, :, 1],
                # 'attr_local_two_2': local_in_two[:, :, 2],
                # 'attr_local_two_3': local_in_two[:, :, 3],
                # 'attr_local_two_4': local_in_two[:, :, 4],
                # 'attr_local_two_5': local_in_two[:, :, 5],
                #'local_in_two_attributes': local_in_two[:, :, 1:],

                # Stats
                'agent_stats': stats[:16],
                'target_stats': stats[16:],
                'prev_action': action_vector
            }

        if self.with_local and self.with_stats and self.with_previous and self.with_class:
            action_vector = np.zeros(17)
            if action != None:
                action_vector[action] = 1

            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two,
                'stats': stats,
                'agent_class': agent_class,
                'enemy_class': enemy_class,
                'action': action_vector
            }

        if self.with_local and self.with_stats and self.with_previous and self.with_hp:
            action_vector = np.zeros(17)
            if action != None:
                action_vector[action] = 1

            observation = {
                'global_in': global_in,
                'local_in': local_in,
                'local_in_two': local_in_two,
                'hp': hp,
                'stats': stats,
                'prev_action': action_vector
            }

        return observation

    def print_observation(self, observation, actions = None, reward = None):
        try:
            # print(observation)
            # print(observation['global_in'])
            print('action = ' + str(actions))
            print('reward = ' + str(reward))
            sum = observation['global_in'][:,:,0]*0
            for i in range(1, 7):
                sum += observation['global_in'][:,:,i]*i
            sum = np.flip(np.transpose(sum), 0)
            print(sum)
            print(" ")
            print(observation['agent_stats'])
            print(observation['target_stats'])
            # print(observation['local_in'])
            # print(observation['local_in_two'])

        except Exception as e:
            pass

    def execute(self, actions, raw_obs=False):

        if self.manual_input:
            input_action = input('...')

            try:
                actions = input_action
            except ValueError:
                pass

        if isinstance(actions, str):
            actions = self.command_to_action(actions)
            print(actions)

        env_info = None
        #signal.alarm(0)
        while env_info == None:
            #signal.signal(signal.SIGALRM, self.handler)
            #signal.alarm(300000)
            try:
                if self.use_double_agent:
                    info = self.unity_env.step({self.default_brain : [actions], self.double_brain : []})
                    env_info = info[self.default_brain]
                else:
                    env_info = self.unity_env.step([actions])[self.default_brain]
            except Exception as exc:
                self.close()
                self.unity_env = self.open_unity_environment(self.game_name, self.no_graphics, seed = int(time.time()),
                                                             worker_id=self.worker_id)
                env_info = self.unity_env.reset(train_mode=True, config=self.config)[self.default_brain]
                print("The environment didn't respond, it was necessary to close and reopen it")

        if self.use_double_agent:
            while len(env_info.vector_observations) <= 0:
                double_info = info[self.double_brain]
                obs = self.get_input_observation(double_info)
                if np.random.rand() > self.play_with_random:
                    act = self.double_agent.eval([obs])[0]
                else:
                    act = np.random.randint(0, self.action_size)
                info = self.unity_env.step({self.default_brain: [], self.double_brain: [act]})
                env_info = info[self.default_brain]

        reward = env_info.rewards[0]
        done = env_info.local_done[0]

        observation = self.get_input_observation(env_info, actions)

        self.count += 1

        # Update global timesteps
        self.global_timesteps += 1

        if self.verbose:
            self.print_observation(observation, reward=reward, actions=actions)

        if not raw_obs:
            return [observation, done, reward]
        else:
            return [env_info, done, reward]

    def command_to_action(self, command):

        switcher = {
            "w": 0,
            "d": 1,
            "x": 2,
            "a": 3,
            "e": 4,
            "c": 5,
            "z": 6,
            "q": 7,
            "s": 8,
            "y": 9,
            "j": 10,
            "n": 11,
            "g": 12,
            "u": 13,
            "m": 14,
            "b": 15,
            "t": 16,
            "f": 17,
            "r": 18
        }

        return switcher.get(command, 99)

    def set_config(self, config):
        self.config = config

    def handler(self, signum, frame):
        print("Timeout!")
        raise Exception("end of time")

    def reset(self, raw_obs=False):

        self.count = 0

        env_info = None

        while env_info == None:
            #signal.signal(signal.SIGALRM, self.handler)
            #signal.alarm(60)
            try:
                logging.getLogger("mlagents.envs").setLevel(logging.WARNING)
                env_info = self.unity_env.reset(train_mode=True, config=self.config)[self.default_brain]
            except Exception as exc:
                self.close()
                self.unity_env = self.open_unity_environment(self.game_name, self.no_graphics, seed=int(time.time()),
                                                             worker_id=self.worker_id)
                env_info = self.unity_env.reset(train_mode=True, config=self.config)[self.default_brain]
                print("The environment didn't respond, it was necessary to close and reopen it")

        if self.use_double_agent:
            while len(env_info.vector_observations) <= 0:
                env_info = self.unity_env.step()[self.default_brain]

        observation = self.get_input_observation(env_info)

        if self.verbose:
            self.print_observation(observation)

        if not raw_obs:
            return observation
        else:
            return env_info

    def close(self):
        self.unity_env.close()

    def open_unity_environment(self, game_name, no_graphics, seed, worker_id):
        return UnityEnvironment(game_name, no_graphics=no_graphics, seed=seed, worker_id=worker_id)

    def add_probs(self, probs):
        self.probabilities.append(probs[0])

    def get_last_entropy(self):
        entropy = 0
        for prob in self.probabilities[-1]:
            entropy += (prob + eps)*(math.log(prob + eps) + eps)

        return -entropy

    def entropy(self, probs):
        entropy = 0
        for prob in probs:
            entropy += (prob + eps) * (math.log(prob + eps) + eps)

        return -entropy

    def states(self):
        return {
            # Categorical values: 2 agent, 1 empty tile, 0 obstacle tile
            # 'global_in': {'shape': (10, 10, 52), 'type': 'float'},
            # 'local_in': {'shape': (5, 5, 52), 'type': 'float'},
            # 'local_in_two': {'shape': (3, 3, 52), 'type': 'float'},
            'global_in': {'shape': (10, 10, 52), 'type': 'float'},
            'local_in': {'shape': (5, 5, 52), 'type': 'float'},
            'local_in_two': {'shape': (3, 3, 52), 'type': 'float'},

            # Transformer
            # I have a max of 20 items, 1 agent and 1 target. Each entity is embedded to have dimension 256.
            # All entities are concatenated together to form a 22(number of entities)x256(number of features)
            # matrix, this will be passed to the Transformer layer.
            # The output of the transformer will alwyas be a 22x256 matrix, and this will be average pooled along the
            # feature axis to form a 256 feature vector.
            # 'items': {'shape': (20, 50), 'type': 'float'},
            # 'agent': {'shape': (1, 85), 'type': 'float'},
            # 'target': {'shape': (1, 8), 'type': 'float'},

            # 'agent_stats': {'shape': (14), 'num_values': 128, 'type': 'int'},
            # 'target_stats': {'shape': (4), 'num_values': 32, 'type': 'int'},

            'agent_stats': {'shape': (16), 'num_values': 129, 'type': 'int'},
            'target_stats': {'shape': (15), 'num_values': 125, 'type': 'int'},

            'prev_action': {'shape': (19), 'type': 'float'}
        }

    def actions(self):
        return {
                    'type': 'int',
                    'num_values': 19
                }

    def max_episode_timesteps(self):
        return self._max_episode_timesteps

class Info():
    def __init__(self, string):
        self.item = string

    def items(self):

        return self.item, self.item


text:
""" Analyze MCMC output - chain length, etc. """

# Built-in libraries
import os
import pickle
# External libraries
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import median_abs_deviation
import xarray as xr
# Local libraries
import pygem.pygem_input as pygem_prms

#from oggm import cfg, utils
from pygem.oggm_compat import single_flowline_glacier_directory

#%%
regions = [1]

option_regional_mb = True

modelprms_fp = pygem_prms.output_filepath + 'calibration/'

#%%
if option_regional_mb:
    debug = True

    reg_df_fp = pygem_prms.output_filepath + 'analysis/'
    reg_df_fn = 'mcmc_reg_comparison.csv'

    for reg in regions:

        # regional filepath
        modelprms_fp_reg = modelprms_fp + str(reg).zfill(2) + '/'

        # find calibrated glaciers
        glac_strs = []
        for i in os.listdir(modelprms_fp_reg):
            if i.endswith('-modelprms_dict.pkl'):
                glac_strs.append(i.split('-')[0])
        glac_strs = sorted(glac_strs)

        # Aggregate regional mass balance data
        reg_mb_m3wea_all = None
        count_glac = 0
        reg_area = 0
        for nglac, glacier_str in enumerate(glac_strs):

            if nglac%500 == 0:
                print(glacier_str)

            # Load calibration data
            try:
                gdir = single_flowline_glacier_directory(glacier_str, logging_level='CRITICAL')
                fls = gdir.read_pickle('inversion_flowlines')
                glacier_area = np.sum(fls[0].widths_m * fls[0].dx_meter)

                mbdata_fn = gdir.get_filepath('mb_obs')
                with open(mbdata_fn, 'rb') as f:
                    gdir.mbdata = pickle.load(f)

                # Load data
                mb_mwea_obs = gdir.mbdata['mb_mwea']
                mb_mwea_err_obs = gdir.mbdata['mb_mwea_err']
                mb_m3wea_obs = mb_mwea_obs * glacier_area

            except:
                gdir = None

            # Load model parameters
            if gdir is not None:
                try:
                    modelprms_fn = glacier_str + '-modelprms_dict.pkl'
                    modelprms_fullfn = modelprms_fp_reg + modelprms_fn
                    with open(modelprms_fullfn, 'rb') as f:
                        modelprms_dict = pickle.load(f)

                    modelprms_all = modelprms_dict['MCMC']
                    mb_mwea_all = np.array(modelprms_all['mb_mwea']['chain_0'])
                    # convert to m3wea for regional aggregation
                    mb_m3wea_all = mb_mwea_all * glacier_area

                    # add to regional mass balance
                    if reg_mb_m3wea_all is None:
                        reg_mb_m3wea_all = mb_m3wea_all
                        reg_mb_m3wea_obs = mb_m3wea_obs
                    else:
                        reg_mb_m3wea_all += mb_m3wea_all
                        reg_mb_m3wea_obs += mb_m3wea_obs

                    # Track statistics
                    count_glac += 1
                    reg_area += glacier_area

                except:
                    print('\n\nNo calibration data for ' + glacier_str)

        # convert back to mwea
        reg_mb_mwea = reg_mb_m3wea_all / reg_area
        reg_mb_mwea_obs = reg_mb_m3wea_obs / reg_area

        # Statistics for comparison
        reg_mb_mwea_mean = np.mean(reg_mb_mwea)
        reg_mb_mwea_std = np.std(reg_mb_mwea)
        reg_mb_mwea_med = np.median(reg_mb_mwea)
        reg_mb_mwea_mad = median_abs_deviation(reg_mb_mwea)
        reg_mb_mwea_perc5 = np.percentile(reg_mb_mwea, 5)
        reg_mb_mwea_perc25 = np.percentile(reg_mb_mwea, 25)
        reg_mb_mwea_perc75 = np.percentile(reg_mb_mwea, 75)
        reg_mb_mwea_perc95 = np.percentile(reg_mb_mwea, 95)

        print('reg_mb_mwea_med:', np.round(reg_mb_mwea_med,2), 'mb_mwea_obs:', np.round(reg_mb_mwea_obs,2))

        reg_df_cns = ['O1Region', 'count', 'area_m2', 'mb_mwea_mean', 'mb_mwea_std', 'mb_mwea_med',
                      'mb_mwea_mad', 'mb_mwea_perc5', 'mb_mwea_perc25', 'mb_mwea_perc75', 'mb_mwea_perc95']
        reg_df_single = pd.DataFrame(np.zeros((1,len(reg_df_cns))), columns=reg_df_cns)
        reg_df_single.loc[0,:] = [reg, count_glac, reg_area, reg_mb_mwea_mean, reg_mb_mwea_std,
                                  reg_mb_mwea_med, reg_mb_mwea_mad, reg_mb_mwea_perc5, reg_mb_mwea_perc25,
                                  reg_mb_mwea_perc75, reg_mb_mwea_perc95]

        # Export csv
        if os.path.exists(reg_df_fp + reg_df_fn):
            reg_df = pd.read_csv(reg_df_fp + reg_df_fn)

            # Add or overwrite existing file
            reg_idx = np.where((reg_df.O1Region == reg))[0]
            if len(reg_idx) > 0:
                reg_df.loc[reg_idx,:] = reg_df_single.values
            else:
                reg_df = pd.concat([reg_df, reg_df_single], axis=0)

        else:
            reg_df = reg_df_single

        reg_df = reg_df.sort_values('O1Region', ascending=True)
        reg_df.reset_index(inplace=True, drop=True)
        reg_df.to_csv(reg_df_fp + reg_df_fn, index=False)

max_stars_repo_name: Ruitangtang/PyGEM

max_issues_repo_name: Ruitangtang/PyGEM

max_forks_repo_name: Ruitangtang/PyGEM

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-

# @Time    : 19-7-15 上午11:33
# @Author  : zj

from builtins import range
from knn_classifier import KNN
import pandas as pd
import numpy as np
from sklearn import utils
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import warnings

warnings.filterwarnings("ignore")

def load_iris(iris_path, shuffle=True, tsize=0.8):
    """
    加载iris数据
    """
    data = pd.read_csv(iris_path, header=0, delimiter=',')

    if shuffle:
        data = utils.shuffle(data)

    species_dict = {
        'Iris-setosa': 0,
        'Iris-versicolor': 1,
        'Iris-virginica': 2
    }
    data['Species'] = data['Species'].map(species_dict)

    data_x = np.array(
        [data['SepalLengthCm'], data['SepalWidthCm'], data['PetalLengthCm'], data['PetalWidthCm']]).T
    data_y = np.array(data['Species'])

    x_train, x_test, y_train, y_test = train_test_split(data_x, data_y, train_size=tsize, test_size=(1 - tsize),
                                                        shuffle=False)

    return np.array(x_train), np.array(x_test), np.array(y_train), np.array(y_test)

def load_german_data(data_path, shuffle=True, tsize=0.8):
    data_list = pd.read_csv(data_path, header=None, sep='\s+')

    data_array = data_list.values
    height, width = data_array.shape[:2]
    data_x = data_array[:, :(width - 1)]
    data_y = data_array[:, (width - 1)]

    x_train, x_test, y_train, y_test = train_test_split(data_x, data_y, train_size=tsize, test_size=(1 - tsize),
                                                        shuffle=shuffle)

    y_train = np.array(list(map(lambda x: 1 if x == 2 else 0, y_train)))
    y_test = np.array(list(map(lambda x: 1 if x == 2 else 0, y_test)))

    return x_train, x_test, y_train, y_test

def compute_accuracy(y, y_pred):
    num = y.shape[0]
    num_correct = np.sum(y_pred == y)
    acc = float(num_correct) / num
    return acc

def cross_validation(x_train, y_train, k_choices, num_folds=5, Classifier=KNN):
    X_train_folds = np.array_split(x_train, num_folds)
    y_train_folds = np.array_split(y_train, num_folds)

    # 计算预测标签和验证集标签的精度
    k_to_accuracies = {}
    for k in k_choices:
        k_accuracies = []
        # 随机选取其中一份为验证集，其余为测试集
        for i in range(num_folds):
            x_folds = X_train_folds.copy()
            y_folds = y_train_folds.copy()

            x_vals = x_folds.pop(i)
            x_trains = np.vstack(x_folds)

            y_vals = y_folds.pop(i)
            y_trains = np.hstack(y_folds)

            classifier = Classifier()
            classifier.train(x_trains, y_trains)

            y_val_pred = classifier.predict(x_vals, k=k)
            k_accuracies.append(compute_accuracy(y_vals, y_val_pred))
        k_to_accuracies[k] = k_accuracies

    return k_to_accuracies

def plot(k_choices, k_to_accuracies):
    # plot the raw observations
    for k in k_choices:
        accuracies = k_to_accuracies[k]
        plt.scatter([k] * len(accuracies), accuracies)

    # plot the trend line with error bars that correspond to standard deviation
    accuracies_mean = np.array([np.mean(v) for k, v in sorted(k_to_accuracies.items())])
    accuracies_std = np.array([np.std(v) for k, v in sorted(k_to_accuracies.items())])
    plt.errorbar(k_choices, accuracies_mean, yerr=accuracies_std)
    plt.title('Cross-validation on k')
    plt.xlabel('k')
    plt.ylabel('Cross-validation accuracy')
    plt.show()

if __name__ == '__main__':
    iris_path = '/home/zj/data/iris-species/Iris.csv'
    x_train, x_test, y_train, y_test = load_iris(iris_path, shuffle=True, tsize=0.8)

    # data_path = '/home/zj/data/german/german.data-numeric'
    # x_train, x_test, y_train, y_test = load_german_data(data_path, shuffle=True, tsize=0.8)

    x_train = x_train.astype(np.double)
    x_test = x_test.astype(np.double)
    mu = np.mean(x_train, axis=0)
    var = np.var(x_train, axis=0)
    eps = 1e-8
    x_train = (x_train - mu) / np.sqrt(var + eps)
    x_test = (x_test - mu) / np.sqrt(var + eps)

    k_choices = [1, 3, 5, 8, 10, 12, 15, 20, 30, 50, 100]
    k_to_accuracies = cross_validation(x_train, y_train, k_choices)

    # print(k_to_accuracies)
    # Print out the computed accuracies
    for k in sorted(k_to_accuracies):
        for accuracy in k_to_accuracies[k]:
            print('k = %d, accuracy = %f' % (k, accuracy))

    plot(k_choices, k_to_accuracies)

    accuracies_mean = np.array([np.mean(v) for k, v in sorted(k_to_accuracies.items())])
    k = k_choices[np.argmax(accuracies_mean)]
    print('最好的k值是：%d' % k)

    # 测试集测试
    classifier = KNN()
    classifier.train(x_train, y_train)

    y_test_pred = classifier.predict(x_test, k=k)
    y_test_acc = compute_accuracy(y_test, y_test_pred)
    print('测试集精度为：%f' % y_test_acc)


text:
"""
 Experimental sample made with data from K. H. Kim et alii report
"""

import numpy as np
import pandas as pd
import apl104lib as apl104

# Settings
expsample=apl104.Sample1D(5,400,256)
expsample.specienames=['La','Zr','Co','O','Li']
padding=-4

# La
Ladata=pd.read_csv('Datasets/La',header=None)
Ladata=Ladata[(Ladata[1]<1020)&(Ladata[0]>0)&(Ladata[1]>0)] # Filtering
Ladata=Ladata.sort_values(0,axis=0) # Sorting
Ladata[0]=Ladata[0]+padding # Padding
Ladata[2]=Ladata[1]/Ladata[Ladata[0]>150][1].mean()*3/24 # Normalization
Ladata[2]=Ladata[2]-Ladata[2].min() # Flooring

# Zr
Zrdata=pd.read_csv('Datasets/Zr',header=None)
Zrdata=Zrdata[(Zrdata[1]<430)&(Zrdata[0]>0)&(Zrdata[1]>0)] # Filtering
Zrdata=Zrdata.sort_values(0,axis=0) # Sorting
Zrdata[0]=Zrdata[0]+padding # Padding
Zrdata[2]=Zrdata[1]/Zrdata[Zrdata[0]>150][1].mean()*2/24 # Normalization
Zrdata[2]=Zrdata[2]-Zrdata[2].min() # Flooring

# Co
Codata=pd.read_csv('Datasets/Co',header=None)
Codata=Codata[(Codata[1]<1100)&(Codata[0]>0)&(Codata[1]>0)] # Filtering
Codata=Codata.sort_values(0,axis=0) # Sorting
Codata[0]=Codata[0]+padding # Padding
Codata[2]=Codata[1]/Codata[Codata[0]<50][1].mean()*1/4 # Normalization
Codata[2]=Codata[2]-Codata[2].min() # Flooring

# O
Odata=pd.read_csv('Datasets/O',header=None)
Odata=Odata[(Odata[1]<390)&(Odata[0]>0)&(Odata[1]>0)] # Filtering
Odata=Odata.sort_values(0,axis=0) # Sorting
Odata[0]=Odata[0]+padding # Padding
Odata[2]=Odata[1]/Odata[1].mean()*2/4 # Normalization
# ~ Odata[2]=Odata[1]-Odata[1].min() # Flooring

# Define experimental sample
Lapoints=np.interp(np.linspace(0,int(expsample.dim/2),int(expsample.res/2)),Ladata[0],Ladata[2])
Zrpoints=np.interp(np.linspace(0,int(expsample.dim/2),int(expsample.res/2)),Zrdata[0],Zrdata[2])
Copoints=np.interp(np.linspace(0,int(expsample.dim/2),int(expsample.res/2)),Codata[0],Codata[2])
Opoints=np.interp(np.linspace(0,int(expsample.dim/2),int(expsample.res/2)),Odata[0],Odata[2])

expsample.c[:,0]=np.concatenate([Lapoints,np.flip(Lapoints,axis=0)]) #La
expsample.c[:,1]=np.concatenate([Zrpoints,np.flip(Zrpoints,axis=0)]) #Zr
expsample.c[:,2]=np.concatenate([Copoints,np.flip(Copoints,axis=0)]) #Co
expsample.c[:,3]=np.concatenate([Opoints,np.flip(Opoints,axis=0)])	 #O
expsample.c[:,4]=1-np.sum(expsample.c,axis=-1)						 #Li
expsample.c[:,4][expsample.c[:,4]<0]=0								 #Filter negative values


text:
"""
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3d hyperbolic honeycombs in Poincaré's ball model
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Note for ideal and hyperideal cases you must make
sure (p, q, r) are ordered so that (p, q) form a
finite spherical polyhedra, this is due to our
construction of the reflection mirrors. For example
(6, 3, 3) won't work, but (3, 3, 6) does.
"""
from functools import partial
from itertools import combinations
import numpy as np
import tqdm
import helpers
from coxeter import CoxeterGroup

def vround(v):
    """Convert a numpy 1d array to a tuple for hashing.
    """
    return tuple(np.round(v, 8))

def centroid(points):
    """Compute the centroid of a list of points.
    """
    mid = helpers.normalize(np.sum(points, axis=0))
    return vround(mid)

class Cell(object):

    """A cell is a uniform polyhedra (or a "room") in the 3d tiling.
    """

    def __init__(self, cox_mat, v0, active, reflections):
        self.cox_mat = cox_mat
        self.v0 = v0
        self.active = active
        self.reflections = reflections
        self.G = CoxeterGroup(cox_mat)

        self.vertices_coords = []
        self.num_vertices = None
        self.edge_coords = []
        self.num_edges = None

    def build_geometry(self, depth, maxcount):
        self.G.init()
        self.word_generator = partial(self.G.traverse, depth=depth, maxcount=maxcount)
        self.get_vertices()
        self.get_edges()
        return self

    def transform(self, word, v):
        for w in reversed(word):
            v = self.reflections[w](v)
        return vround(v)

    def project(self, v):
        return helpers.project_poincare(v)

    def get_vertices(self):
        for word in self.word_generator():
            v = self.transform(word, self.v0)
            if v not in self.vertices_coords:
                self.vertices_coords.append(v)
        self.num_vertices = len(self.vertices_coords)

    def get_edges(self):
        """
        An edge is uniquely determined by the coordinates of its
        middle point. Here I simply use a set to maintain the middle
        points of known edges and avoid duplicates.
        """
        edgehash = set()
        for i, active in enumerate(self.active):
            if active:
                for word in self.word_generator():
                    p1 = self.transform(word, self.v0)
                    p2 = self.transform(word + (i,), self.v0)
                    q = centroid((p1, p2))
                    if q not in edgehash:
                        self.edge_coords.append((p1, p2))
                        edgehash.add(q)
        self.num_edges = len(self.edge_coords)

class Honeycomb(object):

    def __init__(self, coxeter_diagram, init_dist):
        if len(coxeter_diagram) != 6 or len(init_dist) != 4:
            raise ValueError("Invalid input dimension")

        # Coxeter matrix and its rank
        self.cox_mat = helpers.make_symmetry_matrix(coxeter_diagram)
        self.rank = len(self.cox_mat)

        # generators of the symmetry group
        self.gens = tuple(range(self.rank))

        # symmetry group of this tiling
        self.G = CoxeterGroup(self.cox_mat)

        # a mirror is active iff the initial point is not on it
        self.active = tuple(bool(x) for x in init_dist)

        # reflection mirrors
        self.mirrors = self.get_mirrors(coxeter_diagram)

        # reflections (possibly affine) about the mirrors
        self.reflections = self.get_reflections()

        # coordinates of the initial point
        self.init_v = self.get_init_point(init_dist)

        self.edge_hash_set = set()

        self.num_vertices = 0
        self.num_edges = 0

    def get_reflections(self):
        def reflect(v, normal):
            return v - 2 * np.dot(v, normal) * normal

        return [partial(reflect, normal=n) for n in self.mirrors]

    def transform(self, word, v):
        for w in reversed(word):
            v = self.reflections[w](v)
        return vround(v)

    def project(self, v):
        return helpers.project_poincare(v)

    def get_init_point(self, init_dist):
        return helpers.get_point_from_distance(self.mirrors, init_dist)

    def get_mirrors(self, coxeter_diagram):
        return helpers.get_hyperbolic_honeycomb_mirrors(coxeter_diagram)

    def get_fundamental_cells(self, depth=None, maxcount=20000):
        """
        Generate the fundamental cells of the tiling, these cells
        are centered at the vertices of the fundamental tetrahedron
        and are generated by reflecting the initial point about
        the three mirrors meeting at each vertex.
        """
        result = {}
        for triple in combinations(self.gens, 3):
            cox_mat = self.cox_mat[np.ix_(triple, triple)]
            refs = [self.reflections[k] for k in triple]
            active = [self.active[k] for k in triple]
            if not helpers.is_degenerate(cox_mat, active):
                C = Cell(cox_mat, self.init_v, active, refs).build_geometry(depth, maxcount)
                result[triple] = C
        return result

    def is_new_edge(self, edge):
        mid = centroid(edge)
        if mid not in self.edge_hash_set:
            self.edge_hash_set.add(mid)
            return True
        return False

    def collect_fundamental_cell_edges(self):
        result = []
        for C in self.fundamental_cells.values():
            for edge in C.edge_coords:
                if self.is_new_edge(edge):
                    result.append(edge)
        return result

    def export_edge(self, fobj, p1, p2):
        """Export the data of an edge to POV-Ray .inc file."""
        fobj.write("HyperbolicEdge({}, {})\n".format(
            helpers.pov_vector(p1),
            helpers.pov_vector(p2)))

    def generate_povray_data(self, depth=100, maxcount=50000,
                             cell_depth=None, cell_edges=10000,
                             filename="./povray/honeycomb-data.inc",
                             eye=(0, 0, 0.5),
                             lookat=(0, 0, 0)):
        self.G.init()
        self.word_generator = partial(self.G.traverse, depth=depth, maxcount=maxcount)
        self.fundamental_cells = self.get_fundamental_cells(cell_depth, cell_edges)
        init_edges = self.collect_fundamental_cell_edges()
        bar = tqdm.tqdm(desc="processing edges", total=maxcount)
        vertices = set()
        eye = np.array(eye)
        lookat = np.array(lookat)
        viewdir = helpers.normalize(lookat - eye)

        def add_new_edge(edge):
            p1 = self.project(edge[0])
            p2 = self.project(edge[1])
            if np.dot(p1 - eye, viewdir) > 0.5 or np.dot(p2 - eye, viewdir) > 0.5:
                self.export_edge(f, p1, p2)
                self.num_edges += 1
                for v in [p1, p2]:
                    v = vround(v)
                    if v not in vertices:
                        vertices.add(v)
                        self.num_vertices += 1

        with open(filename, "w") as f:
            f.write("#declare camera_loc = {};\n".format(helpers.pov_vector(eye)))
            f.write("#declare lookat = {};\n".format(helpers.pov_vector(lookat)))
            for edge in init_edges:
                add_new_edge(edge)

            for word in self.word_generator():
                for edge in init_edges:
                    edge = [self.transform(word, v) for v in edge]
                    if self.is_new_edge(edge):
                        add_new_edge(edge)

                bar.update(1)
            bar.close()
            verts = "#declare num_vertices = {};\n"
            verts_coords = "#declare vertices = array[{}]{{{}}};\n"
            print("{} vertices and {} edges generated".format(self.num_vertices, self.num_edges))
            f.write(verts.format(self.num_vertices))
            f.write(verts_coords.format(self.num_vertices, helpers.pov_vector_list(vertices)))

max_stars_repo_name: wlbksy/pywonderland

max_issues_repo_name: wlbksy/pywonderland

max_forks_repo_name: wlbksy/pywonderland


text:
from data_loader.conv_emnist_data_loader import ConvEMnistDataLoader
from data_loader.conv_mnist_data_loader import ConvMnistDataLoader
from data_loader.simple_mnist_data_loader import SimpleMnistDataLoader
from data_visualizer.simple_mnist_data_visualizer import SimpleMnistDataVisualizer
from evaluater.conv_mnist_data_predictor import ConvMnistDataPredictor
from models.conv_emnist_model import ConvEMnistModel
from models.conv_mnist_model import ConvMnistModel
from models.simple_mnist_model import SimpleMnistModel
from trainers.conv_mnist_trainer import ConvMnistModelTrainer
from trainers.simple_mnist_trainer import SimpleMnistModelTrainer
from utils.config import process_config, process_image
from utils.dirs import create_dirs
from utils.utils import get_args
import numpy as np
import sklearn.metrics as metrics

def main():
    # capture the config path from the run arguments
    # then process the json configuration file
    global image
    try:
        args = get_args()
        config = process_config(args.config)
        image = process_image(args.image) if args.image is not 'None' else None
    except:
        print("missing or invalid arguments (check correct config or image paths)")
        exit(0)

    # create the experiments dirs
    create_dirs([config.callbacks.tensorboard_log_dir, config.callbacks.checkpoint_dir, config.callbacks.history_dir])

    print('Create the data generator.')
    data_loader = ConvEMnistDataLoader(config)  # Load data set

    print('Some data visualization')
    X_train, y_train = data_loader.get_train_data()  # Get training data
    print("ytrain")
    print(y_train.shape)
    mapp = data_loader.get_map()  # Get map dictionary (Refer to emnist-balanced-mapping.txt file)
    data_visualizer = SimpleMnistDataVisualizer(X_train, y_train, mapp)
    data_visualizer.plot_first_digit()  # Plot first character of training set
    data_visualizer.plot_range()  # Plot several characters

    print('Create the model.')
    model = ConvEMnistModel(config)  # Create the model based on configuration file

    print("Model Summary")
    model.model.summary()  # Print a summary of the model with the respective parameters

    # Custom weight to use instead of training the model
    weight = config.evaluator.weight

    print('Create the trainer')
    trainer = ConvMnistModelTrainer(model.model, data_loader.get_train_data(), config)

    if not config.evaluator.custom_weight:
        print('Start training the model.')
        # if not config.evaluator.custom_weight:
        trainer.train()

        print("Plot loss and accuracy in training model")
        data_visualizer.plot_loss_acc()

        print("Finish training")

    print("Predict")

    predictor = ConvMnistDataPredictor(model.model, data_loader.get_test_data(), mapp, config, weight)
    predict_image = image if image is not None else './test_images/data_representation/0.png'
    predicted_values = predictor.ocr(predict_image)
    print("Predicted values")
    print(predicted_values)
    # predictor.predict3('./test_images/h/1.png')
    # predictor.predict_from_data_set()

    """
    Evaluate model with test set
    """
    predictor.evaluate_model()

    predictor.confusion_matrix()

if __name__ == '__main__':
    main()
    # app.run()