From 5330423b6dd7741a68c2abf5fd9efd0a50fd32f0 Mon Sep 17 00:00:00 2001 From: "Zucker, Jeremy D" Date: Mon, 28 Oct 2024 16:03:54 -0700 Subject: [PATCH] Getting more interesting compilation errors now --- docs/source/hierarchical_sir_model.ipynb | 326 ++++++++++-------- .../mira_integration/compiled_dynamics.py | 9 +- pyciemss/mira_integration/distributions.py | 78 ++++- tests/fixtures.py | 20 +- tests/test_compiled_dynamics.py | 7 +- 5 files changed, 278 insertions(+), 162 deletions(-) diff --git a/docs/source/hierarchical_sir_model.ipynb b/docs/source/hierarchical_sir_model.ipynb index 13c90449..da9d8622 100644 --- a/docs/source/hierarchical_sir_model.ipynb +++ b/docs/source/hierarchical_sir_model.ipynb @@ -9,36 +9,21 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 15, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model parameter beta has a sympy expression for distribution parameter shape which has value beta_mean*gamma_mean\n", - " and free symbol gamma_mean\n", - " and free symbol beta_mean\n", - "Model parameter beta has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter gamma has a sympy expression for distribution parameter shape which has value gamma_mean\n", - " and free symbol gamma_mean\n", - "Model parameter gamma has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter beta_mean has a sympy expression for distribution parameter alpha which has value gamma_mean\n", - " and free symbol gamma_mean\n", - "Model parameter beta_mean has a sympy expression for distribution parameter beta which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter alpha which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter beta which has value 10\n" + "ename": "ValidationError", + "evalue": "2 validation errors for Initial\nconcept\n Field required [type=missing, input_value={'name': 'S', 'value': 1.0}, input_type=dict]\n For further information visit https://errors.pydantic.dev/2.9/v/missing\nexpression\n Field required [type=missing, input_value={'name': 'S', 'value': 1.0}, input_type=dict]\n For further information visit https://errors.pydantic.dev/2.9/v/missing", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mValidationError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[15], line 82\u001b[0m\n\u001b[1;32m 78\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \\\n\u001b[1;32m 79\u001b[0m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m*\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 80\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[0;32m---> 82\u001b[0m sort_mira_dependencies(\u001b[43mtest_acyclic_distribution_expressions\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m)\n", + "Cell \u001b[0;32mIn[15], line 56\u001b[0m, in \u001b[0;36mtest_acyclic_distribution_expressions\u001b[0;34m()\u001b[0m\n\u001b[1;32m 30\u001b[0m gamma \u001b[38;5;241m=\u001b[39m Parameter(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m 31\u001b[0m distribution\u001b[38;5;241m=\u001b[39mDistribution(\u001b[38;5;28mtype\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mInverseGamma1\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 32\u001b[0m parameters\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m: sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m),\n\u001b[1;32m 33\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mscale\u001b[39m\u001b[38;5;124m'\u001b[39m: sympy\u001b[38;5;241m.\u001b[39mFloat(\u001b[38;5;241m0.01\u001b[39m)}))\n\u001b[1;32m 35\u001b[0m \u001b[38;5;66;03m# Make an SIR model with beta and gamma in rate laws\u001b[39;00m\n\u001b[1;32m 36\u001b[0m sir_model \u001b[38;5;241m=\u001b[39m TemplateModel(\n\u001b[1;32m 37\u001b[0m templates\u001b[38;5;241m=\u001b[39m[\n\u001b[1;32m 38\u001b[0m ControlledConversion(\n\u001b[1;32m 39\u001b[0m subject\u001b[38;5;241m=\u001b[39mConcept(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS\u001b[39m\u001b[38;5;124m'\u001b[39m),\n\u001b[1;32m 40\u001b[0m outcome\u001b[38;5;241m=\u001b[39mConcept(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m),\n\u001b[1;32m 41\u001b[0m controller\u001b[38;5;241m=\u001b[39mConcept(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m),\n\u001b[1;32m 42\u001b[0m rate_law\u001b[38;5;241m=\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;241m*\u001b[39m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;241m*\u001b[39m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 43\u001b[0m ),\n\u001b[1;32m 44\u001b[0m NaturalConversion(\n\u001b[1;32m 45\u001b[0m subject\u001b[38;5;241m=\u001b[39mConcept(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m),\n\u001b[1;32m 46\u001b[0m outcome\u001b[38;5;241m=\u001b[39mConcept(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mR\u001b[39m\u001b[38;5;124m'\u001b[39m),\n\u001b[1;32m 47\u001b[0m rate_law\u001b[38;5;241m=\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;241m*\u001b[39m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 48\u001b[0m ),\n\u001b[1;32m 49\u001b[0m ],\n\u001b[1;32m 50\u001b[0m parameters\u001b[38;5;241m=\u001b[39m{\n\u001b[1;32m 51\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m: beta,\n\u001b[1;32m 52\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma\u001b[39m\u001b[38;5;124m'\u001b[39m: gamma,\n\u001b[1;32m 53\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m: beta_mean,\n\u001b[1;32m 54\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m: gamma_mean,\n\u001b[1;32m 55\u001b[0m },\n\u001b[0;32m---> 56\u001b[0m initials\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mS\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[43mInitial\u001b[49m\u001b[43m(\u001b[49m\u001b[43mname\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mS\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mvalue\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1.0\u001b[39;49m\u001b[43m)\u001b[49m, \n\u001b[1;32m 57\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m: Initial(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mI\u001b[39m\u001b[38;5;124m'\u001b[39m, value\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.0\u001b[39m),\n\u001b[1;32m 58\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mR\u001b[39m\u001b[38;5;124m'\u001b[39m: Initial(name\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mR\u001b[39m\u001b[38;5;124m'\u001b[39m, value\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.0\u001b[39m)}\n\u001b[1;32m 59\u001b[0m )\n\u001b[1;32m 61\u001b[0m model \u001b[38;5;241m=\u001b[39m Model(sir_model)\n\u001b[1;32m 62\u001b[0m pn_json \u001b[38;5;241m=\u001b[39m template_model_to_petrinet_json(sir_model)\n", + "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/pydantic/main.py:212\u001b[0m, in \u001b[0;36mBaseModel.__init__\u001b[0;34m(self, **data)\u001b[0m\n\u001b[1;32m 210\u001b[0m \u001b[38;5;66;03m# `__tracebackhide__` tells pytest and some other tools to omit this function from tracebacks\u001b[39;00m\n\u001b[1;32m 211\u001b[0m __tracebackhide__ \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[0;32m--> 212\u001b[0m validated_self \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m__pydantic_validator__\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvalidate_python\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mself_instance\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 213\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m validated_self:\n\u001b[1;32m 214\u001b[0m warnings\u001b[38;5;241m.\u001b[39mwarn(\n\u001b[1;32m 215\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mA custom validator is returning a value other than `self`.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 216\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mReturning anything other than `self` from a top level model validator isn\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mt supported when validating via `__init__`.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 217\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mSee the `model_validator` docs (https://docs.pydantic.dev/latest/concepts/validators/#model-validators) for more details.\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m 218\u001b[0m category\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 219\u001b[0m )\n", + "\u001b[0;31mValidationError\u001b[0m: 2 validation errors for Initial\nconcept\n Field required [type=missing, input_value={'name': 'S', 'value': 1.0}, input_type=dict]\n For further information visit https://errors.pydantic.dev/2.9/v/missing\nexpression\n Field required [type=missing, input_value={'name': 'S', 'value': 1.0}, input_type=dict]\n For further information visit https://errors.pydantic.dev/2.9/v/missing" ] - }, - { - "data": { - "text/plain": [ - "['gamma_mean', 'gamma', 'beta_mean', 'beta']" - ] - }, - "execution_count": 14, - "metadata": {}, - "output_type": "execute_result" } ], "source": [ @@ -50,6 +35,7 @@ "from mira.sources.amr import model_from_json, model_from_url\n", "from pyro.infer.inspect import get_dependencies\n", "from pyciemss.compiled_dynamics import CompiledDynamics\n", + "from pyciemss.mira_integration.distributions import sort_mira_dependencies\n", "import torch\n", "import json\n", "from chirho.dynamical.handlers.solver import TorchDiffEq\n", @@ -74,6 +60,93 @@ " distribution=Distribution(type=\"InverseGamma1\",\n", " parameters={'shape': sympy.Symbol('gamma_mean'),\n", " 'scale': sympy.Float(0.01)}))\n", + " c = {\n", + " \"S\": Concept(name=\"S\", units=person_units(), identifiers={\"ido\": \"0000514\"}), # susceptible\n", + " \"I\": Concept(name=\"I\", units=person_units(), identifiers={\"ido\": \"0000511\"}), # infectious\n", + " \"R\": Concept(name=\"R\", units=person_units(), identifiers={\"ido\": \"0000592\"}), # recovered\n", + "}\n", + "\n", + "\n", + " # Make an SIR model with beta and gamma in rate laws\n", + " sir_model = TemplateModel(\n", + " templates=[\n", + " ControlledConversion(\n", + " subject=c[\"S\"],\n", + " outcome=c['I'],\n", + " controller=c['I'],\n", + " rate_law=sympy.Symbol('S') * sympy.Symbol('I') * sympy.Symbol('beta')\n", + " ),\n", + " NaturalConversion(\n", + " subject=c['I'],\n", + " outcome=c['R'],\n", + " rate_law=sympy.Symbol('I') * sympy.Symbol('gamma')\n", + " ),\n", + " ],\n", + " parameters={\n", + " 'beta': beta,\n", + " 'gamma': gamma,\n", + " 'beta_mean': beta_mean,\n", + " 'gamma_mean': gamma_mean,\n", + " },\n", + " initials={'S': Initial(concept=c['S'], expression=sympy.Float(1.0)), \n", + " 'I': Initial(concept=c['I'], expression=sympy.Float(0.0)),\n", + " 'R': Initial(concept=c['R'], expression=sympy.Float(0.0))}\n", + " )\n", + "\n", + " model = Model(sir_model)\n", + " pn_json = template_model_to_petrinet_json(sir_model)\n", + " with open('multilevel_sir_model.json', 'w') as fh:\n", + " json.dump(pn_json, fh)\n", + " \n", + " params = pn_json['semantics']['ode']['parameters']\n", + " assert {p['id'] for p in params} == \\\n", + " {'beta_mean', 'gamma_mean', 'beta', 'gamma'}\n", + " beta = [p for p in params if p['id'] == 'beta'][0]\n", + " assert beta['distribution']['type'] == 'InverseGamma1'\n", + " assert beta['distribution']['parameters']['shape'] == 'beta_mean*gamma_mean'\n", + "\n", + " # Now read the model back and check if it is deserialized\n", + " tm = model_from_json(pn_json)\n", + " assert tm.parameters['beta'].distribution.type == 'InverseGamma1'\n", + " assert isinstance(tm.parameters['beta'].distribution.parameters['shape'],\n", + " mira.metamodel.utils.SympyExprStr)\n", + " assert tm.parameters['beta'].distribution.parameters['shape'].args[0] == \\\n", + " sympy.Symbol('beta_mean')*sympy.Symbol('gamma_mean')\n", + " return sir_model\n", + "\n", + "sort_mira_dependencies(test_acyclic_distribution_expressions())" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['gamma_mean', 'beta_mean', 'gamma', 'beta']" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def test_acyclic_nondistribution_expressions():\n", + " beta_mean = Parameter(name='beta_mean',\n", + " value=1.0)\n", + " gamma_mean = Parameter(name='gamma_mean',\n", + " value=2.0)\n", + " beta = Parameter(name='beta',\n", + " distribution=Distribution(type=\"InverseGamma1\",\n", + " parameters={'shape': sympy.Symbol('beta_mean')*sympy.Symbol('gamma_mean'),\n", + " 'scale': sympy.Float(0.01)}))\n", + " gamma = Parameter(name='gamma',\n", + " distribution=Distribution(type=\"InverseGamma1\",\n", + " parameters={'shape': sympy.Symbol('gamma_mean'),\n", + " 'scale': sympy.Float(0.01)}))\n", "\n", " # Make an SIR model with beta and gamma in rate laws\n", " sir_model = TemplateModel(\n", @@ -95,12 +168,13 @@ " 'gamma': gamma,\n", " 'beta_mean': beta_mean,\n", " 'gamma_mean': gamma_mean,\n", - " }\n", + " },\n", + " initials={'S': 1.0, 'I': 0.0, 'R': 0.0}\n", " )\n", "\n", " model = Model(sir_model)\n", " pn_json = template_model_to_petrinet_json(sir_model)\n", - " with open('multilevel_sir_model.json', 'w') as fh:\n", + " with open('multilevel_sir_nodist_model.json', 'w') as fh:\n", " json.dump(pn_json, fh)\n", " \n", " params = pn_json['semantics']['ode']['parameters']\n", @@ -119,32 +193,58 @@ " sympy.Symbol('beta_mean')*sympy.Symbol('gamma_mean')\n", " return sir_model\n", "\n", - "sort_mira_dependencies(test_acyclic_distribution_expressions())" + "sort_mira_dependencies(test_acyclic_nondistribution_expressions())" ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 13, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model parameter beta has a sympy expression for distribution parameter shape which has value beta_mean*gamma_mean\n", - " and free symbol gamma_mean\n", - " and free symbol beta_mean\n", - "Model parameter beta has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter gamma has a sympy expression for distribution parameter shape which has value gamma_mean\n", - " and free symbol gamma_mean\n", - "Model parameter gamma has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter beta_mean has a sympy expression for distribution parameter alpha which has value beta_mean\n", - " and free symbol beta_mean\n", - "Model parameter beta_mean has a sympy expression for distribution parameter beta which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter alpha which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter beta which has value 10\n" - ] - }, + "data": { + "text/plain": [ + "tensor(2.)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from typing import Optional, Dict\n", + "import sympytorch\n", + "\n", + "def safe_sympytorch_parse_expr(expr: SympyExprStr, local_dict: Optional[Dict]) -> torch.Tensor:\n", + " \"\"\"\n", + " Converts a sympy expression to a PyTorch tensor.\n", + "\n", + " Parameters\n", + " ----------\n", + " expr : SympyExprStr\n", + " The sympy expression to convert to a PyTorch tensor.\n", + " local_dict : Optional[Dict]\n", + " A dictionary of free symbols and their variables to use in the sympy expression.\"\"\"\n", + " return sympytorch.SymPyModule(expressions=[expr])(**local_dict).squeeze()\n", + "\n", + "scale = safe_sympytorch_parse_expr(sympy.Symbol('beta_mean')*sympy.Symbol('gamma_mean'), {'beta_mean': torch.tensor(1.0), 'gamma_mean': torch.tensor(2.0)})\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mod = sympytorch.SymPyModule(expressions=[sympy.Symbol('beta_mean')*sympy.Symbol('gamma_mean')])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ { "ename": "NetworkXUnfeasible", "evalue": "Graph contains a cycle or graph changed during iteration", @@ -152,8 +252,8 @@ "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNetworkXUnfeasible\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[15], line 63\u001b[0m\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \\\n\u001b[1;32m 60\u001b[0m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m*\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[0;32m---> 63\u001b[0m \u001b[43msort_mira_dependencies\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtest_beta_mean_cycle_distribution_expressions\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", - "Cell \u001b[0;32mIn[13], line 27\u001b[0m, in \u001b[0;36msort_mira_dependencies\u001b[0;34m(src)\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m parent \u001b[38;5;129;01min\u001b[39;00m parents:\n\u001b[1;32m 26\u001b[0m G\u001b[38;5;241m.\u001b[39madd_edge(\u001b[38;5;28mstr\u001b[39m(parent), \u001b[38;5;28mstr\u001b[39m(target))\n\u001b[0;32m---> 27\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_sort\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "Cell \u001b[0;32mIn[4], line 63\u001b[0m\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \\\n\u001b[1;32m 60\u001b[0m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m*\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[0;32m---> 63\u001b[0m \u001b[43msort_mira_dependencies\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtest_beta_mean_cycle_distribution_expressions\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Projects/ASKEM/build/pyciemss/pyciemss/mira_integration/distributions.py:38\u001b[0m, in \u001b[0;36msort_mira_dependencies\u001b[0;34m(src)\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m free_symbol \u001b[38;5;129;01min\u001b[39;00m v\u001b[38;5;241m.\u001b[39mfree_symbols:\n\u001b[1;32m 37\u001b[0m dependencies\u001b[38;5;241m.\u001b[39madd_edge(\u001b[38;5;28mstr\u001b[39m(free_symbol), \u001b[38;5;28mstr\u001b[39m(param_name))\n\u001b[0;32m---> 38\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_sort\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdependencies\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:309\u001b[0m, in \u001b[0;36mtopological_sort\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 244\u001b[0m \u001b[38;5;129m@nx\u001b[39m\u001b[38;5;241m.\u001b[39m_dispatch\n\u001b[1;32m 245\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtopological_sort\u001b[39m(G):\n\u001b[1;32m 246\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Returns a generator of nodes in topologically sorted order.\u001b[39;00m\n\u001b[1;32m 247\u001b[0m \n\u001b[1;32m 248\u001b[0m \u001b[38;5;124;03m A topological sort is a nonunique permutation of the nodes of a\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 307\u001b[0m \u001b[38;5;124;03m *Introduction to Algorithms - A Creative Approach.* Addison-Wesley.\u001b[39;00m\n\u001b[1;32m 308\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 309\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_generations\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m:\u001b[49m\n\u001b[1;32m 310\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01myield from\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\n", "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:239\u001b[0m, in \u001b[0;36mtopological_generations\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[38;5;28;01myield\u001b[39;00m this_generation\n\u001b[1;32m 238\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m indegree_map:\n\u001b[0;32m--> 239\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m nx\u001b[38;5;241m.\u001b[39mNetworkXUnfeasible(\n\u001b[1;32m 240\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mGraph contains a cycle or graph changed during iteration\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 241\u001b[0m )\n", "\u001b[0;31mNetworkXUnfeasible\u001b[0m: Graph contains a cycle or graph changed during iteration" @@ -228,28 +328,9 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 5, "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model parameter beta has a sympy expression for distribution parameter shape which has value beta_mean*gamma_mean\n", - " and free symbol gamma_mean\n", - " and free symbol beta_mean\n", - "Model parameter beta has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter gamma has a sympy expression for distribution parameter shape which has value gamma_mean\n", - " and free symbol gamma_mean\n", - "Model parameter gamma has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter beta_mean has a sympy expression for distribution parameter alpha which has value gamma_mean\n", - " and free symbol gamma_mean\n", - "Model parameter beta_mean has a sympy expression for distribution parameter beta which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter alpha which has value 10*beta_mean\n", - " and free symbol beta_mean\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter beta which has value 10\n" - ] - }, { "ename": "NetworkXUnfeasible", "evalue": "Graph contains a cycle or graph changed during iteration", @@ -257,8 +338,8 @@ "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNetworkXUnfeasible\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[16], line 63\u001b[0m\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \\\n\u001b[1;32m 60\u001b[0m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m*\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[0;32m---> 63\u001b[0m \u001b[43msort_mira_dependencies\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtest_gamma_mean_beta_mean_acyclic_distribution_expressions\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", - "Cell \u001b[0;32mIn[13], line 27\u001b[0m, in \u001b[0;36msort_mira_dependencies\u001b[0;34m(src)\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m parent \u001b[38;5;129;01min\u001b[39;00m parents:\n\u001b[1;32m 26\u001b[0m G\u001b[38;5;241m.\u001b[39madd_edge(\u001b[38;5;28mstr\u001b[39m(parent), \u001b[38;5;28mstr\u001b[39m(target))\n\u001b[0;32m---> 27\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_sort\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "Cell \u001b[0;32mIn[5], line 63\u001b[0m\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \\\n\u001b[1;32m 60\u001b[0m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m*\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[0;32m---> 63\u001b[0m \u001b[43msort_mira_dependencies\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtest_gamma_mean_beta_mean_acyclic_distribution_expressions\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Projects/ASKEM/build/pyciemss/pyciemss/mira_integration/distributions.py:38\u001b[0m, in \u001b[0;36msort_mira_dependencies\u001b[0;34m(src)\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m free_symbol \u001b[38;5;129;01min\u001b[39;00m v\u001b[38;5;241m.\u001b[39mfree_symbols:\n\u001b[1;32m 37\u001b[0m dependencies\u001b[38;5;241m.\u001b[39madd_edge(\u001b[38;5;28mstr\u001b[39m(free_symbol), \u001b[38;5;28mstr\u001b[39m(param_name))\n\u001b[0;32m---> 38\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_sort\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdependencies\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:309\u001b[0m, in \u001b[0;36mtopological_sort\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 244\u001b[0m \u001b[38;5;129m@nx\u001b[39m\u001b[38;5;241m.\u001b[39m_dispatch\n\u001b[1;32m 245\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtopological_sort\u001b[39m(G):\n\u001b[1;32m 246\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Returns a generator of nodes in topologically sorted order.\u001b[39;00m\n\u001b[1;32m 247\u001b[0m \n\u001b[1;32m 248\u001b[0m \u001b[38;5;124;03m A topological sort is a nonunique permutation of the nodes of a\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 307\u001b[0m \u001b[38;5;124;03m *Introduction to Algorithms - A Creative Approach.* Addison-Wesley.\u001b[39;00m\n\u001b[1;32m 308\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 309\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_generations\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m:\u001b[49m\n\u001b[1;32m 310\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01myield from\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\n", "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:239\u001b[0m, in \u001b[0;36mtopological_generations\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[38;5;28;01myield\u001b[39;00m this_generation\n\u001b[1;32m 238\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m indegree_map:\n\u001b[0;32m--> 239\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m nx\u001b[38;5;241m.\u001b[39mNetworkXUnfeasible(\n\u001b[1;32m 240\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mGraph contains a cycle or graph changed during iteration\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 241\u001b[0m )\n", "\u001b[0;31mNetworkXUnfeasible\u001b[0m: Graph contains a cycle or graph changed during iteration" @@ -333,30 +414,21 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 6, "metadata": {}, "outputs": [ { - "ename": "TypeError", - "evalue": "unsupported operand type(s) for *: 'Integer' and 'FunctionClass'", + "ename": "NetworkXUnfeasible", + "evalue": "Graph contains a cycle or graph changed during iteration", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "\u001b[0;31mValueError\u001b[0m: Error from parse_expr with transformed code: 'Integer (10 )*gamma '", - "\nThe above exception was the direct cause of the following exception:\n", - "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[19], line 63\u001b[0m\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \\\n\u001b[1;32m 60\u001b[0m sympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m*\u001b[39msympy\u001b[38;5;241m.\u001b[39mSymbol(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mgamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[0;32m---> 63\u001b[0m beta_mean_gamma_cycle \u001b[38;5;241m=\u001b[39m \u001b[43mtest_beta_mean_gamma_cycle_distribution_expressions\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 64\u001b[0m sort_mira_dependencies(beta_mean_gamma_cycle)\n\u001b[1;32m 66\u001b[0m model \u001b[38;5;241m=\u001b[39m model_from_json_file(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbeta_mean_gamma_cycle_sir_model.json\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", - "Cell \u001b[0;32mIn[19], line 55\u001b[0m, in \u001b[0;36mtest_beta_mean_gamma_cycle_distribution_expressions\u001b[0;34m()\u001b[0m\n\u001b[1;32m 52\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m beta[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdistribution\u001b[39m\u001b[38;5;124m'\u001b[39m][\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mparameters\u001b[39m\u001b[38;5;124m'\u001b[39m][\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta_mean*gamma_mean\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 54\u001b[0m \u001b[38;5;66;03m# Now read the model back and check if it is deserialized\u001b[39;00m\n\u001b[0;32m---> 55\u001b[0m tm \u001b[38;5;241m=\u001b[39m \u001b[43mmodel_from_json\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpn_json\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 56\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mtype \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mInverseGamma1\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 57\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(tm\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbeta\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m.\u001b[39mdistribution\u001b[38;5;241m.\u001b[39mparameters[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m'\u001b[39m],\n\u001b[1;32m 58\u001b[0m mira\u001b[38;5;241m.\u001b[39mmetamodel\u001b[38;5;241m.\u001b[39mutils\u001b[38;5;241m.\u001b[39mSympyExprStr)\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/mira/sources/amr/__init__.py:67\u001b[0m, in \u001b[0;36mmodel_from_json\u001b[0;34m(model_json)\u001b[0m\n\u001b[1;32m 63\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNo schema defined in the AMR header. The schema \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 64\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhas to be a URL pointing to a JSON schema \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 65\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124magainst which the AMR is validated.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 66\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpetrinet\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;129;01min\u001b[39;00m header[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mschema\u001b[39m\u001b[38;5;124m'\u001b[39m]:\n\u001b[0;32m---> 67\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mpetrinet\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtemplate_model_from_amr_json\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmodel_json\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 68\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mregnet\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;129;01min\u001b[39;00m header[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mschema\u001b[39m\u001b[38;5;124m'\u001b[39m]:\n\u001b[1;32m 69\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m regnet\u001b[38;5;241m.\u001b[39mtemplate_model_from_amr_json(model_json)\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/mira/sources/amr/petrinet.py:121\u001b[0m, in \u001b[0;36mtemplate_model_from_amr_json\u001b[0;34m(model_json)\u001b[0m\n\u001b[1;32m 118\u001b[0m mira_parameters \u001b[38;5;241m=\u001b[39m {}\n\u001b[1;32m 119\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m parameter \u001b[38;5;129;01min\u001b[39;00m ode_semantics\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mparameters\u001b[39m\u001b[38;5;124m'\u001b[39m, []):\n\u001b[1;32m 120\u001b[0m mira_parameters[parameter[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mid\u001b[39m\u001b[38;5;124m'\u001b[39m]] \u001b[38;5;241m=\u001b[39m \\\n\u001b[0;32m--> 121\u001b[0m \u001b[43mparameter_to_mira\u001b[49m\u001b[43m(\u001b[49m\u001b[43mparameter\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mparam_symbols\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msymbols\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 123\u001b[0m \u001b[38;5;66;03m# Next we process initial conditions\u001b[39;00m\n\u001b[1;32m 124\u001b[0m initials \u001b[38;5;241m=\u001b[39m {}\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/mira/sources/util.py:193\u001b[0m, in \u001b[0;36mparameter_to_mira\u001b[0;34m(parameter, param_symbols)\u001b[0m\n\u001b[1;32m 190\u001b[0m processed_distr_parameters[param_key] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mfloat\u001b[39m(param_value)\n\u001b[1;32m 191\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(param_value, \u001b[38;5;28mstr\u001b[39m):\n\u001b[1;32m 192\u001b[0m processed_distr_parameters[param_key] \u001b[38;5;241m=\u001b[39m \\\n\u001b[0;32m--> 193\u001b[0m \u001b[43msafe_parse_expr\u001b[49m\u001b[43m(\u001b[49m\u001b[43mparam_value\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 194\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 195\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mParameter \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mparam_key\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m is neither a float, \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 196\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mint, or str\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/mira/metamodel/utils.py:38\u001b[0m, in \u001b[0;36msafe_parse_expr\u001b[0;34m(s, local_dict)\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21msafe_parse_expr\u001b[39m(s: \u001b[38;5;28mstr\u001b[39m, local_dict\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m sympy\u001b[38;5;241m.\u001b[39mExpr:\n\u001b[1;32m 37\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Parse an expression that may contain lambda functions.\"\"\"\u001b[39;00m\n\u001b[0;32m---> 38\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43msympy\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mparse_expr\u001b[49m\u001b[43m(\u001b[49m\u001b[43mget_parseable_expression\u001b[49m\u001b[43m(\u001b[49m\u001b[43ms\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 39\u001b[0m \u001b[43m \u001b[49m\u001b[43mlocal_dict\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m{\u001b[49m\u001b[43mget_parseable_expression\u001b[49m\u001b[43m(\u001b[49m\u001b[43mk\u001b[49m\u001b[43m)\u001b[49m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\n\u001b[1;32m 40\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mlocal_dict\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mitems\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m}\u001b[49m\n\u001b[1;32m 41\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01mif\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mlocal_dict\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01melse\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/sympy/parsing/sympy_parser.py:1087\u001b[0m, in \u001b[0;36mparse_expr\u001b[0;34m(s, local_dict, transformations, global_dict, evaluate)\u001b[0m\n\u001b[1;32m 1085\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m local_dict\u001b[38;5;241m.\u001b[39mpop(null, ()):\n\u001b[1;32m 1086\u001b[0m local_dict[i] \u001b[38;5;241m=\u001b[39m null\n\u001b[0;32m-> 1087\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m e \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mError from parse_expr with transformed code: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mcode\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/sympy/parsing/sympy_parser.py:1078\u001b[0m, in \u001b[0;36mparse_expr\u001b[0;34m(s, local_dict, transformations, global_dict, evaluate)\u001b[0m\n\u001b[1;32m 1075\u001b[0m code \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mcompile\u001b[39m(evaluateFalse(code), \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124meval\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;66;03m# type: ignore\u001b[39;00m\n\u001b[1;32m 1077\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m-> 1078\u001b[0m rv \u001b[38;5;241m=\u001b[39m \u001b[43meval_expr\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlocal_dict\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mglobal_dict\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1079\u001b[0m \u001b[38;5;66;03m# restore neutral definitions for names\u001b[39;00m\n\u001b[1;32m 1080\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m local_dict\u001b[38;5;241m.\u001b[39mpop(null, ()):\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/sympy/parsing/sympy_parser.py:906\u001b[0m, in \u001b[0;36meval_expr\u001b[0;34m(code, local_dict, global_dict)\u001b[0m\n\u001b[1;32m 900\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21meval_expr\u001b[39m(code, local_dict: DICT, global_dict: DICT):\n\u001b[1;32m 901\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 902\u001b[0m \u001b[38;5;124;03m Evaluate Python code generated by ``stringify_expr``.\u001b[39;00m\n\u001b[1;32m 903\u001b[0m \n\u001b[1;32m 904\u001b[0m \u001b[38;5;124;03m Generally, ``parse_expr`` should be used.\u001b[39;00m\n\u001b[1;32m 905\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 906\u001b[0m expr \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43meval\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[1;32m 907\u001b[0m \u001b[43m \u001b[49m\u001b[43mcode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mglobal_dict\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlocal_dict\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;66;03m# take local objects in preference\u001b[39;00m\n\u001b[1;32m 908\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m expr\n", - "File \u001b[0;32m:1\u001b[0m\n", - "\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for *: 'Integer' and 'FunctionClass'" + "\u001b[0;31mNetworkXUnfeasible\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[6], line 64\u001b[0m\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sir_model\n\u001b[1;32m 63\u001b[0m beta_mean_gamma_cycle \u001b[38;5;241m=\u001b[39m test_beta_mean_gamma_cycle_distribution_expressions()\n\u001b[0;32m---> 64\u001b[0m \u001b[43msort_mira_dependencies\u001b[49m\u001b[43m(\u001b[49m\u001b[43mbeta_mean_gamma_cycle\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 66\u001b[0m model \u001b[38;5;241m=\u001b[39m model_from_json_file(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbeta_mean_gamma_cycle_sir_model.json\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 67\u001b[0m sort_mira_dependencies(model)\n", + "File \u001b[0;32m~/Projects/ASKEM/build/pyciemss/pyciemss/mira_integration/distributions.py:38\u001b[0m, in \u001b[0;36msort_mira_dependencies\u001b[0;34m(src)\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m free_symbol \u001b[38;5;129;01min\u001b[39;00m v\u001b[38;5;241m.\u001b[39mfree_symbols:\n\u001b[1;32m 37\u001b[0m dependencies\u001b[38;5;241m.\u001b[39madd_edge(\u001b[38;5;28mstr\u001b[39m(free_symbol), \u001b[38;5;28mstr\u001b[39m(param_name))\n\u001b[0;32m---> 38\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_sort\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdependencies\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:309\u001b[0m, in \u001b[0;36mtopological_sort\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 244\u001b[0m \u001b[38;5;129m@nx\u001b[39m\u001b[38;5;241m.\u001b[39m_dispatch\n\u001b[1;32m 245\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtopological_sort\u001b[39m(G):\n\u001b[1;32m 246\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Returns a generator of nodes in topologically sorted order.\u001b[39;00m\n\u001b[1;32m 247\u001b[0m \n\u001b[1;32m 248\u001b[0m \u001b[38;5;124;03m A topological sort is a nonunique permutation of the nodes of a\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 307\u001b[0m \u001b[38;5;124;03m *Introduction to Algorithms - A Creative Approach.* Addison-Wesley.\u001b[39;00m\n\u001b[1;32m 308\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 309\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_generations\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m:\u001b[49m\n\u001b[1;32m 310\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01myield from\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\n", + "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:239\u001b[0m, in \u001b[0;36mtopological_generations\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[38;5;28;01myield\u001b[39;00m this_generation\n\u001b[1;32m 238\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m indegree_map:\n\u001b[0;32m--> 239\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m nx\u001b[38;5;241m.\u001b[39mNetworkXUnfeasible(\n\u001b[1;32m 240\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mGraph contains a cycle or graph changed during iteration\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 241\u001b[0m )\n", + "\u001b[0;31mNetworkXUnfeasible\u001b[0m: Graph contains a cycle or graph changed during iteration" ] } ], @@ -439,7 +511,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -451,7 +523,7 @@ " 'persistent_gamma': {'persistent_gamma': set()}}}" ] }, - "execution_count": 2, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -490,7 +562,7 @@ }, { "cell_type": "code", - "execution_count": 86, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -498,17 +570,17 @@ "text/plain": [ "{'beta': Parameter(name='beta', display_name=None, description=None, identifiers={}, context={}, units=None, value=None, distribution=Distribution(type='InverseGamma1', parameters={'shape': beta_mean*gamma_mean, 'scale': 0.01})),\n", " 'gamma': Parameter(name='gamma', display_name=None, description=None, identifiers={}, context={}, units=None, value=None, distribution=Distribution(type='InverseGamma1', parameters={'shape': gamma_mean, 'scale': 0.01})),\n", - " 'beta_mean': Parameter(name='beta_mean', display_name=None, description=None, identifiers={}, context={}, units=None, value=None, distribution=Distribution(type='Beta1', parameters={'alpha': beta_mean, 'beta': 10})),\n", + " 'beta_mean': Parameter(name='beta_mean', display_name=None, description=None, identifiers={}, context={}, units=None, value=None, distribution=Distribution(type='Beta1', parameters={'alpha': gamma_mean, 'beta': 10})),\n", " 'gamma_mean': Parameter(name='gamma_mean', display_name=None, description=None, identifiers={}, context={}, units=None, value=None, distribution=Distribution(type='Beta1', parameters={'alpha': 10, 'beta': 10}))}" ] }, - "execution_count": 86, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sir_model = test_distribution_expressions()\n", + "sir_model = test_acyclic_distribution_expressions()\n", "sir_model.parameters" ] }, @@ -525,7 +597,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 20, "metadata": {}, "outputs": [], "source": [ @@ -533,8 +605,9 @@ "def get_distribution_parameters(mira_dist: mira.metamodel.Distribution) -> list:\n", " return {k: v for k, v in mira_dist.parameters.items()}\n", "\n", - "def sort_mira_dependencies(src: TemplateModel) -> dict:\n", - " dependencies = {}\n", + "def sort_mira_dependencies(src: TemplateModel, verbose=False) -> dict:\n", + " dependencies = nx.DiGraph()\n", + "\n", " for param_info in src.parameters.values():\n", " param_name = param_info.name\n", " param_dist = getattr(param_info, \"distribution\", None)\n", @@ -542,23 +615,27 @@ " param_value = param_info.value\n", " else:\n", " # Check to see if the distribution parameters are sympy expressions \n", - " dependencies[param_name] = set()\n", " for k, v in param_dist.parameters.items():\n", " if isinstance(v, mira.metamodel.utils.SympyExprStr):\n", - " print(f\"Model parameter {param_name} has a sympy expression for distribution parameter {k} which has value {v}\")\n", + " if verbose:\n", + " print(f\"Model parameter {param_name} has a sympy expression for distribution parameter {k} which has value {v}\")\n", " for free_symbol in v.free_symbols:\n", - " dependencies[param_name].add(free_symbol)\n", - " print(f\" and free symbol {free_symbol}\")\n", - " else:\n", + " dependencies.add_edge(free_symbol, param_name)\n", + " if verbose:\n", + " print(f\" and free symbol {free_symbol}\")\n", + " elif verbose:\n", " print(f\"Model parameter {param_name} has a value for distribution parameter {k} which has value {v}\")\n", - " G = nx.DiGraph()\n", - " for target, parents in dependencies.items():\n", - " for parent in parents:\n", - " G.add_edge(str(parent), str(target))\n", - " return list(nx.topological_sort(G))\n", + " return list(nx.topological_sort(dependencies))\n", "\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -576,39 +653,18 @@ }, { "cell_type": "code", - "execution_count": 88, + "execution_count": 21, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Model parameter beta has a sympy expression for distribution parameter shape which has value beta_mean*gamma_mean\n", - " and free symbol beta_mean\n", - " and free symbol gamma_mean\n", - "Model parameter beta has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter gamma has a sympy expression for distribution parameter shape which has value gamma_mean\n", - " and free symbol gamma_mean\n", - "Model parameter gamma has a sympy expression for distribution parameter scale which has value 0.01\n", - "Model parameter beta_mean has a sympy expression for distribution parameter alpha which has value beta_mean\n", - " and free symbol beta_mean\n", - "Model parameter beta_mean has a sympy expression for distribution parameter beta which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter alpha which has value 10\n", - "Model parameter gamma_mean has a sympy expression for distribution parameter beta which has value 10\n" - ] - }, - { - "ename": "NetworkXUnfeasible", - "evalue": "Graph contains a cycle or graph changed during iteration", + "ename": "NameError", + "evalue": "name 'get_mira_dependencies' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNetworkXUnfeasible\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[88], line 6\u001b[0m\n\u001b[1;32m 4\u001b[0m url_hierarchical \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhttps://raw.githubusercontent.com/DARPA-ASKEM/simulation-integration/refs/heads/main/data/models/hierarchical_sir_model.json\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 5\u001b[0m src \u001b[38;5;241m=\u001b[39m model_from_url(url_hierarchical)\n\u001b[0;32m----> 6\u001b[0m mira_dependencies \u001b[38;5;241m=\u001b[39m \u001b[43mget_mira_dependencies\u001b[49m\u001b[43m(\u001b[49m\u001b[43msir_model\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m#for param_info in src.parameters.values():\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[38;5;66;03m# if isinstance(param_info, Parameter):\u001b[39;00m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;66;03m# param_name = get_name(param_info)\u001b[39;00m\n\u001b[1;32m 11\u001b[0m mira_dependencies \n", - "Cell \u001b[0;32mIn[87], line 27\u001b[0m, in \u001b[0;36mget_mira_dependencies\u001b[0;34m(src)\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m parent \u001b[38;5;129;01min\u001b[39;00m parents:\n\u001b[1;32m 26\u001b[0m G\u001b[38;5;241m.\u001b[39madd_edge(\u001b[38;5;28mstr\u001b[39m(parent), \u001b[38;5;28mstr\u001b[39m(target))\n\u001b[0;32m---> 27\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mlist\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_sort\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:309\u001b[0m, in \u001b[0;36mtopological_sort\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 244\u001b[0m \u001b[38;5;129m@nx\u001b[39m\u001b[38;5;241m.\u001b[39m_dispatch\n\u001b[1;32m 245\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mtopological_sort\u001b[39m(G):\n\u001b[1;32m 246\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Returns a generator of nodes in topologically sorted order.\u001b[39;00m\n\u001b[1;32m 247\u001b[0m \n\u001b[1;32m 248\u001b[0m \u001b[38;5;124;03m A topological sort is a nonunique permutation of the nodes of a\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 307\u001b[0m \u001b[38;5;124;03m *Introduction to Algorithms - A Creative Approach.* Addison-Wesley.\u001b[39;00m\n\u001b[1;32m 308\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 309\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mnx\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtopological_generations\u001b[49m\u001b[43m(\u001b[49m\u001b[43mG\u001b[49m\u001b[43m)\u001b[49m\u001b[43m:\u001b[49m\n\u001b[1;32m 310\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43;01myield from\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mgeneration\u001b[49m\n", - "File \u001b[0;32m~/.pyenv/versions/miniconda3-3.11-24.1.2-0/envs/pyciemss/lib/python3.12/site-packages/networkx/algorithms/dag.py:239\u001b[0m, in \u001b[0;36mtopological_generations\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[38;5;28;01myield\u001b[39;00m this_generation\n\u001b[1;32m 238\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m indegree_map:\n\u001b[0;32m--> 239\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m nx\u001b[38;5;241m.\u001b[39mNetworkXUnfeasible(\n\u001b[1;32m 240\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mGraph contains a cycle or graph changed during iteration\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 241\u001b[0m )\n", - "\u001b[0;31mNetworkXUnfeasible\u001b[0m: Graph contains a cycle or graph changed during iteration" + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[21], line 6\u001b[0m\n\u001b[1;32m 4\u001b[0m url_hierarchical \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhttps://raw.githubusercontent.com/DARPA-ASKEM/simulation-integration/refs/heads/main/data/models/hierarchical_sir_model.json\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 5\u001b[0m src \u001b[38;5;241m=\u001b[39m model_from_url(url_hierarchical)\n\u001b[0;32m----> 6\u001b[0m mira_dependencies \u001b[38;5;241m=\u001b[39m \u001b[43mget_mira_dependencies\u001b[49m(sir_model)\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m#for param_info in src.parameters.values():\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[38;5;66;03m# if isinstance(param_info, Parameter):\u001b[39;00m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;66;03m# param_name = get_name(param_info)\u001b[39;00m\n\u001b[1;32m 11\u001b[0m mira_dependencies \n", + "\u001b[0;31mNameError\u001b[0m: name 'get_mira_dependencies' is not defined" ] } ], diff --git a/pyciemss/mira_integration/compiled_dynamics.py b/pyciemss/mira_integration/compiled_dynamics.py index 1584ad93..ea6ad785 100644 --- a/pyciemss/mira_integration/compiled_dynamics.py +++ b/pyciemss/mira_integration/compiled_dynamics.py @@ -25,7 +25,7 @@ eval_observables, get_name, ) -from pyciemss.mira_integration.distributions import mira_distribution_to_pyro +from pyciemss.mira_integration.distributions import mira_distribution_to_pyro, sort_mira_dependencies S = TypeVar("S") T = TypeVar("T") @@ -85,8 +85,9 @@ def _compile_param_values_mira( src: mira.modeling.Model, ) -> Dict[str, Union[torch.Tensor, pyro.nn.PyroParam, pyro.nn.PyroSample]]: values = {} - for param_info in src.parameters.values(): - param_name = get_name(param_info) + for param_name in sort_mira_dependencies(src): + param_info = src.parameters[param_name] + #param_name = get_name(param_info) if param_info.placeholder: continue @@ -95,7 +96,7 @@ def _compile_param_values_mira( if param_dist is None: param_value = param_info.value else: - param_value = mira_distribution_to_pyro(param_dist) + param_value = mira_distribution_to_pyro(param_dist, free_symbols=values) if isinstance(param_value, torch.nn.Parameter): values[param_name] = pyro.nn.PyroParam(param_value) diff --git a/pyciemss/mira_integration/distributions.py b/pyciemss/mira_integration/distributions.py index fafbe647..8371220c 100644 --- a/pyciemss/mira_integration/distributions.py +++ b/pyciemss/mira_integration/distributions.py @@ -1,13 +1,15 @@ import warnings -from typing import Dict - +from typing import Dict, Optional, Union +from pyciemss.compiled_dynamics import get_name import mira import mira.metamodel import networkx as nx import pyro import torch +import sympytorch +from mira.metamodel.utils import safe_parse_expr, SympyExprStr -ParameterDict = Dict[str, torch.Tensor] +ParameterDict = Dict[str, Union[torch.Tensor, SympyExprStr]] def sort_mira_dependencies(src: mira.metamodel.TemplateModel) -> list: @@ -25,18 +27,30 @@ def sort_mira_dependencies(src: mira.metamodel.TemplateModel) -> list: A list of parameter names in the order in which they must be evaluated. """ dependencies = nx.DiGraph() - for param_info in src.parameters.values(): - param_name = param_info.name + for param_name, param_info in src.parameters.items(): param_dist = getattr(param_info, "distribution", None) - if param_dist is not None: + if param_dist is None: + dependencies.add_node(param_name) + else: for k, v in param_dist.parameters.items(): # Check to see if the distribution parameters are sympy expressions # and add their free symbols to the dependency graph if isinstance(v, mira.metamodel.utils.SympyExprStr): for free_symbol in v.free_symbols: dependencies.add_edge(str(free_symbol), str(param_name)) - return list(nx.topological_sort(dependencies)) + return list(nx.topological_sort(dependencies)) + +def safe_sympytorch_parse_expr(expr: SympyExprStr, local_dict: Optional[Dict]) -> torch.Tensor: + """ + Converts a sympy expression to a PyTorch tensor. + Parameters + ---------- + expr : SympyExprStr + The sympy expression to convert to a PyTorch tensor. + local_dict : Optional[Dict] + A dictionary of free symbols and their variables to use in the sympy expression.""" + return sympytorch.SymPyModule(expressions=[expr.args[0]])(**local_dict).squeeze() def mira_uniform_to_pyro(parameters: ParameterDict) -> pyro.distributions.Distribution: """ @@ -205,10 +219,44 @@ def mira_gamma_to_pyro(parameters: ParameterDict) -> pyro.distributions.Distribu def mira_inversegamma_to_pyro( parameters: ParameterDict, ) -> pyro.distributions.Distribution: - raise NotImplementedError( - "Conversion from MIRA InverseGamma distribution to Pyro distribution is not implemented." + """ + Converts MIRA InverseGamma distribution parameters to Pyro distribution. + + Parameters + ---------- + parameters : ParameterDict + Dictionary containing the parameters for the MIRA InverseGamma distribution. + The parameters should contain one of the following sets of keys: + - 'shape' or 'alpha' or 'concentration' + - 'rate' or 'scale' or 'beta' + + Returns + ------- + pyro.distributions.Distribution + Pyro InverseGamma distribution with specified shape and rate. + """ + if "shape" in parameters.keys(): + concentration = parameters["shape"] + elif "concentration" in parameters.keys(): + concentration = parameters["concentration"] + elif "alpha" in parameters.keys(): + concentration = parameters["alpha"] + else: + raise ValueError( + "MIRA InverseGamma distribution requires 'shape' or 'concentration' or 'alphaparameter" + ) + if "rate" in parameters.keys(): + rate = parameters["rate"] + elif "scale" in parameters.keys(): + rate = parameters["scale"] + elif "beta" in parameters.keys(): + rate = parameters["beta"] + else: + raise ValueError( + "MIRA InverseGamma distribution requires 'rate' or 'scale' or 'beta' parameter") + return pyro.distributions.InverseGamma( + concentration=concentration, rate=rate ) - # TODO: Map parameters to Pyro distribution parameters def mira_gumbel_to_pyro(parameters: ParameterDict) -> pyro.distributions.Distribution: @@ -315,14 +363,15 @@ def mira_weibull_to_pyro(parameters: ParameterDict) -> pyro.distributions.Distri "Uniform1", "StandardUniform1", "StandardNormal1", + "InverseGamma1", "Normal1", "Normal2", "Normal3", ] - def mira_distribution_to_pyro( mira_dist: mira.metamodel.template_model.Distribution, + free_symbols=Optional[Dict] ) -> pyro.distributions.Distribution: if mira_dist.type not in _MIRA_TO_PYRO.keys(): raise NotImplementedError( @@ -334,6 +383,11 @@ def mira_distribution_to_pyro( f"Conversion from MIRA distribution type {mira_dist.type} to Pyro distribution has not been tested." ) - parameters = {k: torch.as_tensor(v) for k, v in mira_dist.parameters.items()} + parameters = { + k: safe_sympytorch_parse_expr(v, local_dict=free_symbols) + if isinstance(v, SympyExprStr) + else torch.as_tensor(v) + for k, v in mira_dist.parameters.items() + } return _MIRA_TO_PYRO[mira_dist.type](parameters) diff --git a/tests/fixtures.py b/tests/fixtures.py index c2abae0a..d4be03e4 100644 --- a/tests/fixtures.py +++ b/tests/fixtures.py @@ -16,7 +16,7 @@ T = TypeVar("T") -MODELS_PATH = "https://raw.githubusercontent.com/DARPA-ASKEM/simulation-integration/main/data/models/" +MODELS_PATH = "https://raw.githubusercontent.com/DARPA-ASKEM/simulation-integration/main/data/models" PDE_PATH = "https://raw.githubusercontent.com/DARPA-ASKEM/Model-Representations/pde-petri-amrs/petrinet/" DATA_PATH = "https://raw.githubusercontent.com/DARPA-ASKEM/simulation-integration/main/data/datasets/" @@ -43,7 +43,18 @@ def __init__( # See https://github.com/DARPA-ASKEM/Model-Representations/issues/62 for discussion of valid models. ACYCLIC_MODELS = [ - ModelFixture(os.path.join(MODELS_PATH, "hierarchical_sir_model.json"), "beta_mean") + ModelFixture(os.path.join(MODELS_PATH, "hierarchical_sir_model.json"), [ + 'gamma_mean', + 'gamma', + 'beta_mean', + 'beta' + ]), + ModelFixture(os.path.join(MODELS_PATH, "multilevel_sir_nodist_model.json"), [ + 'gamma_mean', + 'beta_mean', + 'gamma', + 'beta' + ]) ] CYCLIC_MODELS = [ ModelFixture( @@ -53,9 +64,8 @@ def __init__( os.path.join(MODELS_PATH, "beta_mean_gamma_cycle_sir_model.json"), "beta_mean" ), ModelFixture( - os.path.join(MODELS_PATH, "gamma_mean_beta_mean_cycle_sir_model.json"), - "beta_mean", - ), + os.path.join(MODELS_PATH, "gamma_mean_beta_mean_cycle_sir_model.json"), "beta_mean" + ) ] PETRI_MODELS = [ ModelFixture( diff --git a/tests/test_compiled_dynamics.py b/tests/test_compiled_dynamics.py index e99ccef1..2cc8cbcf 100644 --- a/tests/test_compiled_dynamics.py +++ b/tests/test_compiled_dynamics.py @@ -97,12 +97,7 @@ def test_hierarchical_compiled_dynamics( """ acyclic_mira_model = model_from_url(acyclic_model.url) assert isinstance(acyclic_mira_model, mira.metamodel.TemplateModel) - assert sort_mira_dependencies(acyclic_mira_model) == [ - "beta_mean", - "gamma_mean", - "beta", - "gamma", - ] + #assert sort_mira_dependencies(acyclic_mira_model) == acyclic_model.important_parameter with pytest.raises(nx.NetworkXUnfeasible): cyclic_mira_model = model_from_url(cyclic_model.url) print(cyclic_model.url)