<jupyter_start><jupyter_text>Simulated data LSTM test  - [1] time series (n=4) of simulated data from Figure 3 of " Robust and Scalable Models of Microbiome Dynamics for Bacteriotherapy Design " of <NAME>, <NAME> - [2] time series (n=13) of simulated data from Figure 4 of "Robust and Scalable Models of Microbiome Dynamics" of <NAME>, <NAME><jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from scipy.interpolate import UnivariateSpline as spline
from scipy.cluster.hierarchy import fcluster
from scipy import stats
import scipy.cluster.hierarchy as hac
from clstm import cLSTM, train_model_gista
import torch<jupyter_output><empty_output><jupyter_text>[1] time series (n=4) of simulated data from Figure 3 of " Robust and Scalable Models of Microbiome Dynamics for Bacteriotherapy Design " of <NAME>, <NAME>[link to image](https://ai2-s2-public.s3.amazonaws.com/figures/2017-08-08/bd0fe950c9b6e25edc80c60685ae0047614b0b93/4-Figure3-1.png)<jupyter_code>from IPython.display import Image
from IPython.core.display import HTML 
Image(url= "https://ai2-s2-public.s3.amazonaws.com/figures/2017-08-08/bd0fe950c9b6e25edc80c60685ae0047614b0b93/4-Figure3-1.png")
df={}
for i in range(1,5):
  df[i]=pd.read_csv('4_series_'+str(i)+'.csv',  names=['x', 'y'], header=None ) 
df[i].head()
color = {
    1:  (  28/255,  130/255,  197/255),
    2:  ( 138/255,   68/255,  153/255),
    3:  ( 238/255,  180/255,  43/255),
    4:  ( 218/255,   87/255,  32/255)
}

for i in range(1,5):
  plt.plot(df[i].x,df[i].y, color=color[i])
plt.show()
i=3
f = spline(df[i].x, df[i].y, k = 5)
x_int = np.linspace(0,601, 601)
y_int = f(x_int)
plt.plot(df[i].x,df[i].y, 'o', color=color[i])
plt.plot(x_int ,y_int , color=color[i])
plt.show()
Y=np.zeros((4,601))
x_int = np.linspace(0,601, 601)

for i in range(1,5):
  f = spline(df[i].x, df[i].y, k = 5)
  Y[i-1,:] = f(x_int)
for i in range(1,5):
  plt.plot(Y[i-1,:], color=color[i])
for i in range(1,5):
  plt.plot(df[i].x,df[i].y, 'o',color=color[i])
plt.show()
def GC_adjust(lstm, threshold=True):
        '''Extract learned Granger causality.
        Args:
          threshold: return norm of weights, or whether norm is nonzero.
        Returns:
          GC: (p x p) matrix. Entry (i, j) indicates whether variable j is
            Granger causal of variable i.
        '''
        GC = [torch.sum(net.lstm.weight_ih_l0, dim=0)
              for net in lstm.networks]
        GC = torch.stack(GC)
        if threshold:
            return (GC > 0).int().cpu().data.numpy()
        else:
            return GC.cpu().data.numpy()
          
def plot_GC(lstm):
  f, ax = plt.subplots(1,2,figsize=(8, 4))
  GC= lstm.GC(threshold=False).cpu().data.numpy()
  im = ax[0].imshow(GC, cmap='Blues');
  f.colorbar(im, ax=ax[0]);

  GC= GC_adjust(lstm, threshold=False)
  im = ax[1].imshow(GC, cmap='coolwarm');
  f.colorbar(im, ax=ax[1]);
  f.tight_layout()
# 13 series
p = 4
n_hidden = 10

gcmodel = cLSTM(p, n_hidden)

top_5_reads = Y.T
X_torch = torch.tensor(top_5_reads[np.newaxis],  dtype=torch.float32)
print(X_torch.shape)

# Train with GISTA
check_every = 50
lam=0.1
lam_ridge=1e-4
lr=0.005
max_iter=50
truncation= 1

train_loss_list, train_mse_list = train_model_gista(
    gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
    check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
plot_GC(gcmodel)
for i in range(0,5):
  train_loss_list, train_mse_list = train_model_gista(
      gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
      check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
  plot_GC(gcmodel)<jupyter_output>----------Iter = 50----------
Total loss = 1.793022
MSE = 1.327402, Ridge = 0.002533, Nonsmooth = 0.463087
Variable usage = 100.00%
----------Iter = 50----------
Total loss = 1.075701
MSE = 0.619965, Ridge = 0.002860, Nonsmooth = 0.452876
Variable usage = 100.00%
----------Iter = 50----------
Total loss = 0.826512
MSE = 0.380146, Ridge = 0.003017, Nonsmooth = 0.443350
Variable usage = 100.00%
----------Iter = 50----------
Total loss = 0.731760
MSE = 0.294599, Ridge = 0.003107, Nonsmooth = 0.434054
Variable usage = 100.00%
----------Iter = 50----------
Total loss = 0.669858
MSE = 0.241832, Ridge = 0.003167, Nonsmooth = 0.424860
Variable usage = 100.00%<jupyter_text>The matrix should look like this:<jupyter_code>x=np.array([
    [-1,  0 , 0 , 1],
    [0 , -4 ,-3 , 0],
    [0 ,  2 ,-4 , 1],
    [1 ,  0 , 2 ,-3],
])
plt.imshow(x, cmap='coolwarm');
plt.colorbar();<jupyter_output><empty_output><jupyter_text>- [2] time series (n=13) of simulated data from Figure 4 of "Robust and Scalable Models of Microbiome Dynamics" of <NAME>, <NAME>[link to image](https://ai2-s2-public.s3.amazonaws.com/figures/2017-08-08/af60e575721d3ef4bf2cf749f1a6024c2a3e4c8c/6-Figure4-1.png)<jupyter_code>Image(url= "https://ai2-s2-public.s3.amazonaws.com/figures/2017-08-08/af60e575721d3ef4bf2cf749f1a6024c2a3e4c8c/6-Figure4-1.png")
df={}
for i in range(1,14):
  df[i]=pd.read_csv(str(i)+'.csv',  names=['x', 'y'], header=None ) 
df[i].head()
color = {
    1:  (  0/255, 114/255, 189/255),
    2:  (126/255,  48/255, 141/255),
    3:  (217/255,  82/255,  25/255),
    4:  (162/255,  18/255,  48/255),
    5:  (119/255, 172/255,  48/255),
    6:  ( 78/255, 190/255, 238/255),
    7:  (238/255, 177/255,  32/255),
    8:  (119/255, 172/255,  48/255),
    9:  (217/255,  82/255,  25/255),
    10: (238/255, 177/255,  32/255),
    11: (  0/255, 114/255, 189/255),
    12: (126/255,  48/255, 141/255),
    13: ( 78/255, 190/255, 238/255)
}

for i in range(1,14):
  plt.plot(df[i].x,df[i].y, color=color[i])
plt.show()
i=6
f = spline(df[i].x, df[i].y, k = 5)
x_int = np.linspace(0,101, 101)
y_int = f(x_int)
plt.plot(df[i].x,df[i].y, color=color[i])
plt.plot(x_int ,y_int )
plt.show()
Y=np.zeros((13,101))
x_int = np.linspace(0,100, 101)

for i in range(1,14):
  f = spline(df[i].x, df[i].y, k = 5)
  Y[i-1,:] = f(x_int)
for i in range(1,14):
  plt.plot(Y[i-1,:], color=color[i])
plt.show()
def plot_dendrogram(Z):
  # Plot dendogram
  plt.figure(figsize=(15, 5))
  plt.title('Hierarchical Clustering Dendrogram')
  plt.xlabel('sample index')
  plt.ylabel('distance')
  hac.dendrogram(
      Z,
      leaf_rotation=90.,  # rotates the x axis labels
      leaf_font_size=8.,  # font size for the x axis labels
  )
  plt.show()
def print_clusters(timeSeries, Z, k, plot=False):
    # k Number of clusters I'd like to extract
    results = fcluster(Z, k, criterion='maxclust')

    # check the results
    s = pd.Series(results)
    clusters = s.unique()

    for c in clusters:
        cluster_indeces = s[s==c].index
        print("Cluster %d number of entries %d" % (c, len(cluster_indeces)))
        if plot:
          for ind in cluster_indeces: 
            plt.plot(timeSeries.iloc[ind,:], color=color[ind+1], label=str(ind+1))
          plt.legend()
          plt.show()
# Do the clustering
timeSeries = pd.DataFrame(Y)
Z = hac.linkage(timeSeries, method='single', metric='correlation')

plot_dendrogram(Z)
print_clusters(timeSeries, Z, 3, plot=False)

print('\nLets print the different clusters:')
print_clusters(timeSeries, Z, 3, plot=True)<jupyter_output><empty_output><jupyter_text>**Changing params method='complete' to get rid of ligne 9 in cluster 2 (to get the same clusters as in paper)**<jupyter_code>timeSeries = pd.DataFrame(Y)
Z = hac.linkage(timeSeries, method='complete')

plot_dendrogram(Z)
print_clusters(timeSeries, Z, 3, plot=False)

print('\nLets print the different clusters:')
print_clusters(timeSeries, Z, 3, plot=True)
def GC_adjust(lstm, threshold=True):
        '''Extract learned Granger causality.
        Args:
          threshold: return norm of weights, or whether norm is nonzero.
        Returns:
          GC: (p x p) matrix. Entry (i, j) indicates whether variable j is
            Granger causal of variable i.
        '''
        GC = [torch.sum(net.lstm.weight_ih_l0, dim=0)
              for net in lstm.networks]
        GC = torch.stack(GC)
        if threshold:
            return (GC > 0).int().cpu().data.numpy()
        else:
            return GC.cpu().data.numpy()
          
def plot_GC(lstm):
  f, ax = plt.subplots(1,2,figsize=(8, 4))
  GC= lstm.GC(threshold=False).cpu().data.numpy()
  im = ax[0].imshow(GC, cmap='Blues');
  f.colorbar(im, ax=ax[0]);

  GC= GC_adjust(lstm, threshold=False)
  im = ax[1].imshow(GC, cmap='coolwarm');
  f.colorbar(im, ax=ax[1]);
  f.tight_layout()
# 13 series
p = 13
n_hidden = 3

gcmodel = cLSTM(p, n_hidden)

top_5_reads = Y.T
X_torch = torch.tensor(top_5_reads[np.newaxis],  dtype=torch.float32)
print(X_torch.shape)

# Train with GISTA
check_every = 50
lam=0.5
lam_ridge=1e-4
lr=0.005
max_iter=100
monotone=True
truncation= 10
switch_tol=10e-3

train_loss_list, train_mse_list = train_model_gista(
    gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
    check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
plot_GC(gcmodel)
for i in range(0,20):
  train_loss_list, train_mse_list = train_model_gista(
      gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
      check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
  plot_GC(gcmodel)
# 13 series
p = 13
n_hidden = 3

gcmodel = cLSTM(p, n_hidden)
top_5_reads = Y.T
X_torch = torch.tensor(top_5_reads[np.newaxis],  dtype=torch.float32)
print(X_torch.shape)

# Train with GISTA
check_every = 25
lam=3
lam_ridge=1e-4
lr=0.005
max_iter=50
monotone=True
truncation= 10
switch_tol=10e-3

train_loss_list, train_mse_list = train_model_gista(
    gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
    check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
plot_GC(gcmodel)
for i in range(0,3):
  train_loss_list, train_mse_list = train_model_gista(
      gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
      check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
  plot_GC(gcmodel)<jupyter_output>----------Iter = 25----------
Total loss = 10.829184
MSE = 4.663805, Ridge = 0.000915, Nonsmooth = 6.164464
Variable usage = 76.92%
----------Iter = 50----------
Total loss = 10.161051
MSE = 3.270681, Ridge = 0.001099, Nonsmooth = 6.889271
Variable usage = 47.93%
----------Iter = 25----------
Total loss = 7.018055
MSE = 2.367607, Ridge = 0.001226, Nonsmooth = 4.649222
Variable usage = 39.05%
----------Iter = 50----------
Total loss = 6.425457
MSE = 2.027633, Ridge = 0.001324, Nonsmooth = 4.396500
Variable usage = 39.05%
----------Iter = 25----------
Total loss = 10.100224
MSE = 4.317703, Ridge = 0.001380, Nonsmooth = 5.781140
Variable usage = 37.87%
Switching to line search<jupyter_text>The matrix should look like one of these two:<jupyter_code>x=np.array([
    [1,1,1,1,1,0,0,0,0,0,0,0,0],
    [1,1,1,1,1,0,0,0,0,0,0,0,0],
    [1,1,1,1,1,0,0,0,0,0,0,0,0],
    [1,1,1,1,1,0,0,0,0,0,0,0,0],
    [1,1,1,1,1,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,0,0,0,0],
    [0,0,0,0,0,1,1,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,1,1,1,1,1,1],
    [0,0,0,0,0,0,0,1,1,1,1,1,1],
    [0,0,0,0,0,0,0,1,1,1,1,1,1],
    [0,0,0,0,0,0,0,1,1,1,1,1,1],
    [0,0,0,0,0,0,0,1,1,1,1,1,1],
    [0,0,0,0,0,0,0,1,1,1,1,1,1],
])
plt.imshow(x, cmap='coolwarm');
plt.colorbar();
x=np.array([
    [0,0,0,0,0,-1,-1,3,3,3,3,3,3],
    [0,0,0,0,0,-1,-1,3,3,3,3,3,3],
    [0,0,0,0,0,-1,-1,3,3,3,3,3,3],
    [0,0,0,0,0,-1,-1,3,3,3,3,3,3],
    [0,0,0,0,0,-1,-1,3,3,3,3,3,3],
    [2,2,2,2,2,0,0,-4,-4,-4,-4,-4,-4],
    [2,2,2,2,2,0,0,-4,-4,-4,-4,-4,-4],
    [0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0],
    [0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0],
    [0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0],
    [0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0],
    [0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0],
    [0,0,0,0,0,0,0, 0, 0, 0, 0, 0, 0],
])
plt.imshow(x, cmap='coolwarm');
plt.colorbar();
# 3 series: collapsing all the 13 into 3 by just taking 1 from every cluster
p = 3
n_hidden = 3

gcmodel = cLSTM(p, n_hidden)

top_5_reads = Y[[0,5,10],:].T
X_torch = torch.tensor(top_5_reads[np.newaxis],  dtype=torch.float32)
print(X_torch.shape)

# Train with GISTA
check_every = 100
lam=0.1
lam_ridge=10e-4
lr=0.005
max_iter=200
monotone=True
truncation= 1
switch_tol=10e-3

for i in range(0,2):
  train_loss_list, train_mse_list = train_model_gista(
      gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
      check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
  plot_GC(gcmodel)
for i in range(0,10):
  train_loss_list, train_mse_list = train_model_gista(
      gcmodel, X_torch, lam=lam, lam_ridge=lam_ridge, lr=lr, max_iter=max_iter, 
      check_every=check_every, monotone=monotone , truncation=truncation, switch_tol=switch_tol)
  plot_GC(gcmodel)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18039245982694685, 'fraction_numerical': 0.10800370828182942, 'mean_word_length': 3.741758241758242, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 33, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '368053', 'license': 'MIT', 'n_tokens_mistral': 6393, 'n_tokens_neox': 5617, 'n_words': 1123, 'path': '3_Simulated_Data.ipynb', 'repo_name': 'dylanrandle/microbiome', 'size': '13020'}