[master] added few options to fasten things up during replay. v 0.0.5

Author: Giacomo

py_replay_bg/environment/__init__.py

@@ -19,6 +19,10 @@ class Environment:
         A string to be attached as suffix to the resulting output files' name.
     save_workspace: bool
             A boolean that specifies whether or not to save the results of ReplayBG in the results/workspaces directory.
+    analyze_results : bool
+            A flag that specifies whether to analyze the resulting trace or not. Setting this flag to False will fasten
+            ReplayBG and it is recommended if ReplayBG will be a component of a bigger framework (e.g., to be used in an
+            iterative process).
     scenario: str, {'single-meal', 'multi-meal'}
         A string that specifies whether the given scenario refers to a single-meal scenario or a multi-meal scenario.
     bolus_source : str, {'data', or 'dss'}
@@ -47,7 +51,8 @@ class Environment:
     None
     """
 
-    def __init__(self, modality, save_name, save_folder, save_suffix='', save_workspace=False, scenario='single_meal',
+    def __init__(self, modality, save_name, save_folder, save_suffix='', save_workspace=False, analyze_results=True,
+                 scenario='single_meal',
                  bolus_source='data', basal_source='data', cho_source='data', seed=1,
                  parallelize=False, n_processes=None, plot_mode=True, verbose=True):
         """
@@ -66,6 +71,10 @@ def __init__(self, modality, save_name, save_folder, save_suffix='', save_worksp
             A string to be attached as suffix to the resulting output files' name.
         save_workspace: bool, optional, default : False
             A boolean that specifies whether or not to save the results of ReplayBG in the results/workspaces directory.
+        analyze_results : bool, optional, default : True
+            A flag that specifies whether to analyze the resulting trace or not. Setting this flag to False will fasten
+            ReplayBG and it is recommended if ReplayBG will be a component of a bigger framework (e.g., to be used in an
+            iterative process).
         scenario: str, {'single-meal', 'multi-meal'}, optional, default : 'single-meal'
             A string that specifies whether the given scenario refers to a single-meal scenario or a multi-meal scenario.
         bolus_source : str, {'data', or 'dss'}, optional, default : 'data'
@@ -120,6 +129,9 @@ def __init__(self, modality, save_name, save_folder, save_suffix='', save_worksp
         # Set the save workspace flag
         self.save_workspace = save_workspace
 
+        # Set the analyze results flag
+        self.analyze_results = analyze_results
+
         # Single-meal or multi-meal scenario?
         self.scenario = scenario
 

py_replay_bg/identification/mcmc.py

@@ -37,8 +37,6 @@ class MCMC:
         Number of steps to use for the main chain.
     thin_factor: int
         Chain thin factor to use.
-    to_sample: int
-        Number of samples to generate via the copula.
     save_chains: bool
         A flag that specifies whether to save the resulting mcmc chains and copula samplers.
     callback_ncheck: int
@@ -52,7 +50,6 @@ class MCMC:
 
     def __init__(self, model,
                  n_steps=10000,
-                 to_sample=1000,
                  save_chains=False,
                  callback_ncheck=100):
         """
@@ -64,8 +61,6 @@ def __init__(self, model,
             An object that represents the physiological model hyperparameters to be used by ReplayBG.
         n_steps: int, optional, default : 10000
             Number of steps to use for the main chain.
-        to_sample: int, optional, default : 1000
-            Number of samples to generate via the copula.
         save_chains: bool, optional, default : False
             A flag that specifies whether to save the resulting mcmc chains and copula samplers.
         callback_ncheck: int, optional, default : 100
@@ -103,9 +98,6 @@ def __init__(self, model,
         # Chain thin factor to use
         self.thin_factor = int(np.ceil(n_steps / 1000))
 
-        # Number of samples to generate via the copula
-        self.to_sample = to_sample
-
         # Save the chains?
         self.save_chains = save_chains
 
@@ -159,7 +151,7 @@ def identify(self, data, rbg_data, rbg):
 
         posterior_func = self.model.log_posterior_single_meal if self.model.is_single_meal else self.model.log_posterior_multi_meal
         sampler = zeus.EnsembleSampler(self.n_walkers, self.n_dim, posterior_func, args=[rbg_data],
-                                       verbose=rbg.environment.verbose, pool=pool, maxsteps=10)
+                                       verbose=rbg.environment.verbose, pool=pool, maxsteps=1000)
         sampler.run_mcmc(start, self.n_steps, callbacks=[cb0, cb1, cb2])
         sampler.summary  # Print summary diagnostics
 
@@ -175,32 +167,35 @@ def identify(self, data, rbg_data, rbg):
         draws = dict()
         for up in range(len(rbg.model.unknown_parameters)):
             draws[rbg.model.unknown_parameters[up]] = dict()
-            draws[rbg.model.unknown_parameters[up]]['samples'] = np.empty(self.to_sample)
+            draws[rbg.model.unknown_parameters[up]]['samples_1000'] = np.empty(1000)
+            draws[rbg.model.unknown_parameters[up]]['samples_100'] = np.empty(100)
+            draws[rbg.model.unknown_parameters[up]]['samples_10'] = np.empty(10)
             draws[rbg.model.unknown_parameters[up]]['chain'] = chain[:, up]
 
-        if rbg.environment.verbose:
-            print('Extracting samples from copula')
-
-        to_be_sampled = [True] * self.to_sample
-        while any(to_be_sampled):
+        to_sample = [1000, 100, 10]
+        for nr in to_sample:
+            if rbg.environment.verbose:
+                print('Extracting samples from copula - ' + str(nr) + ' realizations')
 
-            #Get the idxs of the missing samples
-            tbs = np.where(to_be_sampled)[0]
+            to_be_sampled = [True] * nr
+            while any(to_be_sampled):
 
-            #Get the new samples
-            samples = distributions.sample(len(tbs)).to_numpy()
+                #Get the idxs of the missing samples
+                tbs = np.where(to_be_sampled)[0]
 
-            #For each sample...
-            for i in range(0, len(tbs)):
+                #Get the new samples
+                samples = distributions.sample(len(tbs)).to_numpy()
 
-                #...check if it is ok. If so...
-                if self.model.check_copula_extraction(samples[i]):
+                #For each sample...
+                for i in range(0, len(tbs)):
 
-                    # ...flag it and fill the final vector
-                    to_be_sampled[tbs[i]] = False
-                    for up in range(len(rbg.model.unknown_parameters)):
-                        draws[rbg.model.unknown_parameters[up]]['samples'][tbs[i]] = samples[i,up]
+                    #...check if it is ok. If so...
+                    if self.model.check_copula_extraction(samples[i]):
 
+                        # ...flag it and fill the final vector
+                        to_be_sampled[tbs[i]] = False
+                        for up in range(len(rbg.model.unknown_parameters)):
+                            draws[rbg.model.unknown_parameters[up]]['samples_'+str(nr)][tbs[i]] = samples[i,up]
 
         # Check physiological plausibility
         draws['physiological_plausibility'] = self.__check_physiological_plausibility(draws, data, rbg)
@@ -253,15 +248,15 @@ def __check_physiological_plausibility(self, draws, data, rbg):
         """
 
         if rbg.environment.verbose:
-            print('Running physiological plausibility checks...')
+            print('Running physiological plausibility checks (only on the 1000 samples)...')
 
         # Initialize the return vector
         physiological_plausibility = dict()
 
-        physiological_plausibility['test_1'] = np.full((self.to_sample,), True)
-        physiological_plausibility['test_2'] = np.full((self.to_sample,), True)
-        physiological_plausibility['test_3'] = np.full((self.to_sample,), True)
-        physiological_plausibility['test_4'] = np.full((self.to_sample,), True)
+        physiological_plausibility['test_1'] = np.full((1000,), True)
+        physiological_plausibility['test_2'] = np.full((1000,), True)
+        physiological_plausibility['test_3'] = np.full((1000,), True)
+        physiological_plausibility['test_4'] = np.full((1000,), True)
 
         rbg_fake = copy.copy(rbg)
         rbg_fake.model = copy.copy(rbg.model)
@@ -299,9 +294,9 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
         # Test 1: "if no insulin is injected, BG must go above 300 mg/dl in 1000 min"
         if rbg.environment.verbose:
-            iterations = tqdm(range(self.to_sample), desc='Test 1 of 4')
+            iterations = tqdm(range(1000), desc='Test 1 of 4')
         else:
-            iterations = range(0, self.to_sample)
+            iterations = range(0, 1000)
 
         # Set simulation data
         data_fake_test_1 = copy.copy(data_fake)
@@ -312,7 +307,7 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
             # set the model parameters
             for p in rbg_fake.model.unknown_parameters:
-                setattr(rbg_fake.model.model_parameters,p,draws[p]['samples'][r])
+                setattr(rbg_fake.model.model_parameters,p,draws[p]['samples_1000'][r])
             rbg_fake.model.model_parameters.kgri = rbg_fake.model.model_parameters.kempt
 
             if (rbg_fake.sensors.cgm.model == 'CGM'):
@@ -326,9 +321,9 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
         # Test 2: "if a bolus of 15 U is injected, BG should drop below 100 mg/dl"
         if rbg.environment.verbose:
-            iterations = tqdm(range(self.to_sample), desc='Test 2 of 4')
+            iterations = tqdm(range(1000), desc='Test 2 of 4')
         else:
-            iterations = range(0, self.to_sample)
+            iterations = range(0, 1000)
 
         # Set simulation data
         data_fake_test_2 = copy.copy(data_fake)
@@ -351,7 +346,7 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
             # set the model parameters
             for p in rbg_fake.model.unknown_parameters:
-                setattr(rbg_fake.model.model_parameters, p, draws[p]['samples'][r])
+                setattr(rbg_fake.model.model_parameters, p, draws[p]['samples_1000'][r])
             rbg_fake.model.model_parameters.kgri = rbg_fake.model.model_parameters.kempt
 
             if (rbg_fake.sensors.cgm.model == 'CGM'):
@@ -365,9 +360,9 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
         # Test 3: "it exists a basal insulin value such that glucose stays between 90 and 160 mg/dl", 
         if rbg.environment.verbose:
-            iterations = tqdm(range(self.to_sample), desc='Test 3 of 4')
+            iterations = tqdm(range(1000), desc='Test 3 of 4')
         else:
-            iterations = range(0, self.to_sample)
+            iterations = range(0, 1000)
 
         # Set simulation data
         data_fake_test_3 = copy.copy(data_fake)
@@ -384,7 +379,7 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
             # set the model parameters
             for p in rbg_fake.model.unknown_parameters:
-                setattr(rbg_fake.model.model_parameters, p, draws[p]['samples'][r])
+                setattr(rbg_fake.model.model_parameters, p, draws[p]['samples_1000'][r])
             rbg_fake.model.model_parameters.kgri = rbg_fake.model.model_parameters.kempt
 
             if (rbg_fake.sensors.cgm.model == 'CGM'):
@@ -418,9 +413,9 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
         # Test 4: "a variation of basal insulin of 0.01 U/h does not vary basal glucose more than 20 mg/dl"
         if rbg.environment.verbose:
-            iterations = tqdm(range(self.to_sample), desc='Test 4 of 4')
+            iterations = tqdm(range(1000), desc='Test 4 of 4')
         else:
-            iterations = range(0, self.to_sample)
+            iterations = range(0, 1000)
 
         # Set simulation data
         data_fake_test_4 = copy.copy(data_fake)
@@ -433,7 +428,7 @@ def __check_physiological_plausibility(self, draws, data, rbg):
 
             # set the model parameters
             for p in rbg_fake.model.unknown_parameters:
-                setattr(rbg_fake.model.model_parameters,p,draws[p]['samples'][r])
+                setattr(rbg_fake.model.model_parameters,p,draws[p]['samples_1000'][r])
             rbg_fake.model.model_parameters.kgri = rbg_fake.model.model_parameters.kempt
 
             if (rbg_fake.sensors.cgm.model == 'CGM'):

py_replay_bg/input_validation/__init__.py

@@ -271,30 +271,30 @@ def validate(self):
             raise Exception("'n_steps' input must be an integer.'")
 
 
-class ToSampleValidator:
+class SaveChainsValidator:
     """
-    Class for validating the 'to_sample' input parameter of ReplayBG.
+    Class for validating the 'save_chains' input parameter of ReplayBG.
     """
 
-    def __init__(self, to_sample):
-        self.to_sample = to_sample
+    def __init__(self, save_chains):
+        self.save_chains = save_chains
 
     def validate(self):
-        if not isinstance(self.to_sample, int):
-            raise Exception("'to_sample' input must be an integer.'")
+        if not isinstance(self.save_chains, bool):
+            raise Exception("'save_chains' input must be a boolean.'")
 
 
-class SaveChainsValidator:
+class AnalyzeResultsValidator:
     """
-    Class for validating the 'save_chains' input parameter of ReplayBG.
+    Class for validating the 'analyze_results' input parameter of ReplayBG.
     """
 
-    def __init__(self, save_chains):
-        self.save_chains = save_chains
+    def __init__(self, analyze_results):
+        self.analyze_results = analyze_results
 
     def validate(self):
-        if not isinstance(self.save_chains, bool):
-            raise Exception("'save_chains' input must be a boolean.'")
+        if not isinstance(self.analyze_results, bool):
+            raise Exception("'analyze_results' input must be a boolean.'")
 
 
 class CRValidator:
@@ -550,7 +550,7 @@ def __init__(self, modality, data, bw, scenario, save_name, save_suffix,
                  yts, glucose_model, pathology, exercise, seed,
                  bolus_source, basal_source, cho_source,
                  cgm_model,
-                 n_steps, to_sample, save_chains,
+                 n_steps, save_chains, analyze_results,
                  CR, CF, GT,
                  meal_generator_handler, meal_generator_handler_params,
                  bolus_calculator_handler, bolus_calculator_handler_params,
@@ -578,8 +578,8 @@ def __init__(self, modality, data, bw, scenario, save_name, save_suffix,
         self.cgm_model = cgm_model
 
         self.n_steps = n_steps
-        self.to_sample = to_sample
         self.save_chains = save_chains
+        self.analyze_results = analyze_results
 
         self.CR = CR
         self.CF = CF
@@ -654,12 +654,12 @@ def validate(self):
         # Validate the 'n_steps' input
         NStepsValidator(n_steps=self.n_steps).validate()
 
-        # Validate the 'to_sample' input
-        ToSampleValidator(to_sample=self.to_sample).validate()
-
         # Validate the 'save_chains' input
         SaveChainsValidator(save_chains=self.save_chains).validate()
 
+        # Validate the 'analyze_results' input
+        AnalyzeResultsValidator(analyze_results=self.analyze_results).validate()
+
         # Validate the 'CR' input
         CRValidator(CR=self.CR).validate()