From d427e551e1d0e9a2abda170bb56272394b87091c Mon Sep 17 00:00:00 2001
From: Gaspard Jankowiak <gaspard.jankowiak@uni-graz.at>
Date: Tue, 5 Mar 2024 11:28:57 +0100
Subject: [PATCH] add plotting scripts

---
 plotting/plot_filtered.py |  31 ++++++++++
 plotting/stats.py         | 123 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 154 insertions(+)
 create mode 100644 plotting/plot_filtered.py
 create mode 100644 plotting/stats.py

diff --git a/plotting/plot_filtered.py b/plotting/plot_filtered.py
new file mode 100644
index 0000000..e7faa0e
--- /dev/null
+++ b/plotting/plot_filtered.py
@@ -0,0 +1,31 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+x2 = np.genfromtxt("POPC-test-5k/Results_collection-X2.dat")
+x2_p2p = x2.max() - x2.min()
+params = np.genfromtxt("POPC-test-5k/Results_collection.dat")
+
+def plot_filtered_hist(threshold):
+    threshold_rounded = str(int(threshold*100))
+
+    idx = x2 < x2.min() + threshold * x2_p2p
+    x2_filtered = x2[idx]
+    params_filtered = params[idx,:]
+
+    plt.subplot(221)
+    plt.hist(x2)
+    plt.title(r"$\chi^2$")
+
+    plt.subplot(222)
+    plt.hist(params[:,13])
+    plt.title("area per lipid")
+
+    plt.subplot(223)
+    plt.hist(x2_filtered)
+    plt.title(r"$\chi^2$ (best " + threshold_rounded + "%)")
+
+    plt.subplot(224)
+    plt.hist(params_filtered[:,13])
+    plt.title("area per lipid (best " + threshold_rounded + "%)")
+
+    plt.show()
diff --git a/plotting/stats.py b/plotting/stats.py
new file mode 100644
index 0000000..48d88e2
--- /dev/null
+++ b/plotting/stats.py
@@ -0,0 +1,123 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+import os
+
+labels = [
+    "Norm (Normalization)",
+    "nv (Vesicle concencentration x10^6 (nm^-3))",
+    "Rm (Radius of the vesicle (nm))",
+    "Z (Radius polydispersity)",
+    "n_TR (Tris fraction )",
+    "d_TR (Tris width (nm))",
+    "s_TR (Tris position (nm))",
+    "d_Chol (CholCH3 position (nm))",
+    "s_Chol (CholCH3 width (nm))",
+    "d_PCN (PCN position (nm))",
+    "s_PCN (PCN width (nm))",
+    "d_CG (CG position (nm))",
+    "s_CG (CG width (nm))",
+    "A_L (Area per lipid (nm²))",
+    "s_D_C (HC excess std (nm))",
+    "s_CH2 (HC smearing (nm))",
+    "d_CH (CH position (nm))",
+    "s_CH (CH width (nm))",
+    "s_CH3 (CH3 width (nm))",
+    "r_PCN (V_PCN/V_HL)",
+    "r_CG (V_GG/V_HL)",
+    "r12 (V_CH/V_CH2)",
+    "r32 (V_CH3/V_CH2)",
+    "T (Temperature (°C) )",
+    "V_BW (Bound-water volume (nm³))",
+    "Con (Constant (nm⁻¹))"
+]
+
+
+def print_x2_stats(dir, threshold=1.0):
+    stats = get_x2_stats(dir, threshold=threshold)
+    thres_str = str(int(100*threshold))
+
+    if threshold < 1:
+        print("  x2 (best %s%%)" % thres_str)
+    else:
+        print("  x2 ")
+
+    print(" min: %f" % stats["min"])
+    print(" max: %f" % stats["max"])
+    print("mean: %f" % stats["mean"])
+    print(" std: %f" % stats["std"])
+
+def std_vs_threshold(dir, obs):
+    x2 = np.genfromtxt(os.path.join(dir, "Results_collection-X2.dat"))
+    if obs == "x2":
+        d = x2
+        suptitle= "x2"
+    else:
+        d = np.genfromtxt(os.path.join(dir, "Results_collection.dat"))[:,obs]
+        suptitle = labels[obs]
+
+    x2_min = x2.min()
+    x2_p2p = x2.max() - x2.min()
+
+    d_first_decile = d[x2 < x2_min + 0.1*x2_p2p]
+
+    thres = np.linspace(0.01, 1, 100)
+
+    std = np.zeros_like(thres)
+    card = np.zeros_like(thres)
+
+    for (i,t) in enumerate(thres):
+        d_filtered = d[x2 < x2_min + t*x2_p2p]
+        std[i] = d_filtered.std()
+        card[i] = len(d_filtered)
+
+    fig, axes = plt.subplots(1, 2)
+
+    plt.suptitle(suptitle)
+
+    ax1 = axes[0]
+    ax3 = axes[1]
+
+    ax1.set_xlabel("threshold")
+    ax1.set_ylabel("stddev", color="blue")
+    ax1.tick_params(axis='y', labelcolor="blue")
+
+    ax2 = ax1.twinx()
+    ax2.set_ylabel("cardinality", color="red")
+    ax2.tick_params(axis='y', labelcolor="red")
+
+    ax3.set_ylabel("all samples", color="blue")
+    ax3.tick_params(axis='y', labelcolor="blue")
+
+    ax4 = ax3.twinx()
+    ax4.set_ylabel("1st decile", color="red")
+    ax4.tick_params(axis='y', labelcolor="red")
+
+    ax1.plot(thres, std, color="blue")
+    ax2.semilogy(thres, card, color="red")
+    ax3.hist(d, color="blue")
+    ax4.hist(d_first_decile, color="red")
+    plt.show()
+
+def get_x2_stats(dir, threshold=1.0):
+    x2 = np.genfromtxt(os.path.join(dir, "Results_collection-X2.dat"))
+
+    p2p = x2.max() - x2.min()
+    x2_filtered = x2[x2 < (x2.min() + threshold * p2p)]
+
+    stats = {
+        "min" : x2_filtered.min(),
+        "max" : x2_filtered.max(),
+        "mean" : x2_filtered.mean(),
+        "std" : x2_filtered.std()
+    }
+
+    return stats
+
+if __name__ == "__main__":
+    if len(sys.argv) > 1 and os.path.isdir(sys.argv[1]):
+        dir = sys.argv[1]
+    else:
+        dir = "POPC-test"
+
+    print_x2_stats(dir)