diff --git a/docs/tables/functions.tab b/docs/tables/functions.tab index 915b1d70..7a99a8dc 100644 --- a/docs/tables/functions.tab +++ b/docs/tables/functions.tab @@ -46,3 +46,7 @@ ALLOCATE AVAILABLE "ALLOCATE AVAILABLE(request, pp, avail)" "AllocateAvailable ALLOCATE BY PRIORITY "ALLOCATE BY PRIORITY(request, priority, size, width, supply)" "AllocateByPriorityStructure(request, priority, size, width, supply)" allocate_by_priority(request, priority, width, supply) INITIAL INITIAL(value) init init(value) InitialStructure(value) pysd.statefuls.Initial SAMPLE IF TRUE "SAMPLE IF TRUE(condition, input, initial_value)" "SampleIfTrueStructure(condition, input, initial_value)" pysd.statefuls.SampleIfTrue(...) +RANDOM 0 1 "RANDOM 0 1()" "CallStructure('random_0_1', ())" np.random.uniform(0, 1, size=final_shape) +RANDOM UNIFORM "RANDOM UNIFORM(m, x, s)" "CallStructure('random_uniform', (m, x, s))" np.random.uniform(m, x, size=final_shape) +RANDOM NORMAL "RANDOM NORMAL(m, x, h, r, s)" "CallStructure('random_normal', (m, x, h, r, s))" stats.truncnorm.rvs((m-h)/r, (x-h)/r, loc=h, scale=r, size=final_shape) +RANDOM EXPONENTIAL "RANDOM EXPONENTIAL(m, x, h, r, s)" "CallStructure('random_exponential', (m, x, h, r, s))" stats.truncexpon.rvs((x-np.maximum(m, h))/r, loc=np.maximum(m, h), scale=r, size=final_shape) diff --git a/docs/whats_new.rst b/docs/whats_new.rst index f4e9c675..cd08b503 100644 --- a/docs/whats_new.rst +++ b/docs/whats_new.rst @@ -1,5 +1,33 @@ What's New ========== +v3.14.0 (2024/04/24) +-------------------- +New Features +~~~~~~~~~~~~ +- Support Vensim's `RANDOM EXPONENTIAL `_ function (:issue:`107`). (`@enekomartinmartinez `_) + +Breaking changes +~~~~~~~~~~~~~~~~ + +Deprecations +~~~~~~~~~~~~ + +Bug fixes +~~~~~~~~~ +- Fix truncation in Vensim's `RANDOM NORMAL `_ function translation. (`@enekomartinmartinez `_) + +Documentation +~~~~~~~~~~~~~ +- Add supported random functions to the documentation tables. (`@enekomartinmartinez `_) + +Performance +~~~~~~~~~~~ + +Internal Changes +~~~~~~~~~~~~~~~~ +- Add test for random functions including comparison with Vensim outputs and expected values (:issue:`107`). (`@enekomartinmartinez `_) +- Allow to add multiple imports by the python function call builder. (`@enekomartinmartinez `_) + v3.13.4 (2024/02/29) -------------------- New Features diff --git a/pysd/_version.py b/pysd/_version.py index b15beac2..382c7be0 100644 --- a/pysd/_version.py +++ b/pysd/_version.py @@ -1 +1 @@ -__version__ = "3.13.4" +__version__ = "3.14.0" diff --git a/pysd/builders/python/python_expressions_builder.py b/pysd/builders/python/python_expressions_builder.py index ac7368d4..b207d3b9 100644 --- a/pysd/builders/python/python_expressions_builder.py +++ b/pysd/builders/python/python_expressions_builder.py @@ -606,9 +606,9 @@ def build_function_call(self, arguments: dict) -> BuildAST: """ # Get the function expression from the functionspace expression, modules = functionspace[self.function] - if modules: + for module in modules: # Update module dependencies in imports - self.section.imports.add(*modules) + self.section.imports.add(*module) calls = self.join_calls(arguments) diff --git a/pysd/builders/python/python_functions.py b/pysd/builders/python/python_functions.py index 0635b9fb..f42a91e4 100644 --- a/pysd/builders/python/python_functions.py +++ b/pysd/builders/python/python_functions.py @@ -2,118 +2,120 @@ # functions that can be diretcly applied over an array functionspace = { # directly build functions without dependencies - "elmcount": ("len(%(0)s)", None), + "elmcount": ("len(%(0)s)", ()), # directly build numpy based functions - "pi": ("np.pi", ("numpy",)), - "abs": ("np.abs(%(0)s)", ("numpy",)), - "power": ("np.power(%(0)s,%(1)s)", ("numpy",)), - "min": ("np.minimum(%(0)s, %(1)s)", ("numpy",)), - "max": ("np.maximum(%(0)s, %(1)s)", ("numpy",)), - "exp": ("np.exp(%(0)s)", ("numpy",)), - "sin": ("np.sin(%(0)s)", ("numpy",)), - "cos": ("np.cos(%(0)s)", ("numpy",)), - "tan": ("np.tan(%(0)s)", ("numpy",)), - "arcsin": ("np.arcsin(%(0)s)", ("numpy",)), - "arccos": ("np.arccos(%(0)s)", ("numpy",)), - "arctan": ("np.arctan(%(0)s)", ("numpy",)), - "sinh": ("np.sinh(%(0)s)", ("numpy",)), - "cosh": ("np.cosh(%(0)s)", ("numpy",)), - "tanh": ("np.tanh(%(0)s)", ("numpy",)), - "sqrt": ("np.sqrt(%(0)s)", ("numpy",)), - "ln": ("np.log(%(0)s)", ("numpy",)), - "log": ("(np.log(%(0)s)/np.log(%(1)s))", ("numpy",)), - # NUMPY: "invert_matrix": ("np.linalg.inv(%(0)s)", ("numpy",)), + "pi": ("np.pi", (("numpy",),)), + "abs": ("np.abs(%(0)s)", (("numpy",),)), + "power": ("np.power(%(0)s,%(1)s)", (("numpy",),)), + "min": ("np.minimum(%(0)s, %(1)s)", (("numpy",),)), + "max": ("np.maximum(%(0)s, %(1)s)", (("numpy",),)), + "exp": ("np.exp(%(0)s)", (("numpy",),)), + "sin": ("np.sin(%(0)s)", (("numpy",),)), + "cos": ("np.cos(%(0)s)", (("numpy",),)), + "tan": ("np.tan(%(0)s)", (("numpy",),)), + "arcsin": ("np.arcsin(%(0)s)", (("numpy",),)), + "arccos": ("np.arccos(%(0)s)", (("numpy",),)), + "arctan": ("np.arctan(%(0)s)", (("numpy",),)), + "sinh": ("np.sinh(%(0)s)", (("numpy",),)), + "cosh": ("np.cosh(%(0)s)", (("numpy",),)), + "tanh": ("np.tanh(%(0)s)", (("numpy",),)), + "sqrt": ("np.sqrt(%(0)s)", (("numpy",),)), + "ln": ("np.log(%(0)s)", (("numpy",),)), + "log": ("(np.log(%(0)s)/np.log(%(1)s))", (("numpy",),)), + # NUMPY: "invert_matrix": ("np.linalg.inv(%(0)s)", (("numpy",),)), # vector functions with axis to apply over # NUMPY: - # "prod": "np.prod(%(0)s, axis=%(axis)s)", ("numpy",)), - # "sum": "np.sum(%(0)s, axis=%(axis)s)", ("numpy",)), - # "vmax": "np.max(%(0)s, axis=%(axis)s)", ("numpy", )), - # "vmin": "np.min(%(0)s, axis=%(axis)s)", ("numpy",)) - "prod": ("prod(%(0)s, dim=%(axis)s)", ("functions", "prod")), - "sum": ("sum(%(0)s, dim=%(axis)s)", ("functions", "sum")), - "vmax": ("vmax(%(0)s, dim=%(axis)s)", ("functions", "vmax")), - "vmin": ("vmin(%(0)s, dim=%(axis)s)", ("functions", "vmin")), - "vmax_xmile": ("vmax(%(0)s)", ("functions", "vmax")), - "vmin_xmile": ("vmin(%(0)s)", ("functions", "vmin")), + # "prod": "np.prod(%(0)s, axis=%(axis)s)", (("numpy",),)), + # "sum": "np.sum(%(0)s, axis=%(axis)s)", (("numpy",),)), + # "vmax": "np.max(%(0)s, axis=%(axis)s)", ("numpy",),)), + # "vmin": "np.min(%(0)s, axis=%(axis)s)", (("numpy",),)) + "prod": ("prod(%(0)s, dim=%(axis)s)", (("functions", "prod"),)), + "sum": ("sum(%(0)s, dim=%(axis)s)", (("functions", "sum"),)), + "vmax": ("vmax(%(0)s, dim=%(axis)s)", (("functions", "vmax"),)), + "vmin": ("vmin(%(0)s, dim=%(axis)s)", (("functions", "vmin"),)), + "vmax_xmile": ("vmax(%(0)s)", (("functions", "vmax"),)), + "vmin_xmile": ("vmin(%(0)s)", (("functions", "vmin"),)), "vector_select": ( "vector_select(%(0)s, %(1)s, %(axis)s, %(2)s, %(3)s, %(4)s)", - ("functions", "vector_select") + (("functions", "vector_select"),) ), # functions defined in pysd.py_bakcend.functions "active_initial": ( "active_initial(__data[\"time\"].stage, lambda: %(0)s, %(1)s)", - ("functions", "active_initial")), + (("functions", "active_initial"),)), "if_then_else": ( "if_then_else(%(0)s, lambda: %(1)s, lambda: %(2)s)", - ("functions", "if_then_else")), + (("functions", "if_then_else"),)), "integer": ( "integer(%(0)s)", - ("functions", "integer")), + (("functions", "integer"),)), "invert_matrix": ( # NUMPY: remove "invert_matrix(%(0)s)", - ("functions", "invert_matrix")), # NUMPY: remove + (("functions", "invert_matrix"),)), # NUMPY: remove "modulo": ( "modulo(%(0)s, %(1)s)", - ("functions", "modulo")), + (("functions", "modulo"),)), "pulse": ( "pulse(__data['time'], %(0)s, width=%(1)s)", - ("functions", "pulse")), + (("functions", "pulse"),)), "Xpulse": ( "pulse(__data['time'], %(0)s, magnitude=%(1)s)", - ("functions", "pulse")), + (("functions", "pulse"),)), "pulse_train": ( "pulse(__data['time'], %(0)s, repeat_time=%(1)s, width=%(2)s, "\ "end=%(3)s)", - ("functions", "pulse")), + (("functions", "pulse"),)), "Xpulse_train": ( "pulse(__data['time'], %(0)s, repeat_time=%(1)s, magnitude=%(2)s)", - ("functions", "pulse")), + (("functions", "pulse"),)), "get_time_value": ( "get_time_value(__data['time'], %(0)s, %(1)s, %(2)s)", - ("functions", "get_time_value")), + (("functions", "get_time_value"),)), "quantum": ( "quantum(%(0)s, %(1)s)", - ("functions", "quantum")), + (("functions", "quantum"),)), "Xramp": ( "ramp(__data['time'], %(0)s, %(1)s)", - ("functions", "ramp")), + (("functions", "ramp"),)), "ramp": ( "ramp(__data['time'], %(0)s, %(1)s, %(2)s)", - ("functions", "ramp")), + (("functions", "ramp"),)), "step": ( "step(__data['time'], %(0)s, %(1)s)", - ("functions", "step")), + (("functions", "step"),)), "xidz": ( "xidz(%(0)s, %(1)s, %(2)s)", - ("functions", "xidz")), + (("functions", "xidz"),)), "zidz": ( "zidz(%(0)s, %(1)s)", - ("functions", "zidz")), + (("functions", "zidz"),)), "vector_sort_order": ( "vector_sort_order(%(0)s, %(1)s)", - ("functions", "vector_sort_order")), + (("functions", "vector_sort_order"),)), "vector_reorder": ( "vector_reorder(%(0)s, %(1)s)", - ("functions", "vector_reorder")), + (("functions", "vector_reorder"),)), "vector_rank": ( "vector_rank(%(0)s, %(1)s)", - ("functions", "vector_rank")), + (("functions", "vector_rank"),)), # random functions must have the shape of the component subscripts # most of them are shifted, scaled and truncated - # TODO: it is difficult to find same parametrization in Python, - # maybe build a new model "random_0_1": ( "np.random.uniform(0, 1, size=%(size)s)", - ("numpy",)), + (("numpy",),)), "random_uniform": ( "np.random.uniform(%(0)s, %(1)s, size=%(size)s)", - ("numpy",)), + (("numpy",),)), "random_normal": ( - "stats.truncnorm.rvs(%(0)s, %(1)s, loc=%(2)s, scale=%(3)s," - " size=%(size)s)", - ("scipy", "stats")), + "stats.truncnorm.rvs((%(0)s-%(2)s)/%(3)s, (%(1)s-%(2)s)/%(3)s," + " loc=%(2)s, scale=%(3)s, size=%(size)s)", + (("scipy", "stats"),)), + "random_exponential": ( + "stats.truncexpon.rvs((%(1)s-np.maximum(%(0)s, %(2)s))/%(3)s," + " loc=np.maximum(%(0)s, %(2)s), scale=%(3)s, size=%(size)s)", + (("scipy", "stats"), ("numpy",),)), } diff --git a/tests/conftest.py b/tests/conftest.py index 26542a13..a5f04fc2 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,13 +1,18 @@ import shutil from pathlib import Path +from dataclasses import dataclass + import pytest + from pysd import read_vensim, read_xmile, load from pysd.translators.vensim.vensim_utils import supported_extensions as\ vensim_extensions from pysd.translators.xmile.xmile_utils import supported_extensions as\ xmile_extensions +from pysd.builders.python.imports import ImportsManager + @pytest.fixture(scope="session") def _root(): @@ -21,6 +26,12 @@ def _test_models(_root): return _root.joinpath("test-models/tests") +@pytest.fixture(scope="session") +def _test_random(_root): + # test-models directory + return _root.joinpath("test-models/random") + + @pytest.fixture(scope="class") def shared_tmpdir(tmp_path_factory): # shared temporary directory for each class @@ -58,3 +69,38 @@ def ignore_warns(): "future version. Use timezone-aware objects to represent datetimes " "in UTC.*", ] + + +@pytest.fixture(scope="session") +def random_size(): + # size of generated random samples + return int(1e6) + + +@dataclass +class FakeComponent: + element: str + section: object + subscripts_dict: dict + + +@dataclass +class FakeSection: + namespace: object + macrospace: dict + imports: object + + +@dataclass +class FakeNamespace: + cleanspace: dict + + +@pytest.fixture(scope="function") +def fake_component(): + # fake_component used to translate random functions to python + return FakeComponent( + '', + FakeSection(FakeNamespace({}), {}, ImportsManager()), + {} + ) diff --git a/tests/pytest_pysd/pytest_random.py b/tests/pytest_pysd/pytest_random.py index 37596405..7c83b4d2 100644 --- a/tests/pytest_pysd/pytest_random.py +++ b/tests/pytest_pysd/pytest_random.py @@ -1,11 +1,31 @@ - import pytest import shutil import numpy as np import xarray as xr +import pandas as pd +from scipy import stats import pysd +from pysd.translators.vensim.vensim_element import Component +from pysd.builders.python.python_expressions_builder import\ + CallBuilder, NumericBuilder + + +tolerance = { + 'data': { + 'mean': 1e-2, + 'variance': 1e-2, + 'skewness': 5e-2, + 'kurtosis': 10e-2 + }, + 'expc': { + 'mean': 5e-3, + 'variance': 5e-3, + 'skewness': 2.5e-2, + 'kurtosis': 5e-2 + } +} class TestRandomModel: @@ -56,3 +76,83 @@ def test_translate(self, model_path): values = out[var].values # assert all values are different in each dimension and time step assert len(np.unique(values)) == np.prod(values.shape) + + +class TestRandomVensim(): + @pytest.fixture(scope="function") + def data_raw(self, input_file, _test_random): + file = _test_random.joinpath(input_file) + data = pd.read_table(file, sep='\t', index_col=0) + return data + + @pytest.fixture(scope="function") + def data_python(self, data_raw, fake_component, random_size): + out = {} + for col in data_raw.columns: + component = Component('', ([], []), col) + component.parse() + builder = CallBuilder(component.ast, fake_component) + # TODO get minmax from here based on definition + args = { + i: NumericBuilder(arg, fake_component).build({}) + for i, arg in builder.arguments.items() + } + expr = builder.build(args).expression + expr = expr.replace('()', str(random_size)) + out[col] = eval(expr) + + return pd.DataFrame(out) + + @pytest.mark.parametrize( + "input_file", + [ + ( # uniform distribution + 'random_uniform/uniform_vensim.tab' + ), + ( # truncated normal distribution + 'random_normal/normal_vensim.tab' + ), + ( # truncated exponential distribution + 'random_exponential/exponential_vensim.tab' + ), + ], + ids=["uniform", "normal", "exponential"] + ) + def test_statistics_vensim(self, data_raw, data_python): + raw_desc = stats.describe(data_raw, axis=0, nan_policy='omit') + py_desc = stats.describe(data_python, axis=0, nan_policy='omit') + for stat in ('mean', 'variance', 'skewness', 'kurtosis'): + assert np.allclose( + getattr(py_desc, stat), + getattr(raw_desc, stat), + atol=tolerance['data'][stat], rtol=tolerance['data'][stat] + ), 'Failed when comparing %s:\n\t%s\n\t%s' % ( + stat, getattr(raw_desc, stat), getattr(py_desc, stat)) + + @pytest.mark.parametrize( + "input_file", + [ + ( # uniform distribution + 'random_uniform/uniform_expected.tab' + ), + ( # truncated normal distribution + 'random_normal/normal_expected.tab' + ), + ( # truncated exponential distribution + 'random_exponential/exponential_expected.tab' + ), + ], + ids=["uniform", "normal", "exponential"] + ) + def test_statistics_expected(self, data_raw, data_python): + py_desc = stats.describe(data_python, axis=0, nan_policy='omit') + for stat in ('mean', 'variance', 'skewness', 'kurtosis'): + assert np.allclose( + getattr(py_desc, stat), + data_raw.loc[stat].values, + atol=tolerance['expc'][stat], rtol=tolerance['expc'][stat] + ), 'Failed when comparing %s:\n\t%s\n\t%s' % ( + stat, data_raw.loc[stat], getattr(py_desc, stat)) + + assert np.all(data_raw.loc['min'] < py_desc.minmax[0]) + assert np.all(data_raw.loc['max'] > py_desc.minmax[1]) diff --git a/tests/requirements.txt b/tests/requirements.txt index ec944b8d..a8c583e3 100644 --- a/tests/requirements.txt +++ b/tests/requirements.txt @@ -7,7 +7,7 @@ pytest-xdist coverage coveralls psutil -netCDF4>=1.6.* +netCDF4>=1.6 dask[array] dask[diagnostics] dask[distributed] diff --git a/tests/test-models b/tests/test-models index a752dac6..6a3f4e7f 160000 --- a/tests/test-models +++ b/tests/test-models @@ -1 +1 @@ -Subproject commit a752dac620d1f84cb13b7aa8397ec55b5de623f8 +Subproject commit 6a3f4e7f7427e0a0433d101789ce2b73a56e0385