didnt ad pyproject file

.gitignore

@@ -3,3 +3,6 @@ fit_results_1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_201
 plot_1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_2017.11.27140704.pdf
 plot_1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_2017.11.27140704.png
 naidis_fit.pdf
+__pycache__/Data.cpython-312.pyc
+__pycache__/Model.cpython-312.pyc
+__pycache__/fitter.cpython-312.pyc

Model.py

@@ -33,46 +33,52 @@ class Model:
         self.km_Ca: float = (
             1380  # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange
         )
-        self.k_Na_Ca: float = 292.8  # pA/pF, Scaling factor of Na2+/Ca2+ exchange
+        self.k_Na_Ca: float = 292.8
+        # pA/pF, Scaling factor of Na2+/Ca2+ exchange
 
         self.I_max_pCa: float = 1.0  # pA/pF, Maximum Ca2+ pump current
 
-        # self.Cai = 0.1  # uM, Cytoplasmic Ca2+ concentration fixed at this point
+        # self.Cai = 0.1  # uM, Cytoplasmic Ca2+ concentration fixed at point
 
         self.LTRPN_tot: float = (
-            70.0  # uM, Total myoplasmic troponin low-affinity site concentration
+            70.0  # uM, Total myoplasmic troponin low-affinity site concentr.
         )
         self.HTRPN_tot: float = (
-            140.0  # uM, Total myoplasmic troponin high-affinity site concentration
+            140.0  # uM, Total myoplasmic troponin high-affinity site concentr.
         )
 
-        # 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.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: float = (
-            0.00237  # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affinity sites
+            0.00237
+            # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affin 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
+            0.0327
+            # uM^(-1)/ms, Ca2+ on rate const. for troponin low-affin sites
         )
         self.k_ltrpn_negative: float = (
             0.0196  # ms, Ca2+ off rate const. for troponin low-affinity sites
         )
 
-        self.CMDN_tot: float = 50  # uM, Total myoplasmic calmodulin concentration
-        self.Km_CMDN: float = 0.238  # uM, Ca2 half-saturation constant for calmodulin
+        self.CMDN_tot: float = 50  # uM, Total myoplasmic calmodulin concentr.
+        self.Km_CMDN: float = 0.238  # uM, Ca2 half-sat const for calmodulin
 
         self.nu_1: float = (
-            4.5  # 1/ms, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR
+            4.5  # 1/ms, Max RyR chan Ca2+ perm. Ca2+ leak rate const from NSR
         )
-        self.nu_2: float = 1.74e-5  # ms^(-1), Ca2+ leak rate const. from the NSR
+        self.nu_2: float = 1.74e-5  # ms^(-1), Ca2+ leak rate const. from NSR
         self.nu_3: float = 0.45  # uM/ms, SR Ca2+ -ATPase maximum pupmp rate
 
-        self.Km_up: float = 0.5  # uM, Half-saturation constant for SR Ca2+ -ATPase pump
+        self.Km_up: float = 0.5  # uM, Half-sat const for SR Ca2+ -ATPase pump
         self.km_p_ca: float = (
             0.5  # uM, Ca2+ half-saturation constant for Ca2+ pump current
         )
@@ -84,7 +90,7 @@ class Model:
         )
 
         self.K_pc_max: float = (
-            0.23324  # 1/ms, Maximum time constant for Ca2+-induced inactivation
+            0.23324  # 1/ms, Maximum time constant for Ca2+-induced inact.
         )
         self.K_pc_half: float = (
             20.0  # uM, Half-saturation constant for Ca2+-induced inactivation
@@ -97,7 +103,7 @@ class Model:
 
         self.Cao: float = 1130.0  # uM, Ca2+ outside the cell
         self.gCaL: float = (
-            0.1729  # mS/uF, Specific maximum conductivity for L-type Ca2+ channel
+            0.1729  # mS/uF, Specific max conduct. for L-type Ca2+ channel
         )
         self.ECaL: float = 43.0  # mV, Reversal potential for L-type Ca2
         self.V_ss: float = 1.485e-9  # uL, Dyadic aka subspace volume
@@ -117,9 +123,9 @@ class Model:
             800.0  # uM, Ca2 half-saturation constant for calsequestrin
         )
 
-        self.k_a_positive: float = 0.006075  # uM^(-4)/ms, RyR Pc1 - Po1 rate constant
+        self.k_a_positive: float = 0.006075  # uM^(-4)/ms, RyR Pc1-Po1 rate
         self.k_a_negative: float = 0.07125  # 1/ms, RyR Po1 - Pc1 rate constant
-        self.k_b_positive: float = 0.00405  # uM^(-3)/ms, RyR Po1 - Po2 rate constant
+        self.k_b_positive: float = 0.00405  # uM^(-3)/ms, RyR Po1 - Po2 rate
         self.k_b_negative: float = 0.965  # 1/ms, RyR Po2 - Po1 rate constant
         self.k_c_positive: float = 0.009  # 1/ms, RyR Po1 - Pc2 rate constant
         self.k_c_negative: float = 0.0008  # 1/ms, RyR Pc2 - Po1 rate constant
@@ -160,7 +166,7 @@ class Model:
 
         Kpcf = 13.0 * (1 - np.exp(-((V + 14.5) ** 2) / 100))
 
-        alpha = (  # dimensionless parameter for L-type Ca2+ channel closed states
+        alpha = (  # dimensionless param. for L-type Ca2+ channel closed states
             0.4
             * np.exp((V + 12.0) / 10)
             * (
@@ -267,8 +273,10 @@ class Model:
         P_C1 = 1 - (P_C2 + P_O1 + P_O2)
 
         dP_O1dt = self.k_a_positive * (Cass) ** self.n_ryr * P_C1
-        -self.k_a_negative * P_O1 - self.k_b_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
+        -self.k_a_negative * P_O1
+        -self.k_b_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
 
         dP_O2dt = (
             self.k_b_positive * (Cass) ** self.m_ryr * P_O1 - self.k_b_negative * P_O2
@@ -289,7 +297,8 @@ class Model:
 
         """
         dCaidt = Bi * (J_leak + J_xfer - J_up - J_trpn
-        - (I_Cab - 2 * I_NaCa + I_pCa) * ((self.A_cap * self.C_mem)/(2 * self.V_myo * F)))
+        - (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) #lisatud ette ennetavalt
@@ -308,8 +317,6 @@ class Model:
             * ((self.A_cap * self.C_mem) / (2 * self.V_myo * F))
         )
 
-        # Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Ca_NSR, Ca_JSR, P_RyR, P_O1, P_O2, P_C2, Cai, FCa = states
-
         return [
             dCassdt,
             dOsdt,
@@ -418,7 +425,7 @@ class Model:
         ax1.set_xlabel("time [ms]")
         ax1.set_ylabel(r"$\mu mol/l$")
 
-        ax2.plot(times, self.gCaL * states[1, :] * (V - self.ECaL), label="ICaL ")
+        ax2.plot(times, self.gCaL * states[1, :] * (V - self.ECaL), label="ICaL")
         ax2.legend(frameon=False)
         ax2.set_xlabel("time [ms]")
         ax2.set_ylabel("pA/pF")
@@ -444,6 +451,8 @@ if __name__ == "__main__":
     initial_values = model.get_initial_values()
     times = np.arange(0, 1000, 1.0)
 
-    states, V = model.solve(initial_values=initial_values, tspan=[0, 1000], dt=1.0, times=times)
+    states, V = model.solve(
+        initial_values=initial_values, tspan=[0, 1000], dt=1.0, times=times
+    )
 
     model.plot_results(times, states, V)

experiment_fitter.py

@@ -6,6 +6,45 @@ from Model import Model
 from fitter import Fitter
 from Data import Data
 
+file = "ltcc_current.h5"
+
+
+def print_attrs(name, obj):
+    print(f"\nAttributes for {name}:")
+    for key, val in obj.attrs.items():
+        print(f"  {key}: {val}")
+
+
+with h5py.File(file, "r") as h5_file:
+    h5_file.visititems(print_attrs)
+
+# Dict to hold DFs 'sex', 'tag' & 'spid'
+dfs_by_sex_tag_spid = {}
+
+with h5py.File(file, "r") as h5_file:
+    for eid in h5_file.keys():
+        attributes = h5_file[eid].attrs
+        sex = attributes.get("sex")
+        tag = attributes.get("tag")
+        spid = attributes.get("spid")
+
+        key = f"{sex}_{tag}_{spid}"
+
+        if key not in dfs_by_sex_tag_spid:
+            dfs_by_sex_tag_spid[key] = pd.DataFrame()
+
+        row_data = {"experiment_id": eid, "sex": sex, "tag": tag, "spid": spid}
+        temp_df = pd.DataFrame([row_data])
+
+        dfs_by_sex_tag_spid[key] = pd.concat(
+            [dfs_by_sex_tag_spid[key], temp_df], ignore_index=True
+        )
+
+for key, df in dfs_by_sex_tag_spid.items():
+    print(f"DataFrame for {key}:")
+    print(df)
+    print()
+
 
 def fit_data():
     filename = "ltcc_current.h5"

fitter.py

@@ -61,16 +61,13 @@ class Fitter:
         model.K_pc_half = K_pc_half
         model.tau_xfer = tau_xfer
 
-        model.solve(
+        states, V = model.solve(
             times=self.time_points,
             tspan=self.tspan,
             dt=self.dt,
             initial_values=self.initial_values,
         )
 
-        states = self.initial_values
-        t = self.time_points
-        V = model.mem_potential(t)
         _calc_curr = model.calculated_current(states, V)
         calculated_current = (
             self.convolve_current(_calc_curr, tau=tau_RC)[self.current_time_indecies]
@@ -100,6 +97,7 @@ class Fitter:
     def optimize(self, init_parameters=None):
         t0 = time.time()
         self.iteration = 0
+
         if init_parameters is None:
             m = self.model()
             K_pc_half = m.K_pc_half
@@ -118,7 +116,15 @@ class Fitter:
             (10, 100, 100, 100, 10, 10),
         )
 
-        res = 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,
+            ftol=1e-10,
+            #  verbose=2,
+        )
+
         print()
         print("   Parameters: [gGaL, ECal, K_pc_half, tau_xfer, tau_RC, offset]")
         print("      Initial:", init_parameters.tolist())

multiple_experiment_fitter.py

@@ -0,0 +1,136 @@
+import h5py
+import pandas as pd
+import numpy as np
+import re
+
+import matplotlib as plt
+
+from Model import Model
+from fitter import Fitter
+from Data import Data
+
+file = "ltcc_current.h5"
+
+dfs_by_sex_tag_spid = {}
+
+
+def print_attrs(name, obj):
+    #  print(f"\nAttributes for {name}:")
+    #  for key, val in obj.attrs.items():
+    #  print(f"  {key}: {val}")
+    pass
+
+
+with h5py.File(file, "r") as h5_file:
+    for eid in h5_file.keys():
+        attributes = h5_file[eid].attrs
+        sex = attributes.get("sex")
+        tag = attributes.get("tag")
+        spid = attributes.get("spid")
+
+        key = f"{sex}_{tag}_{spid}"
+
+        if key not in dfs_by_sex_tag_spid:
+            dfs_by_sex_tag_spid[key] = pd.DataFrame()
+
+        row_data = {"experiment_id": eid, "sex": sex, "tag": tag, "spid": spid}
+        temp_df = pd.DataFrame([row_data])
+
+        dfs_by_sex_tag_spid[key] = pd.concat(
+            [dfs_by_sex_tag_spid[key], temp_df], ignore_index=True
+        )
+
+#  for key, df in dfs_by_sex_tag_spid.items():
+#  print(f"DataFrame for {key}:")
+#  print(df)
+#  print()
+
+
+def fit_and_plot_dataframes(dfs_by_sex_tag_spid):
+    for key, df in dfs_by_sex_tag_spid.items():
+        print(f"Fitting and plotting data for {key}...")
+
+        combined_current = []
+        combined_time = []
+
+        for eid in df["experiment_id"].tolist():
+            data = Data(file, group_key=eid)
+            combined_current.append(data.current)
+            combined_time.append(data.current_t)
+
+        combined_current = np.concatenate(combined_current)
+        combined_time = np.concatenate(combined_time)
+
+        # Sort by time for consistency - inspired by data
+        sorted_indices = np.argsort(combined_time)
+        combined_time = combined_time[sorted_indices]
+        combined_current = combined_current[sorted_indices]
+
+        combined_data = Data(file, group_key=eid)
+        combined_data.current = combined_current
+        combined_data.current_t = combined_time
+
+        fit = Fitter(Model, combined_data)
+        fit.optimize()
+        res, fig = fit.optimize()
+
+        plt.figure()
+        plt.title(f"Fit results for {key}")
+
+        for i, eid in enumerate(df["experiment_id"].tolist()):
+            plt.plot(combined_time, combined_current,
+                     label=f"Experiment {eid}")
+
+        plt.plot(
+            combined_data.current_t, combined_data.current, "k-",
+            label="Combined Fit"
+        )
+
+        plt.legend()
+        plt.xlabel("Time")
+        plt.ylabel("Current")
+
+        key_cleaned = re.sub(r"[^\w.-]", "", key)
+        plt.savefig(f"combined_plot_{key_cleaned}.png")
+        plt.savefig(f"combined_plot_{key_cleaned}.pdf")
+        plt.close()
+
+        fit_hist = pd.DataFrame.from_dict(fit.fit_results, orient="index").T
+        fit_hist.index.name = "Iterations"
+
+        res_filename = f"combined_fit_results_{key}.csv"
+        res_filename = res_filename.replace(" ", "_").replace(":", "-")
+        fit_hist.to_csv(res_filename, index=True)
+
+        print(f"Finished fitting for {key}.")
+
+
+fit_and_plot_dataframes(dfs_by_sex_tag_spid)
+
+"""
+def fit_data():
+    filename = "ltcc_current.h5"
+    with h5py.File(filename, "r") as h5:
+        eids = list(h5.keys())
+
+    for eid in eids:
+        data = Data(filename, group_key=eid)
+        fit = Fitter(Model, data)
+        fit.optimize()
+        res, fig = fit.optimize()
+
+        fit_hist = pd.DataFrame.from_dict(fit.fit_results, orient="index").T
+        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)  # Eemaldab kõik eritähed
+        fig.savefig(f"plot_{eid_cleaned}.png")
+        fig.savefig(f"plot_{eid_cleaned}.pdf")
+
+
+fit_data()
+
+"""