Merge pull request #18 from willfinnigan/dev

v1.4.1

docs/conf.py

@@ -24,9 +24,9 @@
 author = 'William Finnigan'
 
 # The short X.Y version
-version = '1.4.0'
+version = '1.4.1'
 # The full version, including alpha/beta/rc tags
-release = '1.4.0'
+release = '1.4.1'
 
 
 # -- General configuration ---------------------------------------------------

kinetics/__init__.py

@@ -1,4 +1,4 @@
-from kinetics.model import Model
+from kinetics.model_module import Model
 
 from kinetics.reaction_classes.general_rate_Law import *
 from kinetics.reaction_classes.irreversible_michaelis_menton import *
@@ -10,11 +10,10 @@
 from kinetics.optimisation.metrics import Metrics, uM_to_mgml
 from kinetics.optimisation.genetic_algorithm import GA_Base_Class
 
-from kinetics.ua_and_sa.sampling import sample_distributions, sample_uniforms, salib_problem, make_saltelli_samples
+from kinetics.ua_and_sa.sampling import sample_distributions, sample_uniforms, salib_problem, make_saltelli_samples, distributions_to_lower_upper_bounds
 from kinetics.ua_and_sa.run_all_models import run_all_models, dataframes_all_runs, dataframes_quartiles
 from kinetics.ua_and_sa.plotting import plot_substrate, plot_ci_intervals, plot_data, remove_st_less_than, plot_sa_total_sensitivity
 from kinetics.ua_and_sa.sensitivity_analysis import get_concentrations_at_timepoint, get_time_to_concentration, analyse_sobal_sensitivity
 
-__version__ = '1.4.0'
 
-#version
+__version__ = '1.4.1'

kinetics/model.py → kinetics/model_module.py

@@ -37,7 +37,7 @@ class Model(list):
     """
 
 
-    def __init__(self, logging=True):
+    def __init__(self, logging=False):
         # Model inherits from list - reaction classes are held in this self list.
         super(Model, self).__init__()
 
@@ -117,7 +117,10 @@ def set_parameters_from_reactions(self):
             # if parameter not set in model, and hasn't been loaded from reaction, take mean of model_distribution
             for name in self.parameter_distributions:
                 if name not in self.parameters:
-                    self.parameters[name] = self.parameter_distributions[name].mean()
+                    if type(self.parameter_distributions[name]) == list or type(self.parameter_distributions[name]) == tuple:
+                        self.parameters[name] = (self.parameter_distributions[name][0] + self.parameter_distributions[name][1]) / 2
+                    else:
+                        self.parameters[name] = self.parameter_distributions[name].mean()
                     if self.logging == True:
                         print(str(name) + ' - ' + str(self.parameters[name]))
 
@@ -159,10 +162,16 @@ def set_species_defaults_to_mean(self):
 
         Called by self.setup_model()
         """
-
+        if self.logging==True:
+            print('-- Setting default species, using means of distributions where undefined: --')
         for name in self.species_distributions:
             if name not in self.species:
-                self.species[name] = self.species_distributions[name].mean()
+                if type(self.species_distributions[name]) == list or type(self.species_distributions[name]) == tuple:
+                    self.species[name] = (self.species_distributions[name][0] + self.species_distributions[name][1])/2
+                else:
+                    self.species[name] = self.species_distributions[name].mean()
+                if self.logging==True:
+                    print(str(name) + ' - ' + str(self.species[name]))
 
     def setup_model(self):
         """
@@ -262,20 +271,23 @@ def results_dataframe(self):
 
         return df
 
-    def plot_substrate(self, substrate, plot=False):
+    def plot_substrate(self, substrate, plot=False, units=['','']):
         """
         Plot a graph of substrate concentration vs time.
 
         Args:
             substrate (str): Name of substrate to plot
             plot (bool): Default False.  If True calls plt.show()
         """
+
         ys_at_t = []
         i = self.run_model_species_names.index(substrate)
         for t in range(len(self.time)):
             ys_at_t.append(self.y[t][i])
 
         plt.plot(self.time, ys_at_t, label=substrate)
+        plt.ylabel(units[0])
+        plt.xlabel(units[1])
         plt.legend()
 
         if plot == True:

kinetics/reaction_classes/reaction_base_class.py

@@ -65,7 +65,10 @@ def __init__(self):
     def set_parameter_defaults_to_mean(self):
         for name in self.parameter_distributions:
             if name not in self.parameters:
-                self.parameters[name] = self.parameter_distributions[name].mean()
+                if type(self.parameter_distributions[name]) == list or type(self.parameter_distributions[name]) == tuple:
+                    self.parameters[name] = (self.parameter_distributions[name][0] + self.parameter_distributions[name][1]) / 2
+                else:
+                    self.parameters[name] = self.parameter_distributions[name].mean()
 
     def get_indexes(self, substrate_names):
         self.substrate_indexes = []
@@ -125,7 +128,8 @@ def reaction(self, y, substrate_names, parameter_dict):
 
         substrates = self.get_substrates(y)
 
-        substrates, parameters = self.calculate_modifiers(substrates, copy.copy(self.run_model_parameters))
+        if len(self.modifiers) != 0:
+            substrates, parameters = self.calculate_modifiers(substrates, copy.copy(self.run_model_parameters))
 
         rate = self.calculate_rate(substrates, parameters)
 

kinetics/ua_and_sa/plotting.py

@@ -2,7 +2,7 @@
 
 """ -- Plotting uncertainty analysis -- """
 def plot_substrate(substrate, dataframes,
-                   colour='blue', xlabel="Time (mins)", ylabel="μM",
+                   colour='blue', units=['',''],
                    alpha=0.1, linewidth=0.1, y_min=True, plot=False):
     """
     Plot every model run for a single substrate.
@@ -20,6 +20,9 @@ def plot_substrate(substrate, dataframes,
 
     """
 
+    ylabel = units[0]
+    xlabel = units[1]
+
     df = dataframes[substrate]
     for i in range(1, len(df.columns)):
         plt.plot(df['Time'], df[str(i)],
@@ -36,7 +39,7 @@ def plot_substrate(substrate, dataframes,
 
 def plot_ci_intervals(substrates_to_add, dataframes, plot=False,
                       colours=['darkred', 'darkblue', 'darkgreen', 'darkorange'],
-                      alpha=0.1, units=['$\mu M$', 'Time (mins)']):
+                      alpha=0.1, units=['','']):
     """
     Plot every model run for a single substrate.
 

kinetics/ua_and_sa/run_all_models.py

@@ -7,7 +7,7 @@ def run_all_models(model, samples, logging=True):
     Run all the models for a set of samples.
 
     Args:
-        model (kinetics.model): A model object
+        model (kinetics.model_module): A model object
         samples (list): A list of samples in the form [(param_dict1, species_dict1), (param_dict2.... ect}
         logging (bool): Show logging and progress bar.  Default = True
 

kinetics/ua_and_sa/sampling.py

@@ -61,21 +61,28 @@ def unpack_sampled_params_and_species(sample_list, parameter_names_list, species
     return parsed_samples
 
 def check_not_neg(sample, name, negative_allowed):
-    def check_sample(sample, name):
-        if (sample < 0) and (name not in negative_allowed):
+    def check_sample(sample_to_check, name_to_check):
+        if (sample_to_check <= 0) and (name_to_check not in negative_allowed):
             return False
 
     if type(sample) == np.ndarray:
         for s in sample:
             print(s)
-            check_sample(s, name)
+            if check_sample(s, name) == False:
+                return False
 
-    elif type(sample) == float:
-        check_sample(sample, name)
+    elif type(sample) == np.float64:
+        if check_sample(sample, name) == False:
+            return False
 
     elif sample == None:
         return False
 
+    else:
+        print('Error no type matches')
+        print(type(sample))
+        return None
+
     return True
 
 def multi_params(distributions):
@@ -102,7 +109,7 @@ def sample_distributions(model, num_samples=1000, negative_allowed=[]):
     Makes a set of samples from the species and parameter distributions in the model.
 
     Args:
-        model (kinetics.model): A model object
+        model (kinetics.model_module): A model object
         num_samples (int): Number of samples to make (default 1000)
         negative_allowed (list): A list of any distributions that can be negative.
 
@@ -118,6 +125,7 @@ def sample_distributions(model, num_samples=1000, negative_allowed=[]):
         species_dict = {}
 
         for name, distribution in model.parameter_distributions.items():
+
             if type(distribution) != list and type(distribution) != tuple:
                 sample = None
                 while not check_not_neg(sample, name, negative_allowed):
@@ -132,7 +140,6 @@ def sample_distributions(model, num_samples=1000, negative_allowed=[]):
                     for multi_name, index in mv_params[name]:
                         parameter_dict[multi_name] = sample[index]
 
-
         for name, distribution in model.species_distributions.items():
             sample=None
             while not check_not_neg(sample, name, negative_allowed):
@@ -143,29 +150,32 @@ def sample_distributions(model, num_samples=1000, negative_allowed=[]):
 
     return samples # samples will = [ (parameter_dict1, species_dict1), (parameter_dict2, species_dict2) ..]
 
-def sample_uniforms(model, num_samples=1000):
+def sample_uniforms(model, num_samples=1000, log=[]):
 
     bounds = []
     names = []
     for name, tuple in model.parameter_distributions.items():
+        names.append(name)
         bounds.append(tuple)
-        names.append(names)
 
-    for name, tuple in  model.species_distributions.items():
+    for name, tuple in model.species_distributions.items():
+        names.append(name)
         bounds.append(tuple)
-        names.append(names)
 
     problem = {'num_vars': len(names),
                'names': names,
                'bounds': bounds}
 
+    problem = salib_problem_to_log_space(problem, log)
+
     lhc_samples = latin.sample(problem, num_samples)
+    lhc_samples = samples_to_normal_space(lhc_samples, problem, log)
 
     samples = parse_samples(lhc_samples, list(model.parameter_distributions.keys()), list(model.species_distributions.keys()))
 
     return samples
 
-def distributions_to_uniforms(model, negative_allowed=[], ppf=(0.05,0.95), save_to_model=False):
+def distributions_to_lower_upper_bounds(model, negative_allowed=[], ppf=(0.05,0.95), save_to_model=False):
     """
     Converts distributions to uniform distributions by taking specified ppf
 
@@ -205,7 +215,7 @@ def distributions_to_uniforms(model, negative_allowed=[], ppf=(0.05,0.95), save_
 
     return bounds
 
-def salib_problem(model, bounds=[]):
+def salib_problem(model, bounds=[], log=[]):
     """
     Make a salib problem by specifying bounds using ppf of scipy distributions
 
@@ -230,9 +240,11 @@ def salib_problem(model, bounds=[]):
                'names': names,
                'bounds': bounds}
 
+    problem = salib_problem_to_log_space(problem, log)
+
     return problem
 
-def make_saltelli_samples(model, salib_problem, num_samples, second_order=False):
+def make_saltelli_samples(model, salib_problem, num_samples, second_order=False, log=[]):
     """
     Use SALib to make saltelli samples
 
@@ -248,7 +260,25 @@ def make_saltelli_samples(model, salib_problem, num_samples, second_order=False)
     """
 
     saltelli_samples = saltelli.sample(salib_problem, num_samples, calc_second_order=second_order)
+    saltelli_samples = samples_to_normal_space(saltelli_samples, salib_problem, log)
+
     samples = parse_samples(saltelli_samples, list(model.parameter_distributions.keys()), list(model.species_distributions.keys()))
 
     return samples
 
+def salib_problem_to_log_space(salib_problem, to_log):
+    for i, name in enumerate(salib_problem['names']):
+        if name in to_log:
+            lower = np.log(salib_problem['bounds'][i][0])
+            upper = np.log(salib_problem['bounds'][i][1])
+            salib_problem['bounds'][i] = [lower, upper]
+
+    return salib_problem
+
+def samples_to_normal_space(salib_samples, salib_problem, to_log):
+    for i, name in enumerate(salib_problem['names']):
+        if name in to_log:
+            for j, sample_set in enumerate(salib_samples):
+                salib_samples[j][i] = float(np.exp(sample_set[i]))
+
+    return salib_samples

setup.py

@@ -7,7 +7,7 @@
               'kinetics.other_analysis',
               'kinetics.reaction_classes',
               'kinetics.ua_and_sa'],
-  version = '1.4.0',
+  version = '1.4.1',
   license='MIT',
   description = 'Python code to run kinetic models of enzyme reactions',
   author = 'William Finnigan',