<jupyter_start><jupyter_text>Tutorial: Estimating Free energy landscapes and clustering with Gaussian mixture modelsThis tutorial shows how to perform density and free energy estimation with Gaussian mixture models (GMM) and model selection. It will also show how to cluster the estimated density with InfleCS clustering. Note that while the tutorial shows an example of a 2D dataset, the estimation is implemented to accommodate datasets of any dimensionality. The visualization part, however, only supports 1D and 2D at the moment. References Free energy estimation with Gaussian mixture modelsInference of Calmodulin’s Ca2+-Dependent Free Energy Landscapes via Gaussian Mixture Model Validation <NAME>, <NAME>, <NAME>, and <NAME>te Journal of Chemical Theory and Computation, 2018 DOI: 10.1021/acs.jctc.7b00346  Clustering with InfleCSInfleCS: Clustering Free Energy Landscapes with Gaussian Mixtures <NAME>, Lucie Delemotte Journal of Chemical Theory and Computation, 2019 DOI: 10.1021/acs.jctc.9b00454 1. Theory 1.1 Gaussian mixture modelsA Gaussian mixture density is a weighted sum of Gaussian components, $$\rho(x) = \sum\limits_{i=1}^{N_{\text{basis}}} \pi_i \mathcal{N}(x|\; \mu_i,\; \Sigma_i),$$where $\pi_i$ is the amplitude, $\mu_i$ is the mean and $\Sigma_i$ is the covariance of the $i$th Gaussian. The amplitudes, means and covariances are here optimized iteratively with Expectation-Maximization. The expectation-maximization algorithm consists of two steps, abbreviated E-step and the M-step. 1.1.1 E-stepIn the E-step, the responsibility of each data point with respect to the Gaussian components is calculated as$$ \gamma_{ni} = \frac{\pi_i \mathcal{N}(x_n|\; \mu_i,\; \Sigma_i)}{\sum\limits_{i=1}^{N_{\text{basis}}} \pi_i \mathcal{N}(x_n|\; \mu_i,\; \Sigma_i)}.$$ 1.1.2 M-stepThe parameters are then updated in the maximization step. The amplitudes are updated with$$\pi_i = \frac{\sum_n \gamma_{ni}}{\sum_{n,i}\gamma_{ni}},$$the means with$$\mu_i = \frac{\sum_n \gamma_{ni}x_n}{\sum_n \gamma_{ni}},$$and covariances with$$\Sigma_i = \frac{\sum_n \gamma_{ni}(x_n-\mu_i)(x_n-\mu_i)^T}{\sum_n \gamma_{ni}} .$$These two steps are carried out until convergence of the log-likelihood, $\mathcal{l} = \frac{1}{N}\sum_{n=1}^N \log (\rho(x_n))$.  1.2 Model selection To select the number of Gaussian components in the density model, we can use either $K$-fold cross-validation or Bayesian information criterion for identifying one model, or Akaike information criterion for weighting multiple models.  1.2.1 $K$-fold cross-validationIn $K$-fold cross-validation, the data is split into $K$ folds. One chunk of the data is used for training model parameters, while the other is used for validating the trained model. In this case, model validation occurs by calculating the log-likelihood of the validation set. In $K$-fold cross-validation, each fold is used as a validation set once per model. This results in $K$ training and validation sets for each model. The total validation log-likelihood score of a model is the summed validation log-likelihoods of all validation sets.  1.2.2 Bayesian information criterionThe second option is to perform model selection with Bayesian information criterion (BIC). In short, the log-likelihood of the fitted data will increase with increased number of parameters. The BIC tries to account for this by adding a penalty for the number of parameters,$$BIC = \log(N) N_{param} - 2\mathcal{l},$$where $\mathcal{l}$ is the log-likelihood of the data. 1.2.3 Weighting multiple GMMs (stacking landscapes)The final density model can also be a weighted average of all estimated models based on the Akaike information criterion. The Bayesian information criterion (BIC) is good for selecting one model (it does so with probability one). However, for weighting models, the Aikake information criterion (AIC) is better suited. AIC follows the same idea as BIC such that complex models are penalized. However, the AIC has a slightly different formula which makes it better suited for weighting models,$$AIC = 2N_{param}-2\mathcal{l},$$where $\mathcal{l}$ is the log-likelihood of the data. The weight of model $m$ with an AIC score of $a_m$ is then set to$$w_m = \frac{e^{-\frac{\Delta a_m}{2}}}{\sum \limits_{m=1}^{N_{models}} e^{-\frac{\Delta a_m}{2}}},$$where $\Delta a_m = a_m - min\big(\{a_{m'}\}_{m'=1}^{N_{models}}\big)$. 1.3 Clustering density maximaUsing the estimated Gaussian mixture density, we can identify core states at density maxima and divide the points into clusters. This is done by first computing the Hessian at each point and determine whether it belongs to a density maximum or not. Each point that belongs to a density maximum is then connected to other points which are part of the same maximum. This creates so called connected components. Once the connected components are constructed, cluster labels can be assigned so that all points within the same connected component have the same cluster label. Transition points (points that do not belong to a density maximum), are not assigned to clusters when identifying core states. However, they can be assigned a cluster label. This is done by sorting the transition points in descending density order and one-by-one assigning them to the same cluster as the closest point with a cluster label.  2. Reading the dataThe data should be on the form [n_samples x n_dimensions]. Each row contains one data sample projected onto a CV of *n_dimensions* dimensions. In this case, the dataset consists of three randomly generated blobs.<jupyter_code>import numpy as np
import free_energy_clustering as FEC
import toy_models

tm = toy_models.Blobs(n_dims=2)
data = tm.sample(4000)<jupyter_output>(4000, 2)<jupyter_text>3. Constructing the object and setting parametersCreate a free energy and clustering object. In principle, *data*, minimum number of components, *min_n_components*, and maximum number of components, *max_n_components* need to be specified. If *max_n_components* is not specified, the free energy estimation and clustering will be done using *min_n_components* components. The settings of the necessary parameters will be printed below the cell. If not supplied by the user when creating the object, the default values will be used. **The following parameters are use to pick a method for model selection:**1. *n_splits*: The number of folds used in cross-validation, $K$. If *n_splits = 1*, model selection will be done with Bayesian information criterion (BIC). 2. *stacking_landscapes*: when set to True, the final density model will be a weighted average of all estimated models, regardless of the value on *n_splits*.**Other parameters are:**1. *shuffle_data*: Whether or not to shuffle the data when performing cross-validation model selection.2. *n_iterations*: The number of times to reestimate the same model, starting from different random initializations. The model with largest log-likelihood is picked as final estimated model.3. *temperature*: Temperature used in the simulation [K].4. *n_grids*: The number of grid cells used along all dimensions when plotting the free energy landscape. Example: for a 2D dataset with n_grids=50, the estimated free energy will be evaluated and visualized on a [50 x 50] grid. Note that the does not affect the density and free energy estimation, just visualization.5. *convergence_tol*: The convergence tolerance (change in log-likelihood during training).6. *x_lims*: Limits on the grid used for visualization. If not supplied, it will default to min/max values.7. *n_components_step*: Example: if set to 3, and *min_n_components=8* and *max_n_components=15*, models with *n_components={8, 11, 14}* will be evaluated.8. *verbose*: Whether or not to print details and estimation progression.9. *data_weights:* Numpy array with the weight of each data point. Default=None (means that all data points have the same weight, as in e.g. unbiased MD).<jupyter_code>fec = FEC.FreeEnergyClustering(data, min_n_components=2, max_n_components=10, temperature=300.0, 
                            n_iterations=5, n_grids=80, n_splits=1,stack_landscapes=False)<jupyter_output>*----------------Gaussian mixture model free energy estimator----------------*
   n_splits = 1
   shuffle_data = False
   n_iterations = 5
   n_grids = 80
   covergence_tol = 0.0001
   stack_landscapes = False
   x_lims (axes limits) = [[-11.64890867921662, -0.6228674363086841], [2.4187153907573857, 12.690299490769473]]
   temperature = 300.0
   min_n_components = 2
   max_n_components = 10
   n_components_step = 1
   Using weighted data: False
*----------------------------------------------------------------------------*<jupyter_text>4. Estimating free energy landscapes> **coords, FE_landscape, FE_points = fec.landscape()**will estimate the density and free energy landscape with Gaussian mixture models. The axes limits of the landscape are printed above. **This returns:**1. *coords*: The coordinates of the grid that is used for visualization.2. *FE_landscape*: The estimated free energy of each grid point.3. *FE_points*: The estimated free energy of each point in *data*.It is also possible to estimate standard error by writing > **standard_error_FE = fec.standard_error()**This returns the standard error of the free energy landscape, evaluated on each grid point.Note: When estimating the standard error for a 1D landscape, the standard error will be shown in the visualization.<jupyter_code>coords, FE_landscape, FE_points = fec.landscape()<jupyter_output>Estimating density with GMM.
# Components = 2
# Components = 3
# Components = 4
# Components = 5
# Components = 6
# Components = 7
# Components = 8
# Components = 9
# Components = 10
Identifying final model with 3 components.
Evaluating density in landscape
Density grid shape: [80, 80]<jupyter_text>5. Visualizing the free energy landscapeIf the dataset is 1D or 2D, the free energy can be visualized. To visualize the estimated free energy landscape, we call the *visualize* function. To visualize with the standard parameter settings, it is just to run > **fec.visualize()** However, you might want to change some of the parameters. **The following parameters can be changed:**1. *title*: The title of the figure2. *fontsize*: Title fontsize, other labels are changed according to this3. *savefig*: True/False (whether to save a .png, .svg and .eps figure)4. *filename*: The figure filename (if saving figure). Default='free_energy_landscape'5. *xlabel*: The x-label6. *ylabel*: The y-label (in case of a 2D dataset)7. *vmax*: Colormap max cutoff value.8. *n_contour_levels*: Determines the number and positions of the contour lines / regions in contourf.9. *show_data*: True/False (whether to scatter-plot the original data and/or clustering)<jupyter_code>fec.visualize(savefig=False, show_data=False, vmax=7, xlabel='X', ylabel='Y',filename='free_energy_landscape',title='Free energy landscape')<jupyter_output><empty_output><jupyter_text>6. Cluster points at free energy minimaUsing the estimated Gaussian mixture density, we can identify core states at density maxima and divide the points into clusters. Clustering the free energy landscape can be done with  > **labels, cluster_centers = fec.cluster(data, FE_points, assign_transition_points=False)** Setting *assign_transition_points*=True will return the full clustering (not just core states). Recall that *FE_points* is an array with the estimated free energy of each point in *data*. To make the clustering slightly more robust, we can do it over the grid that we use for visualization and then return cluster indices of the *data* points,> **labels, cluster_centers = fec.cluster(grid_points, FE_points, data, assign_transition_points=False)**The grid coordinates are kept in *coords*. We also need to parse *data* as a third argument to return the cluster labels and cluster center indices of the *data* points. By default, we will get the core state clustering (core state points have cluster indices > 0, and transition points have cluster indices equal to zero). To assign cluster indices to the non-core state points (transition points), just set *assign_transition_points* to True. **This returns:**1. *labels*: The cluster label/index of each point, denoting which cluster the point belongs to. If zero, the point is a transition point, otherwise it belongs to some intermediate state.2. *cluster_centers*: Contains the indices corresponding to the point of each cluster that has the smallest free energy. In *return_center_coords* is set to True, this contains the coordinates of the cluster centers instead.When the clustering has already been performed, new points can be assigned cluster labels based on the precomputed clustering (assigning points to the nearest neighbor clusters) with> **labels = fec.evaluate_clustering(new_points, assign_transition_points=False)**Note: When "No point in cluster X" is printed, this means that no data point were found within the core state of the identified cluster. This could indicate overfitting of the density model.<jupyter_code>labels, cluster_centers = fec.cluster(coords, FE_points,data, assign_transition_points=False)
print("Cluster center indices: "+str(cluster_centers))

# Computing state populations
state_populations = fec.population_states(n_sampled_points=100000)<jupyter_output>Clustering free energy landscape...
Computing Hessians of density landscape.
Point: 6400/6400
Constructing connected components.
# Graph connections: 23098.0
Clustering data points.
Cluster labels: [0. 1. 2. 3.]
Done clustering.
Cluster center indices: [2466 3143 1848]
Sampling points from density.
Assigning cluster labels based on precomputed density model clustering.
Computing state populations.<jupyter_text>7. Visualizing the free energy landscape with identified clusters<jupyter_code># Visualize free energy landscape with cluster labels
fec.pathways_ = None
fec.visualize(savefig=False, vmax=7, show_data=True, xlabel='X', ylabel='Y', filename='InfleCS_clustering', title='InlfeCS clustering')

# Plotting the state populations
import matplotlib.pyplot as plt
plt.figure(figsize=[15,5]);
plt.plot(np.arange(1,state_populations.shape[0]),state_populations[1::]/state_populations.sum(),linewidth=5,color=[0.7,0.7,0.7],zorder=-1)
plt.scatter(np.arange(1,state_populations.shape[0]),state_populations[1::]/state_populations.sum(),s=500,c=np.arange(1,state_populations.shape[0]),cmap='jet',edgecolor='k')
plt.title('State populations',fontsize=30)
plt.xlabel('State',fontsize=28)
plt.ylabel('Probability',fontsize=28)
plt.xticks(np.arange(1,state_populations.shape[0]))<jupyter_output><empty_output><jupyter_text>Estimating pathways between states<jupyter_code>fec.pathways([cluster_centers[0]],[cluster_centers[1]],max_iter=300,convergence_tol=1e-4,step_size=2e-2,n_points=20)
fec.visualize(savefig=False, vmax=7, show_data=True, xlabel='X', ylabel='Y', filename='pathways', title='')<jupyter_output>Iteration: 300/300
{'fraction_non_alphanumeric': 0.0915356306457009, 'fraction_numerical': 0.021947139511542323, 'mean_word_length': 5.261416003351488, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 36, 'https://': 0, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '125494', 'license': 'MIT', 'n_tokens_mistral': 4335, 'n_tokens_neox': 3974, 'n_words': 1810, 'path': 'tutorial_free_energy_clustering.ipynb', 'repo_name': 'delemottelab/InfleCS-free-energy-clustering-tutorial', 'size': '14996'}