From a22979f5bf4edad94f6292c3c93f4f061f6911a2 Mon Sep 17 00:00:00 2001 From: Eneko Martin-Martinez <61285767+enekomartinmartinez@users.noreply.github.com> Date: Mon, 4 Jul 2022 16:03:56 +0200 Subject: [PATCH] Add support for ALLOCATE AVAILABLE Add support for Vensim's ALLOCATE AVAILABLE --- docs/conf.py | 3 +- docs/installation.rst | 1 + docs/tables/functions.tab | 1 + docs/whats_new.rst | 39 +- pysd/_version.py | 2 +- pysd/builders/python/namespace.py | 7 +- .../python/python_expressions_builder.py | 94 ++- pysd/py_backend/allocation.py | 534 +++++++++++++++++- pysd/tools/benchmarking.py | 6 +- .../structures/abstract_expressions.py | 27 +- pysd/translators/vensim/vensim_structures.py | 1 + requirements.txt | 1 + .../pytest_integration_vensim_pathway.py | 5 + tests/pytest_pysd/pytest_allocation.py | 304 +++++++++- tests/pytest_pysd/pytest_utils.py | 2 +- .../data/pytest_data_with_model.py | 1 - tests/test-models | 2 +- 17 files changed, 997 insertions(+), 33 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 6458aa9d..9cc3267a 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -30,7 +30,7 @@ 'numpy', 'scipy', 'matplotlib', 'matplotlib.pyplot', 'scipy.stats', 'scipy.integrate', 'pandas', 'parsimonious', 'parsimonious.nodes', 'xarray', 'autopep8', 'scipy.linalg', 'parsimonious.exceptions', - 'scipy.stats.distributions', 'progressbar', 'black' + 'scipy.stats.distributions', 'progressbar', 'black', 'scipy.optimize' ] for mod_name in MOCK_MODULES: @@ -152,6 +152,7 @@ 'pandas': ('https://pandas.pydata.org/docs/', None), 'xarray': ('https://docs.xarray.dev/en/stable/', None), 'numpy': ('https://numpy.org/doc/stable/', None), + 'scipy': ('https://docs.scipy.org/doc/scipy/', None), 'pytest': ('https://docs.pytest.org/en/7.1.x/', None) } diff --git a/docs/installation.rst b/docs/installation.rst index 2f9a73c8..f32a5027 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -53,6 +53,7 @@ PySD builds on the core Python data analytics stack, and the following third par * black * openpyxl * progressbar2 +* portion These modules should build automatically if you are installing via `pip`. If you are building from source, or if pip fails to load them, they can be loaded with the same `pip` syntax as diff --git a/docs/tables/functions.tab b/docs/tables/functions.tab index 3cd94775..627125f2 100644 --- a/docs/tables/functions.tab +++ b/docs/tables/functions.tab @@ -39,6 +39,7 @@ VECTOR RANK VECTOR RANK(vec, direction) "CallStructure('vector_rank', (vec, di VECTOR REORDER VECTOR REORDER(vec, svec) "CallStructure('vector_reorder', (vec, svec))" vector_reorder(vec, svec) VECTOR SORT ORDER VECTOR SORT ORDER(vec, direction) "CallStructure('vector_sort_order', (vec, direction))" vector_sort_order(vec, direction) GAME GAME(A) GameStructure(A) A +ALLOCATE AVAILABLE "ALLOCATE AVAILABLE(request, pp, avail)" "AllocateAvailableStructure(request, pp, avail)" allocate_available(request, pp, avail) Not all the priority profiles are included. 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(...) diff --git a/docs/whats_new.rst b/docs/whats_new.rst index ffc8bbb2..6c86db21 100644 --- a/docs/whats_new.rst +++ b/docs/whats_new.rst @@ -1,12 +1,39 @@ What's New ========== +v3.4.0 (2022/06/29) +------------------- + +New Features +~~~~~~~~~~~~ +- Add support for Vensim's `ALLOCATE AVAILABLE `_ (:py:func:`pysd.py_backend.allocation.allocate_available`) function (:issue:`339`). Integer allocation cases have not been implemented neither the fixed quantity and constant elasticity curve priority functions. + +Breaking changes +~~~~~~~~~~~~~~~~ + +Deprecations +~~~~~~~~~~~~ + +Bug fixes +~~~~~~~~~ + +Documentation +~~~~~~~~~~~~~ +- Improve the documentation of the :py:mod:`pysd.py_backend.allocation` module. + +Performance +~~~~~~~~~~~ + +Internal Changes +~~~~~~~~~~~~~~~~ +- Add a class to manage priority profiles so it can be also used by the `many-to-many allocation `_. + v3.3.0 (2022/06/22) ------------------- New Features ~~~~~~~~~~~~ -- Add support for Vensim's `ALLOCATE BY PRIORITY `_ (:func:`pysd.py_backend.allocation.allocate_by_priority`) function (:issue:`263`). +- Add support for Vensim's `ALLOCATE BY PRIORITY `_ (:py:func:`pysd.py_backend.allocation.allocate_by_priority`) function (:issue:`263`). Breaking changes ~~~~~~~~~~~~~~~~ @@ -33,8 +60,8 @@ v3.2.0 (2022/06/10) New Features ~~~~~~~~~~~~ -- Add support for Vensim's `GET TIME VALUE `_ (:func:`pysd.py_backend.functions.get_time_value`) function (:issue:`332`). Not all cases have been implemented. -- Add support for Vensim's `VECTOR SELECT `_ (:func:`pysd.py_backend.functions.vector_select`) function (:issue:`266`). +- Add support for Vensim's `GET TIME VALUE `_ (:py:func:`pysd.py_backend.functions.get_time_value`) function (:issue:`332`). Not all cases have been implemented. +- Add support for Vensim's `VECTOR SELECT `_ (:py:func:`pysd.py_backend.functions.vector_select`) function (:issue:`266`). Breaking changes ~~~~~~~~~~~~~~~~ @@ -61,9 +88,9 @@ v3.1.0 (2022/06/02) New Features ~~~~~~~~~~~~ -- Add support for Vensim's `VECTOR SORT ORDER `_ (:func:`pysd.py_backend.functions.vector_sort_order`) function (:issue:`326`). -- Add support for Vensim's `VECTOR RANK `_ (:func:`pysd.py_backend.functions.vector_rank`) function (:issue:`326`). -- Add support for Vensim's `VECTOR REORDER `_ (:func:`pysd.py_backend.functions.vector_reorder`) function (:issue:`326`). +- Add support for Vensim's `VECTOR SORT ORDER `_ (:py:func:`pysd.py_backend.functions.vector_sort_order`) function (:issue:`326`). +- Add support for Vensim's `VECTOR RANK `_ (:py:func:`pysd.py_backend.functions.vector_rank`) function (:issue:`326`). +- Add support for Vensim's `VECTOR REORDER `_ (:py:func:`pysd.py_backend.functions.vector_reorder`) function (:issue:`326`). Breaking changes ~~~~~~~~~~~~~~~~ diff --git a/pysd/_version.py b/pysd/_version.py index 88c513ea..903a158a 100644 --- a/pysd/_version.py +++ b/pysd/_version.py @@ -1 +1 @@ -__version__ = "3.3.0" +__version__ = "3.4.0" diff --git a/pysd/builders/python/namespace.py b/pysd/builders/python/namespace.py index 5bd98eb4..8246472a 100644 --- a/pysd/builders/python/namespace.py +++ b/pysd/builders/python/namespace.py @@ -6,9 +6,10 @@ # used to create Python safe names with the variable reserved_words from keyword import kwlist from builtins import __dir__ as bidir +from pysd.py_backend.allocation import __dir__ as adir +from pysd.py_backend.cache import __dir__ as cadir from pysd.py_backend.components import __dir__ as cdir from pysd.py_backend.data import __dir__ as ddir -from pysd.py_backend.cache import __dir__ as cadir from pysd.py_backend.external import __dir__ as edir from pysd.py_backend.functions import __dir__ as fdir from pysd.py_backend.statefuls import __dir__ as sdir @@ -29,8 +30,8 @@ class NamespaceManager: """ _reserved_words = set( - dir() + bidir() + cdir() + ddir() + cadir() + edir() + fdir() - + sdir() + udir()).union(kwlist) + dir() + adir() + bidir() + cadir() + cdir() + ddir() + edir() + + fdir() + sdir() + udir()).union(kwlist) def __init__(self, parameters: List[str] = []): self._used_words = self._reserved_words.copy() diff --git a/pysd/builders/python/python_expressions_builder.py b/pysd/builders/python/python_expressions_builder.py index 89e295fc..f77e6709 100644 --- a/pysd/builders/python/python_expressions_builder.py +++ b/pysd/builders/python/python_expressions_builder.py @@ -16,13 +16,14 @@ from pysd.py_backend.utils import compute_shape from pysd.translators.structures.abstract_expressions import\ - AbstractSyntax, AllocateByPriorityStructure, ArithmeticStructure,\ - CallStructure, DataStructure, DelayFixedStructure, DelayStructure,\ - DelayNStructure, ForecastStructure, GameStructure, GetConstantsStructure,\ - GetDataStructure, GetLookupsStructure, InitialStructure,\ - InlineLookupsStructure, IntegStructure, LogicStructure, LookupsStructure,\ - ReferenceStructure, SampleIfTrueStructure, SmoothNStructure,\ - SmoothStructure, SubscriptsReferenceStructure, TrendStructure + AbstractSyntax, AllocateAvailableStructure, AllocateByPriorityStructure,\ + ArithmeticStructure, CallStructure, DataStructure, DelayFixedStructure,\ + DelayStructure, DelayNStructure, ForecastStructure, GameStructure,\ + GetConstantsStructure, GetDataStructure, GetLookupsStructure,\ + InitialStructure, InlineLookupsStructure, IntegStructure,\ + LogicStructure, LookupsStructure, ReferenceStructure,\ + SampleIfTrueStructure, SmoothNStructure, SmoothStructure,\ + SubscriptsReferenceStructure, TrendStructure from .python_functions import functionspace from .subscripts import SubscriptManager @@ -750,6 +751,75 @@ def _compute_axis(self, subscripts: dict) -> tuple: return coords, axis +class AllocateAvailableBuilder(StructureBuilder): + """Builder for allocate_available function.""" + + def __init__(self, allocate_str: AllocateAvailableStructure, + component: object): + super().__init__(None, component) + + pp = allocate_str.pp + pp_sub = self.section.subscripts.elements[pp.reference][-1:] + pp.subscripts.subscripts = pp.subscripts.subscripts[:-1] + pp_sub + self.arguments = { + "request": allocate_str.request, + "pp": pp, + "avail": allocate_str.avail + } + + def build(self, arguments: dict) -> BuildAST: + """ + Build method. + + Parameters + ---------- + arguments: dict + The dictionary of builded arguments. + + Returns + ------- + built_ast: BuildAST + The built object. + + """ + self.section.imports.add("allocation", "allocate_available") + + calls = self.join_calls(arguments) + + # the last sub of the request must be keep last sub of request and + # priority + last_sub = list(arguments["request"].subscripts)[-1] + pp_sub = list(arguments["pp"].subscripts)[-1] + # compute the merged subscripts + final_subscripts = self.get_final_subscripts(arguments) + # remove last sub from request + last_sub_value = final_subscripts[last_sub] + pp_sub_value = final_subscripts[pp_sub] + del final_subscripts[last_sub], final_subscripts[pp_sub] + + # Update the susbcripts of avail + arguments["avail"].reshape( + self.section.subscripts, final_subscripts, True) + + # Include last sub of request in the last position and update + # the subscripts of request + final_subscripts[last_sub] = last_sub_value + arguments["request"].reshape( + self.section.subscripts, final_subscripts, True) + + # Include priority subscripts and update the subscripts of pp + final_subscripts[pp_sub] = pp_sub_value + arguments["pp"].reshape( + self.section.subscripts, final_subscripts, True) + + expression = "allocate_available(%(request)s, %(pp)s, %(avail)s)" + return BuildAST( + expression=expression % arguments, + calls=calls, + subscripts=arguments["request"].subscripts, + order=0) + + class AllocateByPriorityBuilder(StructureBuilder): """Builder for allocate_by_priority function.""" @@ -1731,6 +1801,7 @@ def subscripts(self, subscripts: SubscriptsReferenceStructure): # the reference has a subscripts which is it not # applied (!) and does not appear in the definition # of the variable + not_mapped = True for mapped in self.section.subscripts.mapping[dim]: # check the mapped subscripts # TODO update this and the parser to make it @@ -1741,7 +1812,15 @@ def subscripts(self, subscripts: SubscriptsReferenceStructure): # and it is not already in the variable self.mapping_subscripts[mapped] =\ self.section.subscripts.subscripts[mapped] + not_mapped = False break + if not_mapped: + # manage other not mapped subscripts + # this is necessary for Allocate Available + # where we must force the expression in the + # right to have more subscripts thant the + # expression in the left + self.mapping_subscripts[dim] = coordinates else: # the subscript is in the variable definition, # do not change it @@ -2089,6 +2168,7 @@ class ASTVisitor: ReferenceStructure: ReferenceBuilder, CallStructure: CallBuilder, GameStructure: GameBuilder, + AllocateAvailableStructure: AllocateAvailableBuilder, AllocateByPriorityStructure: AllocateByPriorityBuilder, LogicStructure: OperationBuilder, ArithmeticStructure: OperationBuilder, diff --git a/pysd/py_backend/allocation.py b/pysd/py_backend/allocation.py index 925d83be..61b66910 100644 --- a/pysd/py_backend/allocation.py +++ b/pysd/py_backend/allocation.py @@ -2,16 +2,540 @@ The provided allocation functions have no direct analog in the standard Python data analytics stack. They are provided in a structure that makes it easy for the model elements to call. The functions may be similar to -the original functions given by Vensim or Stella, but sometimes the -number or order of arguments may change. The allocation functions may -call a protected function or class method thatintegrates the algorithm -to compute the allocation. The algorithms are briefly explained in these +the original functions given by Vensim, but sometimes the number or +order of arguments may change. The allocation functions may call a +protected function or class method thatintegrates the algorithm to +compute the allocation. The algorithms are briefly explained in these functions docstring. + +Note +---- +The Allocation functions basis is explained in the Vensim documentation. +https://www.vensim.com/documentation/allocation_overview.html + +Warning +------- +Some allocation function's results may differ from the result given by +Vensim as optimization functions are used to solve the allocation +problems. Those algorithms may not work in the same way or may +have differences in the numerical error propagation. + """ import itertools +from math import erfc -import xarray as xr import numpy as np +import xarray as xr +from scipy.optimize import least_squares +import portion as p + + +class Priorities: + @classmethod + def get_functions(cls, q0, pp, kind): + """ + Get priority functions based on the demand/supply and priority profile. + + Parameters + ---------- + q0: numpy.array + values of maximum demand or supply of each component. + Its shape should be (n,) + pp: numpy.array + pp values array. Its shape should be (n, m). + kind: str ("demand" or "supply") + The kind of priority "demand" or "supply". + + Returns + ------- + functions: list of functions + List of allocation supply or demand function for each element. + + full_allocation: function + Full allocation function. It is the result function of + addying all the functions. + + def_intervals: list of tuples + List of (supply interval, priority interval, mean priority) + where the full_allocation function is extrictly monotonous + (injective). Givin a supply value, this is used to compute + the limits and starting point of the optimization problem. + + """ + if np.any(pp[:, 2] <= 0): + # pwidth values smaller than 0 + raise ValueError("pwidth values must be positive.") + + if kind == "demand": + # Get the list of priority functions and the intervals where + # they are strictly monotonous (injective function) + func_int = [ + cls.get_function_demand(q0[i], pp[i]) + for i in range(pp.shape[0]) + ] + # In order to get the range of the full_allocation function, + # we need to flip the lower and the upper value, as it is a + # decreasing function for demand + int_attr = {"lower": "upper", "upper": "lower"} + elif kind == "supply": # pragma: no cover + # Get the list of priority functions and the intervals where + # they are strictly monotonous (injective function) + func_int = [ + cls.get_function_supply(q0[i], pp[i]) + for i in range(pp.shape[0]) + ] + # In order to get the range of the full_allocation function, + # we need to keep the lower and the upper value, as it is a + # increasing function for supply + int_attr = {"lower": "lower", "upper": "upper"} + else: + raise ValueError( + f"kind='{kind}' is not allowed. kind should be " + "'demand' or 'supply'.") + + functions = [fi[0] for fi in func_int] + intervals = [fi[1] for fi in func_int] + + # Join the intervals of all functions to get the intervals where + # the sum of the functions is strictly monotonous (injective + # function), therefore we can solve the minimization problem in + # strictly monotonous areas of the function, avoiding the crash + # of the algorithm + interval = intervals[0] + for i in intervals[1:]: + interval = interval.union(i) + + # Full allocation function -> function to solve + def full_allocation(x): return np.sum([func(x) for func in functions]) + + def_intervals = [] + for subinterval in interval: + # Iterate over disjoint interval sections + # Each interval section will be converted in supply interval + # and compute the starting point for the supply interval + # as the midpoint in the priority interval + def_intervals.append(( + p.closed( + full_allocation(getattr(subinterval, int_attr["lower"])), + full_allocation(getattr(subinterval, int_attr["upper"])) + ), + subinterval, + .5*(subinterval.upper+subinterval.lower) + )) + + return functions, full_allocation, def_intervals + + @classmethod + def get_function_demand(cls, q0, pp): + """ + Get priority functions for demand based on the priority profile. + + Parameters + ---------- + q0: float [0, +np.inf) + The demand of the target. + pp: numpy.array + pp values array. + + Returns + ------- + priority_func: function + Priority function. + interval: portion.interval + The interval where the priority function is strictly monotonous. + + """ + if q0 == 0: + # No demand is requested return a 0 function with an empty interval + return lambda x: 0, p.empty() + if pp[0] == 0: + # Fixed quantity demand + return cls.fixed_quantity(q0, *pp[1:]) + elif pp[0] == 1: + # Rectangular demand + return cls.rectangular(q0, *pp[1:]) + elif pp[0] == 2: + # Triangular demand + return cls.triangular(q0, *pp[1:]) + elif pp[0] == 3: + # Normal distribution demand + return cls.normal(q0, *pp[1:]) + elif pp[0] == 4: + # Exponential distribution demand + return cls.exponential(q0, *pp[1:]) + elif pp[0] == 5: + # Constant elasticity demand + return cls.constant_elasticity_demand(q0, *pp[1:]) + else: + raise ValueError( + f"The priority function for pprofile={pp[0]} is not valid.") + + @classmethod + def get_function_supply(cls, q0, pp): + """ + Get priority functions for supply based on the priority profile. + + Parameters + ---------- + q0: float [0, +np.inf) + The supply of the producer. + pp: numpy.array + pp values array. + + Returns + ------- + priority_func: function + Priority function. + interval: portion.interval + The interval where the priority function is strictly monotonous. + + """ + # TODO: This function should be similar to the demand function + # it is neccessary for the many-to-many allocation given by + # the set FIND MARKET PLACE, DEMAND AT PRICE, SUPPLY AT PRICE + raise NotImplementedError("get_function_supply is not implemented.") + + @staticmethod + def fixed_quantity(q0, ppriority, pwidth, pextra): + raise NotImplementedError( + "fixed_quantity priority profile is not implemented.") + + @staticmethod + def rectangular(q0, ppriority, pwidth, pextra): + """ + Demand curve for rectangular shape. + The supply curve will be shaped as the integral of a rectangle. + + Parameters + ---------- + q0: float + The total demand/supply of the element. + ppriority: float + Specifies the midpoint of the curve. + pwidth: float + Determines the speed with which the curve goes from 0 to + the specified quantity. + pextra: float + Ignore. + + Returns + ------- + priority_func: function + The priority function. + + """ + def priority_func(x): + if x <= ppriority - pwidth*.5: + return q0 + elif x < ppriority + pwidth*.5: + return q0*(1-(x-ppriority+pwidth*.5)/pwidth) + else: + return 0 + + return ( + priority_func, + p.open(ppriority - pwidth*.5, ppriority + pwidth*.5) + ) + + @staticmethod + def triangular(q0, ppriority, pwidth, pextra): + """ + Demand curve for triangular shape. + The supply curve will be shaped as the integral of a triangle. + + Parameters + ---------- + q0: float + The total demand/supply of the element. + ppriority: float + Specifies the midpoint of the curve. + pwidth: float + Determines the speed with which the curve goes from 0 to the + specified quantity. + pextra: float + Ignore. + + Returns + ------- + priority_func: function + The priority function. + + """ + def priority_func(x): + if x <= ppriority - pwidth*.5: + return q0 + elif x < ppriority: + return q0*(1-2*(x-ppriority+pwidth*.5)**2/pwidth**2) + elif x < ppriority + pwidth*.5: + return 2*q0*(ppriority+pwidth*.5-x)**2/pwidth**2 + else: + return 0 + + return ( + priority_func, + p.open(ppriority - pwidth*.5, ppriority + pwidth*.5) + ) + + @staticmethod + def normal(q0, ppriority, pwidth, pextra): + """ + Demand curve for normal shape. + The supply curve will be shaped as the integral of a normal + distribution. + + Parameters + ---------- + q0: float + The total demand/supply of the element. + ppriority: float + Specifies the midpoint of the curve (the mean of the + underlying distribution). + pwidth: float + Standard deviation of the underlying distribution. + pextra: float + Ignore. + + Returns + ------- + priority_func: function + The priority function. + + """ + def priority_func(x): + return q0*.5*(2-erfc((ppriority-x)/(np.sqrt(2)*pwidth))) + + # Normal distribution CDF is stricty monotonous in (-inf, inf). + # However, numerically it is only in a the range ~ (-8.29*sd, 8.29*sd) + return ( + priority_func, + p.open( + ppriority-8.2923611*pwidth, + ppriority+8.2923611*pwidth + ) + ) + + @staticmethod + def exponential(q0, ppriority, pwidth, pextra): + """ + Demand curve for exponential shape + The supply curve will be shaped as the integral of an + exponential distribution that is symmetric around its mean + (0.5*exp(-ABS(x-ppriority)/pwidth) on -∞ to ∞). + + Parameters + ---------- + q0: float + The total demand/supply of the element. + ppriority: float + Specifies the midpoint of the curve (the mean of the + underlying distribution). + pwidth: float + Multiplier on x in the underlying distribution. + pextra: float + Ignore. + + Returns + ------- + priority_func: function + The priority function. + + """ + def priority_func(x): + if x < ppriority: + return q0*(1-.5*np.exp((x-ppriority)/pwidth)) + else: + return q0*.5*np.exp((ppriority-x)/pwidth) + + # Exponential distribution CDF is stricty monotonous in (-inf, inf). + # However, numerically it is only in a the range ~ (-36.7*sd, 36.7*sd) + return ( + priority_func, + p.open( + ppriority-36.7368005696*pwidth, + ppriority+36.7368005696*pwidth + ) + ) + + @staticmethod + def constant_elasticity_demand(q0, ppriority, pwidth, pextra): + """ + Demand constant elasticity curve. + The curve will be a constant elasticity curve. + + Parameters + ---------- + q0: float + The total demand/supply of the element. + ppriority: float + Specifies the midpoint of the curve (the mean of the + underlying distribution). + pwidth: float + Standard deviation of the underlying distribution. + pextra: positive float + Elasticity exponent. + + Returns + ------- + priority_func: function + The priority function. + + """ + raise NotImplementedError( + "Some results for Vensim showed some bugs when using this " + "priority curve. Therefore, the curve is not implemented in " + "PySD as it cannot be properly tested." + ) + + @staticmethod + def constant_elasticity_supply(ppriority, pwidth, + pextra): # pragma: no cover + """ + Supply constant elasticity curve. + The curve will be a constant elasticity curve. + + Parameters + ---------- + q0: float + The total demand/supply of the element. + ppriority: float + Specifies the midpoint of the curve (the mean of the + underlying distribution). + pwidth: float + Standard deviation of the underlying distribution. + pextra: positive float + Elasticity exponent. + + Returns + ------- + priority_func: function + The priority function. + + """ + raise NotImplementedError( + "Some results for Vensim showed some bugs when using this " + "priority curve. Therefore, the curve is not implemented in " + "PySD as it cannot be properly tested." + ) + + +def _allocate_available_1d(request, pp, avail): + """ + This function implements the algorithm for allocate_available + to be passed for 1d numpy.arrays. The algorithm works as follows: + + 0. If supply > sum(request): return request. In the same way, + if supply = 0: return request*0 + 1. Based on the priority profiles and demands, the priority profiles + are computed. This profiles are returned with the interval where + each of them is strictly monotonous (or injective). + 2. Using the intervals of injectivity the initial guess is + selected depending on the available supply. + 3. The initial guess and injectivity interval are used to compute + the value where the sum of all priority functions is equal to + the avilable supply. This porcess is done using a least_squares + optimization function. + 4. The output from the previous step is used to compute the supply + to each target. + + Parameters + ---------- + request: numpy.ndarray (1D) + The request by target. Values must be non-negative. + pp: numpy.ndarray (2D) + The priority profiles of each target. + avail: float + The available supply. Must be non-negative. + + Returns + ------- + out: numpy.ndarray (1D) + The distribution of the supply. + + """ + if avail >= np.sum(request): + return request + if avail == 0: + return np.zeros_like(request) + + priorities, full_allocation, intervals =\ + Priorities.get_functions(request, pp, "demand") + + for interval, x_interval, x0 in intervals: + if avail in interval: + break + priority = least_squares( + lambda x: full_allocation(x) - avail, + x0, + bounds=(x_interval.lower, x_interval.upper), + method='dogbox', + tr_solver='exact', + ).x[0] + + return [allocate(priority) for allocate in priorities] + + +def allocate_available(request, pp, avail): + """ + Implements Vensim's ALLOCATE AVAILABLE function. + https://www.vensim.com/documentation/fn_allocate_available.html + + Parameters + ----------- + request: xarray.DataArray + Request of each target. Its shape should be the one of the + expected output of the function, having the allocation dimension + in the last position. + pp: xarray.DataArray + Priority of each target. Its shape should be the same as + request with an extra dimension for the priority profiles + in the last position. See Vensim's documentation for more + information https://www.vensim.com/documentation/24335.html + avail: float or xarray.DataArray + The total supply available to fulfill all requests. If the + supply exceeds total requests, all requests are filled, but + none are overfilled. If you wish to conserve material you must + compute supply minus total allocations explicitly. Its shape, + should be the same of request without the last dimension. + + Returns + ------- + out: xarray.DataArray + The distribution of the supply. + + Warning + ------- + This function uses an optimization method for resolution and the + given solution could differ from the one from Vensim. Particularly, + when close to the boundaries of the defined priority profiles. + + """ + if np.any(request < 0): + raise ValueError( + "There are some negative request values. Ensure that " + "your request is always non-negative. Allocation requires " + f"all quantities to be positive or 0.\n{request}") + + if np.any(avail < 0): + raise ValueError( + f"avail={avail} is not allowed. avail should be non-negative." + ) + + if len(request.shape) == 1: + # NUMPY: avoid '.values' and return directly the result of the + # function call + return xr.DataArray( + _allocate_available_1d( + request.values, pp.values, avail), + request.coords + ) + + # NUMPY: use np.empty_like and remove '.values' + out = xr.zeros_like(request, dtype=float) + for comb in itertools.product(*[range(i) for i in avail.shape]): + out.values[comb] = _allocate_available_1d( + request.values[comb], pp.values[comb], avail.values[comb]) + + return out def _allocate_by_priority_1d(request, priority, width, supply): diff --git a/pysd/tools/benchmarking.py b/pysd/tools/benchmarking.py index 694476ac..f493115f 100644 --- a/pysd/tools/benchmarking.py +++ b/pysd/tools/benchmarking.py @@ -252,13 +252,11 @@ def assert_frames_close(actual, expected, assertion="raise", + '\n\nActual values:\n\t'\ + np.array2string(actual[col].values, precision=precision, - separator=', ', - suppress_small=True)\ + separator=', ')\ + '\n\nDifference:\n\t'\ + np.array2string(expected[col].values-actual[col].values, precision=precision, - separator=', ', - suppress_small=True) + separator=', ') if assertion == "raise": raise AssertionError(assertion_details) diff --git a/pysd/translators/structures/abstract_expressions.py b/pysd/translators/structures/abstract_expressions.py index b6ba1fa5..e7453187 100644 --- a/pysd/translators/structures/abstract_expressions.py +++ b/pysd/translators/structures/abstract_expressions.py @@ -143,6 +143,31 @@ def __str__(self) -> str: # pragma: no cover return "GameStructure:\n\t%s" % self.expression +@dataclass +class AllocateAvailableStructure(AbstractSyntax): + """ + Dataclass for a Allocate Available structure. + + Parameters + ---------- + request: AbstractSyntax + The reference to the request variable. + pp: AbstractSyntax + The reference to the priority variable. + avail: AbstractSyntax or float + The total available supply. + + """ + request: AbstractSyntax + pp: AbstractSyntax + avail: Union[AbstractSyntax, float] + + def __str__(self) -> str: # pragma: no cover + return "AllocateAvailableStructure:\n\t%s,\n\t%s,\n\t%s" % ( + self.request, self.pp, self.avail + ) + + @dataclass class AllocateByPriorityStructure(AbstractSyntax): """ @@ -170,7 +195,7 @@ class AllocateByPriorityStructure(AbstractSyntax): def __str__(self) -> str: # pragma: no cover return "AllocateByPriorityStructure:"\ - "\"n\t%s,\n\t%s,\n\t%s,\n\t%s,\n\t%s" % ( + "\n\t%s,\n\t%s,\n\t%s,\n\t%s,\n\t%s" % ( self.request, self.priority, self.size, self.width, self.supply ) diff --git a/pysd/translators/vensim/vensim_structures.py b/pysd/translators/vensim/vensim_structures.py index f983fc26..11c79152 100644 --- a/pysd/translators/vensim/vensim_structures.py +++ b/pysd/translators/vensim/vensim_structures.py @@ -15,6 +15,7 @@ "with_lookup": ae.InlineLookupsStructure, "call": ae.CallStructure, "game": ae.GameStructure, + "allocate_available": ae.AllocateAvailableStructure, "allocate_by_priority": ae.AllocateByPriorityStructure, "get_xls_lookups": ae.GetLookupsStructure, "get_direct_lookups": ae.GetLookupsStructure, diff --git a/requirements.txt b/requirements.txt index 06a5c32c..e0947fec 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,3 +9,4 @@ black openpyxl scipy progressbar2 +portion diff --git a/tests/pytest_integration/pytest_integration_vensim_pathway.py b/tests/pytest_integration/pytest_integration_vensim_pathway.py index f8d6c8e3..cfb36b0a 100644 --- a/tests/pytest_integration/pytest_integration_vensim_pathway.py +++ b/tests/pytest_integration/pytest_integration_vensim_pathway.py @@ -19,6 +19,11 @@ "folder": "active_initial_circular", "file": "test_active_initial_circular.mdl" }, + "allocate_available": { + "folder": "allocate_available", + "file": "test_allocate_available.mdl", + "rtol": 2e-2 + }, "allocate_by_priority": { "folder": "allocate_by_priority", "file": "test_allocate_by_priority.mdl" diff --git a/tests/pytest_pysd/pytest_allocation.py b/tests/pytest_pysd/pytest_allocation.py index 93230067..558f8313 100644 --- a/tests/pytest_pysd/pytest_allocation.py +++ b/tests/pytest_pysd/pytest_allocation.py @@ -1,11 +1,313 @@ import pytest import numpy as np import xarray as xr +import portion as p from pysd.py_backend.allocation import\ + Priorities, allocate_available,\ allocate_by_priority, _allocate_by_priority_1d +class TestPriorities(): + @pytest.mark.parametrize( + "q0,pp", + [ + ( # rectangular + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[1, 5, 1, 0], [1, 10, 3, 0], [1, 4, 0.5, 0], + [1, 4, 0.5, 0], [1, 2, 1, 0], [1, 6, 3, 0]] + ) + ), + ( # triangular + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[2, 5, 1, 0], [2, 10, 3, 0], [2, 4, 0.5, 0], + [2, 4, 0.5, 0], [2, 2, 1, 0], [2, 6, 3, 0]] + ) + ), + ( # normal + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[3, 5, 1, 0], [3, 10, 3, 0], [3, 4, 0.5, 0], + [3, 4, 0.5, 0], [3, 2, 1, 0], [3, 6, 3, 0]] + ) + ), + ( # exponential + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[4, 5, 1, 0], [4, 10, 3, 0], [4, 4, 0.5, 0], + [4, 4, 0.5, 0], [4, 2, 1, 0], [4, 6, 3, 0]] + ) + ), + ], + ids=["rectangular", "triangular", "normal", "exponential"] + ) + def test_supply_demand(self, pp, q0): + d_funcs, d_all, d_inter = Priorities.get_functions(q0, pp, "demand") + + pty = np.arange(0, 15, 0.05) + for i in range(1, len(pty)): + # assert the sum of functions + assert d_all(pty[i]) == np.sum([d_f(pty[i]) for d_f in d_funcs]) + for d_f in d_funcs: + # assert functions are monotonous + assert d_f(pty[i-1]) >= d_f(pty[i]) + + non_monotony = True + for s_i, p_i, p_m in d_inter: + if pty[i] in p_i or pty[i] == p_i.upper: + # assert functions are strictly monotonous in the interval + assert d_all(pty[i-1]) > d_all(pty[i]) + non_monotony = False + + if non_monotony: + # assert functions are constant out the interval + assert d_all(pty[i-1]) == d_all(pty[i]) + + for s_i, p_i, p_m in d_inter: + assert p_m == .5 * (p_i.lower+p_i.upper) + assert s_i.lower == d_all(p_i.upper) + assert s_i.upper == d_all(p_i.lower) + + @pytest.mark.parametrize( + "q0,pp", + [ + ( # rectangular + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[1, 5, 1, 0], [1, 10, 3, 0], [1, 4, 0.5, 0], + [1, 4, 0.5, 0], [1, 2, 1, 0], [1, 6, 3, 0]] + ) + ), + ( # triangular + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[2, 5, 1, 0], [2, 10, 3, 0], [2, 4, 0.5, 0], + [2, 4, 0.5, 0], [2, 2, 1, 0], [2, 6, 3, 0]] + ) + ), + ( # normal + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[3, 5, 1, 0], [3, 10, 3, 0], [3, 4, 0.5, 0], + [3, 4, 0.5, 0], [3, 2, 1, 0], [3, 6, 3, 0]] + ) + ), + ( # exponential + np.array([6, 1, 3, 0, 0.5, 0]), + np.array( + [[4, 5, 1, 0], [4, 10, 3, 0], [4, 4, 0.5, 0], + [4, 4, 0.5, 0], [4, 2, 1, 0], [4, 6, 3, 0]] + ) + ), + ], + ids=["rectangular", "triangular", "normal", "exponential"] + ) + def test_supply_supply(self, pp, q0): + error_message = "get_function_supply is not implemented" + with pytest.raises(NotImplementedError, match=error_message): + Priorities.get_functions(q0, pp, "supply") + + @pytest.mark.parametrize( + "q0,pp,distance", + [ + ( # rectangular + np.array([6, 1, 3, 2, 0.5, 4]), + np.array( + [[1, 5, 1, 0], [1, 10, 3, 0], [1, 4, 0.5, 0], + [1, 4, 0.5, 0], [1, 2, 1, 0], [1, 6, 3, 0]] + ), + 0.5 + ), + ( # triangular + np.array([6, 1, 3, 3, 0.5, 4]), + np.array( + [[2, 5, 1, 0], [2, 10, 3, 0], [2, 4, 0.5, 0], + [2, 4, 0.5, 0], [2, 2, 1, 0], [2, 6, 3, 0]] + ), + 0.5 + ), + ( # normal + np.array([6, 1, 3, 3, 0.5, 4]), + np.array( + [[3, 5, 1, 0], [3, 10, 3, 0], [3, 4, 0.5, 0], + [3, 4, 0.5, 0], [3, 2, 1, 0], [3, 6, 3, 0]] + ), + 8.2923611 + ), + ( # exponential + np.array([6, 1, 3, 3, 0.5, 4]), + np.array( + [[4, 5, 1, 0], [4, 10, 3, 0], [4, 4, 0.5, 0], + [4, 4, 0.5, 0], [4, 2, 1, 0], [4, 6, 3, 0]] + ), + 36.7368005696 + ), + ], + ids=["rectangular", "triangular", "normal", "exponential"] + ) + def test_priority_shape_demand(self, pp, q0, distance): + for i in range(len(q0)): + func, interval = Priorities.get_function_demand(q0[i], pp[i]) + xs = np.linspace( + pp[i][1]-pp[i][2]*distance*0.95, + pp[i][1]+pp[i][2]*distance*0.95, + 100 + ) + for j in range(1, len(xs)): + # assert is monotically increasing + assert func(xs[j-1]) > func(xs[j]) + + xs_lower = np.linspace( + pp[i][1]-3*pp[i][2]*distance, + pp[i][1]-1.01*pp[i][2]*distance, + 100 + ) + for x in xs_lower: + assert func(x) == q0[i] + xs_upper = np.linspace( + pp[i][1]+1.01*pp[i][2]*distance, + pp[i][1]+3*pp[i][2]*distance, + 100 + ) + for x in xs_upper: + assert np.isclose(func(x), 0, atol=1e-15) + + assert interval == p.open( + pp[i][1]-pp[i][2]*distance, + pp[i][1]+pp[i][2]*distance) + + @pytest.mark.parametrize( + "pp,distance", + [ + ( # rectangular + np.array([1, 5, 1, 0]), + 0.5 + ), + ( # triangular + np.array([2, 5, 1, 0]), + 0.5 + ), + ( # normal + np.array([3, 5, 1, 0]), + 8.2923611 + ), + ( # exponential + np.array([4, 5, 1, 0]), + 36.7368005696 + ), + ], + ids=["rectangular", "triangular", "normal", "exponential"] + ) + def test_priority_zeros_demand(self, pp, distance): + func, interval = Priorities.get_function_demand(0, pp) + xs = np.linspace( + pp[1]-3*pp[2]*distance, + pp[1]+3*pp[2]*distance, + 500 + ) + for x in xs: + # assert is monotically increasing + assert func(x) == 0 + + assert interval == p.empty() + + @pytest.mark.parametrize( + "pp,ptype,raise_type,error_message", + [ + ( # invalid-kind + np.array([[1, 10, 1, 0]]), + "invalid", + ValueError, + r"kind='invalid' is not allowed\. kind should be " + r"'demand' or 'supply'\." + ), + ( # fixed-quantity + np.array([[0, 10, 1, 0]]), + "demand", + NotImplementedError, + r"fixed_quantity priority profile is not implemented\." + ), + ( # constant-elasticity-demand + np.array([[5, 10, 1, 0]]), + "demand", + NotImplementedError, + r"Some results for Vensim showed some bugs when using " + r"this priority curve\. Therefore, the curve is not " + r"implemented in PySD as it cannot be properly tested\." + ), + ( # invalid-func + np.array([[8, 10, 1, 0]]), + "demand", + ValueError, + r"The priority function for pprofile=8 is not valid\." + ), + ( # supply + np.array([[1, 10, 1, 0]]), + "supply", + NotImplementedError, + r"get_function_supply is not implemented\." + ), + ( # negative-width + np.array([[1, 10, -1, 0]]), + "demand", + ValueError, + r"pwidth values must be positive\." + ), + ( # zero-width + np.array([[1, 10, 0, 0]]), + "demand", + ValueError, + r"pwidth values must be positive\." + ), + ], + ids=["invalid-kind", "fixed-quantity", "constant-elasticity-demand", + "invalid-func", "supply", "negative-width", "zero-width"] + ) + def test_priorities_errors(self, pp, ptype, raise_type, error_message): + with pytest.raises(raise_type, match=error_message): + Priorities.get_functions(np.array([100]), pp, ptype) + + +class TestAllocateAvailable(): + + @pytest.mark.parametrize( + "requests,pp,avail,raise_type,error_message", + [ + ( # negative-request + xr.DataArray([6, -3, 3], {'dim': ["A", "B", "C"]}), + xr.DataArray( + [[1, 10, 1, 0], [1, 10, 1, 0], [1, 10, 1, 0]], + {'dim': ["A", "B", "C"], + 'pprofile': ["ptype", "ppriority", "pwidth", "pextra"]} + ), + 15, + ValueError, + r"There are some negative request values\. Ensure that " + r"your request is always non-negative\. Allocation requires " + r"all quantities to be positive or 0\.\n.*" + ), + ( # negative supply + xr.DataArray([6, 3, 3], {'dim': ["A", "B", "C"]}), + xr.DataArray( + [[1, 10, 1, 0], [1, 10, 1, 0], [1, 10, 1, 0]], + {'dim': ["A", "B", "C"], + 'pprofile': ["ptype", "ppriority", "pwidth", "pextra"]} + ), + -7.5, + ValueError, + r"avail=-7\.5 is not allowed\. " + r"avail should be non-negative\." + ), + ], + ) + def test_allocate_available_errors(self, requests, pp, avail, + raise_type, error_message): + with pytest.raises(raise_type, match=error_message): + allocate_available(requests, pp, avail) + + class TestAllocateByPriority(): @pytest.mark.parametrize( "requests,priority,width,supply,expected", @@ -116,7 +418,5 @@ def test__allocate_by_priority_1d(self, requests, priority, width, def test_allocate_by_priority_errors(self, requests, priority, width, supply, raise_type, error_message): - # Test some simple cases, the complicate cases are tested with - # a full integration test with pytest.raises(raise_type, match=error_message): allocate_by_priority(requests, priority, width, supply) diff --git a/tests/pytest_pysd/pytest_utils.py b/tests/pytest_pysd/pytest_utils.py index edcd8d96..91d6dafb 100644 --- a/tests/pytest_pysd/pytest_utils.py +++ b/tests/pytest_pysd/pytest_utils.py @@ -259,7 +259,7 @@ def test_make_flat_df_flatten_transposed(self): 'Elem2': ('elem2', {})} actual = pysd.utils.make_flat_df(df, return_addresses, flatten=True) - print(actual.columns) + # check all columns are in the DataFrame assert set(actual.columns) == set(expected.columns) # need to assert one by one as they are xarrays diff --git a/tests/pytest_types/data/pytest_data_with_model.py b/tests/pytest_types/data/pytest_data_with_model.py index de87ee3f..9bd65097 100644 --- a/tests/pytest_types/data/pytest_data_with_model.py +++ b/tests/pytest_types/data/pytest_data_with_model.py @@ -148,7 +148,6 @@ def test_run_error(self, data_model, shared_tmpdir): reason=r"bad scape \e") def test_loading_error(self, data_model, data_files, raise_type, error_message, shared_tmpdir): - print(error_message % (data_files)) with pytest.raises(raise_type, match=error_message % (data_files)): self.model( data_model, data_files, shared_tmpdir) diff --git a/tests/test-models b/tests/test-models index cdb82c08..a62270a4 160000 --- a/tests/test-models +++ b/tests/test-models @@ -1 +1 @@ -Subproject commit cdb82c081b55fcc2d90b73c193a5e4413c5737b0 +Subproject commit a62270a4c079de8288daac75ffcc77dc84a769aa