Changes for flexible no. of sweeps

Author: Frankie Patten-Elliott

methods/classes.py

@@ -53,7 +53,8 @@ def __init__(self, model, protocol, times, win, conc, param_dict):
         if model.split("-")[0] == 'sis':
             self.n_pulses = 10
         elif times[-1] > 15e3:
-            self.n_pulses = 5
+            #self.n_pulses = 5
+            self.n_pulses = int(np.floor(250000/times[-1]))
         else:
             self.n_pulses = sweeps
         self._times = times

methods/funcs.py

@@ -13,7 +13,7 @@ def generate_data(herg_model, drug_vals, prot, sd, max_time, bounds, m_sel, conc
         herg_vals = []
 
     if max_time != 15e3 and herg_model != '2024_Joey_sis_25C':
-        swps = 5
+        swps = int(np.floor(250000/max_time))
     else:
         swps = sweeps
 
@@ -147,7 +147,7 @@ def model_outputs(model_pars, herg_model, prot, times, concs, alt_protocol = Non
     elif herg_model == 'kemp':
         herg_vals = []
     model_out = {}
-    swps = 5 #when generating optimal protocol, we use 5 sweeps
+    swps = int(np.floor(250000/(times[-1]+alt_times[-1])))
     win = (times >= wins[0]) & (times < wins[1])
     if alt_protocol is not None:
         win_alt = (alt_times >= wins_alt[0]) & (alt_times < wins_alt[1])
@@ -180,7 +180,7 @@ def model_outputs(model_pars, herg_model, prot, times, concs, alt_protocol = Non
                 after = model.simulate(binding_params, times)[win]
                 model_milnes = np.append(model_milnes, after/before)
                 control = np.append(control, before)
-                for i in range(swps-1):
+                for i in range(swps*2-1):
                     if (alt_protocol is not None) & ((i % 2) == 0):
                         model.change_protocol(alt_protocol)
                         after = model.simulate(binding_params, alt_times, reset=False)[win_alt]
@@ -196,4 +196,4 @@ def model_outputs(model_pars, herg_model, prot, times, concs, alt_protocol = Non
             except:
                 model_out[m][conc] = np.ones(times.shape) * float('inf')
                 save_control = np.ones(times.shape) * float('inf')
-    return model_out, save_control
+    return model_out, save_control, swps

src/fit_spline.py

@@ -16,9 +16,10 @@
 parser.add_argument('-l', type=float, default=5, help='Smoothing parameter lambda')
 parser.add_argument('-s', type=str, default="[10000]", help="List of integers for the length of each voltage step e.g. '[1,2,3,4,5]'")
 parser.add_argument('-m', type=str, default='milnes', help='Milnes, opt, or real')
+parser.add_argument('-t', type=float, default=15e3, help='Max time')
 args = parser.parse_args()
 
-def main(input, output, lambda_, t_steps, data):
+def main(input, output, lambda_, t_steps, data, max_t):
     df_all = pd.read_csv(f"{input}/synth_Y.csv")
     # Parameters
     if data == 'milnes':
@@ -28,7 +29,7 @@ def main(input, output, lambda_, t_steps, data):
         j = 0
         t_swp = t_steps[j]
         t_total = 0
-        swps = sweeps*3
+        swps = int(np.floor(250000/max_t))*6
     else:
         lens = [3340, 3330, 3340, 3330, 3330]
         t_total = 0
@@ -78,4 +79,4 @@ def main(input, output, lambda_, t_steps, data):
 
 if __name__ == "__main__":
     t_steps = ast.literal_eval(args.s)
-    main(args.i, args.o, float(args.l), t_steps, args.m)
+    main(args.i, args.o, float(args.l), t_steps, args.m, args.t)

src/optimise_protocol.py

@@ -41,11 +41,11 @@ def __call__(self, p):
                 self.v_st, self.t_st = funcs.get_steps(v_steps, t_steps, self.pars)
                 prot = funcs.create_protocol(self.v_st, self.t_st)
             if alt_protocol is not None:
-                model_out, cont = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs,
+                model_out, cont, swps = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs,
                                           alt_protocol = prot_alt, alt_times = np.arange(0, np.sum(self.t_st_alt), 10), wins = [1e3, np.sum(self.t_st[1:4])],
                                           wins_alt = [1e3, np.sum(self.t_st_alt[1:4])])
             else:
-                model_out, cont = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs, wins = [1e3, np.sum(self.t_st[1:4])])
+                model_out, cont, swps = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs, wins = [1e3, np.sum(self.t_st[1:4])])
             ### loop through models and concentrations to get traces
             all_traces = []
             for m in model_out:
@@ -66,6 +66,7 @@ def __call__(self, p):
                         lhoods.append(-2*np.sum(np.log(top/bottom)))
             #out = -np.median(ssq)
             out = np.median(lhoods)
+            print(swps)
             return out
 
         def n_parameters(self):
@@ -126,7 +127,7 @@ def n_parameters(self):
 
     opt.optimiser().set_population_size(CMAES_pop)
     opt.set_max_iterations(max_iter)
-    opt.set_max_unchanged_iterations(iterations=100, threshold=1e-2)
+    opt.set_max_unchanged_iterations(iterations=50, threshold=1)
     opt.set_parallel(-1)
 
     try:
@@ -165,7 +166,7 @@ def main(model_nums, max_time, bounds, herg, output_folder):
         drug_fit_pars[m] = parlist
 
     # perform optimisation
-    p_out, cost = get_opt_prot(drug_fit_pars, herg, v_steps, t_steps, p0_1, CMAES_pop = 7, max_iter = 750, alt_protocol = p0_2)
+    p_out, cost = get_opt_prot(drug_fit_pars, herg, v_steps, t_steps, p0_1, CMAES_pop = 7, max_iter = 500, alt_protocol = p0_2)
     print(f'Final objective cost: {cost}')
     print(f'Final optimised params: {p_out}')
 

src/post_fit_plot.py

@@ -35,7 +35,7 @@ def main(model_nums, prot, max_time, bounds, herg, output_folder):
     if max_time == 15e3 and herg != '2024_Joey_sis_25C':
         swps = sweeps
     elif herg != '2024_Joey_sis_25C':
-        swps = 5
+        swps = int(np.floor(250000/max_time))
     else:
         swps = 10