diff --git a/__pycache__/setup.cpython-39.pyc b/__pycache__/setup.cpython-39.pyc new file mode 100644 index 0000000..5c5dafb Binary files /dev/null and b/__pycache__/setup.cpython-39.pyc differ diff --git a/build/lib/fears/data/ogbunugafor_drugless.csv b/build/lib/fears/data/ogbunugafor_drugless.csv index 9e809e8..ef121f6 100644 --- a/build/lib/fears/data/ogbunugafor_drugless.csv +++ b/build/lib/fears/data/ogbunugafor_drugless.csv @@ -1 +1 @@ -1.398,1.275,1.227,0,1.37,1.375,1.397,1.219,1.119,1.184,1.306,1,1.273,1.282,1.45,1.25 +1.398,1.275,1.227,0,1.37,1.375,1.397,1.219,1.119,1.184,1.306,1,1.273,1.282,1.45,1.25 diff --git a/build/lib/fears/data/pyrimethamine_ic50.csv b/build/lib/fears/data/pyrimethamine_ic50.csv index e5f7c2b..7237556 100644 --- a/build/lib/fears/data/pyrimethamine_ic50.csv +++ b/build/lib/fears/data/pyrimethamine_ic50.csv @@ -1 +1 @@ --6.286,-5.812,-4.239,0,-6.046,-5.774,-3.732,-3.55,-5.724,-5.491,-4.015,-4.6,-5.773,-5.624,-3.587,-3.3 +-6.286,-5.812,-4.239,0,-6.046,-5.774,-3.732,-3.55,-5.724,-5.491,-4.015,-4.6,-5.773,-5.624,-3.587,-3.3 diff --git a/build/lib/fears/experiment.py b/build/lib/fears/experiment.py index c80c4c8..4248891 100644 --- a/build/lib/fears/experiment.py +++ b/build/lib/fears/experiment.py @@ -1,586 +1,1161 @@ -from fears.population import Population -from fears.utils import plotter, dir_manager, fitness -import numpy as np -import pandas as pd -import os -import time -import pickle -import lifelines - -class Experiment(): - - # Initializer - def __init__(self, - n_sims = 1, - curve_types = None, - max_doses = None, - first_dose = None, - third_dose = None, - second_dose = None, - ramp = 100, - fitness_data = 'generate', - n_allele = 4, - null_seascape_dose=None, - transition_times = None, - inoculants = None, - experiment_type = None, - prob_drops = None, - n_impulse=1, - population_options = {}, - slopes=None, - passage=False, - passage_time = 48, - debug = True): # debug = True -> no save - - self.root_path = str(dir_manager.get_project_root()) - - # list of allowed drug curve types - allowed_types = ['linear', - 'constant', - 'heaviside', - 'pharm', - 'pulsed'] - - # list of defined experiments - allowed_experiments = ['inoculant-survival', - 'dose-survival', - 'drug-regimen', - 'dose-entropy', - 'rate-survival', - 'bottleneck', - 'ramp_up_down'] - - if curve_types is not None: - if not all(elem in allowed_types for elem in curve_types): - raise Exception('One or more curve types is not recognized.\nAllowable types are: linear, constant, heaviside, pharm, pulsed.') - - if experiment_type is not None: - if experiment_type not in allowed_experiments: - raise Exception('Experiment type not recognized.\nAllowable types are inoculant-survival, dose-survival, drug-regimen, dose-entropy, and bottleneck.') - - # Curve type: linear, constant, heaviside, pharm, pulsed - if curve_types is None: - self.curve_types = ['constant'] - else: - self.curve_types = curve_types - - if max_doses is None: - self.max_doses = [1] - else: - self.max_doses = max_doses - - if inoculants is None: - self.inoculants = [0] - else: - self.inoculants = inoculants - - self.n_sims = n_sims - - # Common options that are applied to each population - self.population_options = population_options - - # initialize all populations - self.populations = [] - - # initialize a list of figures for saving - self.figures = [] - - if experiment_type is None: - self.experiment_type = 'dose-survival' - # warnings.warn('No experiment type given - set to dose-survival by default.') - else: - self.experiment_type = experiment_type - - # if experiment_type == 'dose-survival' and len(inoculants) > 1: - # raise Exception('The experiment type is set to dose-survival (default), but more than one inoculant is given.') - # elif experiment_type == 'inoculant-survival' and len(self.max_doses) > 1: - # # print('here') - # raise Exception('The experiment type is set to inoculant-survival, but more than one max dose is given.') - - # if experiment_type == 'inoculant-survival' and inoculants is None: - # raise Exception('The experiment type is set to inoculant-survival, but no inoculants are given.') - self.ramp = ramp - if self.experiment_type == 'ramp_up_down': - self.p_landscape = Population(constant_pop = True, - carrying_cap = False, - fitness_data=fitness_data, - n_allele=n_allele, - **self.population_options) - self.p_seascape = Population(constant_pop = True, - carrying_cap = False, - fitness_data=fitness_data, - n_allele=n_allele, - **self.population_options) - - if fitness_data=='random': - self.p_seascape.ic50 = self.p_landscape.ic50 - self.p_seascape.drugless_rates = self.p_landscape.drugless_rates - self.p_landscape.set_null_seascape(null_seascape_dose) - self.null_seascape_dose = null_seascape_dose - - if second_dose is None: - self.second_dose = 10**5 - else: - self.second_dose=second_dose - if third_dose is None: - self.third_dose = 10**1 - else: - self.third_dose=third_dose - if first_dose is None: - self.first_dose = 10**-2 - else: - self.first_dose = first_dose - if transition_times is None: - n_timestep = self.p_landscape.n_timestep - self.transition_times = [int(n_timestep/4),int(3*n_timestep/4)] - else: - self.transition_times = transition_times - - self.set_ramp_ud(self.p_landscape) - self.set_ramp_ud(self.p_seascape) - - - if self.experiment_type == 'bottleneck': - for dc in self.duty_cycles: - self.populations.append(Population(curve_type='on_off', - duty_cycle=dc, - n_sims = self.n_sims, - **self.population_options)) - # p_bottleneck = Population(max_dose = max_doses, - # n_sims = self.n_sims, - # curve_type = 'bottleneck' - # **self.population_options) - # self.populations.append(p_bottleneck) - - # p_no_bottleneck = Population(max_dose = max_doses, - # n_sims = self.n_sims, - # curve_type = 'constant' - # **self.population_options) - # self.populations.append(p_no_bottleneck) - - if self.experiment_type == 'dose-survival': - for curve_type in self.curve_types: - for max_dose in self.max_doses: - fig_title = 'Max dose = ' + str(max_dose) + ', curve type = ' + curve_type - self.populations.append(Population(curve_type=curve_type, - n_sims = self.n_sims, - max_dose = max_dose, - fig_title = fig_title, - **self.population_options)) - - self.n_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) - self.perc_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) - - elif self.experiment_type == 'inoculant-survival': - for curve_type in self.curve_types: - for inoculant in self.inoculants: - fig_title = 'Inoculant = ' + str(inoculant) + ', curve type = ' + curve_type - - init_counts = np.zeros(16) - init_counts[0] = inoculant - - self.populations.append(Population(curve_type=curve_type, - n_sims = self.n_sims, - fig_title = fig_title, - init_counts=init_counts, - **self.population_options)) - - self.n_survive = np.zeros([len(self.curve_types),len(self.inoculants)]) - self.perc_survive = np.zeros([len(self.curve_types),len(self.inoculants)]) - - elif self.experiment_type == 'drug-regimen': - - self.prob_drops = prob_drops - - for prob_drop in self.prob_drops: - curve_type = 'pulsed' - self.populations.append(Population(curve_type=curve_type, - prob_drop=prob_drop, - # n_impulse = 1, - n_sims = 1, - # fig_title = fig_title, - # init_counts=init_counts, - - **self.population_options)) - self.n_survive = np.zeros([len(self.populations)]) - - elif self.experiment_type == 'dose-entropy': - for dose in self.max_doses: - self.populations.append(Population(max_dose = dose, - curve_type=self.curve_types[0])) - self.entropy_results = pd.DataFrame(columns=[]) # will become a dataframe later - - elif self.experiment_type == 'rate-survival': - # if the curve type is 'pharm' then slope will be interpreted as k_abs - self.slopes = slopes - for slope in self.slopes: - if curve_types[0] == 'pharm': - self.populations.append(Population(max_dose=self.max_doses[0], - k_abs=slope, - curve_type='pharm', - n_sims=1, - passage=passage, - passage_time=passage_time, - **self.population_options)) - elif curve_types[0] == 'pulsed': - self.populations.append(Population(max_dose=self.max_doses[0], - k_abs=slope, - curve_type='pulsed', - n_sims=1, - passage=passage, - passage_time=passage_time, - **self.population_options)) - else: - self.populations.append(Population(max_dose=self.max_doses[0], - slope=slope, - curve_type='linear', - n_sims=1, - passage=passage, - passage_time=passage_time, - **self.population_options)) - - # self.rate_survival_results = pd.DataFrame(columns=[]) - - # generate new save folder - - self.debug=debug - - if not debug: - - num = 0 - num_str = str(num).zfill(4) - - date_str = time.strftime('%m%d%Y',time.localtime()) - - save_folder = self.root_path + os.sep + 'results' + os.sep + 'results_' + date_str + '_' + num_str - - # save_folder = os.getcwd() + '//results_' + date_str + '_' + num_str - if not (os.path.exists(self.root_path + os.sep + 'results')): - os.mkdir(self.root_path + os.sep + 'results') - - while(os.path.exists(save_folder)): - num += 1 - num_str = str(num).zfill(4) - save_folder = self.root_path + os.sep + 'results' + os.sep + 'results_' + date_str + '_' + num_str - os.mkdir(save_folder) - - self.results_path = save_folder - # self.experiment_info_path = self.root_path + os.sep + 'results' + os.sep + 'experiment_info_' + date_str + '_' + num_str + '.p' - self.experiment_info_path = self.results_path + os.sep + 'experiment_info_' + date_str + '_' + num_str + '.p' - self.exp_folders = [] - - # self.savename = None - - # self.n_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) - # self.perc_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) -############################################################################### - # Methods for running experiments - - # run experiment and save results - def run_experiment(self): - - n_doses = len(self.max_doses) - n_curves = len(self.curve_types) - n_inoc = len(self.inoculants) - - # pbar = tqdm(total = n_curves*n_doses) # progress bar - - # Loop through each population, execute simulations, and store survival statistics - - if self.experiment_type == 'dose-survival': - # pbar = tqdm(total = n_curves*n_doses) # progress bar - for curve_number in range(n_curves): - for dose_number in range(n_doses): - - exp_num = curve_number*n_doses + dose_number - pop = self.populations[exp_num] # extract population in list of population - c,n_survive_t = pop.simulate() - pop.plot_timecourse() - self.n_survive[curve_number,dose_number] = n_survive_t - # pbar.update() - self.perc_survive = 100*self.n_survive/self.n_sims - - elif self.experiment_type == 'ramp_up_down': - counts_landscape, ft = self.p_landscape.simulate() - drug_curve = self.p_landscape.drug_curve - drug_curve= np.array([drug_curve]) - drug_curve = np.transpose(drug_curve) - counts_seascape, ft = self.p_seascape.simulate() - - if not self.debug: - savedata = np.concatenate((counts_landscape,drug_curve),axis=1) - self.save_counts(savedata, num=None, save_folder=None,prefix = 'landscape_counts') - - savedata = np.concatenate((counts_seascape,drug_curve),axis=1) - self.save_counts(savedata, num=None, save_folder=None,prefix = 'seascape_counts') - - elif self.experiment_type == 'inoculant-survival': - # pbar = tqdm(total = n_curves*n_inoc) # progress bar - for curve_number in range(n_curves): - for inoc_num in range(n_inoc): - - exp_num = curve_number*n_inoc + inoc_num - pop = self.populations[exp_num] # extract population in list of population - c,n_survive_t = pop.simulate() - pop.plot_timecourse() - self.n_survive[curve_number,inoc_num] = n_survive_t - # pbar.update() - self.perc_survive = 100*self.n_survive/self.n_sims - - elif self.experiment_type == 'drug-regimen': - # pbar = tqdm(total=len(self.populations)) - # kk=0 - for p in self.populations: - save_folder = 'p_drop=' + str(p.prob_drop) - save_folder = save_folder.replace('.',',') - # self.exp_folder.append(save_folder) - for i in range(self.n_sims): - # initialize new drug curve - p.drug_curve,u = p.gen_curves() - - counts,n_survive = p.simulate() - drug = p.drug_curve - drug = np.array([drug]) - drug = np.transpose(drug) - - # regimen = self.compute_regimen(p, u) - # regimen = np.array([regimen]) - # regimen = np.transpose(regimen) - # regimen_t = np.zeros(drug.shape) - # regimen_t[0:len(regimen)] = regimen - - u = np.array([u,]) - u = u.transpose() - counts = np.concatenate((counts,drug,u),axis=1) - - if not self.debug: - # self.save_counts(counts,i,save_folder) - data_dict = {'counts':counts, - 'drug_curve':drug, - 'regimen':u} - self.save_dict(data_dict,i,save_folder) - # kk+=1 - # pbar.update() - self.perc_survive = 100*self.n_survive/self.n_sims - - elif self.experiment_type == 'dose-entropy': - # pbar = tqdm(total=len(self.populations)*self.n_sims) - e_survived = [] - e_died = [] - for p in self.populations: - - for i in range(self.n_sims): - c,n_survive = p.simulate() - # e = max(p.entropy()) # compute max entropy - e_t = p.entropy() - e = max(e_t) - # e=1 - # p.plot_timecourse() - - if n_survive == 1: - survive = 'survived' # survived - e_survived.append(e_t) - else: - survive = 'extinct' # died - e_died.append(e_t) - - d = {'dose':[p.max_dose], - 'survive condition':[survive], - 'max entropy':[e]} - - entropy_results_t = pd.DataFrame(d) - self.entropy_results = self.entropy_results.append(entropy_results_t) - # pbar.update() - - elif self.experiment_type == 'rate-survival': - # pbar = tqdm(total=len(self.populations)) - - for p in self.populations: - - for n in range(self.n_sims): - counts,n_survive = p.simulate() - - drug = p.drug_curve - drug = np.array([drug]) - drug = np.transpose(drug) - # counts = np.concatenate((counts,drug),axis=1) - - if self.debug is False: - if (self.curve_types[0] == 'pharm' or - self.curve_types[0] == 'pulsed'): - save_folder = 'k_abs=' + str(p.k_abs) - save_folder.replace('.',',') - else: - save_folder = 'slope=' + str(p.slope) - save_folder.replace('.',',') - self.save_counts(counts,n,save_folder) - data_dict = {'counts':counts, - 'drug_curve':drug} - self.save_dict(data_dict,n,save_folder) - # pickle.dump(self, open(self.experiment_info_path,"wb")) - - - # fig_savename = 'slope = ' + str(p.slope) - # self.figures = self.figures.append(fig) - # pbar.update() - # self.rate_survival_results.index = np.arange(len(self.rate_survival_results)) - - # pbar.close() # close progress bar - - # save counts as a csv in the given subfolder with the label 'num' - def save_counts(self,counts,num,save_folder,prefix='sim_'): - - # check if the desired save folder exists. If not, create it - if save_folder is None: - save_folder = '' - folder_path = self.results_path + os.sep + save_folder - if os.path.exists(folder_path) != True: - os.mkdir(folder_path) - - if num is None: - num = '' - else: - num = str(num).zfill(4) - - savename = self.results_path + os.sep + save_folder + os.sep + prefix + num + '.csv' - np.savetxt(savename, counts, delimiter=",") - # self.savename = savename - return - - def save_dict(self,data_dict,num,save_folder,prefix='sim_'): - # check if the desired save folder exists. If not, create it - if save_folder is None: - save_folder = '' - folder_path = self.results_path + os.sep + save_folder - if os.path.exists(folder_path) != True: - os.mkdir(folder_path) - - if num is None: - num = '' - else: - num = str(num).zfill(4) - - if folder_path not in self.exp_folders: - self.exp_folders.append(folder_path) - - savename = self.results_path + os.sep + save_folder + os.sep + prefix + num + '.p' - pickle.dump(data_dict, open(savename,"wb")) - # np.savetxt(savename, counts, delimiter=",") - return - - def compute_regimen(self,p,u): - gap = int(p.dose_schedule/p.timestep_scale) - n_impulse = int(np.ceil(p.n_timestep/gap)) - regimen = np.zeros(n_impulse) - - for i in range(n_impulse): - if u[(i)*gap] == 1: - regimen[i] = 1 - - return regimen - - def set_ramp_ud(self,p): - - n_timestep = p.n_timestep - drug_curve = np.zeros(n_timestep) - slope = (self.second_dose-self.first_dose)/self.ramp - buffer = self.ramp/2 - - times = [self.transition_times[0]-buffer, - self.transition_times[0]+buffer, - self.transition_times[1]-buffer, - self.transition_times[1]+buffer] - - for t in range(n_timestep): - if t 1: - c = np.sum(counts,axis=1) - else: - c = counts - e = np.argwhere(c 1: - c = counts[:,genotype] - else: - c = counts - - if thresh < 1: - thresh = thresh*pop.max_cells - - e = np.argwhere(c>thresh) - if len(e) == 0: - event_obs = 0 - event_time = len(c) - else: - event_obs = 1 - event_time = e[0] - - timestep_scale = pop.timestep_scale - event_time = event_time*timestep_scale - - return event_obs, event_time - - def log_rank_test(self,durations_A, durations_B, - event_observed_A=None, event_observed_B=None): - - results = lifelines.statistics.logrank_test(durations_A, durations_B, - event_observed_A=event_observed_A, - event_observed_B=event_observed_B) - - return results - \ No newline at end of file +<<<<<<< HEAD:fears/classes/experiment_class.py +from fears.classes import Population +from fears.utils.plotter import * +from fears.utils.dir_manager import * +import numpy as np +import pandas as pd +# import warnings +import os +import time +import pickle +import lifelines + +class Experiment(): + + # Initializer + def __init__(self, + n_sims = 1, + curve_types = None, + max_doses = None, + first_dose = None, + third_dose = None, + second_dose = None, + ramp = 100, + fitness_data = 'generate', + n_allele = 4, + null_seascape_dose=None, + transition_times = None, + inoculants = None, + experiment_type = None, + prob_drops = None, + n_impulse=1, + population_options = {}, + slopes=None, + passage=False, + passage_time = 48, + debug = True): # debug = True -> no save + + self.root_path = str(dir_manager.get_project_root()) + + # list of allowed drug curve types + allowed_types = ['linear', + 'constant', + 'heaviside', + 'pharm', + 'pulsed'] + + # list of defined experiments + allowed_experiments = ['inoculant-survival', + 'dose-survival', + 'drug-regimen', + 'dose-entropy', + 'rate-survival', + 'bottleneck', + 'ramp_up_down'] + + if curve_types is not None: + if not all(elem in allowed_types for elem in curve_types): + raise Exception('One or more curve types is not recognized.\nAllowable types are: linear, constant, heaviside, pharm, pulsed.') + + if experiment_type is not None: + if experiment_type not in allowed_experiments: + raise Exception('Experiment type not recognized.\nAllowable types are inoculant-survival, dose-survival, drug-regimen, dose-entropy, and bottleneck.') + + # Curve type: linear, constant, heaviside, pharm, pulsed + if curve_types is None: + self.curve_types = ['constant'] + else: + self.curve_types = curve_types + + if max_doses is None: + self.max_doses = [1] + else: + self.max_doses = max_doses + + if inoculants is None: + self.inoculants = [0] + else: + self.inoculants = inoculants + + self.n_sims = n_sims + + # Common options that are applied to each population + self.population_options = population_options + + # initialize all populations + self.populations = [] + + # initialize a list of figures for saving + self.figures = [] + + if experiment_type is None: + self.experiment_type = 'dose-survival' + # warnings.warn('No experiment type given - set to dose-survival by default.') + else: + self.experiment_type = experiment_type + + # if experiment_type == 'dose-survival' and len(inoculants) > 1: + # raise Exception('The experiment type is set to dose-survival (default), but more than one inoculant is given.') + # elif experiment_type == 'inoculant-survival' and len(self.max_doses) > 1: + # # print('here') + # raise Exception('The experiment type is set to inoculant-survival, but more than one max dose is given.') + + # if experiment_type == 'inoculant-survival' and inoculants is None: + # raise Exception('The experiment type is set to inoculant-survival, but no inoculants are given.') + self.ramp = ramp + if self.experiment_type == 'ramp_up_down': + self.p_landscape = Population(constant_pop = True, + carrying_cap = False, + fitness_data=fitness_data, + n_allele=n_allele, + **self.population_options) + self.p_seascape = Population(constant_pop = True, + carrying_cap = False, + fitness_data=fitness_data, + n_allele=n_allele, + **self.population_options) + + if fitness_data=='random': + self.p_seascape.ic50 = self.p_landscape.ic50 + self.p_seascape.drugless_rates = self.p_landscape.drugless_rates + self.p_landscape.set_null_seascape(null_seascape_dose) + self.null_seascape_dose = null_seascape_dose + + if second_dose is None: + self.second_dose = 10**5 + else: + self.second_dose=second_dose + if third_dose is None: + self.third_dose = 10**1 + else: + self.third_dose=third_dose + if first_dose is None: + self.first_dose = 10**-2 + else: + self.first_dose = first_dose + if transition_times is None: + n_timestep = self.p_landscape.n_timestep + self.transition_times = [int(n_timestep/4),int(3*n_timestep/4)] + else: + self.transition_times = transition_times + + self.set_ramp_ud(self.p_landscape) + self.set_ramp_ud(self.p_seascape) + + + if self.experiment_type == 'bottleneck': + for dc in self.duty_cycles: + self.populations.append(Population(curve_type='on_off', + duty_cycle=dc, + n_sims = self.n_sims, + **self.population_options)) + # p_bottleneck = Population(max_dose = max_doses, + # n_sims = self.n_sims, + # curve_type = 'bottleneck' + # **self.population_options) + # self.populations.append(p_bottleneck) + + # p_no_bottleneck = Population(max_dose = max_doses, + # n_sims = self.n_sims, + # curve_type = 'constant' + # **self.population_options) + # self.populations.append(p_no_bottleneck) + + if self.experiment_type == 'dose-survival': + for curve_type in self.curve_types: + for max_dose in self.max_doses: + fig_title = 'Max dose = ' + str(max_dose) + ', curve type = ' + curve_type + self.populations.append(Population(curve_type=curve_type, + n_sims = self.n_sims, + max_dose = max_dose, + fig_title = fig_title, + **self.population_options)) + + self.n_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) + self.perc_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) + + elif self.experiment_type == 'inoculant-survival': + for curve_type in self.curve_types: + for inoculant in self.inoculants: + fig_title = 'Inoculant = ' + str(inoculant) + ', curve type = ' + curve_type + + init_counts = np.zeros(16) + init_counts[0] = inoculant + + self.populations.append(Population(curve_type=curve_type, + n_sims = self.n_sims, + fig_title = fig_title, + init_counts=init_counts, + **self.population_options)) + + self.n_survive = np.zeros([len(self.curve_types),len(self.inoculants)]) + self.perc_survive = np.zeros([len(self.curve_types),len(self.inoculants)]) + + elif self.experiment_type == 'drug-regimen': + + self.prob_drops = prob_drops + + for prob_drop in self.prob_drops: + curve_type = 'pulsed' + self.populations.append(Population(curve_type=curve_type, + prob_drop=prob_drop, + # n_impulse = 1, + n_sims = 1, + # fig_title = fig_title, + # init_counts=init_counts, + + **self.population_options)) + self.n_survive = np.zeros([len(self.populations)]) + + elif self.experiment_type == 'dose-entropy': + for dose in self.max_doses: + self.populations.append(Population(max_dose = dose, + curve_type=self.curve_types[0])) + self.entropy_results = pd.DataFrame(columns=[]) # will become a dataframe later + + elif self.experiment_type == 'rate-survival': + # if the curve type is 'pharm' then slope will be interpreted as k_abs + self.slopes = slopes + for slope in self.slopes: + if curve_types[0] == 'pharm': + self.populations.append(Population(max_dose=self.max_doses[0], + k_abs=slope, + curve_type='pharm', + n_sims=1, + passage=passage, + passage_time=passage_time, + **self.population_options)) + elif curve_types[0] == 'pulsed': + self.populations.append(Population(max_dose=self.max_doses[0], + k_abs=slope, + curve_type='pulsed', + n_sims=1, + passage=passage, + passage_time=passage_time, + **self.population_options)) + else: + self.populations.append(Population(max_dose=self.max_doses[0], + slope=slope, + curve_type='linear', + n_sims=1, + passage=passage, + passage_time=passage_time, + **self.population_options)) + + # self.rate_survival_results = pd.DataFrame(columns=[]) + + # generate new save folder + + self.debug=debug + + if not debug: + + num = 0 + num_str = str(num).zfill(4) + + date_str = time.strftime('%m%d%Y',time.localtime()) + + save_folder = self.root_path + os.sep + 'results' + os.sep + 'results_' + date_str + '_' + num_str + + # save_folder = os.getcwd() + '//results_' + date_str + '_' + num_str + if not (os.path.exists(self.root_path + os.sep + 'results')): + os.mkdir(self.root_path + os.sep + 'results') + + while(os.path.exists(save_folder)): + num += 1 + num_str = str(num).zfill(4) + save_folder = self.root_path + os.sep + 'results' + os.sep + 'results_' + date_str + '_' + num_str + os.mkdir(save_folder) + + self.experiment_info_path = self.root_path + os.sep + 'results' + os.sep + 'experiment_info_' + date_str + '_' + num_str + '.p' + + # self.experiment_info_path = experiment_info_path + self.results_path = save_folder + + # save the experiment information + + pickle.dump(self, open(self.experiment_info_path,"wb")) + + # self.n_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) + # self.perc_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) +############################################################################### + # Methods for running experiments + + # run experiment and save results + def run_experiment(self): + + n_doses = len(self.max_doses) + n_curves = len(self.curve_types) + n_inoc = len(self.inoculants) + + # pbar = tqdm(total = n_curves*n_doses) # progress bar + + # Loop through each population, execute simulations, and store survival statistics + + if self.experiment_type == 'dose-survival': + # pbar = tqdm(total = n_curves*n_doses) # progress bar + for curve_number in range(n_curves): + for dose_number in range(n_doses): + + exp_num = curve_number*n_doses + dose_number + pop = self.populations[exp_num] # extract population in list of population + c,n_survive_t = pop.simulate() + pop.plot_timecourse() + self.n_survive[curve_number,dose_number] = n_survive_t + # pbar.update() + self.perc_survive = 100*self.n_survive/self.n_sims + + elif self.experiment_type == 'ramp_up_down': + counts_landscape, ft = self.p_landscape.simulate() + drug_curve = self.p_landscape.drug_curve + drug_curve= np.array([drug_curve]) + drug_curve = np.transpose(drug_curve) + counts_seascape, ft = self.p_seascape.simulate() + + if not self.debug: + savedata = np.concatenate((counts_landscape,drug_curve),axis=1) + self.save_counts(savedata, num=None, save_folder=None,prefix = 'landscape_counts') + + savedata = np.concatenate((counts_seascape,drug_curve),axis=1) + self.save_counts(savedata, num=None, save_folder=None,prefix = 'seascape_counts') + + elif self.experiment_type == 'inoculant-survival': + # pbar = tqdm(total = n_curves*n_inoc) # progress bar + for curve_number in range(n_curves): + for inoc_num in range(n_inoc): + + exp_num = curve_number*n_inoc + inoc_num + pop = self.populations[exp_num] # extract population in list of population + c,n_survive_t = pop.simulate() + pop.plot_timecourse() + self.n_survive[curve_number,inoc_num] = n_survive_t + # pbar.update() + self.perc_survive = 100*self.n_survive/self.n_sims + + elif self.experiment_type == 'drug-regimen': + # pbar = tqdm(total=len(self.populations)) + # kk=0 + for p in self.populations: + for i in range(self.n_sims): + # initialize new drug curve + p.drug_curve,u = p.gen_curves() + + counts,n_survive = p.simulate() + drug = p.drug_curve + drug = np.array([drug]) + drug = np.transpose(drug) + + # regimen = self.compute_regimen(p, u) + # regimen = np.array([regimen]) + # regimen = np.transpose(regimen) + # regimen_t = np.zeros(drug.shape) + # regimen_t[0:len(regimen)] = regimen + + u = np.array([u,]) + u = u.transpose() + counts = np.concatenate((counts,drug,u),axis=1) + + save_folder = 'p_drop=' + str(p.prob_drop) + save_folder = save_folder.replace('.',',') + if not self.debug: + self.save_counts(counts,i,save_folder) + # kk+=1 + # pbar.update() + self.perc_survive = 100*self.n_survive/self.n_sims + + elif self.experiment_type == 'dose-entropy': + # pbar = tqdm(total=len(self.populations)*self.n_sims) + e_survived = [] + e_died = [] + for p in self.populations: + + for i in range(self.n_sims): + c,n_survive = p.simulate() + # e = max(p.entropy()) # compute max entropy + e_t = p.entropy() + e = max(e_t) + # e=1 + # p.plot_timecourse() + + if n_survive == 1: + survive = 'survived' # survived + e_survived.append(e_t) + else: + survive = 'extinct' # died + e_died.append(e_t) + + d = {'dose':[p.max_dose], + 'survive condition':[survive], + 'max entropy':[e]} + + entropy_results_t = pd.DataFrame(d) + self.entropy_results = self.entropy_results.append(entropy_results_t) + # pbar.update() + + elif self.experiment_type == 'rate-survival': + # pbar = tqdm(total=len(self.populations)) + + for p in self.populations: + + for n in range(self.n_sims): + counts,n_survive = p.simulate() + + drug = p.drug_curve + drug = np.array([drug]) + drug = np.transpose(drug) + counts = np.concatenate((counts,drug),axis=1) + + if self.debug is False: + if (self.curve_types[0] == 'pharm' or + self.curve_types[0] == 'pulsed'): + save_folder = 'k_abs=' + str(p.k_abs) + save_folder.replace('.','pnt') + else: + save_folder = 'slope=' + str(p.slope) + save_folder.replace('.','pnt') + self.save_counts(counts,n,save_folder) + + + + # fig_savename = 'slope = ' + str(p.slope) + # self.figures = self.figures.append(fig) + # pbar.update() + # self.rate_survival_results.index = np.arange(len(self.rate_survival_results)) + + # pbar.close() # close progress bar + + # save counts as a csv in the given subfolder with the label 'num' + def save_counts(self,counts,num,save_folder,prefix='sim_'): + + # check if the desired save folder exists. If not, create it + if save_folder is None: + save_folder = '' + folder_path = self.results_path + os.sep + save_folder + if os.path.exists(folder_path) != True: + os.mkdir(folder_path) + + if num is None: + num = '' + else: + num = str(num).zfill(4) + + savename = self.results_path + os.sep + save_folder + os.sep + prefix + num + '.csv' + np.savetxt(savename, counts, delimiter=",") + return + + def compute_regimen(self,p,u): + gap = int(p.dose_schedule/p.timestep_scale) + n_impulse = int(np.ceil(p.n_timestep/gap)) + regimen = np.zeros(n_impulse) + + for i in range(n_impulse): + if u[(i)*gap] == 1: + regimen[i] = 1 + + return regimen + + def set_ramp_ud(self,p): + + n_timestep = p.n_timestep + drug_curve = np.zeros(n_timestep) + slope = (self.second_dose-self.first_dose)/self.ramp + buffer = self.ramp/2 + + times = [self.transition_times[0]-buffer, + self.transition_times[0]+buffer, + self.transition_times[1]-buffer, + self.transition_times[1]+buffer] + + for t in range(n_timestep): + if t 1: + c = np.sum(counts,axis=1) + else: + c = counts + e = np.argwhere(c 1: + c = counts[:,genotype] + else: + c = counts + + if thresh < 1: + thresh = thresh*pop.max_cells + + e = np.argwhere(c>thresh) + if len(e) == 0: + event_obs = 0 + event_time = len(c) + else: + event_obs = 1 + event_time = e[0] + + timestep_scale = pop.timestep_scale + event_time = event_time*timestep_scale + + return event_obs, event_time + + def log_rank_test(self,durations_A, durations_B, + event_observed_A=None, event_observed_B=None): + + results = lifelines.statistics.logrank_test(durations_A, durations_B, + event_observed_A=event_observed_A, + event_observed_B=event_observed_B) + + return results + +======= +from fears.population import Population +from fears.utils import plotter, dir_manager, fitness +import numpy as np +import pandas as pd +import os +import time +import pickle +import lifelines + +class Experiment(): + + # Initializer + def __init__(self, + n_sims = 1, + curve_types = None, + max_doses = None, + first_dose = None, + third_dose = None, + second_dose = None, + ramp = 100, + fitness_data = 'generate', + n_allele = 4, + null_seascape_dose=None, + transition_times = None, + inoculants = None, + experiment_type = None, + prob_drops = None, + n_impulse=1, + population_options = {}, + slopes=None, + passage=False, + passage_time = 48, + debug = True): # debug = True -> no save + + self.root_path = str(dir_manager.get_project_root()) + + # list of allowed drug curve types + allowed_types = ['linear', + 'constant', + 'heaviside', + 'pharm', + 'pulsed'] + + # list of defined experiments + allowed_experiments = ['inoculant-survival', + 'dose-survival', + 'drug-regimen', + 'dose-entropy', + 'rate-survival', + 'bottleneck', + 'ramp_up_down'] + + if curve_types is not None: + if not all(elem in allowed_types for elem in curve_types): + raise Exception('One or more curve types is not recognized.\nAllowable types are: linear, constant, heaviside, pharm, pulsed.') + + if experiment_type is not None: + if experiment_type not in allowed_experiments: + raise Exception('Experiment type not recognized.\nAllowable types are inoculant-survival, dose-survival, drug-regimen, dose-entropy, and bottleneck.') + + # Curve type: linear, constant, heaviside, pharm, pulsed + if curve_types is None: + self.curve_types = ['constant'] + else: + self.curve_types = curve_types + + if max_doses is None: + self.max_doses = [1] + else: + self.max_doses = max_doses + + if inoculants is None: + self.inoculants = [0] + else: + self.inoculants = inoculants + + self.n_sims = n_sims + + # Common options that are applied to each population + self.population_options = population_options + + # initialize all populations + self.populations = [] + + # initialize a list of figures for saving + self.figures = [] + + if experiment_type is None: + self.experiment_type = 'dose-survival' + # warnings.warn('No experiment type given - set to dose-survival by default.') + else: + self.experiment_type = experiment_type + + # if experiment_type == 'dose-survival' and len(inoculants) > 1: + # raise Exception('The experiment type is set to dose-survival (default), but more than one inoculant is given.') + # elif experiment_type == 'inoculant-survival' and len(self.max_doses) > 1: + # # print('here') + # raise Exception('The experiment type is set to inoculant-survival, but more than one max dose is given.') + + # if experiment_type == 'inoculant-survival' and inoculants is None: + # raise Exception('The experiment type is set to inoculant-survival, but no inoculants are given.') + self.ramp = ramp + if self.experiment_type == 'ramp_up_down': + self.p_landscape = Population(constant_pop = True, + carrying_cap = False, + fitness_data=fitness_data, + n_allele=n_allele, + **self.population_options) + self.p_seascape = Population(constant_pop = True, + carrying_cap = False, + fitness_data=fitness_data, + n_allele=n_allele, + **self.population_options) + + if fitness_data=='random': + self.p_seascape.ic50 = self.p_landscape.ic50 + self.p_seascape.drugless_rates = self.p_landscape.drugless_rates + self.p_landscape.set_null_seascape(null_seascape_dose) + self.null_seascape_dose = null_seascape_dose + + if second_dose is None: + self.second_dose = 10**5 + else: + self.second_dose=second_dose + if third_dose is None: + self.third_dose = 10**1 + else: + self.third_dose=third_dose + if first_dose is None: + self.first_dose = 10**-2 + else: + self.first_dose = first_dose + if transition_times is None: + n_timestep = self.p_landscape.n_timestep + self.transition_times = [int(n_timestep/4),int(3*n_timestep/4)] + else: + self.transition_times = transition_times + + self.set_ramp_ud(self.p_landscape) + self.set_ramp_ud(self.p_seascape) + + + if self.experiment_type == 'bottleneck': + for dc in self.duty_cycles: + self.populations.append(Population(curve_type='on_off', + duty_cycle=dc, + n_sims = self.n_sims, + **self.population_options)) + # p_bottleneck = Population(max_dose = max_doses, + # n_sims = self.n_sims, + # curve_type = 'bottleneck' + # **self.population_options) + # self.populations.append(p_bottleneck) + + # p_no_bottleneck = Population(max_dose = max_doses, + # n_sims = self.n_sims, + # curve_type = 'constant' + # **self.population_options) + # self.populations.append(p_no_bottleneck) + + if self.experiment_type == 'dose-survival': + for curve_type in self.curve_types: + for max_dose in self.max_doses: + fig_title = 'Max dose = ' + str(max_dose) + ', curve type = ' + curve_type + self.populations.append(Population(curve_type=curve_type, + n_sims = self.n_sims, + max_dose = max_dose, + fig_title = fig_title, + **self.population_options)) + + self.n_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) + self.perc_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) + + elif self.experiment_type == 'inoculant-survival': + for curve_type in self.curve_types: + for inoculant in self.inoculants: + fig_title = 'Inoculant = ' + str(inoculant) + ', curve type = ' + curve_type + + init_counts = np.zeros(16) + init_counts[0] = inoculant + + self.populations.append(Population(curve_type=curve_type, + n_sims = self.n_sims, + fig_title = fig_title, + init_counts=init_counts, + **self.population_options)) + + self.n_survive = np.zeros([len(self.curve_types),len(self.inoculants)]) + self.perc_survive = np.zeros([len(self.curve_types),len(self.inoculants)]) + + elif self.experiment_type == 'drug-regimen': + + self.prob_drops = prob_drops + + for prob_drop in self.prob_drops: + curve_type = 'pulsed' + self.populations.append(Population(curve_type=curve_type, + prob_drop=prob_drop, + # n_impulse = 1, + n_sims = 1, + # fig_title = fig_title, + # init_counts=init_counts, + + **self.population_options)) + self.n_survive = np.zeros([len(self.populations)]) + + elif self.experiment_type == 'dose-entropy': + for dose in self.max_doses: + self.populations.append(Population(max_dose = dose, + curve_type=self.curve_types[0])) + self.entropy_results = pd.DataFrame(columns=[]) # will become a dataframe later + + elif self.experiment_type == 'rate-survival': + # if the curve type is 'pharm' then slope will be interpreted as k_abs + self.slopes = slopes + for slope in self.slopes: + if curve_types[0] == 'pharm': + self.populations.append(Population(max_dose=self.max_doses[0], + k_abs=slope, + curve_type='pharm', + n_sims=1, + passage=passage, + passage_time=passage_time, + **self.population_options)) + elif curve_types[0] == 'pulsed': + self.populations.append(Population(max_dose=self.max_doses[0], + k_abs=slope, + curve_type='pulsed', + n_sims=1, + passage=passage, + passage_time=passage_time, + **self.population_options)) + else: + self.populations.append(Population(max_dose=self.max_doses[0], + slope=slope, + curve_type='linear', + n_sims=1, + passage=passage, + passage_time=passage_time, + **self.population_options)) + + # self.rate_survival_results = pd.DataFrame(columns=[]) + + # generate new save folder + + self.debug=debug + + if not debug: + + num = 0 + num_str = str(num).zfill(4) + + date_str = time.strftime('%m%d%Y',time.localtime()) + + save_folder = self.root_path + os.sep + 'results' + os.sep + 'results_' + date_str + '_' + num_str + + # save_folder = os.getcwd() + '//results_' + date_str + '_' + num_str + if not (os.path.exists(self.root_path + os.sep + 'results')): + os.mkdir(self.root_path + os.sep + 'results') + + while(os.path.exists(save_folder)): + num += 1 + num_str = str(num).zfill(4) + save_folder = self.root_path + os.sep + 'results' + os.sep + 'results_' + date_str + '_' + num_str + os.mkdir(save_folder) + + self.results_path = save_folder + # self.experiment_info_path = self.root_path + os.sep + 'results' + os.sep + 'experiment_info_' + date_str + '_' + num_str + '.p' + self.experiment_info_path = self.results_path + os.sep + 'experiment_info_' + date_str + '_' + num_str + '.p' + self.exp_folders = [] + + # self.savename = None + + # self.n_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) + # self.perc_survive = np.zeros([len(self.curve_types),len(self.max_doses)]) +############################################################################### + # Methods for running experiments + + # run experiment and save results + def run_experiment(self): + + n_doses = len(self.max_doses) + n_curves = len(self.curve_types) + n_inoc = len(self.inoculants) + + # pbar = tqdm(total = n_curves*n_doses) # progress bar + + # Loop through each population, execute simulations, and store survival statistics + + if self.experiment_type == 'dose-survival': + # pbar = tqdm(total = n_curves*n_doses) # progress bar + for curve_number in range(n_curves): + for dose_number in range(n_doses): + + exp_num = curve_number*n_doses + dose_number + pop = self.populations[exp_num] # extract population in list of population + c,n_survive_t = pop.simulate() + pop.plot_timecourse() + self.n_survive[curve_number,dose_number] = n_survive_t + # pbar.update() + self.perc_survive = 100*self.n_survive/self.n_sims + + elif self.experiment_type == 'ramp_up_down': + counts_landscape, ft = self.p_landscape.simulate() + drug_curve = self.p_landscape.drug_curve + drug_curve= np.array([drug_curve]) + drug_curve = np.transpose(drug_curve) + counts_seascape, ft = self.p_seascape.simulate() + + if not self.debug: + savedata = np.concatenate((counts_landscape,drug_curve),axis=1) + self.save_counts(savedata, num=None, save_folder=None,prefix = 'landscape_counts') + + savedata = np.concatenate((counts_seascape,drug_curve),axis=1) + self.save_counts(savedata, num=None, save_folder=None,prefix = 'seascape_counts') + + elif self.experiment_type == 'inoculant-survival': + # pbar = tqdm(total = n_curves*n_inoc) # progress bar + for curve_number in range(n_curves): + for inoc_num in range(n_inoc): + + exp_num = curve_number*n_inoc + inoc_num + pop = self.populations[exp_num] # extract population in list of population + c,n_survive_t = pop.simulate() + pop.plot_timecourse() + self.n_survive[curve_number,inoc_num] = n_survive_t + # pbar.update() + self.perc_survive = 100*self.n_survive/self.n_sims + + elif self.experiment_type == 'drug-regimen': + # pbar = tqdm(total=len(self.populations)) + # kk=0 + for p in self.populations: + save_folder = 'p_drop=' + str(p.prob_drop) + save_folder = save_folder.replace('.',',') + # self.exp_folder.append(save_folder) + for i in range(self.n_sims): + # initialize new drug curve + p.drug_curve,u = p.gen_curves() + + counts,n_survive = p.simulate() + drug = p.drug_curve + drug = np.array([drug]) + drug = np.transpose(drug) + + # regimen = self.compute_regimen(p, u) + # regimen = np.array([regimen]) + # regimen = np.transpose(regimen) + # regimen_t = np.zeros(drug.shape) + # regimen_t[0:len(regimen)] = regimen + + u = np.array([u,]) + u = u.transpose() + counts = np.concatenate((counts,drug,u),axis=1) + + if not self.debug: + # self.save_counts(counts,i,save_folder) + data_dict = {'counts':counts, + 'drug_curve':drug, + 'regimen':u} + self.save_dict(data_dict,i,save_folder) + # kk+=1 + # pbar.update() + self.perc_survive = 100*self.n_survive/self.n_sims + + elif self.experiment_type == 'dose-entropy': + # pbar = tqdm(total=len(self.populations)*self.n_sims) + e_survived = [] + e_died = [] + for p in self.populations: + + for i in range(self.n_sims): + c,n_survive = p.simulate() + # e = max(p.entropy()) # compute max entropy + e_t = p.entropy() + e = max(e_t) + # e=1 + # p.plot_timecourse() + + if n_survive == 1: + survive = 'survived' # survived + e_survived.append(e_t) + else: + survive = 'extinct' # died + e_died.append(e_t) + + d = {'dose':[p.max_dose], + 'survive condition':[survive], + 'max entropy':[e]} + + entropy_results_t = pd.DataFrame(d) + self.entropy_results = self.entropy_results.append(entropy_results_t) + # pbar.update() + + elif self.experiment_type == 'rate-survival': + # pbar = tqdm(total=len(self.populations)) + + for p in self.populations: + + for n in range(self.n_sims): + counts,n_survive = p.simulate() + + drug = p.drug_curve + drug = np.array([drug]) + drug = np.transpose(drug) + # counts = np.concatenate((counts,drug),axis=1) + + if self.debug is False: + if (self.curve_types[0] == 'pharm' or + self.curve_types[0] == 'pulsed'): + save_folder = 'k_abs=' + str(p.k_abs) + save_folder.replace('.',',') + else: + save_folder = 'slope=' + str(p.slope) + save_folder.replace('.',',') + self.save_counts(counts,n,save_folder) + data_dict = {'counts':counts, + 'drug_curve':drug} + self.save_dict(data_dict,n,save_folder) + # pickle.dump(self, open(self.experiment_info_path,"wb")) + + + # fig_savename = 'slope = ' + str(p.slope) + # self.figures = self.figures.append(fig) + # pbar.update() + # self.rate_survival_results.index = np.arange(len(self.rate_survival_results)) + + # pbar.close() # close progress bar + + # save counts as a csv in the given subfolder with the label 'num' + def save_counts(self,counts,num,save_folder,prefix='sim_'): + + # check if the desired save folder exists. If not, create it + if save_folder is None: + save_folder = '' + folder_path = self.results_path + os.sep + save_folder + if os.path.exists(folder_path) != True: + os.mkdir(folder_path) + + if num is None: + num = '' + else: + num = str(num).zfill(4) + + savename = self.results_path + os.sep + save_folder + os.sep + prefix + num + '.csv' + np.savetxt(savename, counts, delimiter=",") + # self.savename = savename + return + + def save_dict(self,data_dict,num,save_folder,prefix='sim_'): + # check if the desired save folder exists. If not, create it + if save_folder is None: + save_folder = '' + folder_path = self.results_path + os.sep + save_folder + if os.path.exists(folder_path) != True: + os.mkdir(folder_path) + + if num is None: + num = '' + else: + num = str(num).zfill(4) + + if folder_path not in self.exp_folders: + self.exp_folders.append(folder_path) + + savename = self.results_path + os.sep + save_folder + os.sep + prefix + num + '.p' + pickle.dump(data_dict, open(savename,"wb")) + # np.savetxt(savename, counts, delimiter=",") + return + + def compute_regimen(self,p,u): + gap = int(p.dose_schedule/p.timestep_scale) + n_impulse = int(np.ceil(p.n_timestep/gap)) + regimen = np.zeros(n_impulse) + + for i in range(n_impulse): + if u[(i)*gap] == 1: + regimen[i] = 1 + + return regimen + + def set_ramp_ud(self,p): + + n_timestep = p.n_timestep + drug_curve = np.zeros(n_timestep) + slope = (self.second_dose-self.first_dose)/self.ramp + buffer = self.ramp/2 + + times = [self.transition_times[0]-buffer, + self.transition_times[0]+buffer, + self.transition_times[1]-buffer, + self.transition_times[1]+buffer] + + for t in range(n_timestep): + if t 1: + c = np.sum(counts,axis=1) + else: + c = counts + e = np.argwhere(c 1: + c = counts[:,genotype] + else: + c = counts + + if thresh < 1: + thresh = thresh*pop.max_cells + + e = np.argwhere(c>thresh) + if len(e) == 0: + event_obs = 0 + event_time = len(c) + else: + event_obs = 1 + event_time = e[0] + + timestep_scale = pop.timestep_scale + event_time = event_time*timestep_scale + + return event_obs, event_time + + def log_rank_test(self,durations_A, durations_B, + event_observed_A=None, event_observed_B=None): + + results = lifelines.statistics.logrank_test(durations_A, durations_B, + event_observed_A=event_observed_A, + event_observed_B=event_observed_B) + + return results + +>>>>>>> 5293da446901ef348c8e17b1212b19f4df7cd71c:build/lib/fears/experiment.py diff --git a/build/lib/fears/population.py b/build/lib/fears/population.py index 047423e..d958f83 100644 --- a/build/lib/fears/population.py +++ b/build/lib/fears/population.py @@ -2,7 +2,7 @@ import math import random from importlib_resources import files -from fears.utils import dir_manager, pharm, fitness, plotter +from fears.utils import dir_manager, pharm, fitness, plotter, AutoRate class PopParams: """Population parameters class @@ -47,15 +47,39 @@ def __init__(self,**kwargs): p = files('fears.data').joinpath('ogbunugafor_drugless.csv') self.drugless_data_path = str(p) + plate_paths = ['20210929_plate1.csv','20210929_plate2.csv','20210929_plate3.csv'] + plate_paths = [files('fears.data').joinpath(p) for p in plate_paths] + self.plate_paths = [str(p) for p in plate_paths] + self.seascape_drug_conc = [0,0.003,0.0179,0.1072,0.643,3.858,23.1481,138.8889,833.3333,5000] #ug/mL + + # self.growth_rate_data = [] + # for plate_path in self.plate_paths: + # self.growth_rate_data.append(dir_manager.get_growth_rate_data(plate_path)) + self.constant_pop, self.use_carrying_cap = False, True self.carrying_cap = 10**10 + self.init_counts = None self.n_allele, self.n_genotype = None, None self.doubling_time = 1 self.fitness_data = 'two-point' self.seascape_type = 'natural' self.drug_units = '$\u03BC$M' + self.fig_title = None + self.plot_entropy = False + self.plot_drug_curve = True + self.x_lim = None + self.y_lim = None + self.counts_log_scale = None + self.drug_log_scale = False + + self.n_timestep = 1000 + self.timestep_scale = 1 + self.passage = False + self.passage_time = 24 + self.max_cells = 10**9 self.curve_type = 'pharm' + self.prob_drop = 0 self.k_elim = 0.001 self.k_abs = 0.01 self.pad_right = True @@ -63,25 +87,46 @@ def __init__(self,**kwargs): self.dose_schedule = 24 self.p_forget = 0 + self.stop_condition = None + self.state = {} + self.plot = True + self.n_sims = 10 + self.debug = False + + self.landscape_type = 'natural' + + for paramkey in self.__dict__.keys(): for optkey in kwargs.keys(): if paramkey == optkey: td = {paramkey:kwargs.get(paramkey)} self.__dict__.update(td) - self.ic50 = dir_manager.load_fitness(self.ic50_data_path) - self.drugless_rates = dir_manager.load_fitness(self.drugless_data_path) + # self.ic50 = dir_manager.load_fitness(self.ic50_data_path) + # self.drugless_rates = dir_manager.load_fitness(self.drugless_data_path) + # print(self.ic50) + # dr, ic50 = self.initialize_fitness() + # self.drugless_rates = dr + # self.ic50 = ic50 + # if self.n_genotype is None: + # self.n_genotype = int(len(self.ic50)) + # if self.n_allele is None: + # self.n_allele = int(np.log2(self.n_genotype)) + # if int(self.n_allele) != int(np.log2(self.n_genotype)): + # raise Warning('Genotype/allele number mismatch') - if self.n_genotype is None: - self.n_genotype = int(len(self.ic50)) - if self.n_allele is None: - self.n_allele = int(np.log2(self.n_genotype)) - if int(self.n_allele) != int(np.log2(self.n_genotype)): - raise Warning('Genotype/allele number mismatch') - - self.init_counts = np.zeros(self.n_genotype) - self.init_counts[0] = 10**6 + # self.init_counts = np.zeros(self.n_genotype) + # self.init_counts[0] = 10**6 + + # def initialize_fitness(self): + # if self.fitness_data == 'two-point': + # drugless_rates = dir_manager.load_fitness(self.drugless_data_path) + # ic50 = dir_manager.load_fitness(self.ic50_data_path) + # else: + # drugless_rates = None + # ic50 = None + # return drugless_rates, ic50 class Population(PopParams): @@ -89,13 +134,6 @@ class Population(PopParams): def __init__(self,**kwargs): super().__init__(**kwargs) - # initialize fitness data - self.drugless_rates = None - self.ic50 = None - # load data - - self.initialize_fitness() - # initialize constant population condition if self.constant_pop: self.init_counts = self.init_counts*self.max_cells/sum(self.init_counts) @@ -106,18 +144,56 @@ def __init__(self,**kwargs): self.drug_curve = None self.impulses = None self.initialize_drug_curve() + + # initialize fitness + self.initialize_fitness() + + if self.n_genotype is None: + self.n_genotype = int(len(self.ic50)) + if self.n_allele is None: + self.n_allele = int(np.log2(self.n_genotype)) + if int(self.n_allele) != int(np.log2(self.n_genotype)): + raise Warning('Genotype/allele number mismatch') + + # initialize counts + self.counts = np.zeros([self.n_timestep,self.n_genotype]) + + if self.init_counts is None: + self.init_counts = np.zeros(self.n_genotype) + self.init_counts[0] = 10**6 def initialize_fitness(self): - if self.fitness_data == 'two_point': + if self.fitness_data == 'two-point': self.drugless_rates = dir_manager.load_fitness(self.drugless_data_path) self.ic50 = dir_manager.load_fitness(self.ic50_data_path) + elif self.fitness_data == 'estimate': + + # self.growth_rate_library = fitness.gen_growth_rate_library() + # self.seascape_library = fitness.gen_seascape_library() + + f = str(files('fears.data').joinpath('plates')) + e = AutoRate.Experiment(f,drug_conc=self.seascape_drug_conc,moat=True) + e.execute() + self.growth_rate_lib = e.growth_rate_lib + self.seascape_lib = e.seascape_lib + + self.ic50 = np.zeros(self.n_genotype) + self.drugless_rates = np.zeros(self.n_genotype) + + i = 0 + print(self.seascape_lib.keys()) + + for key in self.seascape_lib.keys(): + self.ic50[i] = self.seascape_lib[key]['ic50'] + self.drugless_rates[i] = self.seascape_lib[key]['g_drugless'] + i+=1 def initialize_drug_curve(self): curve,u = pharm.gen_curves(self) self.drug_curve = curve self.impulses = u -############################################################################### + ############################################################################### # ABM helper methods def gen_neighbors(self,genotype): mut = range(self.n_allele) @@ -164,6 +240,22 @@ def random_mutations(self,N): trans_mat = trans_mat/trans_mat.sum(axis=1) return trans_mat + def check_stop_cond(self,counts,mm): + final_landscape = self.gen_fit_land(self.max_dose) + fittest_genotype = final_landscape.argmax() + + most_frequent_genotype = counts.argmax() + stop_cond = False + + if fittest_genotype == most_frequent_genotype: + stop_cond = True + + if mm >= self.n_timestep: + raise Warning('Stop condition not reached. Increase n_timestep or adjust model parameters.') + stop_cond = True + + return stop_cond + def passage_cells(self,mm,counts): """ If self.passage is true, dilute cells according to self.dilution when @@ -189,7 +281,7 @@ def passage_cells(self,mm,counts): counts = np.divide(counts,self.dilution) return counts -############################################################################## + ############################################################################## # core evolutionary model def abm(self,mm,n_genotype,P,counts): @@ -301,19 +393,95 @@ def simulate(self): # n_survive = 0 for i in range(self.n_sims): - # if self.prob_drop > 0: - # self.drug_curve,u = self.gen_curves() - counts, mm = self.run_abm() avg_counts += counts fixation_time.append(mm) if self.plot is True: - self.plot_timecourse(counts_t = counts) + # print(type(counts)) + self.plot_timecourse(counts_t=counts) avg_counts = avg_counts/self.n_sims self.counts = avg_counts return avg_counts, fixation_time + + ############################################################################## + # wrapper methods for plotting + def plot_timecourse(self,**kwargs): + fig,ax = plotter.plot_timecourse(self,**kwargs) + return fig,ax + + def plot_fitness_curves(self,**kwargs): + fig,ax = plotter.plot_fitness_curves(self,**kwargs) + return fig,ax + + def plot_landscape(self,**kwargs): + fig,ax = plotter.plot_landscape(self,**kwargs) + return fig,ax + + ############################################################################## + # wrapper methods for fitness + + def __gen_fl_for_abm(self,conc,counts): + """ + Return the fitness landscape apropriately scaled according to the + population size and carrying capacity + + Parameters + ---------- + conc (float) : drug concentration + counts (list) : vector of genotype population counts + + Returns + ---------- + fit_land (list) : scaled fitness landscape + """ + fit_land = fitness.gen_fl_for_abm(self,conc,counts) + return fit_land + + def gen_fit_land(self,conc,**kwargs): + """ + Returns the fitness landscape at a given drug concentration + + Parameters + ---------- + conc (float) : drug concentration + + Returns + ---------- + fit_land (list) : fitness landscape + + """ + fit_land = fitness.gen_fit_land(self,conc,**kwargs) + return fit_land + + ############################################################################### + # Wrapper methods for generating drug concentration curves + + def pharm_eqn(self,t,k_elim=None,k_abs=None,max_dose=None): + conc = pharm.pharm_eqn(self,t,k_elim=k_elim,k_abs=k_abs,max_dose=max_dose) + return conc + + def convolve_pharm(self,u): + conv = pharm.convolve_pharm(self,u) + return conv + + def gen_impulses(self): + u = pharm.gen_impulses(self) + return u + + def gen_on_off(self,duty_cycle=None): + u = pharm.gen_on_off(self,duty_cycle=duty_cycle) + return u + + def gen_curves(self): + curve, u = pharm.gen_curves(self) + return curve, u + def gen_passage_drug_protocol(self): + drug_curve = pharm.gen_passage_drug_protocol(self) + return drug_curve -# p = Population() \ No newline at end of file + def set_drug_curve(self): + dc = self.gen_curves() + self.drug_curve = dc[0] \ No newline at end of file diff --git a/build/lib/fears/utils/__init__.py b/build/lib/fears/utils/__init__.py index e69de29..3bff057 100644 --- a/build/lib/fears/utils/__init__.py +++ b/build/lib/fears/utils/__init__.py @@ -0,0 +1 @@ +from .plotter import gen_color_cycler, plot_timecourse, plot_fitness_curves, plot_msw, plot_timecourse_to_axes, plot_landscape, add_landscape_to_fitness_curve, get_pos_in_log_space, plot_population_count, plot_kaplan_meier, x_ticks_to_days, shiftx, shifty, shrinky diff --git a/build/lib/fears/utils/fitness.py b/build/lib/fears/utils/fitness.py index a534b6e..048c7ed 100644 --- a/build/lib/fears/utils/fitness.py +++ b/build/lib/fears/utils/fitness.py @@ -1,389 +1,336 @@ import numpy as np +from scipy.optimize import curve_fit # from fears.population import Population -class Fitness: - def __init__(self): - return +def gen_fitness_curves(pop,conc=None): - def gen_fitness_curves(self,pop=None,conc=None): - - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') - - if conc is None: - conc = np.logspace(-3,5,num=1000) - - n_genotype = pop.n_genotype - - fc = {} - for g in range(n_genotype): - f = np.zeros(len(conc)) - i = 0 - for c in conc: - f[i] = self.gen_fitness(g,c,pop=pop) - pop.death_rate - i+=1 - fc[g] = f - - return fc - - # compute fitness given a drug concentration - def gen_fitness(self,genotype,conc,pop=None, - drugless_rate=None,ic50=None): - - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') - - if drugless_rate is None: - drugless_rate = pop.drugless_rates - if ic50 is None: - ic50 = pop.ic50 - - # logistic equation from Ogbunugafor 2016 - conc = conc/10**6 # concentration in uM, convert to M - c = -.6824968 # empirical curve fit - log_eqn = lambda d,i: d/(1+np.exp((i-np.log10(conc))/c)) - if conc <= 0: - fitness = drugless_rate[genotype] - else: - fitness = log_eqn(drugless_rate[genotype],ic50[genotype]) - - return fitness - - def logistic_equation(self,conc,drugless_rate,ic50): - """ - Logistic equation from ogbunugafor et al, PLOS CB, 2016 - - Parameters - ---------- - dugless_rate : float - Drugless growth rate of genotype. - ic50 : float - ic50 of genotype. - conc : float - Drug concentration (in Molarity (M)). - c : float, optional - Logistic curve steepness parameter. The default is -0.6824968. - - Returns - ------- - f : float - Replication rate. - - """ - c=-0.6824968 - conc = conc/10**6 - f = drugless_rate/(1+np.exp((ic50-np.log10(conc))/c)) - - return f - def gen_static_landscape(self,conc,pop=None): - - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') - - # get final landscape and seascape - landscape = np.zeros(pop.n_genotype) - for kk in range(pop.n_genotype): - landscape[kk] = self.gen_fitness(pop, - kk, - pop.static_topo_dose, - pop.drugless_rates, - pop.ic50) - - if min(landscape) == 0: - zero_indx_land = np.argwhere(landscape==0) - landscape_t = np.delete(landscape,zero_indx_land) - min_landscape = min(landscape_t) - else: - min_landscape = min(landscape) - zero_indx_land = [] - - seascape = np.zeros(pop.n_genotype) - for gen in range(pop.n_genotype): - seascape[gen] = self.gen_fitness(gen,conc,pop.drugless_rates,pop.ic50,pop=pop) - - if min(seascape) == 0: - zero_indx_sea = np.argwhere(seascape==0) - seascape_t = np.delete(seascape,zero_indx_sea) - min_seascape = min(seascape_t) - else: - min_seascape = min(seascape) - - landscape = landscape - min_landscape - landscape = landscape/max(landscape) - - rng = max(seascape) - min_seascape - - landscape = landscape*rng + min_seascape - - landscape[zero_indx_land] = 0 - return landscape - - def gen_digital_seascape(self,conc,gen,pop=None,min_fitness=0): - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') - - if pop.mic_estimate is not None: - mic = self.est_mic(gen,Kmic=pop.mic_estimate,pop=pop) - else: - mic = self.est_mic(gen,growth_rate=pop.death_rate,pop=pop) - - if conc >= mic: - fitness = min_fitness - else: - fitness = pop.drugless_rates[gen] - return fitness + if conc is None: + conc = np.logspace(-3,5,num=1000) + + n_genotype = pop.n_genotype + + fc = {} + for g in range(n_genotype): + f = np.zeros(len(conc)) + i = 0 + for c in conc: + f[i] = gen_fitness(g,c,pop=pop) - pop.death_rate + i+=1 + fc[g] = f + + return fc + +# compute fitness given a drug concentration +def gen_fitness(pop,genotype,conc, + drugless_rate=None,ic50=None): + + if drugless_rate is None: + drugless_rate = pop.drugless_rates + if ic50 is None: + ic50 = pop.ic50 + + # logistic equation from Ogbunugafor 2016 + conc = conc/10**6 # concentration in uM, convert to M + c = -.6824968 # empirical curve fit + log_eqn = lambda d,i: d/(1+np.exp((i-np.log10(conc))/c)) + if conc <= 0: + fitness = drugless_rate[genotype] + else: + fitness = log_eqn(drugless_rate[genotype],ic50[genotype]) + + return fitness + +def logistic_equation(conc,drugless_rate,ic50): + """ + Logistic equation from ogbunugafor et al, PLOS CB, 2016 + + Parameters + ---------- + dugless_rate : float + Drugless growth rate of genotype. + ic50 : float + ic50 of genotype. + conc : float + Drug concentration (in Molarity (M)). + c : float, optional + Logistic curve steepness parameter. The default is -0.6824968. + + Returns + ------- + f : float + Replication rate. + + """ + c=-0.6824968 + conc = conc/10**6 + f = drugless_rate/(1+np.exp((ic50-np.log10(conc))/c)) + + return f - def gen_fit_land(self,conc,pop=None,mode=None): - - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') +def gen_static_landscape(pop,conc): + + # get final landscape and seascape + landscape = np.zeros(pop.n_genotype) + for kk in range(pop.n_genotype): + landscape[kk] = gen_fitness(pop, + kk, + pop.static_topo_dose, + pop.drugless_rates, + pop.ic50) + + if min(landscape) == 0: + zero_indx_land = np.argwhere(landscape==0) + landscape_t = np.delete(landscape,zero_indx_land) + min_landscape = min(landscape_t) + else: + min_landscape = min(landscape) + zero_indx_land = [] + + seascape = np.zeros(pop.n_genotype) + for gen in range(pop.n_genotype): + seascape[gen] = gen_fitness(gen,conc,pop.drugless_rates,pop.ic50,pop=pop) + + if min(seascape) == 0: + zero_indx_sea = np.argwhere(seascape==0) + seascape_t = np.delete(seascape,zero_indx_sea) + min_seascape = min(seascape_t) + else: + min_seascape = min(seascape) + + landscape = landscape - min_landscape + landscape = landscape/max(landscape) + + rng = max(seascape) - min_seascape + + landscape = landscape*rng + min_seascape + + landscape[zero_indx_land] = 0 + return landscape - fit_land = np.zeros(pop.n_genotype) - - if pop.fitness_data == 'manual' or mode=='manual': - fit_land = pop.landscape_data/pop.doubling_time +def gen_digital_seascape(pop,conc,gen,min_fitness=0): - else: - - if pop.landscape_type == 'static': - fit_land = self.gen_static_landscape(pop,conc) - - if pop.landscape_type == 'digital': - for kk in range(pop.n_genotype): - fit_land[kk] = self.gen_digital_seascape(pop, conc, kk) - - elif pop.landscape_type == 'natural': - for kk in range(pop.n_genotype): - fit_land[kk] = self.gen_fitness(pop, - kk, - conc, - pop.drugless_rates, - pop.ic50)/pop.doubling_time - - return fit_land + if pop.mic_estimate is not None: + mic = est_mic(gen,Kmic=pop.mic_estimate,pop=pop) + else: + mic = est_mic(gen,growth_rate=pop.death_rate,pop=pop) + + if conc >= mic: + fitness = min_fitness + else: + fitness = pop.drugless_rates[gen] + return fitness - # Generate fitness landscape for use in the abm method - def gen_fl_for_abm(self,conc,counts,pop=None): +def gen_fit_land(pop,conc,mode=None): - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') + fit_land = np.zeros(pop.n_genotype) + + if pop.fitness_data == 'manual' or mode=='manual': + fit_land = pop.landscape_data/pop.doubling_time - fit_land = self.gen_fit_land(pop,conc) - - # # takes the landscape at the max dose and scales the replication rate - # # according to drug concentration - # if pop.static_landscape: - # # max_fitness = max(fit_land) - # # fit_land = pop.gen_fit_land(pop.max_dose) - # # fit_land = fit_land*max_fitness/max(fit_land) - # fit_land = gen_fit_land(pop,conc) - - # if pop.static_topology: - # fit_land = gen_fit_land(pop,conc) - - # Scale division rates based on carrying capacity - if pop.carrying_cap: - division_scale = 1-np.sum(counts)/pop.max_cells - if counts.sum()>pop.max_cells: - division_scale = 0 - else: - division_scale = 1 + else: - fit_land = fit_land*division_scale - - return fit_land + if pop.landscape_type == 'static': + fit_land = gen_static_landscape(pop,conc) + + if pop.landscape_type == 'digital': + for kk in range(pop.n_genotype): + fit_land[kk] = gen_digital_seascape(pop, conc, kk) + + elif pop.landscape_type == 'natural': + for kk in range(pop.n_genotype): + fit_land[kk] = gen_fitness(pop, + kk, + conc, + pop.drugless_rates, + pop.ic50)/pop.doubling_time + + return fit_land - def gen_random_seascape(self,n_allele, - drugless_limits=[1,1.5], - ic50_limits=[-6.5,-1.5]): - - n_genotype = 2**n_allele - - drugless_rates = np.random.uniform(min(drugless_limits), - max(drugless_limits), - n_genotype) - - ic50 = np.random.uniform(min(ic50_limits), - max(ic50_limits), - n_genotype) - - return drugless_rates,ic50 +# Generate fitness landscape for use in the abm method +def gen_fl_for_abm(pop,conc,counts): + + fit_land = gen_fit_land(pop,conc) + + # # takes the landscape at the max dose and scales the replication rate + # # according to drug concentration + # if pop.static_landscape: + # # max_fitness = max(fit_land) + # # fit_land = pop.gen_fit_land(pop.max_dose) + # # fit_land = fit_land*max_fitness/max(fit_land) + # fit_land = gen_fit_land(pop,conc) + + # if pop.static_topology: + # fit_land = gen_fit_land(pop,conc) + + # Scale division rates based on carrying capacity + if pop.use_carrying_cap: + division_scale = 1-np.sum(counts)/pop.carrying_cap + if counts.sum()>pop.max_cells: + division_scale = 0 + else: + division_scale = 1 + + fit_land = fit_land*division_scale + + return fit_land - def randomize_seascape(self,pop=None, +def gen_random_seascape(n_allele, drugless_limits=[1,1.5], ic50_limits=[-6.5,-1.5]): - - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') - - n_genotype = pop.n_genotype - - pop.drugless_rates = np.random.uniform(min(drugless_limits), - max(drugless_limits), - n_genotype) - - pop.ic50 = np.random.uniform(min(ic50_limits), - max(ic50_limits), + + n_genotype = 2**n_allele + + drugless_rates = np.random.uniform(min(drugless_limits), + max(drugless_limits), n_genotype) - - def fit_logistic_curve(self,xdata,ydata): - from scipy.optimize import curve_fit - - popt,var = curve_fit(self.logistic_equation,xdata,ydata) - - return popt + + ic50 = np.random.uniform(min(ic50_limits), + max(ic50_limits), + n_genotype) + + return drugless_rates,ic50 - def gen_null_seascape(self,conc,pop=None): +def randomize_seascape(pop, + drugless_limits=[1,1.5], + ic50_limits=[-6.5,-1.5]): - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') + n_genotype = pop.n_genotype + + pop.drugless_rates = np.random.uniform(min(drugless_limits), + max(drugless_limits), + n_genotype) - landscape = self.gen_fit_land(conc,pop=pop) - start_rates = self.gen_fit_land(10**-3,pop=pop) - final_rates = self.gen_fit_land(10**5,pop=pop) - # mid_rates = gen_fit_land(pop,10**1) - - start_points = self.scale_and_ignore_zeros(landscape,start_rates) - end_points = self.scale_and_ignore_zeros(landscape,final_rates) - # mid_points = scale_and_ignore_zeros(landscape,mid_rates) - mid_points = landscape - - xdata = [10**-3,conc,10**5] - - ic50_new = [] - drugless_rates_new = [] - - for genotype in range(len(landscape)): - ydata = [start_points[genotype], - mid_points[genotype], - end_points[genotype]] - params = self.fit_logistic_curve(xdata,ydata) - ic50_new.append(params[1]) - drugless_rates_new.append(params[0]) - # find the null landscape drugless rates - - drugless_rates_new = self.scale_and_ignore_zeros(drugless_rates_new, - pop.drugless_rates) - - return drugless_rates_new,ic50_new - - def scale_and_ignore_zeros(self,data,target): - """ - Scale data to range of target while ignoring the zero values in data and - target. - - Parameters - ---------- - data : numpy array - Data to be scaled to the range of target. - target : numpy array - Target data range. - - Returns - ------- - scaled_data : numpy array - Scaled data to range of target. Zero values in data are set to zero - in scaled_data and zero values in target are not used to calculate - range. - - """ - # make sure inputs are numpy arrays - - if not isinstance(data,np.ndarray): - data=np.array(data) - if not isinstance(target,np.ndarray): - target=np.array(target) - - if min(data) == 0: - zero_indx_data = np.argwhere(data==0) - data_t = np.delete(data,zero_indx_data) - min_data = min(data_t) - else: - min_data = min(data) - zero_indx_data = [] - - if min(target) == 0: - zero_indx_target = np.argwhere(target==0) - target_t = np.delete(target,zero_indx_target) - min_target = min(target_t) - else: - min_target = min(target) - - data = data - min_data - data = data/max(data) + pop.ic50 = np.random.uniform(min(ic50_limits), + max(ic50_limits), + n_genotype) + +def fit_logistic_curve(xdata,ydata): + + popt,var = curve_fit(logistic_equation,xdata,ydata) + + return popt - rng = max(target) - min_target - - scaled_data = data*rng + min_target +def gen_null_seascape(pop,conc): - scaled_data[zero_indx_data] = 0 - - return scaled_data + landscape = gen_fit_land(conc,pop=pop) + start_rates = gen_fit_land(10**-3,pop=pop) + final_rates = gen_fit_land(10**5,pop=pop) + # mid_rates = gen_fit_land(pop,10**1) + + start_points = scale_and_ignore_zeros(landscape,start_rates) + end_points = scale_and_ignore_zeros(landscape,final_rates) + # mid_points = scale_and_ignore_zeros(landscape,mid_rates) + mid_points = landscape + + xdata = [10**-3,conc,10**5] + + ic50_new = [] + drugless_rates_new = [] + + for genotype in range(len(landscape)): + ydata = [start_points[genotype], + mid_points[genotype], + end_points[genotype]] + params = fit_logistic_curve(xdata,ydata) + ic50_new.append(params[1]) + drugless_rates_new.append(params[0]) + # find the null landscape drugless rates + + drugless_rates_new = scale_and_ignore_zeros(drugless_rates_new, + pop.drugless_rates) + + return drugless_rates_new,ic50_new + +def scale_and_ignore_zeros(data,target): + """ + Scale data to range of target while ignoring the zero values in data and + target. + + Parameters + ---------- + data : numpy array + Data to be scaled to the range of target. + target : numpy array + Target data range. + + Returns + ------- + scaled_data : numpy array + Scaled data to range of target. Zero values in data are set to zero + in scaled_data and zero values in target are not used to calculate + range. + + """ + # make sure inputs are numpy arrays + + if not isinstance(data,np.ndarray): + data=np.array(data) + if not isinstance(target,np.ndarray): + target=np.array(target) + + if min(data) == 0: + zero_indx_data = np.argwhere(data==0) + data_t = np.delete(data,zero_indx_data) + min_data = min(data_t) + else: + min_data = min(data) + zero_indx_data = [] + + if min(target) == 0: + zero_indx_target = np.argwhere(target==0) + target_t = np.delete(target,zero_indx_target) + min_target = min(target_t) + else: + min_target = min(target) + + data = data - min_data + data = data/max(data) + + rng = max(target) - min_target + + scaled_data = data*rng + min_target - def est_mic(self,gen,pop=None,Kmic=None,growth_rate=None): - """ - est_mic: estimates the mic based on a given Kmic (ratio of growth rate to - max growth rate at MIC) or based on a given growth rate. + scaled_data[zero_indx_data] = 0 + + return scaled_data + +def est_mic(pop,gen,Kmic=None,growth_rate=None): + """ + est_mic: estimates the mic based on a given Kmic (ratio of growth rate to + max growth rate at MIC) or based on a given growth rate. + + Parameters + ---------- + pop : population class object + + gen : int + Genotype under consideration. + Kmic : float, optional + Ratio of growth rate to max growth rate at MIC. The default is None. + growth_rate : float, optional + Growth rate at MIC. The default is None. + + Raises + ------ + Exception + Function requires Kmic OR growth_rate to calculate MIC. + + Returns + ------- + mic : float + MIC at a given growth rate or Kmic. + + """ + + if Kmic is None: + if growth_rate is None: + raise Exception('Need a growth rate or Kmic threshold to estimate mic.') + else: + Kmic = growth_rate/pop.drugless_rates[gen] + c=-0.6824968 + mic = 10**(pop.ic50[gen]+6 - c*np.log((1/Kmic)-1)) + return mic - Parameters - ---------- - pop : population class object - - gen : int - Genotype under consideration. - Kmic : float, optional - Ratio of growth rate to max growth rate at MIC. The default is None. - growth_rate : float, optional - Growth rate at MIC. The default is None. - - Raises - ------ - Exception - Function requires Kmic OR growth_rate to calculate MIC. - - Returns - ------- - mic : float - MIC at a given growth rate or Kmic. - - """ - if pop is None: - if type(self) is Population: - pop = self - else: - raise TypeError('Population object required') - - if Kmic is None: - if growth_rate is None: - raise Exception('Need a growth rate or Kmic threshold to estimate mic.') - else: - Kmic = growth_rate/pop.drugless_rates[gen] - c=-0.6824968 - mic = 10**(pop.ic50[gen]+6 - c*np.log((1/Kmic)-1)) - return mic \ No newline at end of file + \ No newline at end of file diff --git a/fears.egg-info/SOURCES.txt b/fears.egg-info/SOURCES.txt index 45bae93..750fc76 100644 --- a/fears.egg-info/SOURCES.txt +++ b/fears.egg-info/SOURCES.txt @@ -10,6 +10,9 @@ fears.egg-info/SOURCES.txt fears.egg-info/dependency_links.txt fears.egg-info/requires.txt fears.egg-info/top_level.txt +fears/classes/__init__.py +fears/classes/experiment_class.py +fears/classes/population_class.py fears/data/cycloguanil_ic50.csv fears/data/ogbunugafor_drugless.csv fears/data/pyrimethamine_ic50.csv diff --git a/fears/utils/__init__.py b/fears/utils/__init__.py index e69de29..3bff057 100644 --- a/fears/utils/__init__.py +++ b/fears/utils/__init__.py @@ -0,0 +1 @@ +from .plotter import gen_color_cycler, plot_timecourse, plot_fitness_curves, plot_msw, plot_timecourse_to_axes, plot_landscape, add_landscape_to_fitness_curve, get_pos_in_log_space, plot_population_count, plot_kaplan_meier, x_ticks_to_days, shiftx, shifty, shrinky diff --git a/fears/utils/dir_manager.py b/fears/utils/dir_manager.py old mode 100755 new mode 100644 index 14edde7..907d68b --- a/fears/utils/dir_manager.py +++ b/fears/utils/dir_manager.py @@ -54,4 +54,4 @@ def load_growth_rate_data(data_path): def load_experiment(exp_path): e = pickle.load(exp_path) - return e \ No newline at end of file + return e diff --git a/fears/utils/pharm.py b/fears/utils/pharm.py index 02b9aa9..217ceaf 100644 --- a/fears/utils/pharm.py +++ b/fears/utils/pharm.py @@ -130,7 +130,11 @@ def gen_curves(pop): curve = gen_passage_drug_protocol(pop) else: for i in range(pop.n_timestep): +<<<<<<< HEAD + curve[i] = pop.pharm_eqn(i) +======= curve[i] = pharm_eqn(pop,i) +>>>>>>> 5293da446901ef348c8e17b1212b19f4df7cd71c # Pulsed convolves an impulse train with the 1-compartment model elif pop.curve_type == 'pulsed': diff --git a/fears/utils/plotter.py b/fears/utils/plotter.py index eeff0bb..0fbf5eb 100644 --- a/fears/utils/plotter.py +++ b/fears/utils/plotter.py @@ -29,7 +29,11 @@ def gen_color_cycler(style=None,palette='bright',n_colors=16): cc = cycler(color=colors) return cc +<<<<<<< HEAD +def plot_timecourse(pop,counts_t=None,title_t=None): +======= def plot_timecourse(pop,counts_t=None,title_t=None,**kwargs): +>>>>>>> 5293da446901ef348c8e17b1212b19f4df7cd71c if (pop.counts == 0).all() and counts_t is None: print('No data to plot!') @@ -99,6 +103,10 @@ def plot_timecourse(pop,counts_t=None,title_t=None,**kwargs): ax1.plot(counts[:,allele],linewidth=3.0,label=None) ax1.legend(loc=(1.25,-.12),frameon=False,fontsize=15) +<<<<<<< HEAD + + ax1.set_xlim(0,pop.x_lim) +======= if pop.x_lim is None: xl = ax1.get_xlim() @@ -106,6 +114,7 @@ def plot_timecourse(pop,counts_t=None,title_t=None,**kwargs): xl = pop.x_lim ax1.set_xlim(xl) +>>>>>>> 5293da446901ef348c8e17b1212b19f4df7cd71c ax1.set_facecolor(color='w') ax1.grid(False) @@ -141,7 +150,11 @@ def plot_timecourse(pop,counts_t=None,title_t=None,**kwargs): ax1.set_xlabel('Days',fontsize=20) plt.show() +<<<<<<< HEAD + return fig +======= return fig,ax1 +>>>>>>> 5293da446901ef348c8e17b1212b19f4df7cd71c def plot_fitness_curves(pop, fig_title='', @@ -732,7 +745,11 @@ def plot_landscape(p,conc=10**0, ax.set_ylim(yl) ax.set_axis_off() +<<<<<<< HEAD + return ax +======= return fig,ax +>>>>>>> 5293da446901ef348c8e17b1212b19f4df7cd71c def add_landscape_to_fitness_curve(c,ax,pop, textcolor='gray', diff --git a/setup.py b/setup.py index cece451..be16df6 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,5 @@ + from setuptools import setup, find_packages setup( @@ -16,7 +17,9 @@ "scipy", "matplotlib", "numpy", - "importlib_resources" + "importlib_resources", + "lifelines", + "pickle" ], include_package_data=True, package_data={'': ['data/*.csv']}