<jupyter_start><jupyter_text>Nystrom Tica aka Spectral-OASIS tutorial Nystrom Tica is based on the Nyström matrix operation theory, which can approximately reconstruct the time-lagged covariance matrix of all input features while using only a subset of features as input. Given an initial input feature set, Spectral-oASIS samples a subset of these features that best reconstructs the leading eigenfunctions of the time lagged covariance matrix obtained from MD simulationsThe end result of Nystrom tica is the same as Sparse Tica but nystrom tica uses different math to get thereNystrom tica is implemented in PyEmma but it is missing from the documentationsee origonal paper here:<NAME>.; <NAME>.; <NAME>.; <NAME>.; <NAME>.; <NAME>.; <NAME>.; <NAME>   Rapid Calculation of Molecular Kinetics Using Compressed Sensing   J. Chem. Theory Comput. 2018, 14 (5), 2771−2783   https://pubs.acs.org/doi/10.1021/acs.jctc.8b00089?ref=PDFX<NAME> spoke very highly of Nystrom tica. See how his lab uses it in their MSM workflow here:<NAME>, <NAME>, <NAME>, <NAME>, and X<NAME>Markov State Models to Study the Functional Dynamics of Proteins in the Wake of Machine Learning   JACS Au 2021   https://doi.org/10.1021/jacsau.1c00254 For this tutorial I will use the D.E Shaw main protease trajectory as a test system. All pairwize alpha carbon distances would be ~ 48,000 distances. Far to many to fit into memory.   I will first load a set of all phi and psi dihedrals (1220 total) and see if I can reduce the number of distances while still keeping a high implied timescale   Second I will see how many distances I can through at Nystrom tica before it crashes due to memory proplems   For more details on the test system see:   https://github.com/jgpattis/Desres-sars-cov-2-apo-mpro<jupyter_code>import pyemma.coordinates as coor
import numpy as np
import pandas as pd
import time 
import matplotlib.pyplot as plt
%matplotlib inline
# load distance data I had already calculated and saved

data_path = '/media/tun37047/HD/tun37047/covid19/desres-protease/monomer_global_apo_msm'

data_a = coor.load(f'{data_path}/feature_data_02/backbone_chain_0.h5')
data_b = coor.load(f'{data_path}/feature_data_02/backbone_chain_1.h5')
data = data_a + data_b<jupyter_output><empty_output><jupyter_text>The benifit of Nystrom tica is you let tica do the hard work of deciding which features are important   I would recomend more general featurizations as input. Depending on your system size things like all pairwize heavy atom distances, all pairwize alpha carbon and beta carbon distances, or all all pairwize alpha carbon distances   for examples of these featurizations see:   https://github.com/vvoelz/msm-best-practices/wiki/03-Featurization   or the PyEmma tutorial<jupyter_code># Here I am workin with 200 trajectories
len(data)
# Each trajectory is 1 microsecond long (1000 frames) with
#   1 ns between frames and 1220 features(sin and cos products of backbone phi and psi dihedral angles)
data[0].shape<jupyter_output><empty_output><jupyter_text>Nystrom tica is not in the PyEmma documentation   so to look at the inputs and outputs look at the PyEmma github:   https://github.com/markovmodel/PyEMMA/blob/devel/pyemma/coordinates/transform/nystroem_tica.pyL51   or alternativly print the doc string<jupyter_code># The docstring is long so use prnt to format it nice
print(coor.tica_nystroem.__doc__)<jupyter_output>Sparse sampling implementation [1]_ of time-lagged independent component analysis (TICA) [2]_, [3]_, [4]_.

    Parameters
    ----------
    max_columns : int
        Maximum number of columns (features) to use in the approximation.
    data : ndarray (T, d) or list of ndarray (T_i, d) or a reader created by
        source function array with the data. With it, the TICA
        transformation is immediately computed and can be used to transform data.
    lag : int, optional, default 10
        lag time
    dim : int, optional, default -1
        Maximum number of significant independent components to use to reduce dimension of input data. -1 means
        all numerically available dimensions (see epsilon) will be used unless reduced by var_cutoff.
        Setting dim to a positive value is exclusive with var_cutoff.
    var_cutoff : float in the range [0,1], optional, default 0.95
        Determines the number of output dimensions by including dimensions until their cumulative kineti[...]<jupyter_text>These options are very similar to regular tica except:   max_columns: How many features to use. I would recomend making an informed decision by plotting different options vs implied timescale (shown below)   nsel: How many columns to scan through at once. The origonal paper uses 20. The lower nsel is the better the tica approximation should be but the longer the calculation takes. As seen below larger numbers here work well<jupyter_code>start = time.time()
Ntica = coor.tica_nystroem(max_columns=200, data=data, lag=15, nsel=10)
end = time.time()
print(f'This took {np.round(end-start)} seconds')
start = time.time()
Ntica_20 = coor.tica_nystroem(max_columns=200, data=data, lag=15, nsel=20)
end = time.time()
print(f'This took {np.round(end-start)} seconds')
start = time.time()
Ntica_30 = coor.tica_nystroem(max_columns=200, data=data, lag=15, nsel=30)
end = time.time()
print(f'This took {np.round(end-start)} seconds')
fig, ax = plt.subplots()
nsel_list = [10, 20, 30]
tica_trials = [Ntica, Ntica_20, Ntica_30]
for i in range(5): # lets look at top 5 timescales
    temp = [k.timescales[i] for k in tica_trials]
    ax.plot(nsel_list, temp, label=f'ts {i + 1}')
ax.set_xlabel('nsel option', fontsize=16)
ax.set_ylabel('Timescale (ns)', fontsize=16)<jupyter_output><empty_output><jupyter_text>in the above, implied timescales go up with increasing nsel, which is the opposite of what I was expecting, and I believe the opposite of what it should be.   I believe this means that would could keep increasing nsel for a nice speedup of the calculation without loss of accuracyI will learn later that larger values of max_columns also slows down the calculation<jupyter_code># Larger numbers of nsel affects the resolution of max_columns
print(f'Features in Ntica, Ntica_20, Ntica_30: {len(Ntica.column_indices)} {len(Ntica_20.column_indices)} {len(Ntica_30.column_indices)}')
# Ok now lets look at the tica components
Ntica_30.get_output()<jupyter_output><empty_output><jupyter_text>So Nystrom tica is refusing to give me the tica components it calculated. I tried lots of different options and was unable to figure out what was going wrong. This is really bad but I will show a work around later   Next lets make a plot of implied timescales vs max_columns to make an informed decision on how sparce we can go<jupyter_code># I would recommend exponetial scaling for column list
column_list = [10, 25, 50, 75, 100, 150, 200, 250, 300, 400, 500, 700, 900]
out_list = []
true_column_list = []
start = time.time()

for i in column_list:
    temp_ntica = coor.tica_nystroem(max_columns=i, data=data, lag=15, nsel=30)
    out_list.append(temp_ntica.timescales)
    true_column_list.append(len(temp_ntica.column_indices))
end = time.time()
print(f'This took {np.round(end-start)} seconds')
reg_tica = coor.tica(data, lag=15)
out_list.append(reg_tica.timescales)
column_list.append(1220) # full number of features
true_column_list.append(1220)
# The true number of columns is not what I input for max_columns
# instead it is around max_columns
true_column_list
fig, ax = plt.subplots()
for i in range(5):  # plot the top i number of timescales
    temp = [out_list[k][i] for k in range(len(out_list))]
    ax.plot(true_column_list, temp, label=f'ts {i +1}')
ax.legend()
ax.set_xlabel('Cloumns (features) in tica output', fontsize=16)
ax.set_ylabel('Timescale (ns)', fontsize=16)<jupyter_output><empty_output><jupyter_text>In the apove plot implied timescales increase with increasing numbers of features. From this we can make a decison on how much 'information' we are willing to loose to make tica simpliar and easier to interpret. These two things area tradeoff one for the other.   In this plot I personally see two good choices:721 features: (aka the free lunch) at 721 features there is a very tiny decrease in timescales from regular (full) tica. With this option we would reduce the number of features by 43% while loosing very little kinetic information (we still describe the motions in the main protease very well). This is what the theory paper suggests a optimal, that you look for the asymptote of the curve   421 features: (aka the compromize) at 421 features there is a signifigent decrease in timescae of the first three tICs, but it's not quite as dramatic as lower than 400 features. This is a happy medium where we are not describing the motions quite as well as we could be but that allows for a 67% reduction in featuresLower that 421 features in my opinion is a bad option as too much kenitic information is lostSee the two publications I cite at the top of this tutorial for more examplesYour results may vary, but this plot helps you make an informed decision<jupyter_code># ok now lets figure out how to actually get tica components
# run nystrom tica with the optimum output of 700 features
start = time.time()
Ntica_700 = coor.tica_nystroem(max_columns=700, data=data, lag=15, nsel=30)
end = time.time()
print(f'This took {np.round(end-start)} seconds')
print(f'This took {np.round(end-start)} seconds')<jupyter_output>This took 82.0 seconds<jupyter_text>notice that moving from a max_columns of 200 to 700 slowed down the calculation from 22 seconds to 82 seconds<jupyter_code># this object contains a column_indices attribute which gives
#  the index of the columns retained
Ntica_700.column_indices.shape
# interesting that it isn't exactly 700??
# Maybe I should be plotting against column_indices
# but anyway, lets now filter my origonal dataset
Ntica_700.column_indices
# go through your list of trajectories and filter each one by our new set of features
filtered_data = [i[:,Ntica_700.column_indices] for i in data]
data[0].shape
# the first dimentions is frames
# the second dimention is features
# check to make sure you filtered along the correct dimention
filtered_data[0].shape
# I have not loaded my PyEmma featurizer object at the biginning because I did this earlier and saved it
# but you could filter your feature discription so you know what these features are

pdb = f'{data_path}/DESRES_protease_chainid.pdb'
''' Bachbone Phi and Psi torsion angles '''
featurizer = coor.featurizer(pdb)
featurizer.add_backbone_torsions(cossin=True, selstr=f'chainid == 0')
featurizer.describe()[:20]
new_discriptions = [featurizer.describe()[i] for i in Ntica_700.column_indices]
# now we know which features were included
new_discriptions[:20]
# now lets put the filtered dataset into regular tica

regular_tica_700 = coor.tica(filtered_data, lag=15, kinetic_map=False, commute_map=True, var_cutoff=0.85)
tica_output = regular_tica_700.get_output()
# Success
# you can now plot this and cluster this
tica_output
# let's make a nice plot because we can
from pyemma.plots import plot_density

cat_output = np.concatenate(tica_output)

fig, ax, cb = plot_density(cat_output[:,0], cat_output[:,1])
ax.set_xlabel('IC 1', fontsize=16)
ax.set_ylabel('IC 2', fontsize=16)<jupyter_output><empty_output><jupyter_text>Everything done so far could also be done with sparse tica (avialable in MSMBuilder)One downside of Tica and Sparse Tica is that they require N_features^3 memory. Sparse Tica makes the full covariance matrix then kicks features out to find a sparse solution. Nystrom tica builds the covariance matrix a little bit at a time and should in theory be able to handle more features.Regular tica will crash on this computer (with 128 GB of RAM) at ~9000 features. Let's see how high Nystrom tica can go  Warning!! The next section may crash your computer, save anything important, and procede at your own risk ... or better yet just look at my results and don't do this<jupyter_code># Let's use pairwize CA distances
traj_num = [f'{i:04d}' for i in range(100)] 
traj_path = '/../DESRES-Trajectory_sarscov2-10880334-no-water-no-ion-glueCA/sarscov2-10880334-no-water-no-ion-glueCA/sarscov2-10880334-no-water-no-ion-glueCA-'
traj_list = [ data_path + traj_path + str(i) + '.dcd' for i in traj_num]

feat = coor.featurizer(pdb)
feat.add_distances(feat.pairs(feat.select('name == CA and chainid == 0'), excluded_neighbors=3))
traj = coor.load(traj_list, feat, stride=5)
traj_cat = np.concatenate(traj)
# working with source also helps minimize RAM usage
# NOTE: I later find out that source does not imporve anyting
#  since Nystrom tica needs all data in memory.
traj_source = coor.source(traj_list, feat, stride=5)
# 46,000 distances !!!
traj[0].shape
# Let's start with 15,000
data_15 = [i[:,:15000] for i in traj]
# Let's see what happens
start = time.time()
Ntica_big = coor.tica_nystroem(max_columns=700, data=data_15, lag=15, nsel=100)
end = time.time()
print(f'This took {np.round(end-start)} seconds')
print(f'This took {np.round((end-start)/60)} minutes')
# Let's go for broke
# use source instead of load to save memory
start = time.time()
Ntica__really_big = coor.tica_nystroem(max_columns=2000, data=traj_source, lag=15, nsel=400)
end = time.time()
print(f'This took {np.round(end-start)} seconds')<jupyter_output>02-09-21 10:46:21 pyemma.coordinates.transform.nystroem_tica.NystroemTICA[245] WARNING  Every iteration of the selection process involves streaming of all data and featurization. Depending on your setup, this might be inefficient.<jupyter_text>Ok so the above ran for 24 hours without crashing so I think it will work I just didn't want to wain for it.Also note that PyEmma through a warning that Nystrom Tica needs to load all data into memory so using source didn't help save memeory<jupyter_code>len(Ntica_big.column_indices)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07695068552149882, 'fraction_numerical': 0.028856507070562056, 'mean_word_length': 5.046875, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 52, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '493729', 'license': 'MIT', 'n_tokens_mistral': 4642, 'n_tokens_neox': 4218, 'n_words': 1752, 'path': 'Nystrom_tica_tutorial.ipynb', 'repo_name': 'vvoelz/msm-best-practices', 'size': '14011'}