Merge pull request #1 from QuanticPony/Add-my-tfg-model

Add my tfg model

.gitignore

@@ -128,4 +128,5 @@ dmypy.json
 # Pyre type checker
 .pyre/
 
-.*
+.*
+*.data

README.md

@@ -27,20 +27,52 @@ limitations under the License. -->
 compartmental
 </h1>
 <h2 align="center">
-.
-</h2><br></br>
-<h3 align="center">
-.
-</h3><br></br>
-
-
+Utility tools for Approximate Bayesian computation on compartmental models 
+</h2>
 
+<br>
 <div align="center">
 
 <a href="https://quanticpony.github.io/compartmental/">
 <img src=https://img.shields.io/github/deployments/QuanticPony/compartmental/github-pages?label=documentation>
 </a>
-<br></br>
+<br></br></br>
+<h3 align="center">
+
+</h3>
+
+
+
+
 
 </div>
-<br></br>
+
+# Instalation
+**compartmental** releases are available as wheel packages on [PyPI](https://pypi.org/project/compartmental/). You can install the last version using `pip`:
+```
+pip install compartmental
+```
+
+
+# Documentation
+Documentations is automatically generated from code on main push and hosted in github-pages [here](https://quanticpony.github.io/compartmental/).
+
+# Help
+Just open an issue with the `question` tag ([or clic here](https://github.com/QuanticPony/compartmental/issues/new?assignees=QuanticPony&labels=question&template=question.md&title=)), I would love to help!
+
+# Contributing
+You can contribute with:
+
+* Examples
+* Documentation
+* [Bug report/fix](https://github.com/QuanticPony/compartmental/issues/new?assignees=QuanticPony&labels=bug&template=bug_report.md&title=)
+* [Features](https://github.com/QuanticPony/compartmental/issues/new?assignees=QuanticPony&labels=new-feature&template=feature_request.md&title=)
+* Code
+
+Even only feedback is greatly apreciated. 
+
+Just create an issue and let me know you want to help! 
+
+
+# Licensing
+**compartmental** is released under the **Apache License Version 2.0**.

compartmental/__init__.py

@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-__version__ = "0.0.3"
+__version__ = "0.1.0"
 _CUPY_MODE_: bool
 
 from . import generic_model
@@ -38,4 +38,4 @@ def update_packages(CNP):
 
 use_numpy()
 
-GenericModel = generic_model.GenericModel
+from .generic_model import GenericModel

compartmental/generic_model.py

@@ -12,19 +12,20 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import __future__
-from io import TextIOWrapper
+from __future__ import annotations
 from typing import TYPE_CHECKING, Any
 
 if TYPE_CHECKING:
     import numpy as CNP
 
-from .parameters import ParametersManager
+
 from .util import *
+from .parameters import ParametersManager
+
 import copy
 
 class GenericModel:
-    """Creates a compartimental model from a dictionary and setting an `evolve` method.
+    """Creates a compartmental model from a dictionary and setting an `evolve` method.
     """
 
     def get_all_params_names(self):
@@ -46,7 +47,7 @@ def __init__(self, configuration: dict[str, Any]):
 
         self.param_to_index: dict[str, int] = { k:i for i,k in enumerate(self.configuration["params"].keys()) }
         self.fixed_param_to_index: dict[str, int] = { k:i for i,k in enumerate(self.configuration["fixed_params"].keys()) }
-        self.compartiment_name_to_index: dict[str, int] = { k:i for i,k in enumerate(self.configuration["compartiments"].keys()) }
+        self.compartment_name_to_index: dict[str, int] = { k:i for i,k in enumerate(self.configuration["compartments"].keys()) }
 
 
     def populate_model_parameters(self, **kargs):
@@ -58,60 +59,54 @@ def populate_model_parameters(self, **kargs):
         # Set offset value if it is a str reference
         if isinstance(REFERENCE_OFFSET, str):
             self.configuration["params"][REFERENCE_OFFSET].update({"type":"int"})
-
-        N_SIMULATIONS = self.configuration["simulation"]["n_simulations"]
-        self.params = CNP.zeros(
-            (len(self.configuration["params"]), N_SIMULATIONS), dtype=CNP.float64
-        )
-        self.fixed_params = CNP.zeros(
-            (len(self.configuration["fixed_params"]), 1), dtype=CNP.float64
-        )
+
+        parameter_manager.populate_params(**kargs)
+        parameter_manager.populate_fixed_params()
 
         for param in self.configuration["params"].keys():
-            setattr(self, param, self.params[self.param_to_index[param]])
+            setattr(self, param, self.params[param])
 
         for fparam in self.configuration["fixed_params"].keys():
             setattr(self, fparam, self.fixed_params[self.fixed_param_to_index[fparam]])
 
         if isinstance(REFERENCE_OFFSET, str):
-            self.reference_offset = self.params[self.param_to_index[REFERENCE_OFFSET]]
+            self.reference_offset = self.params[REFERENCE_OFFSET]
         else:
             self.reference_offset = 0
 
-        parameter_manager.populate_params(self.params, **kargs)
-        parameter_manager.populate_fixed_params(self.fixed_params)
+
 
 
-    def populate_model_compartiments(self, **kargs):
-        """Populates compartiments array. Assigns shortcuts to call them by their name as an attribute.
+    def populate_model_compartments(self, **kargs):
+        """Populates compartments array. Assigns shortcuts to call them by their name as an attribute.
         """
         N_SIMULATIONS = self.configuration["simulation"]["n_simulations"]
         self.state = CNP.zeros(
-            (len(self.configuration["compartiments"]), N_SIMULATIONS), dtype=CNP.float64
+            (len(self.configuration["compartments"]), N_SIMULATIONS), dtype=CNP.float64
         )
         self.log_diff = CNP.zeros((N_SIMULATIONS, 1), dtype=CNP.float64)
 
-        for c,i in self.compartiment_name_to_index.items():
-            C = self.configuration["compartiments"][c]
+        for c,i in self.compartment_name_to_index.items():
+            C = self.configuration["compartments"][c]
             initial_value = C["initial_value"]
             if isinstance(initial_value, str):
                 if initial_value in self.param_to_index.keys():
-                    self.state[i,:] = self.params[self.param_to_index[initial_value]]
+                    self.state[i,:] = self.params[initial_value]
                 continue
             self.state[i,:] = initial_value
 
-        for c,i in self.compartiment_name_to_index.items():
-            C = self.configuration["compartiments"][c]
-            minus = C.get("minus_compartiments", False)
+        for c,i in self.compartment_name_to_index.items():
+            C = self.configuration["compartments"][c]
+            minus = C.get("minus_compartments", False)
             if not minus:
                 continue
             if not isinstance(minus, list):
                 minus = [minus]
             for m in minus:
-                self.state[i,:] -= self.state[self.compartiment_name_to_index[m],:]
+                self.state[i,:] -= self.state[self.compartment_name_to_index[m],:]
 
-        for comp in self.configuration["compartiments"].keys():
-            setattr(self, comp, self.state[self.compartiment_name_to_index[comp]])
+        for comp in self.configuration["compartments"].keys():
+            setattr(self, comp, self.state[self.compartment_name_to_index[comp]])
 
 
     def evolve(self, step, *args, **kargs):
@@ -135,12 +130,15 @@ def get_diff(self, step, reference, reference_mask):
         Returns:
             (list[float]): Distance from simulations to reference(s).
         """
-        index = CNP.clip(step + CNP.int64(self.reference_offset), 0, self.N_STEPS-1)
-        diff = CNP.absolute(CNP.take(self.state, reference_mask, 0)[0].T-reference[index])
+        index = step + self.reference_offset
+        # To only take the diff on the same range for all simulations
+        diff = CNP.absolute(CNP.take(self.state, reference_mask, 0)[0].T-reference[CNP.clip(index, 0, self.N_STEPS-1)]) * \
+               ((self.reference_offset.max()<=index) * (index<=self.N_STEPS))
+
         return CNP.log(diff + 1)
 
 
-    def _internal_run_(self, inner, inner_args: list, outer, outer_args:list,  reference, save_file:str, *args, **kargs):
+    def _internal_run_(self, inner, inner_args: list, outer, outer_args:list,  reference, save_file:str, *args, exclude_pupulate:bool=False, **kargs):
         """Internal function that executes the model.
 
         Args:
@@ -150,19 +148,26 @@ def _internal_run_(self, inner, inner_args: list, outer, outer_args:list,  refer
             outer_args (list): Args given to `outer`.
             reference (list[list[float]]): Reference values used to compare with the simulation.
             save_file (str): Filename of path to file.
+            exclude_populate (bool, optional): If `False` params and compartments are populated with random values. Defaults to False.
         """
         N_EXECUTIONS = self.configuration["simulation"]["n_executions"]
-        REFERENCE_OFFSET = self.configuration["reference"].get("offset", 0)
         self.N_STEPS = self.configuration["simulation"]["n_steps"]
-        
+
         for execution in range(N_EXECUTIONS):
             progress_bar(f"Model running: ", execution, N_EXECUTIONS, len=min(20, max(N_EXECUTIONS,5)))
-            self.log_diff[:] = 0
-            self.populate_model_parameters(**kargs)
-            self.populate_model_compartiments(**kargs)
 
-            for step in range(self.N_STEPS):
+            if not exclude_pupulate:
+                self.populate_model_parameters(**kargs)
+                self.populate_model_compartments(**kargs)
+
+            self._min_offset_: int = self.reference_offset.min()
+            self.log_diff[:] = 0
+
+            # for step in range(self.N_STEPS):
+            step = CNP.int64(0)
+            while (self.reference_offset + step < self.N_STEPS).any():
                 inner(self, step, reference, *inner_args, **kargs)
+                step += 1
             outer(self, *outer_args, execution_number=execution, **kargs)
 
         progress_bar(f"Model running: ", N_EXECUTIONS, N_EXECUTIONS, len=min(20, max(N_EXECUTIONS,5)), end='\n')
@@ -177,7 +182,7 @@ def run_no_diff(self, save_file: str, *args, **kargs):
         self._internal_run_(
             self.evolve, args, 
             save_parameters_no_diff, (save_file, self.param_to_index.keys(),  self.params), 
-            None, save_file, 
+            None, save_file,
             *args, **kargs
         )
 
@@ -190,18 +195,18 @@ def run(self, reference, save_file: str, *args, **kargs):
             save_file (str): Filename of path to file.
         """
 
-        reference_mask = CNP.array([self.compartiment_name_to_index[c] for c in self.configuration["reference"]["compartiments"]])
+        reference_mask = CNP.array([self.compartment_name_to_index[c] for c in self.configuration["reference"]["compartments"]])
 
         def inner(model, step, reference, reference_mask, *args, **kargs):
             model.evolve(model, step, *args, **kargs)
             self.log_diff[:,0] += model.get_diff(step, reference, reference_mask)
 
-        def outer(model, save_file, *args, **kargs):
+        def outer(model, save_file, *args, execution_number, **kargs):
             best_params, best_log_diff = get_best_parameters(model.params, model.log_diff, model.configuration["results"]["save_percentage"])
-            save_parameters(save_file, model.param_to_index.keys(), best_params, best_log_diff)
+            save_parameters(save_file, model.param_to_index.keys(), best_params, best_log_diff, execution_number=execution_number)
 
         self._internal_run_(
-            inner, (reference_mask,), 
+            inner, (reference_mask, *args), 
             outer, (save_file,), 
             reference, save_file, 
             *args, **kargs

compartmental/parameters/parameters_manager.py

@@ -25,58 +25,77 @@ def __init__(self, configuration: dict[str, Any], model):
         self.configuration = configuration
         self.model = model
 
-    def populate_params(self, params, **kargs):
+    def populate_params(self, **kargs):
         """Populates GenericModel.params array. If specific values are given in `kargs` those are used.
 
-        Args:
-            params (ndarray): Params array.
-
         Raises:
             Exception: If the model configuration could not fill all the array.
         """
         counter = 0
+        dtype = []
+
+        dtype = [[]]*len(self.model.param_to_index)
+        N_SIMULATIONS = self.model.configuration["simulation"]["n_simulations"]
+
         for param_name, param_config in self.configuration["params"].items():
             param_index = self.model.param_to_index.get(param_name)
             if param_index is not None:
-                counter += 1
                 TYPE = param_config.get("type", "float64")
-
-                if value := kargs.get(param_name, False):
-                    params[param_index] = value
-                    continue
-
                 if TYPE == 'int':
-                    params[param_index] = CNP.random.randint(
-                        param_config["min"], 
-                        param_config["max"]+1, 
-                        self.configuration["simulation"]["n_simulations"])
-
+                    dtype[param_index] = (param_name, CNP.int64)
                 if TYPE == 'float64':
-                    params[param_index] = CNP.random.random(
-                        self.configuration["simulation"]["n_simulations"]) \
-                            * (param_config["max"] - param_config["min"]) \
-                            + param_config["min"]
+                    dtype[param_index] = (param_name, CNP.float64)
+
+
+        # setattr(self.model, "params", CNP.zeros(
+        #     (len(self.configuration["params"]), N_SIMULATIONS), dtype=dtype
+        # ))
+        setattr(self.model, "params", CNP.zeros(
+            (N_SIMULATIONS), dtype=dtype
+        ))
+
+        for param_name, param_config in self.configuration["params"].items():
+            counter += 1
+            TYPE = param_config.get("type", "float64")
+
+            if value := kargs.get(param_name, False):
+                self.model.params[param_name] = value
+
+                continue
+
+            if TYPE == 'int':
+                self.model.params[param_name] = CNP.random.randint(
+                    param_config["min"], 
+                    param_config["max"]+1, 
+                    self.configuration["simulation"]["n_simulations"])
+
+            if TYPE == 'float64':
+                self.model.params[param_name] = CNP.random.random(
+                    self.configuration["simulation"]["n_simulations"]) \
+                        * (param_config["max"] - param_config["min"]) \
+                        + param_config["min"]
 
         if counter < len(self.model.param_to_index.keys()):
             raise Exception("Parameters array could not be correctly created with current options.")
 
 
-    def populate_fixed_params(self, fixed_params, **kargs):
+    def populate_fixed_params(self, **kargs):
         """Populates GenericModel.fixed_params with the configuration values.
         If specific values are given in `kargs` those are used.
 
-        Args:
-            fixed_params (ndarray): Fixed parameters.
-
         Raises:
             Exception: If the model configuration could not fill all the array.
         """
+        setattr(self.model, "fixed_params", CNP.zeros(
+            (len(self.configuration["fixed_params"]), 1), dtype=CNP.float64
+        ))
+
         counter = 0
         for fparam_name, value in self.configuration["fixed_params"].items():
             fparam_index = self.model.fixed_param_to_index.get(fparam_name)
             if fparam_index is not None:
                 counter += 1 
-                fixed_params[fparam_index] = kargs.get(fparam_name, value)
+                self.model.fixed_params[fparam_index] = kargs.get(fparam_name, value)
 
         if counter < len(self.model.fixed_param_to_index.keys()):
             raise Exception("Fixed parameters array could not be correctly created with current options.")