ralf/Calcium_Model
1
0
Fork 0
Code Issues 2 Pull Requests Actions Packages Projects Releases Wiki Activity

debugging of fitter, many issues still persist

Browse Source
This commit is contained in:
ralf 2024-07-11 13:05:32 +03:00
parent 6e7a727b1a
commit 7cf6590f94
1 changed files with 76 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
fitter.py
View File

@ -6,7 +6,8 @@ import re
from scipy.optimize import least_squares
from model import Model
from scipy.linalg import svd
from Model import Model
from Data import Data
@ -81,38 +82,38 @@ class Fitter:
plt.show()
# exit()
return res # , measured_fluo - calculated_fluo)
return res, err # , measured_fluo - calculated_fluo)
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
tau_xfer = m.tau_xfer
tau_RC = 1.5
offset = 0
d = self.data
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),
(10, 100, 100, 1, 100, 10, 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())
print(" Optimized:", res.x.tolist())
print(" Optim status:", res.status)
print("Optim message:", res.message)
gGaL, ECal, K_pc_half, tau_xfer, tau_RC, offset = res.x
self.fit_results.update({
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
tau_xfer = m.tau_xfer
tau_RC = 1.5
offset = 0
d = self.data
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),
(10, 100, 100, 1, 100, 10, 10),
)
res, err = 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())
print(" Optimized:", res.x.tolist())
print(" Optim status:", res.status)
print("Optim message:", res.message)
gGaL, ECal, K_pc_half, tau_xfer, tau_RC, offset = res.x
self.fit_results.update({
'gGaL': gGaL,
'ECal': ECal,
'K_pc_half': K_pc_half,
@ -120,49 +121,50 @@ class Fitter:
'tau_RC': tau_RC,
'offset': offset,
'mean_squared_error': err})
model = self.model()
model.ECaL = ECal
model.gCaL = gGaL
model.K_pc_half = K_pc_half
model.tau_xfer = tau_xfer
model.solve(times=self.time_points)
_calc_curr = model.calculated_current()
calculated_current = self.convolve_current(_calc_curr, tau=tau_RC) + offset
print("Elapsed time:", time.time() - t0)
fig = plt.figure(figsize=(24, 12))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.plot(1000 * self.data.current_t, self.data.current, label="Measured")
ax1.plot(self.time_points, calculated_current, label="Calculated")
ax1.set_xlabel("time, ms")
ax1.set_ylabel("current, 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, calculated_current[self.current_time_indecies], label="Calculated")
ax2.set_xlabel("time, ms")
ax2.set_ylabel("current, pA/pF")
ax2.legend(frameon=False)
return res, fig
model = self.model()
model.ECaL = ECal
model.gCaL = gGaL
model.K_pc_half = K_pc_half
model.tau_xfer = tau_xfer
model.solve(times=self.time_points)
_calc_curr = model.calculated_current()
calculated_current = self.convolve_current(_calc_curr, tau=tau_RC) + offset
print("Elapsed time:", time.time() - t0)
fig = plt.figure(figsize=(24, 12))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.plot(1000 * self.data.current_t, self.data.current, label="Measured")
ax1.plot(self.time_points, calculated_current, label="Calculated")
ax1.set_xlabel("time, ms")
ax1.set_ylabel("current, 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, calculated_current[self.current_time_indecies], label="Calculated")
ax2.set_xlabel("time, ms")
ax2.set_ylabel("current, pA/pF")
ax2.legend(frameon=False)
return res, fig
def covcor_from_lsq(res):
_, s, VT = 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]
cov = np.dot(VT.T / s**2, VT)
std = np.sqrt(np.diag(cov))
cor = cov / np.outer(std, std)
cor[cov == 0] = 0
return cov, cor
def plot_correlation_matrix(cor):
@ -173,30 +175,32 @@ class Fitter:
plt.ylabel('Variables')
plt.show()
if __name__ == "__main__":
if __name__ == "__main__":
filename = "ltcc_current.h5"
eid = "0033635a51b096dc449eb9964e70443a67fc16b9587ae3ff6564eea1fa0e3437_2018.06.18 14:48:40"
data = Data(filename, eid)
fit = Fitter(Model, data)
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) # 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")
# plot_filename = "fit_plot"
# fig.savefig(f"{plot_filename}.png")
# fig.savefig(f"{plot_filename}.pdf")
fig.savefig("naidis_fit.pdf")
plt.show()