cleared up more errors

.gitignore

@@ -0,0 +1 @@
+ltcc_current.h5

Model.py

@@ -1,18 +1,13 @@
-# -*- coding: utf-8 -*-
-
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import ode
 
-# %% Constants
-
 R = 8.314  # ideal gas constant (J/(mol*K))
 F = 96.485  # coulomb_per_millimole, Faraday constant
 T = 298  # room temperature (K)
 RT = R * T
 
 
-# %% model
 class Model:
     def __init__(self):
 
@@ -20,48 +15,64 @@ class Model:
         self.C_mem = 1.0  # uF/cm2 Specific membrane capacitance
         self.V_myo = 25.84e-6  # uL, Myoplasmic volume
         self.current2flux = self.A_cap * self.C_mem / 2 / F  #
+
         self.period = 1000  # ms, pulse period
         self.V_mem_rest = -80.0  # mV, resting membrane potential
+
         self.Nai = 11000  # uM, Myoplasmic Na+ concentration
         self.Nao = 150000  # uM, Extracellular Na+ concentration
+
         self.eta = 0.35  # Controls voltage dependance of Na/Ca2+ exchange
+
         self.km_Na = 87500  # uM, Na+ half-saturation constant for Na+/Ca2+ exchange
         self.k_sat = (
             0.1  # Na+/Ca2+ exchange saturation factor at very negative potentials
         )
         self.km_Ca = (1380,)  # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange
         self.k_Na_Ca = 292.8  # pA/pF, Scaling factor of Na2+/Ca2+ exchange
+
         self.I_max_pCa = 1.0  # pA/pF, Maximum Ca2+ pump current
+
         # self.Cai = 0.1  # uM, Cytoplasmic Ca2+ concentration fixed at this point
+
         self.LTRPN_tot = (
             70.0  # uM, Total myoplasmic troponin low-affinity site concentration
         )
         self.HTRPN_tot = (
             140.0  # uM, Total myoplasmic troponin high-affinity site concentration
         )
+
         # self.LTRPNCa = 11.2684 # uM, Concentration Ca2+ bound low-affinity troponin-binding sites
         # self.HTRPNCa = 125.290 # uM, Concentration Ca2+ bound high-affinity troponin-binding sites
+
         self.k_htrpn_positive = (
             0.00237  # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affinity sites
         )
         self.k_htrpn_negative = (
             3.2e-5  # ms, Ca2+ off rate const. for troponin high-affinity sites
         )
+
         self.k_ltrpn_positive = (
             0.0327  # uM^(-1)/ms, Ca2+ on rate const. for troponin low-affinity sites
         )
         self.k_ltrpn_negative = (
             0.0196  # ms, Ca2+ off rate const. for troponin low-affinity sites
         )
+
         self.CMDN_tot = 50  # uM, Total myoplasmic calmodulin concentration
         self.Km_CMDN = 0.238  # uM, Ca2 half-saturation constant for calmodulin
+
         self.nu_1 = 4.5  # 1/ms, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR
         self.nu_2 = 1.74e-5  # ms^(-1), Ca2+ leak rate const. from the NSR
         self.nu_3 = 0.45  # uM/ms, SR Ca2+ -ATPase maximum pupmp rate
+
         self.Km_up = 0.5  # uM, Half-saturation constant for SR Ca2+ -ATPase pump
         self.km_p_ca = 0.5  # uM, Ca2+ half-saturation constant for Ca2+ pump current
+
         # self.Ca_NSR = 1299.50 # uM,NSR Ca2+ concentration
+
         self.tau_xfer = 8.0  # ms, Time constant for transfer from subspace to myoplasm
+
         self.K_pc_max = (
             0.23324  # 1/ms, Maximum time constant for Ca2+-induced inactivation
         )
@@ -70,29 +81,37 @@ class Model:
         )
         self.Kpcb = 0.0005  # 1/ms, Voltage-insensitive rate constant for inactivation
         self.km_Ca = 1380  # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange
+
         self.Gcab = 0.000367  # mS/uF
+
         self.Cao = 1130.0  # uM, Ca2+ outside the cell
         self.gCaL = (
             0.1729  # mS/uF, Specific maximum conductivity for L-type Ca2+ channel
         )
         self.ECaL = 43.0  # mV, Reversal potential for L-type Ca2⫹ channel, kas arvutame voi jaab const?
         self.V_ss = 1.485e-9  # uL, Dyadic aka subspace volume
+
         # self.Ca_JSR = 1299.50 # uM, JSR Ca2+ concentration
         self.F_tot = 25  # uM, total concentration of Fluo-4
         self.k_on = 0.1  # 1/uM * 1/ms, Fluo-4 reaction rate constant
         self.k_off = 0.11  # 1/ms, Fluo-4 reaction rate constant 2
+
         self.V_NSR = 2.098e-6  # ul, Network SR volume
         self.tau_tr = 20  # ms, Time const for transfer from NSR to JSR
+
         self.CSQN_tot = 15000.0  # uM, total junctional SR calsequestrin concentration
         self.Km_CSQN = 800.0  # uM, Ca2 half-saturation constant for calsequestrin
+
         self.k_a_positive = 0.006075  # uM^(-4)/ms, RyR Pc1 - Po1 rate constant
         self.k_a_negative = 0.07125  # 1/ms, RyR Po1 - Pc1 rate constant
         self.k_b_positive = 0.00405  # uM^(-3)/ms, RyR Po1 - Po2 rate constant
         self.k_b_negative = 0.965  # 1/ms, RyR Po2 - Po1 rate constant
         self.k_c_positive = 0.009  # 1/ms, RyR Po1 - Pc2 rate constant
         self.k_c_negative = 0.0008  # 1/ms, RyR Pc2 - Po1 rate constant
+
         self.n_ryr = 4  # RyR Ca2+ cooperativity parameter Pc1 - Po1
         self.m_ryr = 3  # RyR Ca2+ cooperativity parameter Po1 - Po2
+
         self.I_CaL_max = 7.0  # pA/pF, normalization constant for L-type Ca2+ current
         self.V_JSR = 0.12e-6  # ul, Junctional SR volume
 
@@ -180,11 +199,11 @@ class Model:
 
         Bi = 1 / (
             1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cai) ** 2
-        )  # A6 !!!!!!! korras V8
+        )
 
         Bss = 1 / (
             1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2
-        )  # A7 !!!!!!!!! korras V8
+        )
 
         J_xfer = (Cass - Cai) / self.tau_xfer
 
@@ -199,7 +218,7 @@ class Model:
             - self.k_htrpn_negative * HTRPNCa
             + self.k_ltrpn_positive
             * Cai
-            * (self.LTRPN_tot - LTRPNCa)  # !!!!!! korras v8
+            * (self.LTRPN_tot - LTRPNCa)
             - self.k_ltrpn_negative * LTRPNCa
         )
 
@@ -207,7 +226,7 @@ class Model:
 
         J_tr = (Ca_NSR - Ca_JSR) / self.tau_tr
 
-        J_leak = self.nu_2 * (Ca_NSR - Cai)  # puudu CaNSR diff vorrand
+        J_leak = self.nu_2 * (Ca_NSR - Cai)
 
         dCa_NSRdt = (J_up - J_leak) * (self.V_myo / self.V_NSR)
         -J_tr * (self.V_JSR / self.V_NSR)
@@ -239,7 +258,9 @@ class Model:
 
         P_C1 = 1 - (P_C2 + P_O1 + P_O2)
 
-        dP_O1dt = self.k_a_positive * (Cass) ** self.n_ryr * P_C1
+        dP_O1dt = (
+            self.k_a_positive * (Cass) ** self.n_ryr * P_C1
+        )
         -self.k_a_negative * P_O1 - self.k_a_positive * (Cass) ** self.m_ryr * P_O1
         +self.k_b_negative * P_O2 - self.k_c_positive * P_O1 + self.k_c_negative * P_C2
 
@@ -265,7 +286,7 @@ class Model:
         - (I_Cab - 2 * I_NaCa + I_pCa) * ((self.A_cap * self.C_mem)/(2 * self.V_myo * F)))
         """
 
-        # J_rel_caf = self.nu_1 * (self.Ca_JSR - Cass)
+        # J_rel_caf = self.nu_1 * (self.Ca_JSR - Cass) #lisatud ette ennetavalt
 
         JFCa = self.k_on * (self.F_tot - FCa) * Cai - self.k_off * FCa
 

fitter.py

@@ -4,18 +4,19 @@ import pylab as plt
 import pandas as pd
 import re
 
-
 from scipy.optimize import least_squares
-from scipy.linalg import svd
 from Model import Model
 from Data import Data
 
 
 class Fitter:
-    def __init__(self, model: Model, data: Data, current_fit_range: tuple = (107, 341)) -> None:
+    def __init__(
+        self, model: Model, data: Data, current_fit_range: tuple = (107, 341)
+    ) -> None:
         """
         current_fit_range : tuple (t0, t1), t0 nad t1 are start and stop times in between current is fitted
         """
+
         self.model = model
         self.data = data
         self.current_fit_range = current_fit_range
@@ -60,19 +61,21 @@ class Fitter:
         model.solve(times=self.time_points)
 
         _calc_curr = model.calculated_current()
-        calculated_current = self.convolve_current(_calc_curr, tau=tau_RC)[self.current_time_indecies] + offset
+        calculated_current = (
+            self.convolve_current(_calc_curr, tau=tau_RC)[self.current_time_indecies]
+            + offset
+        )
 
         res = self.measured_current - calculated_current
         err = np.mean(res**2)  # mean squared error
         self.iteration += 1
         print(self.iteration, parameters.tolist(), "err", err)
-        # measured_fluo = self.data.fluo
-        # fluo_interplolator = interp1d(self.time_domain, model.calculated_fluo)
-        # calculated_fluo = fluo_interplolator(self.data.fluo_time)
 
         if self.iteration < -0:
             t = self.time_points[self.current_time_indecies]
-            plt.plot(t, _calc_curr[self.current_time_indecies], label="calculated current")
+            plt.plot(
+                t, _calc_curr[self.current_time_indecies], label="calculated current"
+            )
             plt.plot(t, self.measured_current, label="measured current")
             plt.plot(t, calculated_current, label="conv calculated current")
             plt.plot(t, self.measured_current - calculated_current, label="error")
@@ -80,9 +83,8 @@ class Fitter:
             plt.ylabel("current, pA/pF")
             plt.legend(frameon=False)
             plt.show()
-            # exit()
 
-        return res, err  # , measured_fluo - calculated_fluo)
+        return res  # , measured_fluo - calculated_fluo)
 
     def optimize(self, init_parameters=None):
         t0 = time.time()
@@ -95,15 +97,17 @@ class Fitter:
             offset = 0
 
             d = self.data
-            init_parameters = np.array([d.gGaL, d.ECal, K_pc_half, tau_xfer, tau_RC, offset])
+            init_parameters = np.array(
+                [d.gGaL, d.ECal, K_pc_half, tau_xfer, tau_RC, offset]
+            )
             print(init_parameters.tolist())
 
         bounds = (
-            (0.01, 10, 0.1, 0.1, 0.1, -5, -10),
+            (0.01, 10, 0.1, 0.1, -5, -10),
-            (10, 100, 100, 1, 100, 10, 10),
+            (10, 100, 100, 100, 10, 10),
         )
 
-        res, err = least_squares(self.cost_func, init_parameters, bounds=bounds, xtol=1e-10)
+        res = least_squares(self.cost_func, init_parameters, bounds=bounds, xtol=1e-10)
         print()
         print("   Parameters: [gGaL, ECal, K_pc_half, tau_xfer, tau_RC, offset]")
         print("      Initial:", init_parameters.tolist())
@@ -113,14 +117,19 @@ class Fitter:
 
         gGaL, ECal, K_pc_half, tau_xfer, tau_RC, offset = res.x
 
-        self.fit_results.update({
+        err = self.cost_func(res.x)
-            'gGaL': gGaL,
+
-            'ECal': ECal,
+        self.fit_results.update(
-            'K_pc_half': K_pc_half,
+            {
-            'tau_xfer': tau_xfer,
+                "gGaL": gGaL,
-            'tau_RC': tau_RC,
+                "ECal": ECal,
-            'offset': offset,
+                "K_pc_half": K_pc_half,
-            'mean_squared_error': err})
+                "tau_xfer": tau_xfer,
+                "tau_RC": tau_RC,
+                "offset": offset,
+                "mean_squared_error": np.mean((err) ** 2),
+            }
+        )
 
         model = self.model()
 
@@ -135,27 +144,26 @@ class Fitter:
         calculated_current = self.convolve_current(_calc_curr, tau=tau_RC) + offset
         print("Elapsed time:", time.time() - t0)
 
-        fig = plt.figure(figsize=(24, 12))
+        fig = plt.figure(figsize=(6, 3))
         ax1 = fig.add_subplot(121)
         ax2 = fig.add_subplot(122)
-        ax1.plot(1000 * self.data.current_t, self.data.current, label="Measured")
+        ax1.plot(1000 * self.data.current_t, self.data.current, label="Mõõdetud")
-        ax1.plot(self.time_points, calculated_current, label="Calculated")
+        ax1.plot(self.time_points, calculated_current, label="Arvutatud")
-        ax1.set_xlabel("time, ms")
+        ax1.set_xlabel("Aeg [ms]")
-        ax1.set_ylabel("current, pA/pF")
+        ax1.set_ylabel("Vool [pA/pF]")
         ax1.legend(frameon=False)
 
         tp = self.time_points[self.current_time_indecies]
-        ax2.plot(tp, self.measured_current, label="Measured")
+        ax2.plot(tp, self.measured_current, label="Mõõdetud")
-        ax2.plot(tp, calculated_current[self.current_time_indecies], label="Calculated")
+        ax2.plot(tp, calculated_current[self.current_time_indecies], label="Arvutatud")
-        ax2.set_xlabel("time, ms")
+        ax2.set_xlabel("Aeg [ms]")
-        ax2.set_ylabel("current, pA/pF")
+        ax2.set_ylabel("Vool [pA/pF]")
         ax2.legend(frameon=False)
-
         return res, fig
 
     def covcor_from_lsq(res):
 
-        _, s, VT = svd(res.jac, full_matrices=False)
+        _, s, VT = np.linalg.svd(res.jac, full_matrices=False)
         threshold = np.finfo(float).eps * max(res.jac.shape) * s[0]
         s = s[s > threshold]
         VT = VT[: s.size]
@@ -168,18 +176,18 @@ class Fitter:
         return cov, cor
 
     def plot_correlation_matrix(cor):
-        plt.imshow(cor, cmap='viridis', interpolation='nearest')
+        plt.imshow(cor, cmap="viridis", interpolation="nearest")
-        plt.colorbar(label='Correlation')
+        plt.colorbar(label="Correlation")
-        plt.title('Correlation Matrix')
+        plt.title("Correlation Matrix")
-        plt.xlabel('Variables')
+        plt.xlabel("Variables")
-        plt.ylabel('Variables')
+        plt.ylabel("Variables")
         plt.show()
 
 
 if __name__ == "__main__":
 
     filename = "ltcc_current.h5"
-    eid = "0033635a51b096dc449eb9964e70443a67fc16b9587ae3ff6564eea1fa0e3437_2018.06.18 14:48:40"
+    eid = "1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_2017.11.27 14:07:04"
 
     data = Data(filename, eid)
 
@@ -187,14 +195,16 @@ if __name__ == "__main__":
 
     res, fig = fit.optimize()
 
-    fit_hist = pd.DataFrame.from_dict(fit.fit_results, orient='index').T
+    fit_hist = pd.DataFrame.from_dict(fit.fit_results, orient="index").T
-    fit_hist.index.name = 'Iterations'
+    fit_hist.index.name = "Iterations"
 
     res_filename = f"fit_results_{eid}.csv"
     res_filename = res_filename.replace(" ", "_").replace(":", "-")
     fit_hist.to_csv(res_filename, index=True)
 
-    eid_cleaned = re.sub(r'[^\w.-]', '', eid)  # Eemalda kõik eritähed ja jääb alles alphanumbrilised tähed, sidekriipsud ja punktid
+    eid_cleaned = re.sub(
+        r"[^\w.-]", "", eid
+    )  # Eemalda kõik eritähed ja jääb alles alphanumbrilised tähed, sidekriipsud ja punktid
     fig.savefig(f"plot_{eid_cleaned}.png")
     fig.savefig(f"plot_{eid_cleaned}.pdf")
 
@@ -202,5 +212,4 @@ if __name__ == "__main__":
     # fig.savefig(f"{plot_filename}.png")
     # fig.savefig(f"{plot_filename}.pdf")
     fig.savefig("naidis_fit.pdf")
-
     plt.show()