# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np
seed = 14
np.random.seed(seed)
import keras 
from keras import backend as K, activations
from keras.models import Sequential, Model
from sklearn.model_selection import GridSearchCV
from keras.wrappers.scikit_learn import KerasClassifier
from preprocessing import*
import numpy as np
from keras.layers import Input, Dense, Activation, Dropout
from numpy.random import seed as set_seed
from sklearn.model_selection import train_test_split

features.head(10)

GE.head(10)


# # Building the Neural Network

def tune_shallow(optimizer,activation,n_dimension):
    model = Sequential() 
    model.add(Dense(n_dimension,input_dim=414, activation="relu",kernel_regularizer=regularizers.l2(0.01), name="hidd_layer_1"))
    model.add(Dense(n_dimension, activation="relu",kernel_regularizer=regularizers.l2(0.01),name="hidd_layer_2"))
    model.add(Dense(4, activation="softmax",name="out_layer"))
    model.compile(loss='categorical_crossentropy', optimizer=optimizer, metrics=['accuracy'])
    return model


# # Hyper Paramters Tuning

param_grid = {
             "n_dimension":[8,10,11,12,16],
              "optimizer": ['adam'],
              "activation": ['relu'],
              "batch_size": [16],
              "epochs":[25]
             }
classifier = KerasClassifier(tune_shallow)
validator = GridSearchCV(classifier,param_grid, cv = 10, n_jobs=-1,verbose=0)
validator.fit(X, dummyY,verbose =0 )

validator.best_params_

validator.best_score_

# # Model Fitting

x_train, x_test, y_train, y_test = train_test_split(X, dummyY, test_size=0.33,random_state=123)

myModel = None
myModel = tune_shallow(optimizer='adam' ,activation='relu',n_dimension=8)
history = myModel.fit(x_train, y_train, epochs=50, batch_size=16,verbose = 0)
scores = myModel.evaluate(x_test, y_test, verbose=0)
print("Model accuracy: " + str(scores[1]*100))

print(myModel.summary())

# # Get the weights

get_0rd_layer_output = K.function([myModel.layers[0].input],
                                  [myModel.layers[0].output])
layer_output = get_0rd_layer_output([X])
layer_1=[]
for i in layer_output[0]:
    layer_1.append(np.argmax(i))

get_1rd_layer_output = K.function([myModel.layers[0].input],
                                  [myModel.layers[1].output])
layer_output = get_1rd_layer_output([X])
layer_2=[]
for i in layer_output[0]:
    layer_2.append(np.argmax(i))

get_2rd_layer_output = K.function([myModel.layers[0].input],
                                  [myModel.layers[2].output])
layer_output = get_2rd_layer_output([X])
layer_3=[]
for i in layer_output[0]:
    layer_3.append(np.argmax(i))

overall = []
for i in range(0,2441):
    overall.append(str(layer_1[i]) + str(layer_2[i]) + str(layer_3[i]))
    #print(str(layer_1[i]) + str(layer_2[i]) + str(layer_3[i]))

from collections import Counter
Counter(overall)

classes
# 310: IFNb
# 401: LPS
# 542: dNS1
# 733: unstim

overall_weights = []
gene_list = list(GE.columns)
dict_aux = {}
for i in range(0,len(myModel.layers[0].get_weights()[0])):
    l1 = np.argmax(myModel.layers[0].get_weights()[0][i])
    l2 = np.argmax(myModel.layers[1].get_weights()[0][l1])
    l_out = np.argmax(myModel.layers[2].get_weights()[0][l2])
    overall_weights.append(str(l1)+str(l2)+str(l_out))
    path = str(l1)+str(l2)+str(l_out)
    if path not in dict_aux.keys():
        dict_aux[path] = []
        dict_aux[path].append(gene_list[i])
    else:
       
        dict_aux[path].append(gene_list[i])

for i in range(0,len(myModel.layers[0].get_weights()[0])):
        l1_f_w_p = np.argsort(myModel.layers[0].get_weights()[0][i])[-1]
        l1_s_w_p = np.argsort(myModel.layers[0].get_weights()[0][i])[-2]
        l2_f_w_p = np.argsort(myModel.layers[1].get_weights()[0][l1_f_w_p])[-1]
        l2_s_w_p = np.argsort(myModel.layers[1].get_weights()[0][l1_s_w_p])[-2]
        l_out_f_w_p = np.argmax(myModel.layers[2].get_weights()[0][l2_f_w_p])
        l_out_s_w_p = np.argmax(myModel.layers[2].get_weights()[0][l2_s_w_p])


overall_weights = []
gene_list = list(GE.columns)
dict_aux = {}
for i in range(0,len(myModel.layers[0].get_weights()[0])):
    l1 = np.argmax(myModel.layers[0].get_weights()[0][i])
    l2 = np.argmax(myModel.layers[1].get_weights()[0][l1])
    l_out = np.argmax(myModel.layers[2].get_weights()[0][l2])
    overall_weights.append(str(l1)+str(l2)+str(l_out))
    path = str(l1)+str(l2)+str(l_out)
    #print(path)


# # Part B

# +
def get_output(ln,X=X):
    get_layer_output = K.function([myModel.layers[0].input],
                                  [myModel.layers[ln].output])
    return get_layer_output([X])

def get_paths(nclasses,ylabels,nlayers,model):
    paths_classes = []
    for ci in range(nclasses):
        curr_class = np.where(ylabels[:,ci] == 1)[0]
        all_neur = []    
        
        for li in range(nlayers):
            all_neur.append([np.argsort(l)[::-1][0] for l in get_output(li)[0][curr_class]])

        count_curr_class = Counter(all_neur[0])
        count_i = list(count_curr_class.keys())[0]
        paths_classes.append([count_curr_class])
        
        for ni in range(1,nlayers):
                
            buff = Counter([all_neur[ni][i] for i in np.where(np.array(all_neur[ni-1])== count_i)[0]])
            paths_classes[ci].append(buff)
            count_i = list(buff.keys())[0]
    return paths_classes


# +
from collections import Counter
nreal = 10

edf = pd.DataFrame(columns=list(GE.columns),index=None)
paths_weighted_in_output = []

for I in range(nreal):

    x_train, x_test, y_train, y_test = train_test_split(X, dummyY, test_size=0.33,random_state=123)

    myModel = None
    myModel = tune_shallow(optimizer='adam' ,activation='relu',n_dimension=8)
    history = myModel.fit(x_train, y_train, epochs=50, batch_size=16,verbose = 0)
    scores = myModel.evaluate(x_test, y_test, verbose=0)
    print("Model accuracy: " + str(scores[1]*100))
    
    #most likely paths made on the basis of layer outputs
    cpaths = get_paths(4,dummyY,3,myModel)
    cpaths_most_likely = [[list(c.keys())[0] for c in c_] for c_ in cpaths]
    cpaths_0 = [''.join([str(i) for i in a]) for a in cpaths_most_likely]

    #find the most weighted neuron, for each of the 414 GE, in the first 'input' layer (414x8)
    tm_0 = [''.join([str(ii) for ii in np.argsort(i)[::-1][:1]]) for i in myModel.layers[0].get_weights()[0]]
    tm_0_counts = Counter(sorted(tm_0)) # make a count of the neuron vs how many times it was the heaviest

    #for each neuron above, find the first most and second most important neuron in the second layer (8x8)
    tm_01_0_ = [myModel.layers[1].get_weights()[0][int(j)] for j in tm_0_counts  ] # storing the (1x8) 8 times
    tm_01_0 = [np.argsort(i)[::-1][0] for i in tm_01_0_] # getting the most heavy neuron in the second layer
    tm_01_1 = [np.argsort(i)[::-1][1] for i in tm_01_0_] # getting the second most heavy neuron in the second layer

    #for each neuron above (i.e. in the heaviest in the first layer), find the first most and second most 
    # important neuron in the third layer (8x4)
    tm_12_0_ =[myModel.layers[2].get_weights()[0][int(j)] for j in tm_01_0]
    tm_12_0 = [np.argsort(i)[::-1][0] for i in tm_12_0_]
    tm_12_1 = [np.argsort(i)[::-1][1] for i in tm_12_0_]
    
    # create the first most likely path
    paths_0 = [str(list(tm_0_counts.keys())[i]) + str(tm_01_0[i]) + str(tm_12_0[i]) for i in range(len(tm_0_counts))]
    # create the second most likely path
    paths_1 = [str(list(tm_0_counts.keys())[i]) + str(tm_01_1[i]) + str(tm_12_1[i]) for i in range(len(tm_0_counts))]

#     tm_1 = [''.join([str(ii) for ii in np.argsort(i)[::-1][:2]]) for i in myModel.layers[1].get_weights()[0]]
#     tm_2 = [''.join([str(ii) for ii in np.argsort(i)[::-1][:2]]) for i in myModel.layers[2].get_weights()[0]]

#     tm_23_0 = [str(np.argsort(myModel.layers[2].get_weights()[0][int(i[0])])[::-1][0]) for i in tm_1]
#     tm_23_1 = [str(np.argsort(myModel.layers[2].get_weights()[0][int(i[0])])[::-1][1]) for i in tm_1]

#     tm_12_0 = [str(i)+tm_1[i][0]+tm_23_0[i] for i in range(len(tm_1))]
#     tm_12_1 = [str(i)+tm_1[i][1]+tm_23_1[i] for i in range(len(tm_1))]

    paths_weighted_in_output.append([(len(cpaths_0) - len(set(cpaths_0) - set(paths_0)))/len(cpaths_0),
                                     (len(cpaths_0) - len(set(cpaths_0) - set(paths_1)))/len(cpaths_0)])
    
    paths = []
    for neuron in myModel.layers[0].get_weights()[0]:
        paths.append([(np.argmax(neuron),np.argmax(layer.get_weights()[0][np.argmax(neuron)])) for layer in myModel.layers[1:-1]][0])

    ge_paths_binary = {GE.columns[i]:paths[i] for i in range(len(GE.columns))}

    ge_feature_dict = {GE.columns[i]:np.argmax(myModel.layers[0].get_weights()[0][i]) for i in range(len(GE.columns))}

    ge_neuron_dict = {str(i):[key for key,val in ge_feature_dict.items() if val == i] for i in range(8)}

    ge_neuron_dict_out = {str(i):ge_neuron_dict[cpaths_0[i][0]] for i in range(4)}
    
    edf = edf.append(pd.DataFrame(-1,columns=list(GE.columns),index=[I]))
    for oi in list(ge_neuron_dict_out.keys()):
        edf.loc[I,ge_neuron_dict_out[oi]] = int(oi) 
# -

plt.plot(np.array(paths_weighted_in_output)[:,0],label='O(0)')
plt.plot(np.array(paths_weighted_in_output)[:,1],label='O(1)')
plt.plot(np.array(paths_weighted_in_output)[:,0]+np.array(paths_weighted_in_output)[:,1],label='O(0+1)')
plt.legend()

yh,xh=np.histogram(np.array(paths_weighted_in_output)[:,0]+np.array(paths_weighted_in_output)[:,1])
plt.bar([(xh[i]+xh[i+1])*0.5*100 for i in range(len(xh)-1)],yh*100/nreal,width=5)
plt.title('% of paths using layer weights O(1,2)\nrecovered in paths using layer output')
plt.ylabel('% of cases')
plt.xlabel('% of recovery')

# +
sns.set(style='whitegrid')
fig,ax = plt.subplots(1,1,figsize=(9,35))
col_class =[ 'royalblue','tomato','gold','green']
cc= 0
yl = []
threshold = 2
gene_dict = {0:[],1:[],2:[],3:[]}

for ci in range(4):
    buff_df = edf[edf==ci].dropna(how='all',axis=1).count()
    buff_df=buff_df[buff_df>=threshold].sort_values()
    
    for bdi,buff_df_i in enumerate(buff_df):
        gene_dict[ci].append(buff_df.index[bdi])
        plt.plot([0,buff_df_i*10],[cc,cc],color=col_class[ci])
        plt.annotate(buff_df.index[bdi],xy=[buff_df_i*10 +2*(cc%2),cc],size=5,xycoords='data')

        cc+=1
    plt.annotate(classes[ci],xy=[-3,cc-20],xycoords='data',rotation=90,size=15)
    cc+=2
    yl.extend(buff_df.index)
ax.set_yticks([])
ax.set_yticklabels('')
plt.xlim(-5,110)
plt.ylim(-5,len(yl)+15)
plt.xlabel('% recovery at a neuron\nbelonging to a given class')
plt.savefig('GE_dist.png',dpi=300);

# +
mat = [8,8,4]

mat_imp = [[],[],[]]

sns.set(style='dark',context='talk')
plt.figure(figsize=(14,9))

x_pts = np.arange(-0.5,len(mat)+0.5)

for mi in range(len(mat)):
    print(mi)
    if len(mat_imp[mi])>0:
        plt.scatter(np.zeros(mat[mi])+mi, np.arange(mat[mi]) + (np.max(mat) - mat[mi])*0.5,s=1200*np.array(mat_imp[mi]))
    else:
        plt.scatter(np.zeros(mat[mi])+mi, np.arange(mat[mi]) + (np.max(mat) - mat[mi])*0.5,s=250)

# for i,pi in enumerate(np.unique(tm_12_0)):
#     for mi in range(len(pi)-1):
#         plt.plot([mi,mi+1],[int(pi[mi])+ (np.max(mat) - mat[mi])*0.5,int(pi[mi+1]) + (np.max(mat) - mat[mi+1])*0.5],c='k',zorder=-1)

# for i,pi in enumerate(np.unique(tm_12_1)):
#     for mi in range(len(pi)-1):
#         plt.plot([mi,mi+1],[int(pi[mi])+ (np.max(mat) - mat[mi])*0.5,int(pi[mi+1]) + (np.max(mat) - mat[mi+1])*0.5],c='k',zorder=-1)

for i,pi in enumerate(paths_0):
    for mi in range(len(pi)-1):
        plt.plot([mi,mi+1],[int(pi[mi])+ (np.max(mat) - mat[mi])*0.5,int(pi[mi+1]) + (np.max(mat) - mat[mi+1])*0.5],c='k',zorder=-1)

for i,pi in enumerate(paths_1):
    for mi in range(len(pi)-1):
        plt.plot([mi,mi+1],[int(pi[mi])+ (np.max(mat) - mat[mi])*0.5,int(pi[mi+1]) + (np.max(mat) - mat[mi+1])*0.5],c='grey',linestyle='dotted',zorder=-1)


for i,pi in enumerate(cpaths_most_likely):
    print(pi)
    for mi in range(len(pi)-1):
        plt.plot([mi,mi+1],[pi[mi]+ (np.max(mat) - mat[mi])*0.5,pi[mi+1]+ (np.max(mat) - mat[mi+1])*0.5],c='gold',lw=5,alpha=0.5,zorder=-1)

for i in range(4):
    plt.annotate('class '+str(i),xy=[2+.2,i+ (np.max(mat) - np.min(mat))*0.5],
                 xycoords='data')
        

plt.xlim(x_pts[0],x_pts[-1])
plt.ylim(-1,np.max(mat))

plt.title('two most preferred paths')
# plt.annotate(str(cmat),xy=(1.5,7), xycoords='data')
plt.axis('off')
plt.tight_layout()
plt.savefig('GE_most_pref_paths.png',dpi=250)


StarcoderdataJupyterScriptsDedupFiltered