Model ND energy Gaussians¶

import os
import numpy as np
import pandas as pd
from scipy.stats import pearsonr, norm, normaltest, shapiro
import matplotlib.pyplot as plt

import NDutils as nd # my lib

Paci¶

10 protein chains that fold in a two-state manner (Paci dataset)
ND energy profile. EDV = 'a', NDV = 'cxary'
Normality of ND energy distribution at each $Q$.
- 'a' is the only ND variant with 0 not-normals $Q < 0.6$ at $\alpha = 0.005$.

fn = os.path.join('P_ND3', 'P_frates.txt')
df = nd.load_frates(fn)
df.head()

dirname = "P_ND3"
PIDs = df.pid.tolist(); num_pids = len(PIDs)
fold_rates = df.k_f.values
print(num_pids, PIDs)
print("Range of fold rates:", np.ptp(fold_rates))

10 ['1imq', '1lmb4', '1bf4', '2ptl', '2ci2', '1aps', '1fmk', '1bk2', '1shf', '1ten']
Range of fold rates: 4.328

# Get ND energy data
edv = 'a'
dirn = os.path.join(dirname, "edv_" + edv)
c_ndstats = nd.get_NDstats(dirn, PIDs, 'c', edv)
c_mje_params = nd.get_params(c_ndstats, 'sc_mje')
c_peakQs, c_peakEs = nd.get_peaks(c_mje_params, PIDs)
print("'c'", pearsonr(fold_rates, c_peakEs), np.mean(c_peakQs), np.std(c_peakQs), 
      nd.mjes_range(c_peakEs))

x_ndstats = nd.get_NDstats(dirn, PIDs, 'x', edv)
x_mje_params = nd.get_params(x_ndstats, 'sc_mje')
x_peakQs, x_peakEs = nd.get_peaks(x_mje_params, PIDs)
print("'x'", pearsonr(fold_rates, x_peakEs), np.mean(x_peakQs), np.std(x_peakQs), 
      nd.mjes_range(x_peakEs))

r_ndstats = nd.get_NDstats(dirn, PIDs, 'r', edv)
r_mje_params = nd.get_params(r_ndstats, 'sc_mje')
r_peakQs, r_peakEs = nd.get_peaks(r_mje_params, PIDs)
print("'r'", pearsonr(fold_rates, r_peakEs), np.mean(r_peakQs), np.std(r_peakQs),
      nd.mjes_range(r_peakEs))

a_ndstats = nd.get_NDstats(dirn, PIDs, 'a', edv)
a_mje_params = nd.get_params(a_ndstats, 'sc_mje')
a_peakQs, a_peakEs = nd.get_peaks(a_mje_params, PIDs)
print("'a'", pearsonr(fold_rates, a_peakEs), np.mean(a_peakQs), np.std(a_peakQs),
      nd.mjes_range(a_peakEs))

y_ndstats = nd.get_NDstats(dirn, PIDs, 'y', edv)
y_mje_params = nd.get_params(y_ndstats, 'sc_mje')
y_peakQs, y_peakEs = nd.get_peaks(y_mje_params, PIDs)
print("'y'", pearsonr(fold_rates, y_peakEs), np.mean(y_peakQs), np.std(y_peakQs),
      nd.mjes_range(y_peakEs))

'c' (-0.8326901471750437, 0.0027867362980185874) 0.26 0.04898979485566354 1.8711453949497967
'x' (-0.8955232013098765, 0.0004587141873736154) 0.25000000000000006 0.04999999999999999 2.2487040312412834
'r' (-0.05563325121177174, 0.8786782256482984) 0.19 0.03 1.4174191962184999
'a' (-0.8479417111836289, 0.0019386037169391614) 0.25000000000000006 0.04999999999999999 1.84833757867314
'y' (0.2567066594223849, 0.47401892205094687) 0.15000000000000005 0.05000000000000001 1.621999990733038

# Example content of dicts
c_mje_params['2ci2', 0.1]

(26.3204, 24.15196936236894)

c_ndstats['2ci2'].loc[0.1].sc_mje

z
0.1    36.43
0.1    47.72
0.1    54.19
0.1     7.26
0.1    49.53
       ...  
0.1    25.30
0.1    31.79
0.1    40.00
0.1   -10.48
0.1    10.47
Name: sc_mje, Length: 100, dtype: float64

Test normality of ND energy distribution at each Q.¶

Check sample data normality:

With Shapiro-Wilk scipy.stats.shapiro():
- Null hypothesis: the sample comes from a normal distribution
- Uses quantiles to compare two distributions, more specific than Kolmogorov-Smirnoff.
- Q-Q plot to look at the correlation between the percentile points from two distributions (one of which is normal).
- If p-value < $\alpha$, the null hypothesis can be rejected (data are unlikely normal).
With scipy.stats.normaltest():
- Null hypothesis: the sample comes from a normal distribution.
- Test statistic combines skew and kurtosis. p-value is 2-sided $\chi^2$ probability.
- If p-value < $\alpha$, the null hypothesis can be rejected (data are unlikely normal).

def check_normality(ndstats: dict, col: str, alpha: float, method='shapiro'):
    """Method is either 'shapiro' (default) or 'normaltest'."""
    pvals_dict = {}
    for pid in PIDs:
        if method == 'normaltest': # pvals are two-tailed, want the tail area upto statistic
            pvals = [ (normaltest(ndstats[pid].loc[q][[col]].values)[1]/2).item() 
                     for q in nd.Qs ]
        else:
            pvals = [ shapiro(ndstats[pid][[col]].loc[q].values)[1] for q in nd.Qs ]
        pvals_dict[pid] = pvals # np.round(pvals, 3)
    df = pd.DataFrame(pvals_dict, index = nd.Qs)
    db = df.apply(lambda x: x < alpha)
    null_rejects = np.sum(db.values, axis=0) # times null hypothesis rejected
    return df, db, null_rejects

# shapiro-wilk, p-value < 0.005, reject normality (True)
a_df, a_db, a_nr = check_normality(a_ndstats.copy(), "sc_mje", 0.005, "shapiro")
print(a_db, '\n')
print(a_nr, np.sum(a_nr))
print("not-normals in Q < 0.6", np.count_nonzero(a_db.loc[0.1:0.5].values), '\n')
print(np.round(a_df, 3))

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  False  False  False  False  False  False  False  False  False  False
0.2  False  False  False  False  False  False  False  False  False  False
0.3  False  False  False  False  False  False  False  False  False  False
0.4  False  False  False  False  False  False  False  False  False  False
0.5  False  False  False  False  False  False  False  False  False  False
0.6  False  False  False  False  False  False  False  False  False  False
0.7  False  False  False  False  False  False  False  False  False  False
0.8  False  False   True  False  False  False  False  False  False  False
0.9  False   True   True   True  False  False  False  False  False  False
1.0   True   True   True   True   True   True   True   True   True   True 

[1 2 3 2 1 1 1 1 1 1] 14
not-normals in Q < 0.6 0 

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  0.473  0.070  0.211  0.380  0.456  0.152  0.341  0.590  0.692  0.152
0.2  0.351  0.416  0.233  0.313  0.836  0.120  0.814  0.391  0.243  0.683
0.3  0.727  0.740  0.324  0.370  0.165  0.208  0.419  0.015  0.662  0.641
0.4  0.444  0.142  0.045  0.569  0.220  0.039  0.610  0.418  0.506  0.630
0.5  0.007  0.103  0.541  0.671  0.735  0.046  0.022  0.791  0.231  0.419
0.6  0.118  0.597  0.304  0.006  0.603  0.776  0.320  0.079  0.307  0.467
0.7  0.497  0.162  0.915  0.357  0.670  0.447  0.747  0.343  0.556  0.293
0.8  0.196  0.026  0.001  0.908  0.974  0.064  0.153  0.717  0.506  0.170
0.9  0.005  0.000  0.000  0.000  0.750  0.333  0.757  0.339  0.328  0.136
1.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000

# shapiro-wilk, p-value < 0.005, reject normality (True)
c_df, c_db, c_nr = check_normality(c_ndstats.copy(), "sc_mje", 0.005, "shapiro")
print(c_db, '\n')
print(c_nr, np.sum(c_nr))
print("not-normals in Q < 0.6", np.count_nonzero(c_db.loc[0.1:0.5].values), '\n')
print(np.round(c_df, 3))

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  False  False  False  False  False  False  False  False  False  False
0.2  False  False  False  False   True  False  False  False  False  False
0.3  False  False  False  False  False  False  False  False  False  False
0.4  False  False  False  False  False  False  False  False  False  False
0.5  False  False  False  False  False  False  False  False  False  False
0.6  False  False  False  False  False  False  False  False  False  False
0.7  False  False  False  False  False  False  False  False  False  False
0.8  False  False   True  False  False  False  False  False  False   True
0.9  False   True   True  False  False  False  False  False  False  False
1.0   True   True   True   True   True   True   True   True   True   True 

[1 2 3 1 2 1 1 1 1 2] 15
not-normals in Q < 0.6 1 

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  0.396  0.031  0.006  0.809  0.650  0.150  0.039  0.010  0.593  0.125
0.2  0.728  0.372  0.662  0.987  0.000  0.524  0.713  0.473  0.099  0.263
0.3  0.582  0.772  0.861  0.773  0.006  0.459  0.229  0.551  0.594  0.829
0.4  0.175  0.581  0.151  0.662  0.669  0.345  0.983  0.620  0.996  0.850
0.5  0.443  0.262  0.183  0.465  0.210  0.066  0.050  0.816  0.734  0.032
0.6  0.089  0.970  0.463  0.187  0.136  0.868  0.075  0.033  0.590  0.607
0.7  0.726  0.573  0.035  0.761  0.562  0.141  0.844  0.258  0.549  0.716
0.8  0.005  0.077  0.001  0.054  0.042  0.096  0.098  0.925  0.108  0.003
0.9  0.088  0.001  0.001  0.048  0.113  0.583  0.551  0.039  0.104  0.959
1.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000

print(c_df.loc[0.2][['2ci2']].item())

0.0001698280539130792

x_df, x_db, x_nr = check_normality(x_ndstats.copy(), "sc_mje", 0.005, "shapiro")
print(x_db, '\n')
print(x_nr, np.sum(x_nr))
print("not-normals in Q < 0.6", np.count_nonzero(x_db.loc[0.1:0.5].values), '\n')
print(np.round(x_df, 3))

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  False  False  False   True  False  False  False  False  False  False
0.2  False  False  False  False  False  False  False  False  False  False
0.3  False  False  False  False  False  False  False  False  False  False
0.4  False  False  False  False   True  False  False  False  False  False
0.5  False  False  False  False  False  False  False  False  False   True
0.6  False  False  False  False  False  False  False  False  False  False
0.7  False  False   True  False  False  False  False  False  False  False
0.8  False  False   True  False  False  False  False  False  False  False
0.9  False  False  False  False  False  False  False  False  False  False
1.0   True   True   True   True   True   True   True   True   True   True 

[1 1 3 2 2 1 1 1 1 2] 15
not-normals in Q < 0.6 3 

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  0.046  0.040  0.038  0.000  0.451  0.919  0.874  0.204  0.097  0.090
0.2  0.349  0.591  0.169  0.123  0.774  0.674  0.109  0.311  0.370  0.392
0.3  0.849  0.455  0.370  0.222  0.571  0.940  0.065  0.652  0.304  0.462
0.4  0.374  0.663  0.342  0.268  0.004  0.291  0.665  0.568  0.713  0.037
0.5  0.019  0.277  0.095  0.429  0.045  0.028  0.133  0.505  0.516  0.003
0.6  0.534  0.284  0.328  0.615  0.017  0.836  0.287  0.294  0.740  0.957
0.7  0.212  0.014  0.004  0.467  0.618  0.247  0.331  0.502  0.790  0.757
0.8  0.038  0.075  0.001  0.276  0.052  0.554  0.196  0.083  0.367  0.425
0.9  0.294  0.054  0.007  0.209  0.217  0.578  0.416  0.342  0.045  0.006
1.0  0.000  0.000  0.000  0.002  0.000  0.000  0.000  0.000  0.000  0.000

print(x_df.loc[0.1][['2ptl']].item())
print(x_df.loc[0.4][['2ci2']].item())
print(x_df.loc[0.5][['1ten']].item())

0.00043271438335068524
0.004232669249176979
0.0033112899400293827

r_df, r_db, r_nr = check_normality(r_ndstats.copy(), "sc_mje", 0.005, "shapiro")
print(r_db, '\n')
print(r_nr, np.sum(r_nr))
print("not-normals in Q < 0.6", np.count_nonzero(r_db.loc[0.1:0.5].values), '\n')
print(np.round(r_df, 3))

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  False  False  False  False  False  False  False  False  False  False
0.2  False  False  False  False  False  False  False  False  False  False
0.3  False  False  False  False  False  False  False  False  False  False
0.4   True  False  False  False  False  False  False  False  False  False
0.5  False  False  False  False  False  False  False  False  False  False
0.6  False  False  False  False  False  False  False  False  False  False
0.7  False  False  False  False  False  False  False  False  False  False
0.8  False  False  False  False  False  False  False  False  False  False
0.9  False  False  False   True  False  False  False  False  False  False
1.0   True   True   True   True   True   True   True   True   True   True 

[2 1 1 2 1 1 1 1 1 1] 12
not-normals in Q < 0.6 1 

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  0.554  0.250  0.008  0.026  0.134  0.080  0.022  0.552  0.746  0.623
0.2  0.840  0.110  0.879  0.016  0.038  0.851  0.048  0.471  0.955  0.054
0.3  0.078  0.101  0.416  0.089  0.007  0.807  0.471  0.183  0.635  0.796
0.4  0.004  0.047  0.532  0.952  0.072  0.255  0.730  0.581  0.722  0.723
0.5  0.110  0.632  0.469  0.601  0.871  0.118  0.893  0.882  0.787  0.518
0.6  0.039  0.263  0.332  0.574  0.166  0.335  0.915  0.310  0.874  0.661
0.7  0.434  0.930  0.243  0.870  0.257  0.840  0.126  0.957  0.681  0.139
0.8  0.018  0.009  0.691  0.097  0.603  0.605  0.006  0.272  0.339  0.255
0.9  0.403  0.423  0.008  0.003  0.268  0.829  0.365  0.182  0.169  0.921
1.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000

print(r_df.loc[0.4][['1imq']].item())

0.0037239319644868374

y_df, y_db, y_nr = check_normality(y_ndstats.copy(), "sc_mje", 0.005, "shapiro")
print(y_db, '\n')
print(y_nr, np.sum(y_nr))
print("not-normals in Q < 0.6", np.count_nonzero(y_db.loc[0.1:0.5].values), '\n')
print(np.round(y_df, 3))

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  False  False  False  False  False  False  False  False  False  False
0.2  False  False  False  False  False  False  False  False  False  False
0.3  False  False  False  False  False  False  False  False  False  False
0.4  False   True  False  False  False  False  False  False  False  False
0.5   True  False   True   True  False   True   True   True  False   True
0.6  False  False  False   True  False  False  False  False   True  False
0.7  False  False  False  False  False  False  False  False  False  False
0.8  False  False  False  False  False  False  False  False  False  False
0.9  False  False  False  False  False  False  False  False  False  False
1.0   True   True   True   True   True   True   True   True   True   True 

[2 2 2 3 1 2 2 2 2 2] 20
not-normals in Q < 0.6 8 

      1imq  1lmb4   1bf4   2ptl   2ci2   1aps   1fmk   1bk2   1shf   1ten
0.1  0.977  0.096  0.406  0.265  0.569  0.033  0.172  0.219  0.377  0.081
0.2  0.682  0.467  0.581  0.418  0.352  0.086  0.269  0.049  0.806  0.663
0.3  0.674  0.970  0.220  0.193  0.439  0.331  0.118  0.011  0.669  0.777
0.4  0.248  0.000  0.648  0.672  0.055  0.528  0.282  0.720  0.379  0.127
0.5  0.000  0.302  0.000  0.000  0.505  0.000  0.000  0.000  0.960  0.001
0.6  0.440  0.137  0.474  0.003  0.053  0.056  0.013  0.840  0.002  0.028
0.7  0.097  0.237  0.347  0.844  0.933  0.755  0.247  0.957  0.700  0.331
0.8  0.542  0.549  0.090  0.983  0.408  0.451  0.144  0.815  0.430  0.280
0.9  0.149  0.316  0.547  0.913  0.160  0.109  0.026  0.322  0.123  0.466
1.0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000

ND energy distributions of native-state ND nets, are not normal but skewed right.
$Q$ = 1.0 means all native shortcuts have been identified; the PRN0 may not be 100% complete and there may still be a few non-native shortcuts. Hence the distributions are not identical.

fig, axes = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True, sharey=True)
plt.suptitle("ND energy distributions at $Q$=1.0")

for pid in PIDs:
    mjes = c_ndstats[pid].sc_mje.loc[1.0].values
    axes[0].hist(mjes, label=pid)
    mjes = x_ndstats[pid].sc_mje.loc[1.0].values
    axes[1].hist(mjes, label=pid)
    mjes = r_ndstats[pid].sc_mje.loc[1.0].values
    axes[2].hist(mjes, label=pid)
    
axes[0].set_title("'c' ND variant")
axes[1].set_title("'x' ND variant")
axes[2].set_title("'r' ND variant")
_ = axes[0].legend(loc='upper left', bbox_to_anchor=(0.3, 1.0), title='PTS pids')

# Q-Q plots for distributions identified as not-normal
# c 2ci2 0.2
# x 2ci2 0.4
# x 2ptl 0.1
# x 1ten 0.5
# r 1imq 0.4
# may want to check out statsmodels.graphics.gofplots.qqplot
# https://www.statsmodels.org/stable/generated/statsmodels.graphics.gofplots.qqplot.html

def qq_data(pid, q, ndstats: dict, mje_params: dict):
    pp = np.arange(0, 101, 5)[1:] # percentile points
    vals = ndstats[pid].loc[q].sc_mje.values
    mu, stdev = mje_params[pid, q]
    x = np.percentile(norm.rvs(size=100), pp) # independent variable
    y = np.percentile((vals-mu)/stdev, pp)
    return x, y, pearsonr(x, y)

fig, axes = plt.subplots(2, 3, figsize=(12, 8), constrained_layout=True)
plt.suptitle('Percentile-percentile plots')
z = np.linspace(-2, 4)

ax = axes[0][0]
x, y, p = qq_data('2ci2', 0.3, c_ndstats.copy(), c_mje_params.copy())
print("'c' 2ci2 0.3 Example of a normal ", np.round(p, 5))
ax.scatter(x, y, alpha=0.5)
ax.set_title(" 'c' 2ci2 Q=0.3 Example of a normal")
ax.set_ylabel("Standardized 'c' ND energy")
ax.plot(z, z, color='red'); ax.grid(True)

ax = axes[0][1]
x, y, p = qq_data('2ci2', 0.2, c_ndstats.copy(), c_mje_params.copy())
print("'c' 2ci2 0.2 ", np.round(p, 5))
ax.scatter(x, y, alpha=0.5)
ax.set_title(" 'c' 2ci2 Q=0.2 ")
ax.set_ylabel("Standardized 'c' ND energy")
ax.plot(z, z, color='red'); ax.grid(True)

ax = axes[0][2]
x, y, p = qq_data('2ci2', 0.4, x_ndstats.copy(), x_mje_params.copy())
print("'x' 2ci2 0.4 ", np.round(p, 5))
ax.scatter(x, y, alpha=0.5)
ax.set_title(" 'x' 2ci2 Q=0.4 ")
ax.set_ylabel("Standardized 'x' ND energy")
ax.plot(z, z, color='red'); ax.grid(True)

ax = axes[1][0]
x, y, p = qq_data('2ptl', 0.1, x_ndstats.copy(), x_mje_params.copy())
print("'x' 2ptl 0.1 ", np.round(p, 5))
ax.scatter(x, y, alpha=0.5)
ax.set_title(" 'x' 2ptl Q=0.1 ")
ax.set_ylabel("Standardized 'x' ND energy")
ax.plot(z, z, color='red'); ax.grid(True)

ax = axes[1][1]
x, y, p = qq_data('1ten', 0.5, x_ndstats.copy(), x_mje_params.copy())
print("'x' 1ten 0.5 ", np.round(p, 5))
ax.scatter(x, y, alpha=0.5)
ax.set_title(" 'x' 1ten Q=0.5 ")
ax.set_ylabel("Standardized 'x' ND energy")
ax.plot(z, z, color='red'); ax.grid(True)

ax = axes[1][2]
x, y, p = qq_data('1imq', 0.4, r_ndstats.copy(), r_mje_params.copy())
print("'r' 1imq 0.4 ", np.round(p, 5))
ax.scatter(x, y, alpha=0.5)
ax.set_title(" 'r' 1imq Q=0.4 ")
ax.set_ylabel("Standardized 'r' ND energy")
ax.plot(z, z, color='red'); ax.grid(True)

# plots may vary with different runs, but correlations are > 0.9.
# closer the dots are to the red line, "more normal".

'c' 2ci2 0.3 Example of a normal  [0.96687 0.     ]
'c' 2ci2 0.2  [0.92506 0.     ]
'x' 2ci2 0.4  [0.98147 0.     ]
'x' 2ptl 0.1  [0.96097 0.     ]
'x' 1ten 0.5  [0.98605 0.     ]
'r' 1imq 0.4  [0.97909 0.     ]

Visualize some ND energy histograms and their Gaussians.¶

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', \
          'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan']

# adapted from UCSD Fundamentals of Machine Learning edx course
# plot Gaussians at Qs, for a pid
def plot_gaussians(ax, ndv, pid, mje_params):    
    for i, q in enumerate(nd.Qs[:-1]):
        mu, stdev = mje_params[pid, q]
        x_axis = np.linspace(mu - 3*stdev, mu + 3*stdev, 1000)
        ax.plot(x_axis, norm.pdf(x_axis, mu, stdev), colors[i], lw=2, label=str(q))    
    ax.set_title("'" + ndv + "' ND energy Gaussians " + pid, fontsize=12)
    ax.set_xlabel("ND energy (E)", fontsize=12)
    ax.set_ylabel("Density", fontsize=12)
    ax.legend(title='Q')
    return

def plot_energy_profile(ax, nd_type, pid, mje_params, nd_peakQs, nd_peakEs):
    ax.plot(nd.Qs, [ mje_params[pid, q][0] for q in nd.Qs ], linewidth=1)
    p = PIDs.index(pid)
    q = nd_peakQs[p]; i = int(q*10-1)
    ax.scatter(q, nd_peakEs[p], marker='^', color=colors[i], label='peak')      
    ax.set_title("'" + nd_type + "' ND energy profile " + pid)
    ax.set_xlabel('Q', fontsize=12); ax.set_xticks(nd.Qs)
    ax.set_ylabel('ND energy (E)', fontsize=12);     
    ax.legend(); ax.grid()
    return

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,4), constrained_layout=True)

pid = '2ci2'; ndv = 'a'
plot_gaussians(ax1, ndv, pid, a_mje_params)
plot_energy_profile(ax2, ndv, pid, a_mje_params, a_peakQs, a_peakEs)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,4), constrained_layout=True)

pid = '2ci2'; ndv = 'y'
plot_gaussians(ax1, ndv, pid, y_mje_params)
plot_energy_profile(ax2, ndv, pid, y_mje_params, y_peakQs, y_peakEs)

fig, axes = plt.subplots(2, 5, figsize=(20, 8), constrained_layout=True)
plt.suptitle("Gaussians of 'c' ND energy by Q")
for i, p in enumerate(PIDs):
    r = i // 5; c = i % 5 
    plot_gaussians(axes[r, c], 'c', p, c_mje_params)

fig, axes = plt.subplots(2, 5, figsize=(20, 8), constrained_layout=True)
plt.suptitle("Gaussians of 'a' ND energy by Q")
for i, p in enumerate(PIDs):
    r = i // 5; c = i % 5 
    plot_gaussians(axes[r, c], 'a', p, a_mje_params)

fig, axes = plt.subplots(2, 5, figsize=(20, 8), constrained_layout=True)
plt.suptitle("Gaussians of 'x' ND energy by Q")
for i, p in enumerate(PIDs):
    r = i // 5; c = i % 5    
    plot_gaussians(axes[r, c], 'x', p, x_mje_params)

fig, axes = plt.subplots(2, 5, figsize=(20, 8), constrained_layout=True)
plt.suptitle("Gaussians of 'r' ND energy by Q")
for i, p in enumerate(PIDs):
    r = i // 5; c = i % 5    
    plot_gaussians(axes[r, c], 'r', p, r_mje_params)

fig, axes = plt.subplots(2, 5, figsize=(20, 8), constrained_layout=True)
plt.suptitle("Gaussians of 'y' ND energy by Q")
for i, p in enumerate(PIDs):
    r = i // 5; c = i % 5    
    plot_gaussians(axes[r, c], 'y', p, y_mje_params)

End of notebook¶

SKLC August 2020

	pid	N	k_f
0	1imq	86	3.160
1	1lmb4	92	3.690
2	1bf4	63	3.018
3	2ptl	64	1.778
4	2ci2	65	1.681