diff --git a/fitter.py b/fitter.py index 0e337e2..b771b90 100644 --- a/fitter.py +++ b/fitter.py @@ -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()