Test

Backend/ModelCreation/CreateModel.ipynb

@@ -250,7 +250,7 @@
  "name": "python",
  "nbconvert_exporter": "python",
  "pygments_lexer": "ipython3",
- "version": "3.12.0"
+ "version": "3.10.6"
  }
  },
  "nbformat": 4,

Backend/ModelCreation/Example regional SIR.ipynb

@@ -22,7 +22,7 @@
  },
  {
  "cell_type": "code",
- "execution_count": 6,
+ "execution_count": 1,
  "metadata": {},
  "outputs": [],
  "source": [
@@ -33,7 +33,7 @@
  },
  {
  "cell_type": "code",
- "execution_count": 7,
+ "execution_count": 2,
  "metadata": {},
  "outputs": [],
  "source": [
@@ -42,43 +42,43 @@
  },
  {
  "cell_type": "code",
- "execution_count": 8,
+ "execution_count": 3,
  "metadata": {},
  "outputs": [],
  "source": [
  "model_constants = {\n",
- " \"Exposure_rate\":0.11,\n",
- " \"Infection_rate\":0.2,\n",
- " \"Recovery_rate\":0.3,\n",
+ " \"Exposure_rate\":2,\n",
+ " \"Infection_rate\":1,\n",
+ " \"Recovery_rate\":0.2,\n",
  " \"Death_proportion\":0.01\n",
  "}"
  ]
  },
  {
  "cell_type": "code",
- "execution_count": 9,
+ "execution_count": 4,
  "metadata": {},
  "outputs": [],
  "source": [
  "def epidemologyRules():\n",
  " # Frequency vs prevalance - /N but can also use a frequency propensity (fits some dieases better than others).\n",
  " exposure = SupplyChainRules.SingleLocationProductionRule(\"Region\", [\"S\"], [1], [\"E\"], [1], \n",
- " f\"{model_constants[\"Exposure_rate\"]}*I*S/(S+E+I+R)\",\n",
+ " f\"{model_constants['Exposure_rate']}*I*S/(S+E+I+R)*model_month_jan + {model_constants['Exposure_rate']}*I*S/(S+E+I+R)*model_month_feb\",\n",
  " [\"S\",\"E\",\"I\",\"R\"], \"Exposure\")\n",
  " infection = SupplyChainRules.SingleLocationProductionRule(\"Region\", [\"E\"], [1], [\"I\"], [1], \n",
- " f\"{model_constants[\"Infection_rate\"]}*E\",\n",
+ " f\"{model_constants['Infection_rate']}*E\",\n",
  " [\"E\"], \"Infection\") \n",
- " recovery = SupplyChainRules.SingleLocationProductionRule(\"Region\", [\"I\"], [1], [\"R\"], [1], \n",
- " f\"{model_constants[\"Recovery_rate\"]}*I\",\n",
+ " recovery = SupplyChainRules.SingleLocationProductionRule(\"Region\", [\"I\"], [1], [\"R\"], [1], pp9\n",
+ " f\"{model_constants['Recovery_rate']}*I\",\n",
  " [\"I\"], \"Recovery\")\n",
  " death = SupplyChainRules.ExitEntranceRule(\"Region\", \"I\", 1, \n",
- " f\"{model_constants[\"Death_proportion\"]}*{model_constants[\"Recovery_rate\"]}*I\", [\"I\"], \"Death\")\n",
+ " f\"{model_constants['Death_proportion']}*{model_constants['Recovery_rate']}*I\", [\"I\"], \"Death\")\n",
  " return [exposure, infection, recovery, death]"
  ]
  },
  {
  "cell_type": "code",
- "execution_count": 10,
+ "execution_count": 5,
  "metadata": {},
  "outputs": [],
  "source": [
@@ -94,49 +94,53 @@
  "def sirLocations():\n",
  " all_locations = []\n",
  " # Midpoints from wikipedia, South East and London uses South East midpoint (combined to match DEFRA reporting)\n",
- " region_infos = [[54.075, -2.75, \"North East\"], [55, -1.87, \"North West\"], [53.566667, -1.2, \"Yorkshire & The Humber\"], [52.98, -0.75, \"East Midlands\"], [52.478861, -2.256306, \"West Midlands\"], \n",
+ " region_infos = [[0, 0, \"Test Region\"], [54.075, -2.75, \"North East\"], [55, -1.87, \"North West\"], [53.566667, -1.2, \"Yorkshire & The Humber\"], [52.98, -0.75, \"East Midlands\"], [52.478861, -2.256306, \"West Midlands\"], \n",
  " [52.24, 0.41, \"East of England\"], [51.515447, -0.09214, \"London\"], [51.3, -0.8, \"South East\"], [50.96, -3.22, \"South West\"], [56.816738, -4.183963, \"Scotland\"], \n",
  " [52.33022, -3.766409,\"Wales\"]]\n",
  " # ONS mid 2022 estimates\n",
- " region_populations = [2683040, 7516113, 5541262, 4934939, 6021653, 6398497, 8866180, 9379833, 5764881, 5447700, 3131640]\n",
+ " region_populations = [10000, 2683040, 7516113, 5541262, 4934939, 6021653, 6398497, 8866180, 9379833, 5764881, 5447700, 3131640]\n",
  "\n",
  " for region_index, region_info in enumerate(region_infos):\n",
  " region = Region(*region_info)\n",
- " region.setInitialConditions({\"S\":region_populations[region_index]})\n",
+ " region.setInitialConditions({\"S\":region_populations[region_index],\"I\":100})\n",
  " all_locations.append(region)\n",
  "\n",
  " return all_locations"
  ]
  },
  {
  "cell_type": "code",
- "execution_count": 14,
+ "execution_count": 6,
  "metadata": {},
  "outputs": [
  {
  "name": "stdout",
  "output_type": "stream",
  "text": [
- "0.11*I*S/(S+E+I+R)\n",
- "0.0275000000000000\n",
- "0.2*E\n",
- "0.200000000000000\n",
- "0.3*I\n",
- "0.300000000000000\n",
- "0.01*0.3*I\n",
- "0.00300000000000000\n",
+ "2*I*S/(S+E+I+R)*model_month_jan + 2*I*S/(S+E+I+R)*model_month_feb\n",
+ "1*E\n",
+ "0.2*I\n",
+ "0.01*0.2*I\n",
+ "1.00000000000000\n",
  "['S', 'E']\n",
+ "1.00000000000000\n",
  "['E', 'I']\n",
+ "0.200000000000000\n",
  "['I', 'R']\n",
+ "0.00200000000000000\n",
  "['I']\n",
- "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]}\n",
- "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]}\n",
- "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]}\n",
- "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]}\n",
- "['0.11*x1*x3/(x3+x0+x1+x2)']\n",
- "['0.2*x0']\n",
- "['0.3*x1']\n",
- "['0.01*0.3*x1']\n"
+ "1.00000000000000\n",
+ "1.00000000000000\n",
+ "0.200000000000000\n",
+ "0.00200000000000000\n",
+ "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]}\n",
+ "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]}\n",
+ "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]}\n",
+ "MATCHED {0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]}\n",
+ "['2*x1*x3/(x3+x0+x1+x2)*x10 + 2*x1*x3/(x3+x0+x1+x2)*x9']\n",
+ "['1*x0']\n",
+ "['0.2*x1']\n",
+ "['0.01*0.2*x1']\n"
  ]
  }
  ],
@@ -176,7 +180,7 @@
  "name": "python",
  "nbconvert_exporter": "python",
  "pygments_lexer": "ipython3",
- "version": "3.12.0"
+ "version": "3.10.6"
  }
  },
  "nbformat": 4,

Backend/ModelCreation/ModelClasses.py

@@ -4,7 +4,7 @@ class Classes:
  def __init__(self, allow_base_model_classes = True):
  self.class_def_dict = {}
  self.class_names = set([])
- self.base_model_string = "model."
+ self.base_model_string = "model_"
  if allow_base_model_classes:
  for built_in_class in self.returnBuiltInClasses():
  self.addClass(*built_in_class)
@@ -16,17 +16,20 @@ def addClass(self, class_name, class_measurement_unit, class_restriction = "None
  self.class_def_dict[class_name] = {"class_measurement_unit":class_measurement_unit, "class_restriction":class_restriction}
  else:
  raise(ValueError(f"Duplicate type {class_name} provided."))
+
  def returnBuiltInClasses(self):
  builtin_classes = [[f"{self.base_model_string}day", "Day of year"], 
  [f"{self.base_model_string}hour","Hour of day", "integer"], 
- [f"{self.base_model_string}month", "Month of year", "integer"]]
+ [f"{self.base_model_string}months", "Month of year", "integer"]]
  months = [["jan", "January"], ["feb", "February"], ["mar", "March"], ["apr", "April"],
  ["may", "May"], ["jun", "June"], ["jul", "July"], ["aug", "August"], 
  ["sept", "September"], ["oct", "October"], ["nov", "November"], ["dec", "December"]
  ]
  for month in months:
- builtin_classes.append([f"{self.base_model_string}month.{month[0]}", month[1], "indictator"])
+ builtin_classes.append([f"{self.base_model_string}month_{month[0]}", month[1], "indictator"])
+ builtin_classes.sort()
  return builtin_classes
+
  def writeClassJSON(self, filepath):
 
  json_classes = json.dumps(self.class_def_dict, indent=4, sort_keys=True)

Backend/ModelCreation/ModelDefinition.py

@@ -50,36 +50,11 @@ def matchRules(self):
  additional_classes = []
  if self.builtin_classes:
  additional_classes = ModelClasses.Classes().returnBuiltInClasses()
- RuleMatching.writeMatchedRuleJSON(self.rules, self.locations, f"{self.model_folder}{self.matched_rules_filename}", 
+ RuleMatching.writeMatchedRuleJSON(self.rules, self.locations, f"{self.model_folder}{self.matched_rules_filename}",
  additional_classes)
 
  def build(self):
  self.defineClasses()
  self.createLocations()
  self.createRules()
  self.matchRules()
-
-
-def supplyChainRules():
- transport_nitrogen = SupplyChainRules.TransportRule(source="ChemicalPlant NitrogenPlant",
- target="ChemicalPlant AmmoniaPlant", 
- transport_class="N2", propensities=["1","1"], transport_amount=2,
- propensity_classes=[["N2"], ["N2"]], rule_name="Transport Nitrogen")
- manufacture_ammonia = SupplyChainRules.SingleLocationProductionRule(target="ChemicalPlant AmmoniaPlant",
- reactant_classes=["N2","H2"], reactant_amount=[1,3], 
- product_classes=["NH4"],product_amount=[1],propensity="N2*H2",
- propensity_classes=["N2", "H2"], rule_name="Make Ammonia")
- return [transport_nitrogen, manufacture_ammonia]
-
-def supplyChainLocations():
- ammonia_plant = SupplyChainLocations.ChemicalPlant(1,1,"Example plant")
- nitrogen_plant = SupplyChainLocations.ChemicalPlant(1,1,"Example plant 2")
-
- ammonia_plant.addAmmoniaManufacturing()
- nitrogen_plant.addNitrogenManufacturing()
- return [ammonia_plant, nitrogen_plant]
-
-#supplyChainClasses = [["NH4", "tonnes"], ["N2", "m^3"], ["H2", "m^3"], ["CH4", "m^3"]]
-
-#model = ModelDefinition(supplyChainClasses, supplyChainLocations, supplyChainRules, model_folder="Backend/ModelFiles/")
-#model.build()

Backend/ModelCreation/ModelRules.py

@@ -2,6 +2,17 @@
 import sympy
 import json
 
+def returnSympyClassVarsDict(classes):
+ symbols_str = ""
+ for class_label in classes:
+ symbols_str += class_label+" "
+ symbols = sympy.symbols(symbols_str)
+ if not isinstance(symbols, (list, tuple)):
+ symbols = [symbols]
+ symbols_dict = {class_label:symbols[index] for index, class_label in enumerate(classes)}
+ return symbols_dict
+
+
 class Rule:
  def __init__(self, rule_name, targets):
 
@@ -27,16 +38,6 @@ def addLinearStoichiomety(self, target_indices, stoichiometies, required_target_
  self.stoichiometry_classes[index] = required_target_classes[i]
  else:
  raise(ValueError(f"Unrecognised stoichiomety of type {type(stoichiometies[i])}, for target index {index}"))
-
- def returnSympyClassVarsDict(self, classes):
- symbols_str = ""
- for class_label in classes:
- symbols_str += class_label+" "
- symbols = sympy.symbols(symbols_str)
- if not isinstance(symbols, (list, tuple)):
- symbols = [symbols]
- symbols_dict = {class_label:symbols[index] for index, class_label in enumerate(classes)}
- return symbols_dict
 
  def addSimplePropensityFunction(self, target_indices, values, required_target_classes):
  # Accepts matrix or constant values at the moment
@@ -47,14 +48,9 @@ def addSimplePropensityFunction(self, target_indices, values, required_target_cl
  if not self.propensities[index] is None:
  raise(ValueError(f"Overwriting already set propensity is forbidden. Target location {self.targets[index]} at position {str(index+1)}"))
 
- if isinstance(value, str):
- # TODO validate formula here
- # We expect that the array has a single entry.
- symbols = self.returnSympyClassVarsDict(required_target_classes[index])
- sympy_formula = sympy.parse_expr(value, local_dict=symbols)
- self.validateFormula(sympy_formula, symbols)
- else:
+ if not isinstance(value, str):
  raise(ValueError(f"Unrecognised propensity function of type {type(values[i])}, for target index {index}"))
+
  self.propensity_classes[index] = required_target_classes[i]
  self.propensities[index] = value
 
@@ -69,7 +65,7 @@ def validateFormula(self, formula, class_symbols, safe_num = 1):
  assert(isinstance(res, (sympy.core.numbers.Float, sympy.core.numbers.Zero)))
  return True
 
- def checkRuleDefinition(self):
+ def checkRuleDefinition(self, builtin_class_symbols):
  for i in range(len(self.targets)):
  if self.stoichiometies[i] is None:
  raise(ValueError(f"Incorrect rule definition for {self.rule_name}\nThe Location type {self.targets[i]} at rule position {str(i+1)} has no defined stochiometry."))
@@ -80,6 +76,11 @@ def checkRuleDefinition(self):
  elif self.stoichiometry_classes[i] is None:
  raise(ValueError(f"Incorrect rule definition for {self.rule_name}\nThe Location type {self.targets[i]} at rule position {str(i+1)} has no defined stochiometry class requirement."))
 
+ for index in range(len(self.propensities)):
+ symbols = returnSympyClassVarsDict(self.propensity_classes[index]) | builtin_class_symbols
+ sympy_formula = sympy.parse_expr(self.propensities[index], local_dict=symbols)
+ self.validateFormula(sympy_formula, symbols)
+
  def mergeClassLists(self):
  # Maps new index to class label
  self.rule_classes = []
@@ -96,8 +97,8 @@ def mergeClassLists(self):
  self.stoichiometies[i] = list(new_stoichiometry)
  self.rule_classes.append(tmp_rule_class_dict)
 
- def returnRuleDict(self):
- self.checkRuleDefinition()
+ def returnRuleDict(self, builtin_class_symbols):
+ self.checkRuleDefinition(builtin_class_symbols)
  self.mergeClassLists()
 
  # UNPACK DEF TO SPECIFIC LOCATION
@@ -109,6 +110,15 @@ class Rules:
  def __init__(self, defined_classes):
  self.rules = []
  self.defined_classes = defined_classes
+ self.model_prefix = "model_"
+
+ builtin_classes = []
+
+ for classes in defined_classes:
+ if self.model_prefix in classes:
+ builtin_classes.append(classes)
+
+ self.builtin_symbols = returnSympyClassVarsDict(builtin_classes)
 
  def addRule(self, rule:Rule):
  if isinstance(rule, Rule):
@@ -122,7 +132,7 @@ def addRules(self, rules):
 
  def checkRules(self):
  for rule in self.rules:
- rule.checkRuleDefinition()
+ rule.checkRuleDefinition(self.builtin_symbols)
 
  for propensity_class_index in range(len(rule.propensity_classes)):
  location_propensity_class = rule.propensity_classes[propensity_class_index]
@@ -138,12 +148,11 @@ def checkRules(self):
  raise ValueError(f"Class required for the stoichiometry, {loc_stoich_class} not defined (Rule name: {rule.rule_name})")
  return True
 
-
  def writeJSON(self, filename):
  self.checkRules()
  rules_dict = {}
  for i, rule in enumerate(self.rules):
- rule_dict = rule.returnRuleDict()
+ rule_dict = rule.returnRuleDict(self.builtin_symbols)
  rules_dict[i] = rule_dict
 
  json_rules = json.dumps(rules_dict, indent=4, sort_keys=True)

Backend/ModelCreation/RuleMatching.py

@@ -81,8 +81,8 @@ def obtainPropensity(rule, locations, builtin_classes):
  for label_i in list(new_label_mapping.keys()):
  new_propensity = new_propensity.replace(new_label_mapping[label_i], f"x{label_i}")
  # Add the built in class at the end of the location classes
- for builtin_class in builtin_classes:
- new_propensity = new_propensity.replace(builtin_class[0], f"x{label_i+len(new_label_mapping)}")
+ for built_in_i, builtin_class in enumerate(builtin_classes):
+ new_propensity = new_propensity.replace(builtin_class[0], f"x{built_in_i+len(new_label_mapping)}")
  new_propensities.append(new_propensity)
  print(new_propensities)
  return new_propensities