diff --git a/.gitignore b/.gitignore index 9aa6e8c..001b32a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,7 @@ +typings/* boario/archive.py +archive.py +*.el *.lprof *.prof *.zip @@ -162,3 +165,5 @@ workflow/config_full.json DOC_BUILDING_REMINDER new_docs +pyproject.toml.bak +.vscode/settings.json diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index 474d270..10da363 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -3,3 +3,102 @@ Contributing ################ .. _contributions: + +Contributions are welcome, and they are greatly appreciated! +Every little bit helps, and credit will always be given. +If you are interested in contributing, please read all this section +before doing anything and in case of doubt, `contact the developer`_! + +.. _contact the developer: pro@sjuhel.org + +---------------------- +Types of Contributions +---------------------- + +You can contribute in many ways: + +Report Bugs +=========== + +Before reporting a bug, **Ensure the bug was not already reported** by searching on GitHub under `Issues `_. + +If you are reporting a bug, please use the "Bug report template" `here `_, and include: + +* Your operating system name and version. +* Any details about your local setup that might be helpful in troubleshooting. +* Detailed steps to reproduce the bug. + + +Fix Bugs +======== + +Look through the Github `issues `_ for bugs. +If you want to start working on a bug then please write short message on the issue tracker to prevent duplicate work. + + +Implement Features +================== + +Look through the Github issues for features or add your own request +using the Code feature request template `here `_ +and assign yourself. + +Write Documentation +=================== + +BoARIO could always use more documentation, whether as part of the official documentation website, in docstrings, in Jupyter Notebook examples, +or even on the web in blog posts, articles, and such. + +BoARIO uses `Sphinx `_ for the user manual (that you are currently reading). + + +Submit Feedback +=============== + +The best way to send feedback is to file an issue at `here `_ + +If you are proposing a feature: + +* Explain in detail how it would work. +* Keep the scope as narrow as possible, to make it easier to implement. +* Remember that this is a volunteer-driven project, and that contributions are welcome :) + +----------------------- +Pull Request Guidelines +----------------------- + +To update the documentation, fix bugs or add new features you need to create a Pull Request. +A PR is a change you make to your local copy of the code for us to review and potentially integrate into the code base. + +To create a Pull Request you need to do these steps: + +1. Create a Github account. +2. Fork the repository. +3. Clone your fork locally. +4. Go to the created BoARIO folder with :code:`cd boario`. +5. Create a new branch with :code:`git checkout -b `. +6. Make your changes to the code or documentation. +7. Run :code:`git add .` to add all the changed files to the commit (to see what files will be added you can run :code:`git add . --dry-run`). +8. To commit the added files use :code:`git commit`. (This will open a command line editor to write a commit message. These should have a descriptive 80 line header, followed by an empty line, and then a description of what you did and why. To use your command line text editor of choice use (for example) :code:`export GIT_EDITOR=vim` before running :code:`git commit`). +9. Now you can push your changes to your Github copy of BoARIO by running :code:`git push origin `. +10. If you now go to the webpage for your Github copy of BoARIO you should see a link in the sidebar called "Create Pull Request". +11. Now you need to choose your PR from the menu and click the "Create pull request" button. Be sure to change the pull request target branch to ! + +If you want to create more pull requests, first run :code:`git checkout main` and then start at step 5. with a new branch name. + +Feel free to ask questions about this if you want to contribute to BoARIO :) + +------------------ +Testing Guidelines +------------------ + +To ensure that you do not introduce bugs into BoARIO, you should test your code thouroughly. + +From the base BoARIO folder you call :code:`pytest` to run all the tests, or choose one specific test: + +.. code-block:: console + + $ pytest + $ pytest tests/tests.py::test_log_input + +If you introduce a new feature you should add a new test to the tests directory. See the folder for examples. \ No newline at end of file diff --git a/README.rst b/README.rst index 7c5d0be..e72bd8f 100644 --- a/README.rst +++ b/README.rst @@ -7,6 +7,14 @@ BoARIO BoARIO : The Adaptative Regional Input Output model in python. +.. _`Documentation Website`: https://spjuhel.github.io/BoARIO/boario-what-is.html + +Disclaimer +=========== + +This work is in progress. The code has not been extensively tested. Any results produced with the model should be interpreted with great care. +Do not hesitate to contact the author when using the model ! + What is BoARIO ? ================= @@ -23,18 +31,22 @@ It is still an ongoing project (in parallel of a PhD project). .. _`Hallegatte 2013`: https://doi.org/10.1111/j.1539-6924.2008.01046.x .. _`Pymrio`: https://pymrio.readthedocs.io/en/latest/intro.html -Here is a non-exhaustive chronology of academic works with or about the ARIO model : +You can find most academic literature using ARIO or related models :ref:`here ` -.. image:: https://raw.githubusercontent.com/spjuhel/BoARIO/master/imgs/chronology.svg?sanitize=true - :width: 900 - :alt: ARIO academic work chronology What is ARIO ? =============== -ARIO stands for Adaptive Regional Input-Output. It is an hybrid input-output / agent-based economic model, designed to compute indirect costs from economic shocks. Its first version dates back to 2008 and has originally been developed to assess the indirect costs of natural disasters (Hallegatte 2008). +ARIO stands for Adaptive Regional Input-Output. It is an hybrid input-output / agent-based economic model, +designed to compute indirect costs from economic shocks. Its first version dates back to 2008 and has originally +been developed to assess the indirect costs of natural disasters (Hallegatte 2008). -In ARIO, the economy is modelled as a set of economic sectors and a set of regions. Each economic sector produces its generic product and draws inputs from an inventory. Each sector answers to a total demand consisting of a final demand (household consumption, public spending and private investments) of all regions (local demand and exports) and intermediate demand (through inputs inventory resupply). An initial equilibrium state of the economy is built based on multi-regional input-output tables (MRIO tables). +In ARIO, the economy is modelled as a set of economic sectors and a set of regions. +Each economic sector produces its generic product and draws inputs from an inventory. +Each sector answers to a total demand consisting of a final demand (household consumption, +public spending and private investments) of all regions (local demand and exports) and +intermediate demand (through inputs inventory resupply). An initial equilibrium state of +the economy is built based on multi-regional input-output tables (MRIO tables). Where to get BoARIO ? @@ -92,7 +104,7 @@ This work was supported by the French Environment and Energy Management Agency Development ------------ -* Samuel Juhel +* Samuel Juhel (pro@sjuhel.org) Contributions --------------- @@ -102,4 +114,5 @@ All :ref:`contributions` to the project are welcome ! Acknowledgements ------------------ -I would like to thank Vincent Viguie, Fabio D'Andrea my PhD supervisors as well as Célian Colon and Alessio Ciulo for their inputs during the model implementation. +I would like to thank Vincent Viguie, Fabio D'Andrea my PhD supervisors as well as Célian Colon, Alessio Ciulo and Adrien Delahais +for their inputs during the model implementation. diff --git a/boario/__init__.py b/boario/__init__.py index 4f601ba..7311b8e 100644 --- a/boario/__init__.py +++ b/boario/__init__.py @@ -20,16 +20,11 @@ else: _has_coloredlogs = True -try: - import pygit2 -except: - pass - import pathlib import logging from functools import lru_cache -__version__ = "v0.4.1b0" +__version__ = "v0.5.0a" __author__ = "sjuhel " # __minimum_python_version__ = "3.8" @@ -49,13 +44,13 @@ datefmt="%H:%M:%S", ) - try: + import pygit2 + __git_branch__ = pygit2.Repository(__file__).head.name -except Exception: - pass -else: - logger.info(f"You are using boario from branch {__git_branch__}") + logger.info("You are using boario from branch %s", __git_branch__) +except ModuleNotFoundError: + logger.info("Unable to tell git branch as pygit2 was not found.") @lru_cache(10) diff --git a/boario/event.py b/boario/event.py index f95b1bb..9c7a37a 100644 --- a/boario/event.py +++ b/boario/event.py @@ -14,15 +14,23 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . from __future__ import annotations -import abc -from typing import Callable, Optional, Union, List, Tuple -from numpy.typing import ArrayLike +from abc import ABC, abstractmethod +from typing import Callable, Optional, Union, List, Tuple, get_origin, get_args +import numpy.typing as npt import numpy as np import pandas as pd -from boario import logger import math import inspect from functools import partial +import warnings + +from boario import logger +from boario.utils.recovery_functions import ( + concave_recovery, + convexe_recovery, + linear_recovery, + convexe_recovery_scaled, +) __all__ = [ "Event", @@ -36,160 +44,33 @@ "RegionsList", ] -Impact = Union[int, float, list, dict, np.ndarray, pd.DataFrame, pd.Series] -IndustriesList = Union[List[Tuple[str, str]], pd.MultiIndex, np.ndarray] -SectorsList = Union[List[str], pd.Index, np.ndarray] -RegionsList = Union[List[str], pd.Index, np.ndarray] - - -def linear_recovery( - elapsed_temporal_unit: int, init_kapital_destroyed: np.ndarray, recovery_time: int -): - r"""Linear Kapital recovery function - - Kapital is entirely recovered when `recovery_time` has passed since event - started recovering - - Parameters - ---------- - - init_kapital_destroyed : float - Initial kapital destroyed - - elapsed_temporal_unit : int - Elapsed time since event started recovering - - recovery_time : int - Total time it takes the event to fully recover - - """ - - return init_kapital_destroyed * (1 - (elapsed_temporal_unit / recovery_time)) - - -def convexe_recovery( - elapsed_temporal_unit: int, init_kapital_destroyed: np.ndarray, recovery_time: int -): - r"""Convexe Kapital recovery function - - Kapital is recovered with characteristic time `recovery_time`. (This doesn't mean Kapital is fully recovered after this time !) - This function models a recovery similar as the one happening in the rebuilding case, for the same characteristic time. - - Parameters - ---------- - - init_kapital_destroyed : float - Initial kapital destroyed - - elapsed_temporal_unit : int - Elapsed time since event started recovering - - recovery_time : int - Total time it takes the event to fully recover - - """ - - return init_kapital_destroyed * (1 - (1 / recovery_time)) ** elapsed_temporal_unit +VectorImpact = Union[list, dict, np.ndarray, pd.DataFrame, pd.Series] +ScalarImpact = Union[int, float, np.integer] +Impact = Union[VectorImpact, ScalarImpact] +IndustriesList = Union[List[Tuple[str, str]], pd.MultiIndex] +SectorsList = Union[List[str], pd.Index, str] +RegionsList = Union[List[str], pd.Index, str] +FinalCatList = Union[List[str], pd.Index, str] +REBUILDING_FINALDEMAND_CAT_REGEX = ( + r"(?i)(?=.*household)(?=.*final)(?!.*NPISH|profit).*|HFCE" +) -def convexe_recovery_scaled( - elapsed_temporal_unit: int, - init_kapital_destroyed: np.ndarray, - recovery_time: int, - scaling_factor: float = 4, -): - r"""Convexe Kapital recovery function (scaled to match other recovery duration) +LOW_DEMAND_THRESH = 10 - Kapital is mostly recovered (>95% by default for most cases) when `recovery_time` has passed since event - started recovering. - Parameters - ---------- - - init_kapital_destroyed : float - Initial kapital destroyed - - elapsed_temporal_unit : int - Elapsed time since event started recovering - - recovery_time : int - Total time it takes the event to fully recover - - scaling_factor: float - Used to scale the exponent in the function so that kapital is mostly rebuilt after `recovery_time`. A value of 4 insure >95% of kapital is recovered for a reasonable range of `recovery_time` values. - - """ - - return init_kapital_destroyed * (1 - (1 / recovery_time)) ** ( - scaling_factor * elapsed_temporal_unit - ) - - -def concave_recovery( - elapsed_temporal_unit: int, - init_kapital_destroyed: np.ndarray, - recovery_time: int, - steep_factor: float = 0.000001, - half_recovery_time: Optional[int] = None, -): - r"""Concave (s-shaped) Kapital recovery function - - Kapital is mostly (>95% in most cases) recovered when `recovery_time` has passed since event started recovering. - - Parameters - ---------- - - init_kapital_destroyed : float - Initial kapital destroyed - - elapsed_temporal_unit : int - Elapsed time since event started recovering - - recovery_time : int - Total time it takes the event to fully recover - - steep_factor: float - This coefficient governs the slope of the central part of the s-shape, smaller values lead to a steeper slope. As such it also affect the percentage of kapital rebuilt after `recovery_time` has elapsed. A value of 0.000001 should insure 95% of the kapital is rebuild for a reasonable range of recovery duration. - - half_recovery_time : int - This can by use to change the time the inflexion point of the s-shape curve is attained. By default it is set to half the recovery duration. - - """ - - if half_recovery_time is None: - tau_h = 2 - else: - tau_h = recovery_time / half_recovery_time - exponent = (np.log(recovery_time) - np.log(steep_factor)) / ( - np.log(recovery_time) - np.log(tau_h) - ) - return (init_kapital_destroyed * recovery_time) / ( - recovery_time + steep_factor * (elapsed_temporal_unit**exponent) - ) - - -class Event(metaclass=abc.ABCMeta): +class Event(ABC): # Class Attributes - __required_class_attributes = [ - "possible_sectors", - "possible_regions", - "temporal_unit_range", - "z_shape", - "y_shape", - "x_shape", - "regions_idx", - "sectors_idx", - "monetary_factor", - "sectors_gva_shares", - "Z_distrib", - "mrio_name", - ] - possible_sectors: SectorsList = pd.Index([]) + + possible_sectors: pd.Index = pd.Index([]) r"""List of sectors present in the MRIO used by the model""" - possible_regions: RegionsList = pd.Index([]) + possible_regions: pd.Index = pd.Index([]) r"""List of regions present in the MRIO used by the model""" + possible_final_demand_cat: pd.Index = pd.Index([]) + r"""List of final demand categories present in the MRIO used by the model""" + temporal_unit_range: int = 0 r"""Maximum temporal unit simulated""" @@ -202,37 +83,38 @@ class Event(metaclass=abc.ABCMeta): x_shape: tuple[int, int] = (0, 0) r"""Shape of the x (production) vector in the model""" - regions_idx: np.ndarray = np.array([]) + regions_idx: npt.NDArray = np.array([]) r"""lexicographic region indexes""" - sectors_idx: np.ndarray = np.array([]) + sectors_idx: npt.NDArray = np.array([]) r"""lexicographic sector indexes""" - monetary_factor: int = 0 + model_monetary_factor: int = 1 r"""Amount of unitary currency used in the MRIO (e.g. 1000000 if in € millions)""" - sectors_gva_shares: np.ndarray = np.array([]) + gva_df: pd.Series = pd.Series([], dtype="float64") + r"""GVA per (region,sector)""" + + sectors_gva_shares: npt.NDArray = np.array([]) r"""Fraction of total (regional) GVA for each sectors""" - Z_distrib: np.ndarray = np.array([]) + Z_distrib: npt.NDArray = np.array([]) r"""Normalized intermediate consumption matrix""" + Y_distrib: npt.NDArray = np.array([]) + r"""Normalized final consumption matrix""" + mrio_name: str = "" r"""MRIO identification""" + @abstractmethod def __init__( self, - impact: Impact, - aff_regions: Optional[RegionsList] = None, - aff_sectors: Optional[SectorsList] = None, - aff_industries: Optional[IndustriesList] = None, - impact_industries_distrib=None, - impact_regional_distrib=None, - impact_sectoral_distrib_type="custom", - impact_sectoral_distrib=None, - name="Unnamed", - occurrence=1, - duration=1, + *, + impact: pd.Series, + name: str | None = "Unnamed", + occurrence: int = 1, + duration: int = 1, ) -> None: r"""Create an event shocking the model from a dictionary. @@ -252,206 +134,25 @@ def __init__( FIXME: Add docs. """ - - self._aff_sectors_idx = None - self._aff_sectors = None - self._aff_regions_idx = None - self._aff_regions = None - logger.debug("Initializing new Event") + logger.info("Initializing new Event") logger.debug("Checking required Class attributes are defined") - for v in Event.__required_class_attributes: - if Event.__dict__[v] is None: - raise AttributeError( - "Required Event Class attribute {} is not set yet so instantiating an Event isn't possible".format( - v - ) - ) - - self.name: str = name - r"""An identifying name for the event (for convenience mostly)""" - - self.occurrence = occurrence - self.duration = duration - self.impact_df: pd.Series = pd.Series( - 0, - dtype="float64", - index=pd.MultiIndex.from_product( - [self.possible_regions, self.possible_sectors], - names=["region", "sector"], - ), - ) - r"""A pandas Series with all possible industries as index, holding the impact vector of the event. The impact is defined for each sectors in each region.""" - - ################## DATAFRAME INIT ################# - # CASE VECTOR 1 (everything is there and regrouped) (only df creation) - if isinstance(impact, pd.Series): - logger.debug("Given impact is a pandas Series") - self.impact_df.loc[impact.index] = impact - if self.name == "Unnamed" and not impact.name is None: - self.name = str(impact.name) - elif isinstance(impact, dict): - logger.debug("Given impact is a dict, converting it to pandas Series") - impact = pd.Series(impact) - self.impact_df.loc[impact.index] = impact - elif isinstance(impact, pd.DataFrame): - logger.debug("Given Impact is a pandas DataFrame, squeezing it to a Series") - impact = impact.squeeze() - if not isinstance(impact, pd.Series): - raise ValueError( - "The given impact DataFrame is not a Series after being squeezed" - ) - self.impact_df.loc[impact.index] = impact - # CASE VECTOR 2 (everything is there but not regrouped) AND CASE SCALAR (Only df creation) - elif ( - isinstance(impact, (int, float, list, np.ndarray)) - and aff_industries is not None - ): - logger.debug( - f"Given Impact is a {type(impact)} and list of impacted industries given. Proceeding." - ) - self.impact_df.loc[aff_industries] = impact - elif ( - isinstance(impact, (int, float, list, np.ndarray)) - and aff_regions is not None - and aff_sectors is not None - ): - logger.debug( - f"Given Impact is a {type(impact)} and lists of impacted regions and sectors given. Proceeding." - ) - if isinstance(aff_regions, str): - aff_regions = [aff_regions] - if isinstance(aff_sectors, str): - aff_sectors = [aff_sectors] - - self.impact_df.loc[ - pd.MultiIndex.from_product([aff_regions, aff_sectors]) - ] = impact - else: - raise ValueError("Invalid input format. Could not initiate pandas Series.") - # Check for <0 values and remove 0. - if (self.impact_df < 0).any(): - logger.warning( - "Found negative values in impact vector. This should raise concern" + if np.size(self.possible_regions) == 0 or np.size(self.possible_sectors) == 0: + raise AttributeError( + "It appears that no model has been instantiated as some class attributes are not initialized (possible_regions, possible_sectors). Events require to instantiate a model and a simulation context before they can be instantiated" ) - # SORT DF - # at this point impact_df is built, and can be sorted. Note that if impact was a scalar, impact_df contains copies of this scalar. - logger.debug("Sorting Impact Series") - self.impact_df = self.impact_df.sort_index() - - # Init self.impact_sectoral_distrib_type, - self.impact_sectoral_distrib_type = impact_sectoral_distrib_type - ################################################# - - # SET INDEXES ATTR - # note that the following also sets aff_regions and aff_sectors - assert isinstance(self.impact_df.index, pd.MultiIndex) - # Only look for industries where impact is greater than 0 - self.aff_industries = self.impact_df.loc[self.impact_df > 0].index - - logger.debug( - f"Impact df at the moment:\n {self.impact_df.loc[self.aff_industries]}" - ) - - ###### SCALAR DISTRIBUTION ###################### - # if impact_industries_distrib is given, set it. We assume impact is scalar ! - # CASE SCALAR + INDUS DISTRIB - if impact_industries_distrib is not None and not isinstance( - impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) - ): - logger.debug("Impact is Scalar and impact_industries_distrib was given") - self.impact_industries_distrib = np.array(impact_industries_distrib) - self.impact_df.loc[self.aff_industries] = ( - self.impact_df.loc[self.aff_industries] * self.impact_industries_distrib - ) - # if impact_reg_dis and sec_dis are give, deduce the rest. We also assume impact is scalar ! - # CASE SCALAR + REGION and SECTOR DISTRIB - elif ( - impact_regional_distrib is not None - and impact_sectoral_distrib is not None - and not isinstance( - impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) - ) - ): - logger.debug( - "Impact is Scalar and impact_regional_distrib and impact_sectoral_distrib were given" - ) - if len(impact_regional_distrib) != len(self.aff_regions) or len( - impact_sectoral_distrib - ) != len(self.aff_sectors): - raise ValueError( - "Lengths of `impact_regional_distrib` and/or `impact_sectoral_distrib` are incompatible with `aff_regions` and/or `aff_sectors`." - ) - else: - self.impact_regional_distrib = np.array(impact_regional_distrib) - self.impact_sectoral_distrib = np.array(impact_sectoral_distrib) - self.impact_industries_distrib = ( - self.impact_regional_distrib[:, np.newaxis] - * self.impact_sectoral_distrib - ).flatten() - self.impact_df.loc[self.aff_industries] = ( - self.impact_df.loc[self.aff_industries] - * self.impact_industries_distrib - ) - # CASE SCALAR + 'gdp' distrib - elif ( - impact_regional_distrib is not None - and impact_sectoral_distrib_type is not None - and impact_sectoral_distrib_type == "gdp" - and not isinstance( - impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) + if self.temporal_unit_range == 0: + raise AttributeError( + "It appears that no simulation context has been instantiated as some class attributes are not initialized (temporal_unit_range). Events require to instantiate a model and a simulation context before they can be instantiated" ) - ): - logger.debug("Impact is Scalar and impact_sectoral_distrib_type is 'gdp'") - - self.impact_regional_distrib = np.array(impact_regional_distrib) - shares = self.sectors_gva_shares.reshape( - (len(self.possible_regions), len(self.possible_sectors)) - ) - self.impact_sectoral_distrib = ( - shares[self._aff_regions_idx][:, self._aff_sectors_idx] - / shares[self._aff_regions_idx][:, self._aff_sectors_idx].sum(axis=1)[ - :, np.newaxis - ] - ) - self.impact_industries_distrib = ( - self.impact_regional_distrib[:, np.newaxis] - * self.impact_sectoral_distrib - ).flatten() - self.impact_df.loc[self.aff_industries] = ( - self.impact_df.loc[self.aff_industries] * self.impact_industries_distrib - ) - self.impact_sectoral_distrib_type = "gdp" - # CASE SCALAR + NO DISTRIB + list of industries - # if neither was given, we use default values. Again impact should be scalar here ! - elif isinstance(aff_industries, (list, np.ndarray)) and not isinstance( - impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) - ): - logger.debug( - "Impact is Scalar and no distribution was given but a list of affected industries was given" - ) - self._default_distribute_impact_from_industries_list() - self.impact_sectoral_distrib_type = "default (shared equally between affected regions and then affected sectors)" - # CASE SCALAR + NO DISTRIB + list of region + list of sectors - elif ( - aff_regions is not None - and aff_sectors is not None - and not isinstance( - impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) - ) - ): - logger.debug( - "Impact is Scalar and no distribution was given but lists of regions and sectors affected were given" - ) - self._default_distribute_impact_from_industries_list() - self.impact_sectoral_distrib_type = "default (shared equally between affected regions and then affected sectors)" - elif not isinstance(impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray)): - raise ValueError(f"Invalid input format: Could not compute impact") + self.name: str = name if name else "unnamed" + r"""An identifying name for the event (for convenience mostly)""" - self._finish_init() - ################################################## + self.occurrence = occurrence + self.duration = duration + self.impact_df = impact self.happened: bool = False r"""States if the event happened""" @@ -468,7 +169,6 @@ def __init__( "impact": self.total_impact, "impact_industries_distrib": list(self.impact_industries_distrib), "impact_regional_distrib": list(self.impact_regional_distrib), - "impact_sectoral_distrib_type": self.impact_sectoral_distrib_type, "globals_vars": { "possible_sectors": list(self.possible_sectors), "possible_regions": list(self.possible_regions), @@ -476,53 +176,358 @@ def __init__( "z_shape": self.z_shape, "y_shape": self.y_shape, "x_shape": self.x_shape, - "monetary_factor": self.monetary_factor, + "model_monetary_factor": self.model_monetary_factor, "mrio_used": self.mrio_name, }, } r"""Store relevant information about the event""" - def _default_distribute_impact_from_industries_list(self): - # at this point, impact should still be scalar. - logger.debug("Using default impact distribution to industries") - logger.debug( - f"Impact df at the moment:\n {self.impact_df.loc[self.aff_industries]}" - ) - self.impact_regional_distrib = np.full( - len(self.aff_regions), 1 / len(self.aff_regions) + @classmethod + @abstractmethod + def _instantiate( + cls, + impact: pd.Series, + *, + occurrence: int = 1, + duration: int = 1, + name: Optional[str] = None, + ): + return cls(impact=impact, occurrence=occurrence, duration=duration, name=name) + + @classmethod + def from_series( + cls, + impact: pd.Series, + *, + occurrence: int = 1, + duration: int = 1, + name: Optional[str] = None, + **kwarg, + ) -> Event: + """Create an event from an impact pd.Series. + + Args: + impact (pd.Series): A vector definition of the impact per (region,sector) + occurrence (int, optional): The ordinal of occurrence of the event (requires to be > 0). Defaults to 1. + duration (int, optional): The duration of the event (entire impact applied during this number of steps). Defaults to 1. + name (Optional[str], optional): A possible name for the event, for convenience. Defaults to None. + + Raises: + ValueError: Raised if impact is empty or contains negative values. + + Returns: + Event: an Event object (or one of its subclass). + """ + if impact.size == 0: + raise ValueError( + "Empty impact Series at init, did you not set the impact correctly ?" + ) + impact = impact[impact != 0] + if np.less_equal(impact, 0).any(): + logger.debug(f"Impact has negative values:\n{impact}\n{impact[impact<0]}") + raise ValueError("Impact has negative values") + return cls._instantiate( + impact=impact, + occurrence=occurrence, + duration=duration, + name=name, + **kwarg, ) - logger.debug( - f"self.impact_regional_distrib: {list(self.impact_regional_distrib)}" - ) - logger.debug(f"len aff_regions: {len(self.aff_regions)}") - self.impact_df.loc[self.aff_industries] = ( - self.impact_df.loc[self.aff_industries] * 1 / len(self.aff_regions) - ) - impact_sec_vec = np.array( - [ - 1 / len(self.aff_industries.to_series().loc[reg]) - for reg in self.aff_regions - ] + @classmethod + def from_dataframe( + cls, + impact: pd.DataFrame, + *, + occurrence: int = 1, + duration: int = 1, + name: Optional[str] = None, + **kwarg, + ) -> Event: + """Convenience function for DataFrames. See :meth:`~boario.event.Event.from_series` + + Raises: + ValueError: If impact cannot be squeezed to a Series + + Returns: + Event: an Event object (or one of its subclass). + """ + impact = impact.squeeze() + if not isinstance(impact, pd.Series): + raise ValueError("Could not squeeze impact dataframe to a serie.") + + return cls.from_series( + impact=impact, + occurrence=occurrence, + duration=duration, + name=name, + **kwarg, ) - self.impact_df.loc[self.aff_industries] = ( - self.impact_df.loc[self.aff_industries] * impact_sec_vec + + @classmethod + def distribute_impact_by_gva(cls, impact_vec: pd.Series) -> pd.Series: + """Distribute a vector of impact by the GVA of affected industries. + + Each values of the given impact are mutliplied by the share of the GVA + the industry has over the GVA of all affected industries. + + Args: + impact_vec (pd.Series): The impact values to be reweigthed. + Current use-case assumes all values are the total impact. + + Returns: + pd.Series: The impact where each value was multiplied by the share + of GVA of each affected industry (over total GVA affected). + """ + gva = cls.gva_df.loc[impact_vec.index] + gva = gva.transform(lambda x: x / sum(x)) + return impact_vec * gva + + @classmethod + def distribute_impact_equally(cls, impact_vec: pd.Series) -> pd.Series: + """Distribute an impact equally between all affected regions. + + Assume impact is given as a vector with all value being the + total impact to distribute. + + Args: + impact_vec (pd.Series): The impact to distribute. + + Returns: + pd.Series: The impact vector equally distributed among affected industries. + """ + dfg = impact_vec.groupby("region") + return dfg.transform(lambda x: x / (dfg.ngroups * x.count())) + + @classmethod + def from_scalar_industries( + cls, + impact: ScalarImpact, + *, + industries: IndustriesList, + impact_industries_distrib: Optional[npt.ArrayLike] = None, + gva_distrib: bool = False, + occurrence: int = 1, + duration: int = 1, + name: Optional[str] = None, + **kwarg, + ) -> Event: + """Creates an Event from a scalar and a list of industries affected. + + The scalar impact is distributed evenly by default. Otherwise it can be distributed proportionnaly to the GVA of + affected industries, or to a custom distribution. + + Args: + impact (ScalarImpact): The scalar impact. + industries (IndustriesList): The list of industries affected by the impact + impact_industries_distrib (Optional[npt.ArrayLike], optional): A vector of equal size to the list of industries, stating the + share of the impact each industry should receive. Defaults to None. + gva_distrib (bool, optional): A boolean stating if the impact should be distributed proportionnaly to GVA. Defaults to False. + occurrence (int, optional): The ordinal of occurrence of the event (requires to be > 0). Defaults to 1. + duration (int, optional): The duration of the event (entire impact applied during this number of steps). Defaults to 1. + name (Optional[str], optional): A possible name for the event, for convenience. Defaults to None. + + Raises: + ValueError: Raise if Impact is null, if len(industries) < 1 or if the sum of impact_industries_distrib differs from 1.0. + + Returns: + Event: An Event object or one of its subclass. + """ + if impact <= 0: + raise ValueError(f"Impact is null") + + if len(industries) < 1: + raise ValueError(f"Null sized affected industries ?") + + if isinstance(industries, list): + industries = pd.MultiIndex.from_tuples( + industries, names=["region", "sector"] + ) + + impact_vec = pd.Series(impact, dtype="float64", index=industries) + + if impact_industries_distrib: + if np.sum(impact_industries_distrib) != 1.0: + raise ValueError("Impact distribution doesn't sum up to 1.0") + else: + impact_vec *= impact_industries_distrib # type: ignore + + elif gva_distrib: + impact_vec = cls.distribute_impact_by_gva(impact_vec) + + else: + impact_vec = cls.distribute_impact_equally(impact_vec) + + return cls.from_series( + impact=impact_vec, + occurrence=occurrence, + duration=duration, + name=name, + **kwarg, ) - logger.debug( - f"Impact df after default distrib:\n {self.impact_df.loc[self.aff_industries]}" + + @classmethod + def from_scalar_regions_sectors( + cls, + impact: ScalarImpact, + *, + regions: RegionsList, + sectors: SectorsList, + impact_regional_distrib: Optional[npt.ArrayLike] = None, + impact_sectoral_distrib: Optional[Union[str, npt.ArrayLike]] = None, + occurrence: int = 1, + duration: int = 1, + name: Optional[str] = None, + **kwarg, + ) -> Event: + """Creates an Event from a scalar, a list of regions and a list of sectors affected. + + + + Args: + impact (ScalarImpact): The scalar impact. + regions (RegionsList): The list of regions affected. + sectors (SectorsList): The list of sectors affected in each region. + impact_regional_distrib (Optional[npt.ArrayLike], optional): A vector of equal size to the list of regions affected, stating the + share of the impact each industry should receive. Defaults to None. + impact_sectoral_distrib (Optional[Union[str, npt.ArrayLike]], optional): A vector of equal size to the list of sectors affected, stating the + share of the impact each industry should receive. Defaults to None. + occurrence (int, optional): The ordinal of occurrence of the event (requires to be > 0). Defaults to 1. + duration (int, optional): The duration of the event (entire impact applied during this number of steps). Defaults to 1. + name (Optional[str], optional): A possible name for the event, for convenience. Defaults to None. + + Raises: + ValueError: Raise if Impact is null, if len(regions) or len(sectors) < 1, + + Returns: + Event: An Event object or one of its subclass + """ + if impact <= 0: + raise ValueError(f"Impact is null") + + if isinstance(regions, str): + regions = [regions] + + if isinstance(sectors, str): + sectors = [sectors] + + _regions = pd.Index(regions, name="region") + _sectors = pd.Index(sectors, name="sector") + + if len(_regions) < 1: + raise ValueError(f"Null sized affected regions ?") + + if len(_sectors) < 1: + raise ValueError(f"Null sized affected sectors ?") + + if _sectors.duplicated().any(): + warnings.warn( + UserWarning( + "Multiple presence of the same sector in affected sectors. (Will remove duplicate)" + ) + ) + _sectors = _sectors.drop_duplicates() + + if _regions.duplicated().any(): + warnings.warn( + UserWarning( + "Multiple presence of the same region in affected region. (Will remove duplicate)" + ) + ) + _regions = _regions.drop_duplicates() + + industries = pd.MultiIndex.from_product( + [_regions, _sectors], names=["region", "sector"] ) - def _finish_init(self): - logger.debug("Finishing Event init") - self.impact_vector = self.impact_df.to_numpy() - self.total_impact = self.impact_vector.sum() - self.impact_industries_distrib = ( - self.impact_vector[self.impact_vector > 0] / self.total_impact + impact_vec = pd.Series(impact, dtype="float64", index=industries) + + assert isinstance(_regions, pd.Index) + if impact_regional_distrib is None: + regional_distrib = pd.Series(1.0 / len(_regions), index=_regions) + elif not isinstance(impact_regional_distrib, pd.Series): + impact_regional_distrib = pd.Series(impact_regional_distrib, index=_regions) + regional_distrib = pd.Series(0.0, index=_regions) + regional_distrib.loc[ + impact_regional_distrib.index + ] = impact_regional_distrib + else: + regional_distrib = pd.Series(0.0, index=_regions) + try: + regional_distrib.loc[ + impact_regional_distrib.index + ] = impact_regional_distrib + except KeyError: + regional_distrib.loc[_regions] = impact_regional_distrib.values + + assert isinstance(_sectors, pd.Index) + if impact_sectoral_distrib is None: + sectoral_distrib = pd.Series(1.0 / len(_sectors), index=_sectors) + elif ( + isinstance(impact_sectoral_distrib, str) + and impact_sectoral_distrib == "gdp" + ): + gva = cls.gva_df.loc[(_regions, _sectors)] + sectoral_distrib = gva.groupby("region").transform(lambda x: x / sum(x)) + elif not isinstance(impact_sectoral_distrib, pd.Series): + impact_sectoral_distrib = pd.Series(impact_sectoral_distrib, index=_sectors) + sectoral_distrib = pd.Series(0.0, index=_sectors) + sectoral_distrib.loc[ + impact_sectoral_distrib.index + ] = impact_sectoral_distrib + else: + sectoral_distrib = pd.Series(0.0, index=_sectors) + try: + sectoral_distrib.loc[ + impact_sectoral_distrib.index + ] = impact_sectoral_distrib + except KeyError: + sectoral_distrib.loc[_sectors] = impact_sectoral_distrib.values + + logger.debug(f"{sectoral_distrib}") + logger.debug(f"{regional_distrib}") + if isinstance(sectoral_distrib.index, pd.MultiIndex): + industries_distrib = sectoral_distrib * regional_distrib + else: + industries_distrib = pd.Series( + np.outer(regional_distrib.values, sectoral_distrib.values).flatten(), # type: ignore + index=pd.MultiIndex.from_product( + [regional_distrib.index, sectoral_distrib.index] + ), + ) + + impact_vec *= industries_distrib + + return cls.from_series( + impact=impact_vec, + occurrence=occurrence, + duration=duration, + name=name, + **kwarg, ) - self.impact_regional_distrib = ( - self.impact_df.loc[self.aff_industries].groupby("region").sum().values - / self.total_impact + + @property + def impact_df(self) -> pd.Series: + r"""A pandas Series with all possible industries as index, holding the impact vector of the event. The impact is defined for each sectors in each region.""" + return self._impact_df + + @impact_df.setter + def impact_df(self, value: pd.Series): + self._impact_df = pd.Series( + 0, + dtype="float64", + index=pd.MultiIndex.from_product( + [self.possible_regions, self.possible_sectors], + names=["region", "sector"], + ), ) + self._impact_df[value.index] = value + logger.debug("Sorting impact Series") + self._impact_df.sort_index(inplace=True) + self.aff_industries = self.impact_df.loc[self.impact_df > 0].index # type: ignore + self.impact_industries_distrib = self.impact_df.transform(lambda x: x / sum(x)) + self.total_impact = self.impact_df.sum() + self.impact_vector = self.impact_df.values @property def aff_industries(self) -> pd.MultiIndex: @@ -565,12 +570,12 @@ def aff_industries(self, index: pd.MultiIndex): self.aff_regions = regions self.aff_sectors = sectors logger.debug( - f"Setting _aff_industries. There are {len(index)} affected industries" + f"Setting _aff_industries. There are {np.size(index)} affected industries" ) self._aff_industries = index self._aff_industries_idx = np.array( [ - len(self.possible_sectors) * ri + si + np.size(self.possible_sectors) * ri + si for ri in self._aff_regions_idx for si in self._aff_sectors_idx ] @@ -584,9 +589,9 @@ def occurrence(self) -> int: @occurrence.setter def occurrence(self, value: int): - if not 0 <= value <= self.temporal_unit_range: + if not 0 < value <= self.temporal_unit_range: raise ValueError( - "Occurrence of event is not in the range of simulation steps : {} not in [0-{}]".format( + "Occurrence of event is not in the range of simulation steps (cannot be 0) : {} not in ]0-{}]".format( value, self.temporal_unit_range ) ) @@ -619,16 +624,16 @@ def aff_regions(self) -> pd.Index: return self._aff_regions @property - def aff_regions_idx(self) -> np.ndarray: + def aff_regions_idx(self) -> npt.NDArray: r"""The array of lexicographically ordered affected region indexes""" return self._aff_regions_idx @aff_regions.setter - def aff_regions(self, value: ArrayLike | str): + def aff_regions(self, value: npt.ArrayLike | str): if isinstance(value, str): value = [value] - value = pd.Index(value) + value = pd.Index(value) # type: ignore impossible_regions = np.setdiff1d(value, self.possible_regions) if impossible_regions.size > 0: raise ValueError( @@ -646,16 +651,16 @@ def aff_sectors(self) -> pd.Index: return self._aff_sectors @property - def aff_sectors_idx(self) -> np.ndarray: + def aff_sectors_idx(self) -> npt.NDArray: r"""The array of lexicographically ordered affected sectors indexes""" return self._aff_sectors_idx @aff_sectors.setter - def aff_sectors(self, value: ArrayLike | str): + def aff_sectors(self, value: npt.ArrayLike | str): if isinstance(value, str): value = [value] - value = pd.Index(value, name="sector") + value = pd.Index(value, name="sector") # type: ignore impossible_sectors = np.setdiff1d(value, self.possible_sectors) if impossible_sectors.size > 0: raise ValueError( @@ -667,76 +672,33 @@ def aff_sectors(self, value: ArrayLike | str): self._aff_sectors_idx = np.searchsorted(self.possible_sectors, value) @property - def impact_regional_distrib(self) -> np.ndarray: - r"""The array specifying how damages are distributed among affected regions""" + def impact_regional_distrib(self) -> pd.Series: + r"""The series specifying how damages are distributed among affected regions""" return self._impact_regional_distrib - @impact_regional_distrib.setter - def impact_regional_distrib(self, value: ArrayLike): - if self.aff_regions is None: - raise AttributeError("Affected regions attribute isn't set yet") - value = np.array(value) - if value.size != self.aff_regions.size: - raise ValueError( - "There are {} affected regions by the event and length of given damage distribution is {}".format( - self.aff_regions.size, value.size - ) - ) - s = value.sum() - if not math.isclose(s, 1): - raise ValueError( - "Damage distribution doesn't sum up to 1 but to {}, which is not valid".format( - s - ) - ) - self._impact_regional_distrib = value - - # @property - # def impact_sectoral_distrib(self) -> np.ndarray: - # r"""The array specifying how damages are distributed among affected sectors""" - # return self._impact_sectoral_distrib - - # @impact_sectoral_distrib.setter - # def impact_sectoral_distrib(self, value: ArrayLike): - # if self.aff_sectors is None: - # raise AttributeError("Affected sectors attribute isn't set yet") - # value = np.array(value) - # if value.size != self.aff_sectors.size: - # raise ValueError( - # "There are {} affected sectors by the event and length of given damage distribution is {}".format( - # self.aff_sectors.size, value.size - # ) - # ) - # s = value.sum() - # if not math.isclose(s, 1): - # raise ValueError( - # "Damage distribution doesn't sum up to 1 but to {}, which is not valid".format( - # s - # ) - # ) - # self._impact_sectoral_distrib = value - @property - def impact_sectoral_distrib_type(self) -> str: - r"""The type of damages distribution among sectors (currently only 'gdp')""" + def impact_industries_distrib(self) -> pd.Series: + r"""The series specifying how damages are distributed among affected industries (regions,sectors)""" - return self._impact_sectoral_distrib_type + return self._impact_industries_distrib - @impact_sectoral_distrib_type.setter - def impact_sectoral_distrib_type(self, value: str): - logger.debug(f"Setting _impact_sectoral_distrib_type to {value}") - self._impact_sectoral_distrib_type = value + @impact_industries_distrib.setter + def impact_industries_distrib(self, value: pd.Series): + self._impact_industries_distrib = value + self._impact_regional_distrib = self._impact_industries_distrib.groupby( + "region" + ).sum() def __repr__(self): # TODO: find ways to represent long lists - return f""" [Representation WIP] + return f""" Event( name = {self.name}, occur = {self.occurrence}, duration = {self.duration} - aff_regions = {self.aff_regions}, - aff_sectors = {self.aff_sectors}, + aff_regions = {self.aff_regions.to_list()}, + aff_sectors = {self.aff_sectors.to_list()}, ) """ @@ -744,159 +706,336 @@ def __repr__(self): class EventArbitraryProd(Event): def __init__( self, - impact: Impact, - aff_regions: Optional[RegionsList] = None, - aff_sectors: Optional[SectorsList] = None, - aff_industries: Optional[IndustriesList] = None, - impact_industries_distrib=None, - impact_regional_distrib=None, - impact_sectoral_distrib_type="equally shared", - impact_sectoral_distrib=None, - name="Unnamed", - occurrence=1, - duration=1, + *, + impact: pd.Series, + recovery_time: int = 1, + recovery_function: str = "linear", + name: str | None = "Unnamed", + occurrence: int = 1, + duration: int = 1, ) -> None: + if (impact > 1.0).any(): + raise ValueError( + "Impact is greater than 100% (1.) for at least an industry." + ) + super().__init__( - impact, - aff_regions, - aff_sectors, - aff_industries, - impact_industries_distrib, - impact_regional_distrib, - impact_sectoral_distrib_type, - impact_sectoral_distrib, - name, - occurrence, - duration, + impact=impact, + name=name, + occurrence=occurrence, + duration=duration, + ) + + self._prod_cap_delta_arbitrary_0 = ( + self.impact_vector.copy() + ) # np.zeros(shape=len(self.possible_sectors)) + self.prod_cap_delta_arbitrary = ( + self.impact_vector.copy() + ) # type: ignore # np.zeros(shape=len(self.possible_sectors)) + self.recovery_time = recovery_time + r"""The characteristic recovery duration after the event is over""" + + self.recovery_function = recovery_function + + logger.info("Initialized") + + @classmethod + def _instantiate( + cls, + impact: pd.Series, + *, + recovery_time: int = 1, + recovery_function: str = "linear", + occurrence: int = 1, + duration: int = 1, + name: Optional[str] = None, + ): + return cls( + impact=impact, + occurrence=occurrence, + duration=duration, + name=name, + recovery_time=recovery_time, + recovery_function=recovery_function, ) - self._prod_cap_delta_arbitrary = np.zeros(shape=len(self.possible_sectors)) - self._aff_sectors_idx = None - self._aff_sectors = None - raise NotImplementedError() @property - def prod_cap_delta_arbitrary(self) -> np.ndarray: - r"""The optional array storing arbitrary (as in not related to kapital destroyed) production capacity loss""" + def prod_cap_delta_arbitrary(self) -> npt.NDArray: + r"""The optional array storing arbitrary (as in not related to productive_capital destroyed) production capacity loss""" return self._prod_cap_delta_arbitrary @prod_cap_delta_arbitrary.setter - def prod_cap_delta_arbitrary(self, value: dict[str, float]): - if self.aff_regions is None: - raise AttributeError("Affected regions attribute isn't set yet") - aff_sectors = np.array(list(value.keys())) - aff_shares = np.array(list(value.values())) - impossible_sectors = np.setdiff1d(aff_sectors, self.possible_sectors) - if impossible_sectors.size > 0: - raise ValueError( - "These sectors are not in the model : {}".format(impossible_sectors) + def prod_cap_delta_arbitrary(self, value: npt.NDArray): + self._prod_cap_delta_arbitrary = value + + @property + def recoverable(self) -> bool: + r"""A boolean stating if an event can start recover""" + return self._recoverable + + @recoverable.setter + def recoverable(self, current_temporal_unit: int): + reb = (self.occurrence + self.duration) <= current_temporal_unit + if reb and not self.recoverable: + logger.info( + "Temporal_Unit : {} ~ Event named {} that occured at {} in {} has started recovering (arbitrary production capacity loss)".format( + current_temporal_unit, + self.name, + self.occurrence, + self.aff_regions.to_list(), + ) ) - self._aff_sectors = aff_sectors - self._aff_sectors_idx = np.searchsorted(self.possible_sectors, aff_sectors) - aff_industries_idx = np.array( - [ - len(self.possible_sectors) * ri + si - for ri in self.regions_idx - for si in self._aff_sectors_idx - ] - ) - self._prod_cap_delta_arbitrary[aff_industries_idx] = np.tile( - aff_shares, self._aff_regions.size + self._recoverable = reb + + @property + def recovery_function(self) -> Callable: + r"""The recovery function used for recovery (`Callable`)""" + return self._recovery_fun + + @recovery_function.setter + def recovery_function(self, r_fun: str | Callable | None): + if r_fun is None: + r_fun = "instant" + if self.recovery_time is None: + raise AttributeError( + "Impossible to set recovery function if no recovery time is given." + ) + if isinstance(r_fun, str): + if r_fun == "linear": + fun = linear_recovery + elif r_fun == "convexe": + fun = convexe_recovery_scaled + elif r_fun == "convexe noscale": + fun = convexe_recovery + elif r_fun == "concave": + fun = concave_recovery + else: + raise NotImplementedError( + "No implemented recovery function corresponding to {}".format(r_fun) + ) + elif callable(r_fun): + r_fun_argsspec = inspect.getfullargspec(r_fun) + r_fun_args = r_fun_argsspec.args + r_fun_argsspec.kwonlyargs + if not all( + args in r_fun_args + for args in [ + "init_impact_stock", + "elapsed_temporal_unit", + "recovery_time", + ] + ): + raise ValueError( + "Recovery function has to have at least the following keyword arguments: {}".format( + [ + "init_impact_stock", + "elapsed_temporal_unit", + "recovery_time", + ] + ) + ) + fun = r_fun + + else: + raise ValueError("Given recovery function is not a str or callable") + + r_fun_partial = partial( + fun, + init_impact_stock=self._prod_cap_delta_arbitrary_0, + recovery_time=self.recovery_time, ) + self._recovery_fun = r_fun_partial + + def recovery(self, current_temporal_unit: int): + elapsed = current_temporal_unit - (self.occurrence + self.duration) + if elapsed < 0: + raise RuntimeError("Trying to recover before event is over") + if self.recovery_function is None: + raise RuntimeError( + "Trying to recover event while recovery function isn't set yet" + ) + res = self.recovery_function(elapsed_temporal_unit=elapsed) + precision = int(math.log10(self.model_monetary_factor)) + 1 + res = np.around(res, precision) + if not np.any(res): + self.over = True + self._prod_cap_delta_arbitrary = res -class EventKapitalDestroyed(Event, metaclass=abc.ABCMeta): +class EventKapitalDestroyed(Event, ABC): def __init__( self, - impact: Impact, - aff_regions: Optional[RegionsList] = None, - aff_sectors: Optional[SectorsList] = None, - aff_industries: Optional[IndustriesList] = None, - impact_industries_distrib=None, - impact_regional_distrib=None, - impact_sectoral_distrib_type="equally shared", - impact_sectoral_distrib=None, - name="Unnamed", - occurrence=1, - duration=1, + *, + impact: pd.Series, + households_impact: Optional[Impact] = None, + name: str | None = "Unnamed", + occurrence: int = 1, + duration: int = 1, + event_monetary_factor: Optional[int] = None, ) -> None: + if event_monetary_factor is None: + logger.info( + f"No event monetary factor given. Assuming it is the same as the model ({self.model_monetary_factor})" + ) + self.event_monetary_factor = self.model_monetary_factor + r"""The monetary factor for the impact of the event (e.g. 10**6, 10**3, ...)""" + + else: + self.event_monetary_factor = event_monetary_factor + if self.event_monetary_factor != self.model_monetary_factor: + logger.warning( + f"Event monetary factor ({self.event_monetary_factor}) differs from model monetary factor ({self.model_monetary_factor}). Be careful to define your impact with the correct unit (ie in event monetary factor)." + ) + + if (impact < LOW_DEMAND_THRESH / self.event_monetary_factor).all(): + raise ValueError( + "Impact is too small to be considered by the model. Check you units perhaps ?" + ) + super().__init__( - impact, - aff_regions, - aff_sectors, - aff_industries, - impact_industries_distrib, - impact_regional_distrib, - impact_sectoral_distrib_type, - impact_sectoral_distrib, - name, - occurrence, - duration, + impact=impact, + name=name, + occurrence=occurrence, + duration=duration, ) - # The only thing we have to do is affecting/computing the regional_sectoral_kapital_destroyed - self.total_kapital_destroyed = self.total_impact - self.total_kapital_destroyed /= self.monetary_factor - self.remaining_kapital_destroyed = self.total_kapital_destroyed - self._regional_sectoral_kapital_destroyed_0 = ( - self.impact_vector / self.monetary_factor + + # The only thing we have to do is affecting/computing the regional_sectoral_productive_capital_destroyed + self.impact_df = self.impact_df * ( + self.event_monetary_factor / self.model_monetary_factor + ) + self.total_productive_capital_destroyed = self.total_impact + logger.info( + f"Total impact on productive capital is {self.total_productive_capital_destroyed} (in model unit)" + ) + self.remaining_productive_capital_destroyed = ( + self.total_productive_capital_destroyed ) - self.regional_sectoral_kapital_destroyed = ( - self.impact_vector / self.monetary_factor + self._regional_sectoral_productive_capital_destroyed_0 = ( + self.impact_vector.copy() + ) + self.regional_sectoral_productive_capital_destroyed = self.impact_vector.copy() # type: ignore + self.households_impact_df: pd.Series = pd.Series( + 0, + dtype="float64", + index=pd.MultiIndex.from_product( + [self.possible_regions, self.possible_final_demand_cat], + names=["region", "final_demand_cat"], + ), + ) + r"""A pandas Series with all possible (regions, final_demand_cat) as index, holding the households impacts vector of the event. The impact is defined for each region and each final demand category.""" + + if households_impact: + try: + rebuilding_demand_idx = self.possible_final_demand_cat[ + self.possible_final_demand_cat.str.match( + REBUILDING_FINALDEMAND_CAT_REGEX + ) + ] # .values[0] + if len(rebuilding_demand_idx) > 1: + raise ValueError( + f"The rebuilding demand index ({rebuilding_demand_idx}) contains multiple values which is a problem. Contact the dev to update the matching regexp." + ) + + except IndexError: + logger.warning( + f"No final demand category matched common rebuilding final demand category, hence we will put it in the first available ({self.possible_final_demand_cat[0]})." + ) + rebuilding_demand_idx = self.possible_final_demand_cat[0] + if isinstance(households_impact, pd.Series): + logger.debug("Given household impact is a pandas Series") + self.households_impact_df.loc[ + households_impact.index, rebuilding_demand_idx + ] = households_impact # type: ignore + elif isinstance(households_impact, dict): + logger.debug("Given household impact is a dict") + households_impact = pd.Series(households_impact, dtype="float64") + self.households_impact_df.loc[ + households_impact.index, rebuilding_demand_idx + ] = households_impact # type: ignore + elif isinstance(households_impact, pd.DataFrame): + logger.debug("Given household impact is a dataframe") + households_impact = households_impact.squeeze() + if not isinstance(households_impact, pd.Series): + raise ValueError( + "The given households_impact DataFrame is not a Series after being squeezed" + ) + self.households_impact_df.loc[ + households_impact.index, rebuilding_demand_idx + ] = households_impact # type: ignore + elif isinstance(households_impact, (list, np.ndarray)): + if np.size(households_impact) != np.size(self.aff_regions): + raise ValueError( + f"Length mismatch between households_impact ({np.size(households_impact)} and affected regions ({np.size(self.aff_regions)}))" + ) + else: + self.households_impact_df.loc[ + self.aff_regions, rebuilding_demand_idx + ] = households_impact # type: ignore + elif isinstance(households_impact, (float, int)): + if self.impact_regional_distrib is not None: + logger.warning( + f"households impact ({households_impact}) given as scalar, distributing among region following `impact_regional_distrib` ({self.impact_regional_distrib}) to {self.aff_regions}" + ) + logger.debug(f"{rebuilding_demand_idx}") + self.households_impact_df.loc[ + self.aff_regions, rebuilding_demand_idx + ] = ( + households_impact * self.impact_regional_distrib + ) # type: ignore + self.households_impact_df *= ( + self.event_monetary_factor / self.model_monetary_factor + ) + logger.info( + f"Total impact on households is {self.households_impact_df.sum()} (in model unit)" ) @property - def regional_sectoral_kapital_destroyed(self) -> np.ndarray: - r"""The (optional) array of kapital destroyed per industry (ie region x sector)""" - return self._regional_sectoral_kapital_destroyed + def regional_sectoral_productive_capital_destroyed(self) -> npt.NDArray: + r"""The (optional) array of productive_capital destroyed per industry (ie region x sector)""" + return self._regional_sectoral_productive_capital_destroyed - @regional_sectoral_kapital_destroyed.setter - def regional_sectoral_kapital_destroyed(self, value: ArrayLike): + @regional_sectoral_productive_capital_destroyed.setter + def regional_sectoral_productive_capital_destroyed(self, value: npt.ArrayLike): value = np.array(value) + logger.debug( + f"Changing regional_sectoral_productive_capital_destroyed with {value}\n Sum is {value.sum()}" + ) if value.shape != self.x_shape: raise ValueError( - "Incorrect shape give for regional_sectoral_kapital_destroyed: {} given and {} expected".format( + "Incorrect shape give for regional_sectoral_productive_capital_destroyed: {} given and {} expected".format( value.shape, self.x_shape ) ) - self._regional_sectoral_kapital_destroyed = value + self._regional_sectoral_productive_capital_destroyed = value class EventKapitalRebuild(EventKapitalDestroyed): def __init__( self, - impact: Impact, - rebuilding_sectors: Union[dict[str, float], pd.Series], - rebuild_tau=60, - aff_regions: Optional[RegionsList] = None, - aff_sectors: Optional[SectorsList] = None, - aff_industries: Optional[IndustriesList] = None, - impact_industries_distrib=None, - impact_regional_distrib=None, - impact_sectoral_distrib_type="equally shared", - impact_sectoral_distrib=None, - name="Unnamed", - occurrence=1, - duration=1, - rebuilding_factor=None, + *, + impact: pd.Series, + households_impact: Impact | None = None, + name: str | None = "Unnamed", + occurrence: int = 1, + duration: int = 1, + event_monetary_factor: Optional[int] = None, + rebuild_tau: int, + rebuilding_sectors: dict[str, float] | pd.Series, + rebuilding_factor: float = 1.0, ) -> None: super().__init__( - impact, - aff_regions, - aff_sectors, - aff_industries, - impact_industries_distrib, - impact_regional_distrib, - impact_sectoral_distrib_type, - impact_sectoral_distrib, - name, - occurrence, - duration, + impact=impact, + households_impact=households_impact, + name=name, + occurrence=occurrence, + duration=duration, + event_monetary_factor=event_monetary_factor, ) - - self._rebuildable = None + self._rebuildable = False self.rebuild_tau = rebuild_tau self.rebuilding_sectors = rebuilding_sectors + self.rebuilding_factor = rebuilding_factor - rebuilding_demand = np.zeros(shape=self.z_shape) + industrial_rebuilding_demand = np.zeros(shape=self.z_shape) tmp = np.zeros(self.z_shape, dtype="float") mask = np.ix_( np.union1d( @@ -904,23 +1043,38 @@ def __init__( ), self._aff_industries_idx, ) - rebuilding_demand = np.outer( - self.rebuilding_sectors_shares, self.regional_sectoral_kapital_destroyed + industrial_rebuilding_demand = np.outer( + self.rebuilding_sectors_shares, + self.regional_sectoral_productive_capital_destroyed, + ) + industrial_rebuilding_demand = industrial_rebuilding_demand * rebuilding_factor + tmp[mask] = self.Z_distrib[mask] * industrial_rebuilding_demand[mask] + + precision = int(math.log10(self.model_monetary_factor)) + 1 + np.around(tmp, precision, tmp) + + self.rebuilding_demand_indus = tmp + + households_rebuilding_demand = np.zeros(shape=self.y_shape) + tmp = np.zeros(self.y_shape, dtype="float") + mask = np.ix_( # type: ignore + np.union1d( # type: ignore + self._rebuilding_industries_RoW_idx, self._rebuilding_industries_idx + ), + list(range(self.y_shape[1])), + ) + + households_rebuilding_demand = np.outer( + self.rebuilding_sectors_shares, self.households_impact_df.to_numpy() ) - rebuilding_demand = rebuilding_demand * rebuilding_factor + households_rebuilding_demand = households_rebuilding_demand * rebuilding_factor + tmp[mask] = self.Y_distrib[mask] * households_rebuilding_demand[mask] + np.around(tmp, precision, tmp) + self.rebuilding_demand_house = tmp logger.debug( - f"kapital destroyed vec is {self.regional_sectoral_kapital_destroyed}" + f"house rebuilding demand is: {pd.DataFrame(self.rebuilding_demand_house, index=pd.MultiIndex.from_product([self.possible_regions,self.possible_sectors]))}" ) - logger.debug(f"rebuilding sectors shares are {self.rebuilding_sectors_shares}") - logger.debug(f"Z shape is {self.z_shape}") - logger.debug(f"Z_distrib[mask] has shape {self.Z_distrib[mask].shape}") - logger.debug(f"reb_demand: {rebuilding_demand}") - logger.debug(f"reb_demand: {rebuilding_demand.shape}") - # reb_tiled = np.tile(rebuilding_demand, (len(self.possible_regions), 1)) - # reb_tiled = reb_tiled[mask] - tmp[mask] = self.Z_distrib[mask] * rebuilding_demand[mask] - self.rebuilding_demand_indus = tmp - self.rebuilding_demand_house = np.zeros(shape=self.y_shape) + self.event_dict["rebuilding_sectors"] = { sec: share for sec, share in zip( @@ -928,18 +1082,44 @@ def __init__( ) } + @classmethod + def _instantiate( + cls, + impact: pd.Series, + *, + households_impact: Optional[Impact] = None, + occurrence: int = 1, + duration: int = 1, + name: str | None = None, + event_monetary_factor: Optional[int] = None, + rebuild_tau: int, + rebuilding_sectors: dict[str, float] | pd.Series, + rebuilding_factor: float = 1.0, + ): + return cls( + impact=impact, + households_impact=households_impact, + occurrence=occurrence, + duration=duration, + name=name, + event_monetary_factor=event_monetary_factor, + rebuild_tau=rebuild_tau, + rebuilding_sectors=rebuilding_sectors, + rebuilding_factor=rebuilding_factor, + ) + @property def rebuilding_sectors(self) -> pd.Index: r"""The (optional) array of rebuilding sectors""" return self._rebuilding_sectors @property - def rebuilding_sectors_idx(self) -> np.ndarray: + def rebuilding_sectors_idx(self) -> npt.NDArray: r"""The (optional) array of indexes of the rebuilding sectors (in lexicographic order)""" return self._rebuilding_sectors_idx @property - def rebuilding_sectors_shares(self) -> np.ndarray: + def rebuilding_sectors_shares(self) -> npt.NDArray: r"""The array specifying how rebuilding demand is distributed along the rebuilding sectors""" return self._rebuilding_sectors_shares @@ -949,6 +1129,7 @@ def rebuilding_sectors(self, value: dict[str, float] | pd.Series): reb_sectors = pd.Series(value) else: reb_sectors = value + assert np.isclose(reb_sectors.sum(), 1.0) impossible_sectors = np.setdiff1d(reb_sectors.index, self.possible_sectors) if impossible_sectors.size > 0: raise ValueError( @@ -963,36 +1144,36 @@ def rebuilding_sectors(self, value: dict[str, float] | pd.Series): self._rebuilding_sectors_shares = np.zeros(self.x_shape) self._rebuilding_industries_idx = np.array( [ - len(self.possible_sectors) * ri + si + np.size(self.possible_sectors) * ri + si for ri in self._aff_regions_idx for si in self._rebuilding_sectors_idx ] ) self._rebuilding_industries_RoW_idx = np.array( [ - len(self.possible_sectors) * ri + si - for ri in range(len(self.possible_regions)) + np.size(self.possible_sectors) * ri + si + for ri in range(np.size(self.possible_regions)) if ri not in self._aff_regions_idx for si in self._rebuilding_sectors_idx ] ) self._rebuilding_sectors_shares[self._rebuilding_industries_idx] = np.tile( - np.array(reb_sectors.values), len(self.aff_regions) + np.array(reb_sectors.values), np.size(self.aff_regions) ) self._rebuilding_sectors_shares[ self._rebuilding_industries_RoW_idx ] = np.tile( np.array(reb_sectors.values), - (len(self.possible_regions) - len(self.aff_regions)), + (np.size(self.possible_regions) - np.size(self.aff_regions)), ) @property - def rebuilding_demand_house(self) -> np.ndarray: + def rebuilding_demand_house(self) -> npt.NDArray: r"""The optional array of rebuilding demand from households""" return self._rebuilding_demand_house @rebuilding_demand_house.setter - def rebuilding_demand_house(self, value: ArrayLike): + def rebuilding_demand_house(self, value: npt.ArrayLike): value = np.array(value) if value.shape != self.y_shape: raise ValueError( @@ -1000,15 +1181,16 @@ def rebuilding_demand_house(self, value: ArrayLike): value.shape, self.y_shape ) ) + value[value < LOW_DEMAND_THRESH / self.model_monetary_factor] = 0.0 self._rebuilding_demand_house = value @property - def rebuilding_demand_indus(self) -> np.ndarray: - r"""The optional array of rebuilding demand from industries (to rebuild their kapital)""" + def rebuilding_demand_indus(self) -> npt.NDArray: + r"""The optional array of rebuilding demand from industries (to rebuild their productive_capital)""" return self._rebuilding_demand_indus @rebuilding_demand_indus.setter - def rebuilding_demand_indus(self, value: ArrayLike): + def rebuilding_demand_indus(self, value: npt.ArrayLike): value = np.array(value) if value.shape != self.z_shape: raise ValueError( @@ -1016,9 +1198,12 @@ def rebuilding_demand_indus(self, value: ArrayLike): value.shape, self.z_shape ) ) + value[value < LOW_DEMAND_THRESH / self.model_monetary_factor] = 0.0 self._rebuilding_demand_indus = value - # Also update kapital destroyed - self._regional_sectoral_kapital_destroyed = value.sum(axis=0) + # Also update productive_capital destroyed + self.regional_sectoral_productive_capital_destroyed = value.sum(axis=0) * ( + 1 / self.rebuilding_factor + ) @property def rebuildable(self) -> bool: @@ -1034,8 +1219,8 @@ def rebuildable(self, current_temporal_unit: int): current_temporal_unit, self.name, self.occurrence, - self._aff_regions.values, - self.total_kapital_destroyed, + self.aff_regions.to_list(), + self.total_productive_capital_destroyed, ) ) self._rebuildable = reb @@ -1044,37 +1229,52 @@ def rebuildable(self, current_temporal_unit: int): class EventKapitalRecover(EventKapitalDestroyed): def __init__( self, - impact: Impact, + *, + impact: pd.Series, recovery_time: int, recovery_function: str = "linear", - aff_regions: Optional[RegionsList] = None, - aff_sectors: Optional[SectorsList] = None, - aff_industries: Optional[IndustriesList] = None, - impact_industries_distrib=None, - impact_regional_distrib=None, - impact_sectoral_distrib_type="equally shared", - impact_sectoral_distrib=None, - name="Unnamed", + households_impact: Optional[Impact] = None, + name: str | None = "Unnamed", occurrence=1, duration=1, + event_monetary_factor=None, ) -> None: super().__init__( - impact, - aff_regions, - aff_sectors, - aff_industries, - impact_industries_distrib, - impact_regional_distrib, - impact_sectoral_distrib_type, - impact_sectoral_distrib, - name, - occurrence, - duration, + impact=impact, + households_impact=households_impact, + name=name, + occurrence=occurrence, + duration=duration, + event_monetary_factor=event_monetary_factor, ) self.recovery_time = recovery_time self.recovery_function = recovery_function self._recoverable = False + @classmethod + def _instantiate( + cls, + impact: pd.Series, + *, + households_impact: Optional[Impact] = None, + occurrence: int = 1, + duration: int = 1, + name: str | None = None, + event_monetary_factor: Optional[int] = None, + recovery_time: int, + recovery_function: str = "linear", + ): + return cls( + impact=impact, + households_impact=households_impact, + occurrence=occurrence, + duration=duration, + name=name, + event_monetary_factor=event_monetary_factor, + recovery_time=recovery_time, + recovery_function=recovery_function, + ) + @property def recoverable(self) -> bool: r"""A boolean stating if an event can start recover""" @@ -1089,8 +1289,8 @@ def recoverable(self, current_temporal_unit: int): current_temporal_unit, self.name, self.occurrence, - self._aff_regions, - self.total_kapital_destroyed, + self.aff_regions.to_list(), + self.total_productive_capital_destroyed, ) ) self._recoverable = reb @@ -1127,7 +1327,7 @@ def recovery_function(self, r_fun: str | Callable | None): if not all( args in r_fun_args for args in [ - "init_kapital_destroyed", + "init_impact_stock", "elapsed_temporal_unit", "recovery_time", ] @@ -1135,7 +1335,7 @@ def recovery_function(self, r_fun: str | Callable | None): raise ValueError( "Recovery function has to have at least the following keyword arguments: {}".format( [ - "init_kapital_destroyed", + "init_impact_stock", "elapsed_temporal_unit", "recovery_time", ] @@ -1148,7 +1348,7 @@ def recovery_function(self, r_fun: str | Callable | None): r_fun_partial = partial( fun, - init_kapital_destroyed=self._regional_sectoral_kapital_destroyed_0, + init_impact_stock=self._regional_sectoral_productive_capital_destroyed_0, recovery_time=self.recovery_time, ) self._recovery_fun = r_fun_partial @@ -1162,8 +1362,9 @@ def recovery(self, current_temporal_unit: int): "Trying to recover event while recovery function isn't set yet" ) res = self.recovery_function(elapsed_temporal_unit=elapsed) - precision = int(math.log10(self.monetary_factor)) + 1 + precision = int(math.log10(self.model_monetary_factor)) + 1 res = np.around(res, precision) + res[res < 0] = 0.0 if not np.any(res): self.over = True - self.regional_sectoral_kapital_destroyed = res + self.regional_sectoral_productive_capital_destroyed = res diff --git a/boario/extended_models.py b/boario/extended_models.py index c12a30f..579367f 100644 --- a/boario/extended_models.py +++ b/boario/extended_models.py @@ -15,8 +15,9 @@ # along with this program. If not, see . from __future__ import annotations -from typing import Dict +from typing import Dict, Optional import numpy as np +import pandas as pd from boario import logger from boario.model_base import * from boario.event import * @@ -37,9 +38,16 @@ def __init__( rebuild_tau=60, main_inv_dur=90, monetary_factor=10**6, - psi_param=0.90, + temporal_units_by_step: int = 1, + iotable_year_to_temporal_unit_factor: int = 365, + infinite_inventories_sect: Optional[list] = None, + inventory_dict: Optional[dict] = None, + productive_capital_vector: Optional[ + pd.Series | np.ndarray | pd.DataFrame + ] = None, + productive_capital_to_VA_dict: Optional[dict] = None, + psi_param=0.80, inventory_restoration_tau: int | Dict[str, int] = 60, - **kwargs, ) -> None: """An ARIO3 model with some additional features @@ -55,7 +63,12 @@ def __init__( rebuild_tau=rebuild_tau, main_inv_dur=main_inv_dur, monetary_factor=monetary_factor, - **kwargs, + temporal_units_by_step=temporal_units_by_step, + iotable_year_to_temporal_unit_factor=iotable_year_to_temporal_unit_factor, + infinite_inventories_sect=infinite_inventories_sect, + inventory_dict=inventory_dict, + productive_capital_vector=productive_capital_vector, + productive_capital_to_VA_dict=productive_capital_to_VA_dict, ) logger.debug("Model is an ARIOPsiModel") @@ -66,10 +79,14 @@ def __init__( elif isinstance(psi_param, float): self.psi = psi_param + elif isinstance(psi_param, int): + self.psi = float(psi_param) else: raise ValueError( "'psi_param' parameter is neither a str rep of a float or a float" ) + if self.psi > 1.0: + raise ValueError("'psi_param' parameter must be less or equal than 1.") if isinstance(inventory_restoration_tau, int): restoration_tau = [ @@ -85,11 +102,19 @@ def __init__( ) for _, value in inventory_restoration_tau.items(): - if not isinstance(value, int): + if isinstance(value, float): + if not value.is_integer(): + raise ValueError( + f"Invalid value: {value} ({type(value)}) in inventory_restoration_tau, values should be integer: {inventory_restoration_tau}" + ) + elif not isinstance(value, int): raise ValueError( - "Invalid value in inventory_restoration_tau, values should be integer." + f"Invalid value: {value} ({type(value)}) in inventory_restoration_tau, values should be integer: {inventory_restoration_tau}" ) + inventory_restoration_tau = { + key: int(val) for key, val in inventory_restoration_tau.items() + } inventory_restoration_tau = dict(sorted(inventory_restoration_tau.items())) restoration_tau = [ (self.n_temporal_units_by_step / v) if v >= INV_THRESHOLD else v @@ -100,7 +125,7 @@ def __init__( f"Invalid inventory_restoration_tau: expected dict or int got {type(inventory_restoration_tau)}" ) self.restoration_tau = np.array(restoration_tau) - """numpy.ndarray of int: Array of size :math:`n` setting for each inputs its characteristic restoration time :math:`\tau_{\textrm{INV}}` in ``n_temporal_units_by_step``. (see :ref:`boario-math`).""" + r"""numpy.ndarray of int: Array of size :math:`n` setting for each inputs its characteristic restoration time :math:`\tau_{\textrm{INV}}` in ``n_temporal_units_by_step``. (see :ref:`boario-math`).""" ################################################################# @property @@ -117,50 +142,6 @@ def calc_production(self, current_temporal_unit: int) -> np.ndarray: r"""Computes and updates actual production. The difference with :class:`ARIOBaseModel` is in the way inventory constraints are computed. See :ref:`boario-math-prod`. - 1. Computes ``production_opt`` and ``inventory_constraints`` as : - - .. math:: - :nowrap: - - \begin{alignat*}{4} - \iox^{\textrm{Opt}}(t) &= (x^{\textrm{Opt}}_{f}(t))_{f \in \firmsset} &&= \left ( \min \left ( d^{\textrm{Tot}}_{f}(t), x^{\textrm{Cap}}_{f}(t) \right ) \right )_{f \in \firmsset} && \text{Optimal production}\\ - \ioinv^{\textrm{Cons}}(t) &= (\omega^{\textrm{Cons},f}_p(t))_{\substack{p \in \sectorsset\\f \in \firmsset}} &&= - \begin{bmatrix} - s^{1}_1 & \hdots & s^{p}_1 \\ - \vdots & \ddots & \vdots\\ - s^1_n & \hdots & s^{p}_n - \end{bmatrix} - \odot \begin{bmatrix} \iox^{\textrm{Opt}}(t)\\ - \vdots\\ - \iox^{\textrm{Opt}}(t) \end{bmatrix} \odot \ioa^{\sectorsset} && \text{Inventory constraints} \\ - &&&= \begin{bmatrix} - s^{1}_1 x^{\textrm{Opt}}_{1}(t) a_{11} & \hdots & s^{p}_1 x^{\textrm{Opt}}_{p}(t) a_{1p}\\ - \vdots & \ddots & \vdots\\ - s^1_n x^{\textrm{Opt}}_{1}(t) a_{n1} & \hdots & s^{p}_n x^{\textrm{Opt}}_{p}(t) a_{np} - \end{bmatrix} - \cdot \psi && \\ - &&&= \begin{bmatrix} - \tau^{1}_1 x^{\textrm{Opt}}_{1}(t) a_{11} & \hdots & \tau^{p}_1 x^{\textrm{Opt}}_{p}(t) a_{1p}\\ - \vdots & \ddots & \vdots\\ - \tau^1_n x^{\textrm{Opt}}_{1}(t) a_{n1} & \hdots & \tau^{p}_n x^{\textrm{Opt}}_{p}(t) a_{np} - \end{bmatrix} && \\ - \end{alignat*} - - 2. If stocks do not meet ``inventory_constraints`` for any inputs, then decrease production accordingly : - - .. math:: - :nowrap: - - \begin{alignat*}{4} - \iox^{a}(t) &= (x^{a}_{f}(t))_{f \in \firmsset} &&= \left \{ \begin{aligned} - & x^{\textrm{Opt}}_{f}(t) & \text{if $\omega_{p}^f(t) \geq \omega^{\textrm{Cons},f}_p(t)$}\\ - & x^{\textrm{Opt}}_{f}(t) \cdot \min_{p \in \sectorsset} \left ( \frac{\omega^s_{p}(t)}{\omega^{\textrm{Cons,f}}_p(t)} \right ) & \text{if $\omega_{p}^f(t) < \omega^{\textrm{Cons},f}_p(t)$} - \end{aligned} \right. \quad && - \end{alignat*} - - Also warns in logs if such shortages happen. - - Parameters ---------- current_temporal_unit : int diff --git a/boario/indicators.py b/boario/indicators.py index a4899e9..373724b 100644 --- a/boario/indicators.py +++ b/boario/indicators.py @@ -79,7 +79,7 @@ class Indicators(object): "overproduction", "rebuild_demand", "rebuild_prod", - "kapital_to_recover", + "productive_capital_to_recover", ] params_list = ["simulated_params", "simulated_events"] @@ -91,7 +91,7 @@ def __init__( params, regions, sectors, - kapital, + productive_capital, prod, prodmax, overprod, @@ -123,7 +123,9 @@ def __init__( self.n_rows = n_temporal_units_to_sim self.n_temporal_units_by_step = n_temporal_units_by_step - self.kapital_to_recover_df = create_df(kapital, regions, sectors) + self.productive_capital_to_recover_df = create_df( + productive_capital, regions, sectors + ) self.production_realised_df = create_df(prod, regions, sectors) self.production_capacity_df = create_df(prodmax, regions, sectors) self.overproduction_df = create_df(overprod, regions, sectors) @@ -187,15 +189,16 @@ def __init__( .reset_index() ) - self.limiting_inputs_df = pd.DataFrame( - limiting_stocks.reshape(n_temporal_units_to_sim * len(sectors), -1), - index=pd.MultiIndex.from_product( - [steps, sectors], names=["step", "stock of"] - ), - columns=pd.MultiIndex.from_product( - [regions, sectors], names=["region", "sector"] - ), - ) + self.limiting_inputs_df = limiting_stocks + # pd.DataFrame( + # limiting_stocks.reshape(n_temporal_units_to_sim * len(sectors), -1), + # index=pd.MultiIndex.from_product( + # [steps, sectors], names=["step", "stock of"] + # ), + # columns=pd.MultiIndex.from_product( + # [regions, sectors], names=["region", "sector"] + # ), + # ) self.aff_regions = [] self.aff_sectors = [] for e in events: @@ -259,9 +262,9 @@ def __init__( def from_simulation( cls, sim: Simulation, + results_storage: Optional[str | pathlib.Path] = None, stocks_treatment: bool = False, include_crash: bool = False, - results_storage: Optional[str | pathlib.Path] = None, ) -> Indicators: if isinstance(sim.model, ARIOPsiModel): psi = sim.model.psi @@ -270,27 +273,27 @@ def from_simulation( psi = None inventory_restoration_tau = None - prod = getattr(sim, "production_evolution") - kapital = getattr(sim, "regional_sectoral_kapital_destroyed_evolution") - prodmax = getattr(sim, "production_cap_evolution") - overprod = getattr(sim, "overproduction_evolution") - final_demand = getattr(sim, "final_demand_evolution") - io_demand = getattr(sim, "io_demand_evolution") - r_demand = getattr(sim, "rebuild_demand_evolution") - r_prod = getattr(sim, "rebuild_production_evolution") - fd_unmet = getattr(sim, "final_demand_unmet_evolution") + prod = getattr(sim, "production_realised") + productive_capital = getattr(sim, "productive_capital_to_recover") + prodmax = getattr(sim, "production_capacity") + overprod = getattr(sim, "overproduction") + final_demand = getattr(sim, "final_demand") + io_demand = getattr(sim, "intermediate_demand") + r_demand = getattr(sim, "rebuild_demand") + r_prod = getattr(sim, "rebuild_prod") + fd_unmet = getattr(sim, "final_demand_unmet") if stocks_treatment: - stocks = getattr(sim, "inputs_evolution") + stocks = getattr(sim, "inputs_stocks") else: stocks = None - limiting_stocks = getattr(sim, "limiting_inputs_evolution") + limiting_stocks = getattr(sim, "limiting_inputs") return cls( has_crashed=sim.has_crashed, params=sim.params_dict, regions=sim.model.regions, sectors=sim.model.sectors, - kapital=kapital, + productive_capital=productive_capital, prod=prod, prodmax=prodmax, overprod=overprod, @@ -373,8 +376,8 @@ def from_folder( dtype="float64", shape=(t, indexes["n_industries"]), ) - kapital = np.memmap( - records_path / "kapital_to_recover", + productive_capital = np.memmap( + records_path / "productive_capital_to_recover", mode="r+", dtype="float64", shape=(t, indexes["n_industries"]), @@ -450,7 +453,7 @@ def from_folder( params=params, regions=regions, sectors=sectors, - kapital=kapital, + productive_capital=productive_capital, prod=prod, prodmax=prodmax, overprod=overprod, @@ -706,8 +709,8 @@ def save_dfs(self): self.production_realised_df.to_parquet( self.parquets_path / "production_realised_df.parquet" ) - self.kapital_to_recover_df.to_parquet( - self.parquets_path / "kapital_to_recover_df.parquet" + self.productive_capital_to_recover_df.to_parquet( + self.parquets_path / "productive_capital_to_recover_df.parquet" ) self.production_capacity_df.to_parquet( self.parquets_path / "production_capacity_df.parquet" diff --git a/boario/model_base.py b/boario/model_base.py index 3d24c8f..1b3c63d 100644 --- a/boario/model_base.py +++ b/boario/model_base.py @@ -20,10 +20,17 @@ import typing from typing import Optional import numpy as np +import numpy.typing as npt import pandas as pd -from boario import logger, warn_once -from boario.event import * from pymrio.core.mriosystem import IOSystem + +from boario import logger, warn_once +from boario.event import ( + Event, + EventKapitalDestroyed, + EventKapitalRebuild, + EventArbitraryProd, +) from boario.utils.misc import lexico_reindex @@ -37,7 +44,9 @@ INV_THRESHOLD = 0 # 20 #temporal_units -TECHNOLOGY_THRESHOLD = 0.00001 +TECHNOLOGY_THRESHOLD = ( + 0.00001 # Do not care about input if producing requires less than this value +) VALUE_ADDED_NAMES = [ "VA", @@ -84,8 +93,10 @@ def __init__( iotable_year_to_temporal_unit_factor: int = 365, infinite_inventories_sect: Optional[list] = None, inventory_dict: Optional[dict] = None, - kapital_vector: Optional[pd.Series | np.ndarray | pd.DataFrame] = None, - kapital_to_VA_dict: Optional[dict] = None, + productive_capital_vector: Optional[ + pd.Series | npt.NDArray | pd.DataFrame + ] = None, + productive_capital_to_VA_dict: Optional[dict] = None, ) -> None: r"""The core of an ARIO model. Handles the different arrays containing the mrio tables. @@ -94,7 +105,7 @@ def __init__( An ARIOBaseModel instance based is build on the given IOSystem. Most default parameters are the same as in :cite:`2013:hallegatte` and - default order mechanism comes from :cite:`2020:guan`. By default each + default order mechanisms comes from :cite:`2020:guan`. By default each industry capital stock is 4 times its value added (:cite:`2008:hallegatte`). Parameters @@ -111,30 +122,30 @@ def __init__( alpha_tau : int, default 365 Characteristic time of overproduction :math:`\tau_{\alpha}` in ``n_temporal_units_by_step`` (default should be 365 days). rebuild_tau : int, default 60 - Rebuilding characteristic time :math:`\tau_{\textrm{REBUILD}}` (see :ref:`boario-math`). + Rebuilding characteristic time :math:`\tau_{\textrm{REBUILD}}` (see :ref:`boario-math`). Overwritten by per event value if it exists. main_inv_dur : int, default 90 The default numbers of days for inputs inventory to use if it is not defined for an input. monetary_factor : int, default 10**6 Monetary unit factor (i.e. if the tables unit is 10^6 € instead of 1 €, it should be set to 10^6). temporal_units_by_step: int, default 1 - The number of temporal_units between each step. (Current version of the model was not tested with values other than `1`). + The number of temporal_units between each step. (Current version of the model showed inconsistencies with values other than `1`). iotable_year_to_temporal_unit_factor : int, default 365 The (inverse of the) factor by which MRIO data should be multiplied in order to get "per temporal unit value", assuming IO data is yearly. infinite_inventories_sect: list, optional List of inputs (sector) considered never constraining for production. inventory_dict: dict, optional Dictionary in the form input:initial_inventory_size, (where input is a sector, and initial_inventory_size is in "temporal_unit" (defaults to a day)) - kapital_vector: ArrayLike, optional + productive_capital_vector: ArrayLike, optional Array of capital per industry if you need to give it exogenously. - kapital_to_VA_dict: dict, optional + productive_capital_to_VA_dict: dict, optional Dictionary in the form sector:ratio. Where ratio is used to estimate capital stock based on the value added of the sector. Notes ----- - It is recommended to use ``kapital_to_VA_dict`` if you have a more precise estimation of - the ratio of (Capital Stock / Value Added) per sectors. - You may also feed in directly a ``kapital_vector`` if you did your estimation before-hand. + It is recommended to use ``productive_capital_to_VA_dict`` if you have a more precise estimation of + the ratio of (Capital Stock / Value Added) per sectors than the default 4/1 ratio. + You may also feed in directly a ``productive_capital_vector`` if you did your estimation before-hand. (This is especially useful if you have events based of an exposure layer for instance) Regarding inventories, they default to 90 days for all inputs (ie sectors). @@ -142,18 +153,25 @@ def __init__( in ``infinite_inventories_sect`` or directly feed in a dictionary of the inventory duration for each input. - - Examples - -------- - FIXME: Add docs. - - """ logger.debug("Initiating new ARIOBaseModel instance") super().__init__() ################ Parameters variables ####################### - logger.info("IO system metadata :\n{}".format(str(pym_mrio.meta))) + try: + logger.info( + "IO system metadata :\n{}\n{}\n{}\n{}".format( + str(pym_mrio.meta.description), + str(pym_mrio.meta.name), + str(pym_mrio.meta.system), + str(pym_mrio.meta.version), + ) + ) + except AttributeError: + logger.warning( + "It seems the MRIOT you loaded doesn't have metadata to print." + ) + pym_mrio = lexico_reindex(pym_mrio) self.main_inv_dur: int = main_inv_dur r"""int, default 90: The default numbers of days for inputs inventory to use if it is not defined for an input.""" @@ -174,11 +192,16 @@ def __init__( self.n_sectors = len(sec) r"""int : The number :math:`n` of sectors.""" - self.industries = pd.MultiIndex.from_product([self.regions, self.sectors]) + self.industries = pd.MultiIndex.from_product( + [self.regions, self.sectors], names=["region", "sector"] + ) r"""pandas.MultiIndex : A pandas MultiIndex of the industries (region,sector) of the model.""" + self.n_industries = len(self.industries) + r"""int : The number :math:`m * n` of industries.""" + try: - self.fd_cat = np.array(sorted(list(pym_mrio.get_Y_categories()))) # type: ignore + self.final_demand_cat = np.array(sorted(list(pym_mrio.get_Y_categories()))) # type: ignore r"""numpy.ndarray of str: An array of the final demand categories of the model (``["Final demand"]`` if there is only one)""" self.n_fd_cat = len(pym_mrio.get_Y_categories()) # type: ignore @@ -186,12 +209,18 @@ def __init__( except KeyError: self.n_fd_cat = 1 - self.fd_cat = np.array(["Final demand"]) + self.final_demand_cat = np.array(["Final demand"]) except IndexError: self.n_fd_cat = 1 - self.fd_cat = np.array(["Final demand"]) + self.final_demand_cat = np.array(["Final demand"]) - self.monetary_factor = monetary_factor + if hasattr(pym_mrio, "monetary_factor"): + logger.warning( + f"Custom monetary factor found in the mrio pickle file, continuing with this one ({getattr(pym_mrio,'monetary_factor')})" + ) + self.monetary_factor = getattr(pym_mrio, "monetary_factor") + else: + self.monetary_factor = monetary_factor r"""int, default 10^6: Monetary unit factor (i.e. if the tables unit is 10^6 € instead of 1 €, it should be set to 10^6).""" self.n_temporal_units_by_step = temporal_units_by_step @@ -231,18 +260,24 @@ def __init__( r"""numpy.ndarray of int: Array :math:`\mathbf{s}` of size :math:`n` (sectors), setting for each input the initial number of ``n_temporal_units_by_step`` of stock for the input. (see :ref:`boario-math`).""" else: - inv = inventory_dict + infinite_inventories_sect = ( + [] if infinite_inventories_sect is None else infinite_inventories_sect + ) self.inventories = [ - np.inf if inv[k] == "inf" else inv[k] for k in sorted(inv.keys()) + np.inf + if inventory_dict[k] in ["inf", "Inf", "Infinity", "infinity"] + or k in infinite_inventories_sect + else inventory_dict[k] + for k in sorted(inventory_dict.keys()) ] + logger.debug(f"inventories: {self.inventories}") + logger.debug(f"n_temporal_units_by_step: {self.n_temporal_units_by_step}") self.inv_duration = np.array(self.inventories) / self.n_temporal_units_by_step if (self.inv_duration <= 1).any(): logger.warning( - "At least one product has inventory duration lower than the numbers of temporal units in one step ({}), model will set it to 2 by default, but you should probably check this !".format( - self.n_temporal_units_by_step - ) + f"At least one product has inventory duration lower than the numbers of temporal units in one step ({self.n_temporal_units_by_step}), model will set it to 2 by default, but you should probably check this !" ) self.inv_duration[self.inv_duration <= 1] = 2 ################################################################# @@ -263,7 +298,7 @@ def __init__( self._matrix_id = np.eye(self.n_sectors) self._matrix_I_sum = np.tile(self._matrix_id, self.n_regions) self.Z_0 = pym_mrio.Z.to_numpy() - r"""numpy.ndarray of float: 2-dim square matrix array :math:`\ioz` of size :math:`(n \times m, n \times m)` representing the intermediate (transaction) matrix (see :ref:`boario-math-init`).""" + r"""numpy.ndarray of float: 2-dim square matrix array :math:`\ioz` of size :math:`(n \times m, n \times m)` representing the daily intermediate (transaction) matrix (see :ref:`boario-math-init`).""" self.Z_C = self._matrix_I_sum @ self.Z_0 r"""numpy.ndarray of float: 2-dim matrix array :math:`\ioz^{\sectorsset}` of size :math:`(n, n \times m)` representing the intermediate (transaction) matrix aggregated by inputs (see :ref:`here `).""" @@ -277,17 +312,36 @@ def __init__( self.Z_distrib = np.nan_to_num(self.Z_distrib) self.Z_0 = pym_mrio.Z.to_numpy() * self.steply_factor + self.Y_0 = pym_mrio.Y.to_numpy() + self.Y_C = self._matrix_I_sum @ self.Y_0 + r"""numpy.ndarray of float: 2-dim matrix array :math:`\ioy^{\sectorsset}` of size :math:`(n, m \times \text{number of final demand categories})` representing the final demand matrix aggregated by inputs (see :ref:`here `).""" + + with np.errstate(divide="ignore", invalid="ignore"): + self.Y_distrib = np.divide( + self.Y_0, (np.tile(self.Y_C, (self.n_regions, 1))) + ) + r"""numpy.ndarray of float: math:`\ioy` normalised by :math:`\ioy^{\sectorsset}`, i.e. representing for each input the share of the total ordered transiting from an industry to final demand (see :ref:`here `).""" + + self.Y_distrib = np.nan_to_num(self.Y_distrib) + self.Y_0 = pym_mrio.Y.to_numpy() * self.steply_factor - r"""numpy.ndarray of float: 2-dim array :math:`\ioy` of size :math:`(n \times m, m \times \text{number of final demand categories})` representing the final demand matrix.""" + r"""numpy.ndarray of float: 2-dim array :math:`\ioy` of size :math:`(n \times m, m \times \text{number of final demand categories})` representing the daily final demand matrix.""" tmp = np.array(pym_mrio.x.T) self.X_0 = tmp.flatten() * self.steply_factor - r"""numpy.ndarray of float: Array :math:`\iox(0)` of size :math:`n \times m` representing the initial gross production.""" + r"""numpy.ndarray of float: Array :math:`\iox(0)` of size :math:`n \times m` representing the initial daily gross production.""" del tmp value_added = pym_mrio.x.T - pym_mrio.Z.sum(axis=0) value_added = value_added.reindex(sorted(value_added.index), axis=0) value_added = value_added.reindex(sorted(value_added.columns), axis=1) + if (value_added < 0).any().any(): + logger.warning( + f"""Found negative values in the value added, will set to 0. Note that sectors with null value added are not impacted by capital destruction. + industries with negative VA: {(value_added[value_added < 0].dropna(axis=1)).columns} + """ + ) + value_added[value_added < 0] = 0.0 self.gdp_df = value_added.groupby("region", axis=1).sum() r"""pandas.DataFrame: Dataframe of the total GDP of each region of the model""" @@ -298,18 +352,18 @@ def __init__( self.tech_mat = (self._matrix_I_sum @ pym_mrio.A).to_numpy() # type: ignore #to_numpy is not superfluous ! r"""numpy.ndarray of float: 2-dim array :math:`\ioa` of size :math:`(n \times m, n \times m)` representing the technical coefficients matrix.""" - if kapital_vector is not None: - self.k_stock = kapital_vector + if productive_capital_vector is not None: + self.k_stock = productive_capital_vector """numpy.ndarray of float: Array of size :math:`n \times m` representing the estimated stock of capital of each industry.""" if isinstance(self.k_stock, pd.DataFrame): self.k_stock = self.k_stock.squeeze().sort_index().to_numpy() - elif kapital_to_VA_dict is None: + elif productive_capital_to_VA_dict is None: logger.warning(f"No capital to VA dictionary given, considering 4/1 ratio") self.kstock_ratio_to_VA = 4 self.k_stock = self.VA_0 * self.kstock_ratio_to_VA else: - kratio = kapital_to_VA_dict + kratio = productive_capital_to_VA_dict kratio_ordered = [kratio[k] for k in sorted(kratio.keys())] self.kstock_ratio_to_VA = np.tile(np.array(kratio_ordered), self.n_regions) self.k_stock = self.VA_0 * self.kstock_ratio_to_VA @@ -354,7 +408,7 @@ def __init__( self.intmd_demand = self.Z_0.copy() self.final_demand = self.Y_0.copy() self.rebuilding_demand = None - self.kapital_lost = np.zeros(self.production.shape) + self.productive_capital_lost = np.zeros(self.production.shape) r"""numpy.ndarray of float: Array of size :math:`n \times m` representing the estimated stock of capital currently destroyed for each industry.""" self.order_type = order_type @@ -362,10 +416,18 @@ def __init__( ################## SIMULATION TRACKING VARIABLES ################ self.in_shortage = False + r"""Boolean stating if at least one industry is in shortage (i.e.) if at least one of its inputs inventory is low enough to reduce production.""" + self.had_shortage = False + r"""Boolean stating if at least one industry was in shortage at some point.""" + self.rebuild_prod = np.zeros(shape=self.X_0.shape) + r"""numpy.ndarray of float: Array of size :math:`n \times m` representing the remaining stock of rebuilding demand asked of each industry.""" + self.final_demand_not_met = np.zeros(self.Y_0.shape) - ############################################################################# + r"""numpy.ndarray of float: Array of size :math:`n \times m` representing the final demand that could not be met at current step for each industry.""" + + ################################################################# ### Internals self._prod_delta_type = None @@ -375,59 +437,72 @@ def __init__( logger.debug( f"Setting possible regions (currently: {Event.possible_regions}) to: {self.regions}" ) - Event.possible_regions = self.regions.copy() + Event.possible_regions = pd.CategoricalIndex( + self.regions, name="region", copy=True + ) logger.debug(f"Possible regions is now {Event.possible_regions}") Event.regions_idx = np.arange(self.n_regions) - Event.possible_sectors = self.sectors.copy() + Event.possible_sectors = pd.CategoricalIndex( + self.sectors, name="sector", copy=True + ) + Event.possible_final_demand_cat = pd.CategoricalIndex( + pym_mrio.get_Y_categories(), name="final demand", copy=True + ) Event.sectors_idx = np.arange(self.n_sectors) Event.z_shape = self.Z_0.shape Event.y_shape = self.Y_0.shape Event.x_shape = self.X_0.shape - Event.monetary_factor = monetary_factor + Event.model_monetary_factor = monetary_factor Event.sectors_gva_shares = self.gdp_share_sector.copy() Event.Z_distrib = self.Z_distrib.copy() + Event.Y_distrib = self.Y_distrib.copy() + # logger.warning(f"{value_added}") + Event.gva_df = value_added.iloc[0].copy() meta = pym_mrio.meta.metadata + assert isinstance(meta, dict) + # meta = {"name":"unnamed", "description":"", "system":"unknown","version":"unknown"} try: - Event.mrio_name = ( - meta["name"] - + "_" - + meta["description"] - + "_" - + meta["system"] - + "_" - + meta["version"] + Event.mrio_name = "_".join( + [ + str(meta["name"]), + str(meta["description"]), + str(meta["system"]), + str(meta["version"]), + ] ) except TypeError: Event.mrio_name = "custom - method WIP" # initialize those (it's not very nice, but otherwise python complains) - self._indus_rebuild_demand_tot = None - self._house_rebuild_demand_tot = None - self._indus_rebuild_demand = None - self._house_rebuild_demand = None - self._tot_rebuild_demand = None - self._kapital_lost = np.zeros(self.VA_0.shape) - self._prod_cap_delta_kapital = None - self._prod_cap_delta_arbitrary = None - self._prod_cap_delta_tot = None + self._indus_rebuild_demand_tot: npt.NDArray = np.array([]) + self._house_rebuild_demand_tot: npt.NDArray = np.array([]) + self._indus_rebuild_demand: npt.NDArray = np.array([]) + self._house_rebuild_demand: npt.NDArray = np.array([]) + self._tot_rebuild_demand: npt.NDArray = np.array([]) + self._productive_capital_lost: npt.NDArray = np.zeros(self.VA_0.shape) + self._prod_cap_delta_productive_capital: npt.NDArray = np.zeros( + len(self.industries) + ) + self._prod_cap_delta_arbitrary: npt.NDArray = np.zeros(len(self.industries)) + self._prod_cap_delta_tot: npt.NDArray = np.zeros(len(self.industries)) ## Properties @property - def tot_rebuild_demand(self) -> Optional[np.ndarray]: + def tot_rebuild_demand(self) -> npt.NDArray: r"""Returns current total rebuilding demand (as the sum of rebuilding demand addressed to each industry)""" tmp = [] logger.debug("Trying to return tot_rebuilding demand") - if self._indus_rebuild_demand_tot is not None: + if np.any(self._indus_rebuild_demand_tot > 0): tmp.append(self._indus_rebuild_demand_tot) - if self._house_rebuild_demand_tot is not None: + if np.any(self._house_rebuild_demand_tot > 0): tmp.append(self._house_rebuild_demand_tot) if tmp: ret = np.concatenate(tmp, axis=1).sum(axis=1) self._tot_rebuild_demand = ret else: - self._tot_rebuild_demand = None + self._tot_rebuild_demand = np.array([]) return self._tot_rebuild_demand @tot_rebuild_demand.setter @@ -446,21 +521,18 @@ def tot_rebuild_demand(self, source: list[EventKapitalRebuild]): logger.debug(f"Trying to set tot_rebuilding demand from {source}") if not isinstance(source, list): ValueError( - "Setting tot_rebuild_demand can only be done with a list of events, not a {}".format( - type(source) - ) + f"Setting tot_rebuild_demand can only be done with a list of events, not a {type(source)}" ) self.house_rebuild_demand = source self.indus_rebuild_demand = source @property - def house_rebuild_demand(self) -> Optional[np.ndarray]: - r"""Returns household rebuilding demand matrix or `None` if - there is no such demand. + def house_rebuild_demand(self) -> npt.NDArray: + r"""Returns household rebuilding demand matrix. Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\ioy`, containing the sum of all currently rebuildable final demand stock. """ @@ -468,13 +540,12 @@ def house_rebuild_demand(self) -> Optional[np.ndarray]: return self._house_rebuild_demand @property - def house_rebuild_demand_tot(self) -> Optional[np.ndarray]: - r"""Returns total household rebuilding demand vector or `None` if - there is no such demand. + def house_rebuild_demand_tot(self) -> npt.NDArray: + r"""Returns total household rebuilding demand vector. Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\iox`, containing the sum of all currently rebuildable households demands. """ @@ -498,27 +569,23 @@ def house_rebuild_demand(self, source: list[EventKapitalRebuild]): events : 'list[Event]' A list of Event objects - Notes - ----- - - So far the model wasn't tested with such a rebuilding demand. Only intermediate demand was considered. """ tmp = [] - for ev in source: - if ev.rebuildable: - if ev.rebuild_tau is None: + for evnt in source: + if evnt.rebuildable: + if evnt.rebuild_tau is None: rebuild_tau = self.rebuild_tau else: - rebuild_tau = ev.rebuild_tau + rebuild_tau = evnt.rebuild_tau warn_once(logger, "Event has a custom rebuild_tau") - if ev.rebuilding_demand_house is not None: + if len(evnt.rebuilding_demand_house) > 0: tmp.append( - ev.rebuilding_demand_house + evnt.rebuilding_demand_house * (self.n_temporal_units_by_step / rebuild_tau) ) if not tmp: - self._house_rebuild_demand = None - self._house_rebuild_demand_tot = None + self._house_rebuild_demand = np.array([]) + self._house_rebuild_demand_tot = np.array([]) else: house_reb_dem = np.stack(tmp, axis=-1) tot = house_reb_dem.sum(axis=1) @@ -526,14 +593,12 @@ def house_rebuild_demand(self, source: list[EventKapitalRebuild]): self._house_rebuild_demand_tot = tot @property - def indus_rebuild_demand(self) -> Optional[np.ndarray]: - r"""Returns industrial rebuilding demand matrix or `None` if - there is no such demand. - + def indus_rebuild_demand(self) -> npt.NDArray: + r"""Returns industrial rebuilding demand matrix. Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\ioz`, containing the sum of all currently rebuildable intermediate demand stock. """ @@ -541,13 +606,12 @@ def indus_rebuild_demand(self) -> Optional[np.ndarray]: return self._indus_rebuild_demand @property - def indus_rebuild_demand_tot(self) -> Optional[np.ndarray]: - r"""Returns total industrial rebuilding demand vector or `None` if - there is no such demand. + def indus_rebuild_demand_tot(self) -> npt.NDArray: + r"""Returns total industrial rebuilding demand vector. Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\iox`, containing the sum of all currently rebuildable intermediate demands. """ @@ -570,42 +634,44 @@ def indus_rebuild_demand(self, source: list[EventKapitalRebuild]): """ tmp = [] - for ev in source: - if ev.rebuildable: - if ev.rebuild_tau is None: + for evnt in source: + if evnt.rebuildable: + if evnt.rebuild_tau is None: rebuild_tau = self.rebuild_tau else: - rebuild_tau = ev.rebuild_tau + rebuild_tau = evnt.rebuild_tau warn_once(logger, "Event has a custom rebuild_tau") tmp.append( - ev.rebuilding_demand_indus + evnt.rebuilding_demand_indus * (self.n_temporal_units_by_step / rebuild_tau) ) if not tmp: - self._indus_rebuild_demand = None - self._indus_rebuild_demand_tot = None + self._indus_rebuild_demand = np.array([]) + self._indus_rebuild_demand_tot = np.array([]) else: indus_reb_dem = np.stack(tmp, axis=-1) self._indus_rebuild_demand = indus_reb_dem self._indus_rebuild_demand_tot = indus_reb_dem.sum(axis=1) - logger.debug(f"Setting indus_rebuild_demand_tot to {indus_reb_dem}") + # logger.debug(f"Setting indus_rebuild_demand_tot to {indus_reb_dem}") @property - def kapital_lost(self) -> np.ndarray: + def productive_capital_lost(self) -> npt.NDArray: r"""Returns current stock of destroyed capital Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\iox`, containing the "stock" of capital currently destroyed for each industry. """ - return self._kapital_lost + return self._productive_capital_lost - @kapital_lost.setter - def kapital_lost(self, source: list[EventKapitalDestroyed] | np.ndarray) -> None: + @productive_capital_lost.setter + def productive_capital_lost( + self, source: list[EventKapitalDestroyed] | npt.NDArray + ) -> None: r"""Computes current capital lost and update production delta accordingly. Computes and sets the current stock of capital lost by each industry of @@ -615,58 +681,67 @@ def kapital_lost(self, source: list[EventKapitalDestroyed] | np.ndarray) -> None Parameters ---------- - source : list[EventKapitalDestroyed] | np.ndarray + source : list[EventKapitalDestroyed] | npt.NDArray Either a list of events to consider for the destruction of capital or directly a vector of destroyed capital for each industry. """ - logger.debug("Updating kapital lost from list of events") + logger.debug("Updating productive_capital lost from list of events") if isinstance(source, list): if source: - self._kapital_lost = np.add.reduce( - np.array([e.regional_sectoral_kapital_destroyed for e in source]) + self._productive_capital_lost = np.add.reduce( + np.array( + [ + e.regional_sectoral_productive_capital_destroyed + for e in source + ] + ) ) else: - self._kapital_lost = np.zeros(self.VA_0.shape) + self._productive_capital_lost = np.zeros(self.VA_0.shape) elif isinstance(source, np.ndarray): - self._kapital_lost = source + self._productive_capital_lost = source - productivity_loss_from_K = np.zeros(shape=self._kapital_lost.shape) + productivity_loss_from_capital_destroyed = np.zeros( + shape=self._productive_capital_lost.shape + ) np.divide( - self._kapital_lost, + self._productive_capital_lost, self.k_stock, - out=productivity_loss_from_K, + out=productivity_loss_from_capital_destroyed, where=self.k_stock != 0, ) - logger.debug("Updating production delta from kapital loss") - self._prod_cap_delta_kapital = productivity_loss_from_K - if (self._prod_cap_delta_kapital > 0.0).any(): + logger.debug("Updating production delta from productive_capital loss") + self._prod_cap_delta_productive_capital = ( + productivity_loss_from_capital_destroyed + ) + if (self._prod_cap_delta_productive_capital > 0.0).any(): if self._prod_delta_type is None: - self._prod_delta_type = "from_kapital" + self._prod_delta_type = "from_productive_capital" elif self._prod_delta_type == "from_arbitrary": - self._prod_delta_type = "mixed_from_kapital_from_arbitrary" + self._prod_delta_type = "mixed_from_productive_capital_from_arbitrary" @property - def prod_cap_delta_arbitrary(self) -> Optional[np.ndarray]: + def prod_cap_delta_arbitrary(self) -> npt.NDArray: r"""Return the possible "arbitrary" production capacity lost vector if it was set. Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\iox`, stating the amount of production capacity lost arbitrarily (ie exogenous). """ return self._prod_cap_delta_arbitrary @prod_cap_delta_arbitrary.setter - def prod_cap_delta_arbitrary(self, source: list[Event] | np.ndarray): + def prod_cap_delta_arbitrary(self, source: list[EventArbitraryProd] | npt.NDArray): """Computes and sets the loss of production capacity from "arbitrary" sources. Parameters ---------- - source : list[Event] | np.ndarray + source : list[Event] | npt.NDArray Either a list of Event objects with arbitrary production losses set, or directly a vector of production capacity loss. @@ -674,42 +749,46 @@ def prod_cap_delta_arbitrary(self, source: list[Event] | np.ndarray): if isinstance(source, list): event_arb = np.array( - [ev for ev in source if ev.prod_cap_delta_arbitrary is not None] + [ + ev.prod_cap_delta_arbitrary + for ev in source + if len(ev.prod_cap_delta_arbitrary) > 0 + ] ) if event_arb.size == 0: self._prod_cap_delta_arbitrary = np.zeros(shape=self.X_0.shape) else: - self._prod_cap_delta_arbitrary = np.max.reduce(event_arb) + self._prod_cap_delta_arbitrary = np.maximum.reduce(event_arb) else: self._prod_capt_delta_arbitrary = source assert self._prod_cap_delta_arbitrary is not None if (self._prod_cap_delta_arbitrary > 0.0).any(): if self._prod_delta_type is None: self._prod_delta_type = "from_arbitrary" - elif self._prod_delta_type == "from_kapital": - self._prod_delta_type = "mixed_from_kapital_from_arbitrary" + elif self._prod_delta_type == "from_productive_capital": + self._prod_delta_type = "mixed_from_productive_capital_from_arbitrary" @property - def prod_cap_delta_kapital(self) -> Optional[np.ndarray]: + def prod_cap_delta_productive_capital(self) -> npt.NDArray: r"""Return the possible production capacity lost due to capital destroyed vector if it was set. Returns ------- - np.ndarray + npt.NDArray An array of same shape as math:`\iox`, stating the amount of production capacity lost due to capital destroyed. """ - return self._prod_cap_delta_kapital + return self._prod_cap_delta_productive_capital @property - def prod_cap_delta_tot(self) -> np.ndarray: + def prod_cap_delta_tot(self) -> npt.NDArray: r"""Computes and return total current production delta. Returns ------- - np.ndarray + npt.NDArray The total production delta (ie share of production capacity lost) for each industry. @@ -719,30 +798,30 @@ def prod_cap_delta_tot(self) -> np.ndarray: tmp = [] if self._prod_delta_type is None: raise AttributeError("Production delta doesn't appear to be set yet.") - elif self._prod_delta_type == "from_kapital": - logger.debug("Production delta is only set from kapital destruction") - tmp.append(self._prod_cap_delta_kapital) + elif self._prod_delta_type == "from_productive_capital": + logger.debug( + "Production delta is only set from productive_capital destruction" + ) + tmp.append(self._prod_cap_delta_productive_capital) elif self._prod_delta_type == "from_arbitrary": logger.debug("Production delta is only set from arbitrary delta") tmp.append(self._prod_cap_delta_arbitrary) - elif self._prod_delta_type == "mixed_from_kapital_from_arbitrary": + elif self._prod_delta_type == "mixed_from_productive_capital_from_arbitrary": logger.debug( - "Production delta is a mixed form of kapital destruction and arbitrary delta" + "Production delta is a mixed form of productive_capital destruction and arbitrary delta" ) - tmp.append(self._prod_cap_delta_kapital) + tmp.append(self._prod_cap_delta_productive_capital) tmp.append(self._prod_cap_delta_arbitrary) else: raise NotImplementedError( - "Production delta type {} not recognised".format(self._prod_delta_type) + f"Production delta type {self._prod_delta_type} not recognised" ) - tmp.append(np.ones(shape=self.X_0.shape)) + # tmp.append(np.ones(shape=self.X_0.shape)) # logger.debug("tmp: {}".format(tmp)) - self._prod_cap_delta_tot = np.amin(np.stack(tmp, axis=-1), axis=1) + self._prod_cap_delta_tot = np.amax(np.stack(tmp, axis=-1), axis=1) assert ( self._prod_cap_delta_tot.shape == self.X_0.shape - ), "expected shape {}, received {}".format( - self.X_0.shape, self._prod_cap_delta_tot.shape - ) + ), f"expected shape {self.X_0.shape}, received {self._prod_cap_delta_tot.shape}" return self._prod_cap_delta_tot @prod_cap_delta_tot.setter @@ -760,12 +839,10 @@ def prod_cap_delta_tot(self, source: list[Event]): logger.debug("Updating total production delta from list of events") if not isinstance(source, list): ValueError( - "Setting prod_cap_delta_tot can only be done with a list of events, not a {}".format( - type(source) - ) + f"Setting prod_cap_delta_tot can only be done with a list of events, not a {type(source)}" ) - self.kapital_lost = [ + self.productive_capital_lost = [ event for event in source if isinstance(event, EventKapitalDestroyed) ] self.prod_cap_delta_arbitrary = [ @@ -790,7 +867,7 @@ def update_system_from_events(self, events: list[Event]) -> None: ] @property - def production_cap(self) -> np.ndarray: + def production_cap(self) -> npt.NDArray: r"""Compute and update production capacity. Compute and update production capacity from current total production delta and overproduction. @@ -818,8 +895,8 @@ def production_cap(self) -> np.ndarray: return production_cap @property - def total_demand(self) -> np.ndarray: - r"""Computes and return total demand as the sum of intermediate demand (orders), final demand, and possible rebuilding demand.""" + def total_demand(self) -> npt.NDArray: + r"""Computes and returns total demand as the sum of intermediate demand (orders), final demand, and possible rebuilding demand.""" if (self.matrix_orders < 0).any(): raise RuntimeError("Some matrix orders are negative which shouldn't happen") @@ -827,12 +904,12 @@ def total_demand(self) -> np.ndarray: tot_dem = self.matrix_orders.sum(axis=1) + self.final_demand.sum(axis=1) if (tot_dem < 0).any(): raise RuntimeError("Some total demand are negative which shouldn't happen") - if self.tot_rebuild_demand is not None: + if len(self.tot_rebuild_demand) > 0: tot_dem += self.tot_rebuild_demand return tot_dem @property - def production_opt(self) -> np.ndarray: + def production_opt(self) -> npt.NDArray: r"""Computes and returns "optimal production" :math:`\iox^{textrm{Opt}}`, as the per industry minimum between total demand and production capacity. @@ -841,17 +918,17 @@ def production_opt(self) -> np.ndarray: return np.fmin(self.total_demand, self.production_cap) @property - def inventory_constraints_opt(self) -> np.ndarray: + def inventory_constraints_opt(self) -> npt.NDArray: r"""Computes and returns inventory constraints for "optimal production" (see :meth:`calc_inventory_constraints`)""" return self.calc_inventory_constraints(self.production_opt) @property - def inventory_constraints_act(self) -> np.ndarray: + def inventory_constraints_act(self) -> npt.NDArray: r"""Computes and returns inventory constraints for "actual production" (see :meth:`calc_inventory_constraints`)""" return self.calc_inventory_constraints(self.production) - def calc_production(self, current_temporal_unit: int) -> np.ndarray: + def calc_production(self, current_temporal_unit: int) -> npt.NDArray: r"""Computes and updates actual production. See :ref:`boario-math-prod`. 1. Computes ``production_opt`` and ``inventory_constraints`` as : @@ -906,9 +983,7 @@ def calc_production(self, current_temporal_unit: int) -> np.ndarray: ).any(): if not self.in_shortage: logger.info( - "At least one industry entered shortage regime. (step:{})".format( - current_temporal_unit - ) + f"At least one industry entered shortage regime. (step:{current_temporal_unit})" ) self.in_shortage = True self.had_shortage = True @@ -934,15 +1009,13 @@ def calc_production(self, current_temporal_unit: int) -> np.ndarray: if self.in_shortage: self.in_shortage = False logger.info( - "All industries exited shortage regime. (step:{})".format( - current_temporal_unit - ) + f"All industries exited shortage regime. (step:{current_temporal_unit})" ) assert not (production_opt < 0).any() self.production = production_opt return stock_constraint - def calc_inventory_constraints(self, production: np.ndarray) -> np.ndarray: + def calc_inventory_constraints(self, production: npt.NDArray) -> npt.NDArray: r"""Compute inventory constraints (no psi parameter, for the psi version, the recommended one, see :meth:`~boario.extended_models.ARIOPsiModel.calc_inventory_constraints`) @@ -950,12 +1023,12 @@ def calc_inventory_constraints(self, production: np.ndarray) -> np.ndarray: Parameters ---------- - production : np.ndarray + production : npt.NDArray The production vector to consider. Returns ------- - np.ndarray + npt.NDArray For each input, for each industry, the size of the inventory required to produce at `production` level for the duration goal (`inv_duration`). @@ -1054,16 +1127,23 @@ def distribute_production( """ if scheme != "proportional": - raise ValueError("Scheme %s not implemented" % scheme) + raise ValueError(f"Scheme {scheme} is currently not implemented") # list_of_demands = [self.matrix_orders, self.final_demand] ## 1. Calc demand from rebuilding requirements (with characteristic time rebuild_tau) + house_reb_dem_per_event = ( + house_reb_dem_tot_per_event + ) = indus_reb_dem_per_event = indus_reb_dem_tot_per_event = None + if rebuildable_events: logger.debug("There are rebuildable events") n_events = len(rebuildable_events) tot_rebuilding_demand_summed = self.tot_rebuild_demand.copy() # debugging assert - assert tot_rebuilding_demand_summed.shape == self.X_0.shape + if not tot_rebuilding_demand_summed.shape == self.X_0.shape: + raise RuntimeError( + f"received {tot_rebuilding_demand_summed.shape}, expected {self.X_0.shape}" + ) indus_reb_dem_tot_per_event = self.indus_rebuild_demand_tot.copy() indus_reb_dem_per_event = self.indus_rebuild_demand.copy() @@ -1155,14 +1235,8 @@ def distribute_production( if not np.allclose(stock_add, stock_use): self.matrix_stock = self.matrix_stock - stock_use + stock_add if (self.matrix_stock < 0).any(): - self.matrix_stock.dump(self.records_storage / "matrix_stock_dump.pkl") - logger.error( - "Negative values in the stocks, matrix has been dumped in the results dir : \n {}".format( - self.records_storage / "matrix_stock_dump.pkl" - ) - ) raise RuntimeError( - "stock_add (restocking) contains negative values, this should not happen" + "matrix_stock contains negative values, this should not happen" ) # 6. Compute final demand not met due to rationing @@ -1192,22 +1266,68 @@ def distribute_production( ) self.rebuild_prod = rebuild_prod.copy() + logger.debug( + f"tot_rebuilding_demand_summed: {pd.Series(tot_rebuilding_demand_summed, index=self.industries)}" + ) + if h_reb and ind_reb: - # logger.debug("Entering here") + logger.debug( + "There are both household rebuilding demand and industry rebuilding demand" + ) + logger.debug( + f"indus_reb_dem_tot_per_event: {pd.DataFrame(indus_reb_dem_tot_per_event, index=self.industries)}" + ) + logger.debug( + f"house_reb_dem_tot_per_event: {pd.DataFrame(house_reb_dem_tot_per_event, index=self.industries)}" + ) + logger.debug( + f"dem_tot_per_event: {pd.DataFrame(house_reb_dem_tot_per_event+indus_reb_dem_tot_per_event, index=self.industries)}" # type: ignore + ) + + tot_reb_dem_divider = np.tile( + tot_rebuilding_demand_summed[:, np.newaxis], + (1, indus_reb_dem_tot_per_event.shape[1]), + ) indus_shares = np.divide( indus_reb_dem_tot_per_event, - tot_rebuilding_demand_summed, - where=(tot_rebuilding_demand_summed != 0), + tot_reb_dem_divider, + out=np.zeros_like(indus_reb_dem_tot_per_event), + where=tot_reb_dem_divider != 0, ) house_shares = np.divide( house_reb_dem_tot_per_event, - tot_rebuilding_demand_summed, - where=(tot_rebuilding_demand_summed != 0), + tot_reb_dem_divider, + out=np.zeros_like(house_reb_dem_tot_per_event), + where=tot_reb_dem_divider != 0, + ) + # np.around(house_shares, 6, out=house_shares) + # np.around(indus_shares, 6, out=indus_shares) + logger.debug( + f"indus_shares: {pd.DataFrame(indus_shares, index=self.industries)}" + ) + logger.debug( + f"house_shares: {pd.DataFrame(house_shares, index=self.industries)}" ) - # logger.debug("indus_shares: {}".format(indus_shares)) + + # tot_reb = pd.Series(tot_rebuilding_demand_summed, index=self.industries) + # indus_reb_dem_tot = pd.DataFrame(indus_reb_dem_tot_per_event, index=self.industries)[0] + # house_reb_dem_tot = pd.DataFrame(house_reb_dem_tot_per_event, index=self.industries)[0] + # np.around(indus_shares,10,indus_shares) + # np.around(house_shares,10,house_shares) + # np.around(tmp,10,tmp) + # indus_df = pd.DataFrame(indus_shares,index=self.industries)[0] + # house_df = pd.DataFrame(house_shares,index=self.industries)[0] + # tmp_df=pd.DataFrame(tmp, index=self.industries)[0] + + # debug_df = pd.concat([tot_reb, indus_reb_dem_tot, house_reb_dem_tot, tot_reb-(house_reb_dem_tot+indus_reb_dem_tot), indus_df, house_df, tmp_df], axis=1) + # #logger.info(debug_df[debug_df.iloc[:,6]!=0].to_string()) + # debug_str = f""" + # {debug_df[debug_df.iloc[:,6]!=0].to_string()} + # """ assert np.allclose( - indus_shares + house_shares, np.ones(shape=indus_shares.shape) + (indus_shares + house_shares)[tot_rebuilding_demand_summed != 0], 1.0 ) + elif h_reb: house_shares = np.ones(house_reb_dem_tot_per_event.shape) indus_shares = np.zeros(indus_reb_dem_tot_per_event.shape) @@ -1218,9 +1338,23 @@ def distribute_production( else: return [] + logger.debug(f"rebuild_prod: {rebuild_prod}") + logger.debug(f"indus_shares: {indus_shares}") + indus_rebuild_prod = rebuild_prod[:, np.newaxis] * indus_shares # type:ignore house_rebuild_prod = rebuild_prod[:, np.newaxis] * house_shares # type:ignore + if ind_reb: + logger.debug( + f"indus_rebuild_prod: {pd.DataFrame(indus_rebuild_prod, index=self.industries)}" + ) + if h_reb: + logger.debug( + f"house_rebuild_prod: {pd.DataFrame(house_rebuild_prod, index=self.industries)}" + ) + # logger.debug(f"tot_rebuilding_demand_summed: {pd.DataFrame(tot_rebuilding_demand_summed, index=self.industries)}") + # logger.debug(f"house_rebuild_prod: {pd.Series(house_rebuild_prod.flatten(), index=self.industries)}") + indus_rebuild_prod_distributed = np.zeros(shape=indus_reb_dem_per_event.shape) house_rebuild_prod_distributed = np.zeros(shape=house_reb_dem_per_event.shape) if ind_reb: @@ -1252,7 +1386,7 @@ def distribute_production( if h_reb: house_rebuilding_demand_shares = np.zeros( - shape=indus_reb_dem_per_event.shape + shape=house_reb_dem_per_event.shape ) house_rebuilding_demand_broad = np.broadcast_to( house_reb_dem_tot_per_event[:, np.newaxis], @@ -1275,18 +1409,22 @@ def distribute_production( # update rebuilding demand events_to_remove = [] - for e_id, e in enumerate(rebuildable_events): - if e.rebuilding_demand_indus is not None: - e.rebuilding_demand_indus -= indus_rebuild_prod_distributed[:, :, e_id] - if e.rebuilding_demand_house is not None: - e.rebuilding_demand_house -= house_rebuild_prod_distributed[:, :, e_id] - if (e.rebuilding_demand_indus < (10 / self.monetary_factor)).all() and ( - e.rebuilding_demand_house < (10 / self.monetary_factor) + for e_id, evnt in enumerate(rebuildable_events): + if len(evnt.rebuilding_demand_indus) > 0: + evnt.rebuilding_demand_indus -= indus_rebuild_prod_distributed[ + :, :, e_id + ] + if len(evnt.rebuilding_demand_house) > 0: + evnt.rebuilding_demand_house -= house_rebuild_prod_distributed[ + :, :, e_id + ] + if (evnt.rebuilding_demand_indus < (10 / self.monetary_factor)).all() and ( + evnt.rebuilding_demand_house < (10 / self.monetary_factor) ).all(): - events_to_remove.append(e) + events_to_remove.append(evnt) return events_to_remove - def calc_matrix_stock_gap(self, matrix_stock_goal) -> np.ndarray: + def calc_matrix_stock_gap(self, matrix_stock_goal) -> npt.NDArray: """Computes and returns inputs stock gap matrix The gap is simply the difference between the goal (given as argument) @@ -1294,12 +1432,12 @@ def calc_matrix_stock_gap(self, matrix_stock_goal) -> np.ndarray: Parameters ---------- - matrix_stock_goal : np.ndarray of float + matrix_stock_goal : npt.NDArray of float The target inventories. Returns ------- - np.ndarray + npt.NDArray The (only positive) gap between goal and current inventories. Raises @@ -1340,7 +1478,7 @@ def calc_orders(self) -> None: with np.errstate(invalid="ignore"): matrix_stock_goal *= self.inv_duration[:, np.newaxis] if np.allclose(self.matrix_stock, matrix_stock_goal): - matrix_stock_gap = matrix_stock_goal * 0 + matrix_stock_gap = np.zeros(matrix_stock_goal.shape) else: matrix_stock_gap = self.calc_matrix_stock_gap(matrix_stock_goal) matrix_stock_gap += ( @@ -1363,7 +1501,7 @@ def calc_orders(self) -> None: self.matrix_orders = tmp def calc_overproduction(self) -> None: - r"""Computes and update the overproduction vector. + r"""Computes and updates the overproduction vector. See :ref:`Overproduction module ` @@ -1386,92 +1524,15 @@ def calc_overproduction(self) -> None: self.overprod += overprod_chg self.overprod[self.overprod < 1.0] = 1.0 - # def check_stock_increasing(self, current_temporal_unit: int): - # tmp = np.full(self.matrix_stock.shape, 0.0) - # mask = np.isfinite(self.matrix_stock_0) - # np.subtract(self.matrix_stock, self.matrix_stock_0, out=tmp, where=mask) - # check_1 = tmp > 0.0 - # tmp = np.full(self.matrix_stock.shape, 0.0) - # np.subtract( - # self.stocks_evolution[current_temporal_unit], - # self.stocks_evolution[current_temporal_unit - 1], - # out=tmp, - # where=mask, - # ) - # check_2 = tmp >= 0.0 - # return (check_1 & check_2).all() - - # def check_production_eq_strict(self): - # return ( - # (np.isclose(self.production, self.X_0)) - # | np.greater(self.production, self.X_0) - # ).all() - - # def check_production_eq_soft( - # self, current_temporal_unit: int, period: int = 10 - # ) -> bool: - # return self.check_monotony( - # self.production_evolution, current_temporal_unit, period - # ) - - # def check_stocks_monotony( - # self, current_temporal_unit: int, period: int = 10 - # ) -> bool: - # return self.check_monotony(self.stocks_evolution, current_temporal_unit, period) - - # def check_initial_equilibrium(self) -> bool: - # return np.allclose(self.production, self.X_0) and np.allclose( - # self.matrix_stock, self.matrix_stock_0 - # ) - - # def check_equilibrium_soft(self, current_temporal_unit: int): - # return ( - # self.check_stock_increasing(current_temporal_unit) - # and self.check_production_eq_strict - # ) - - # def check_equilibrium_monotony( - # self, current_temporal_unit: int, period: int = 10 - # ) -> bool: - # return self.check_production_eq_soft( - # current_temporal_unit, period - # ) and self.check_stocks_monotony(current_temporal_unit, period) - - # @staticmethod - # def check_monotony(x, current_temporal_unit: int, period: int = 10) -> bool: - # return np.allclose( - # x[current_temporal_unit], x[current_temporal_unit - period], atol=0.0001 - # ) - - # def check_crash(self, prod_threshold: float = 0.80) -> int: - # """Check for economic crash - - # This method look at the production vector and returns the number of - # industries which production is less than a certain share (default 20%) of the starting - # production. - - # Parameters - # ---------- - # prod_threshold : float, default: 0.8 - # An industry is counted as 'crashed' if its current production is less than its starting production times (1 - `prod_threshold`). - - # """ - # tmp = np.full(self.production.shape, 0.0) - # checker = np.full(self.production.shape, 0.0) - # mask = self.X_0 != 0 - # np.subtract(self.X_0, self.production, out=tmp, where=mask) - # np.divide(tmp, self.X_0, out=checker, where=mask) - # return np.where(checker >= prod_threshold)[0].size - def reset_module( self, ) -> None: """Resets the model to initial state [Deprecated] - This method is currently not functioning. + This method has not been checked extensively. """ - self.kapital_lost = np.zeros(self.production.shape) + self.productive_capital_lost = np.zeros(self.production.shape) self.overprod = np.full( (self.n_regions * self.n_sectors), self.overprod_base, dtype=np.float64 ) @@ -1481,51 +1542,18 @@ def reset_module( self.intmd_demand = self.Z_0.copy() self.final_demand = self.Y_0.copy() self.final_demand_not_met = np.zeros(self.Y_0.shape) - self.rebuilding_demand = None + self.rebuilding_demand = np.array([]) self.in_shortage = False self.had_shortage = False self._prod_delta_type = None - self._indus_rebuild_demand_tot = None - self._house_rebuild_demand_tot = None - self._indus_rebuild_demand = None - self._house_rebuild_demand = None - self._tot_rebuild_demand = None - self._prod_cap_delta_kapital = None - self._prod_cap_delta_arbitrary = None - self._prod_cap_delta_tot = None - - def update_params(self, new_params: dict) -> None: - """Update the parameters of the model. - - Replace each parameters with given new ones. - - .. warning:: - Be aware this method calls :meth:`~boario.model_base.reset_record_files`, which resets the memmap files located in the results directory ! - - Parameters - ---------- - new_params : dict - Dictionary of new parameters to use. - - """ - logger.warning("This method is quite probably deprecated") - self.n_temporal_units_by_step = new_params["n_temporal_units_by_step"] - self.iotable_year_to_temporal_unit_factor = new_params[ - "year_to_temporal_unit_factor" - ] - self.rebuild_tau = new_params["rebuild_tau"] - self.overprod_max = new_params["alpha_max"] - self.overprod_tau = new_params["alpha_tau"] - self.overprod_base = new_params["alpha_base"] - if self.records_storage != pathlib.Path( - new_params["output_dir"] + "/" + new_params["results_storage"] - ): - self.records_storage = pathlib.Path( - new_params["output_dir"] + "/" + new_params["results_storage"] - ) - self.reset_record_files( - new_params["n_temporal_units_to_sim"], new_params["register_stocks"] - ) + self._indus_rebuild_demand_tot = np.array([]) + self._house_rebuild_demand_tot = np.array([]) + self._indus_rebuild_demand = np.array([]) + self._house_rebuild_demand = np.array([]) + self._tot_rebuild_demand = np.array([]) + self._prod_cap_delta_productive_capital = np.array([]) + self._prod_cap_delta_arbitrary = np.array([]) + self._prod_cap_delta_tot = np.array([]) def write_index(self, index_file: str | pathlib.Path) -> None: """Write the indexes of the different dataframes of the model in a json file. @@ -1547,7 +1575,7 @@ def write_index(self, index_file: str | pathlib.Path) -> None: indexes = { "regions": list(self.regions), "sectors": list(self.sectors), - "fd_cat": list(self.fd_cat), + "fd_cat": list(self.final_demand_cat), "n_sectors": self.n_sectors, "n_regions": self.n_regions, "n_industries": self.n_sectors * self.n_regions, @@ -1555,10 +1583,10 @@ def write_index(self, index_file: str | pathlib.Path) -> None: if isinstance(index_file, str): index_file = pathlib.Path(index_file) index_file.parent.mkdir(parents=True, exist_ok=True) - with index_file.open("w") as f: - json.dump(indexes, f) + with index_file.open("w") as ffile: + json.dump(indexes, ffile) - def change_inv_duration(self, new_dur: int, old_dur: Optional[int] = None) -> None: + def change_inv_duration(self, new_dur, old_dur=None) -> None: # replace this method by a property ! if old_dur is None: old_dur = self.main_inv_dur diff --git a/boario/simulation.py b/boario/simulation.py index 0b4e4d7..43b045b 100644 --- a/boario/simulation.py +++ b/boario/simulation.py @@ -32,14 +32,21 @@ from typing import Optional, Union import numpy as np +import pandas as pd import progressbar from boario import DEBUGFORMATTER from boario import logger -from boario.event import * -from boario.extended_models import * +from boario.event import ( + Event, + EventArbitraryProd, + EventKapitalDestroyed, + EventKapitalRebuild, + EventKapitalRecover, +) +from boario.extended_models import ARIOPsiModel from boario.model_base import ARIOBaseModel -from boario.utils.misc import CustomNumpyEncoder, TempMemmap +from boario.utils.misc import CustomNumpyEncoder, TempMemmap, sizeof_fmt, print_summary __all__ = ["Simulation"] @@ -63,43 +70,58 @@ class Simulation: "rebuild_prod", "inputs_stocks", "limiting_inputs", - "kapital_to_recover", + "productive_capital_to_recover", ] __file_save_array_specs = { "production_realised": ( "float64", - "production_evolution", + "_production_evolution", "industries", np.nan, ), "production_capacity": ( "float64", - "production_cap_evolution", + "_production_cap_evolution", + "industries", + np.nan, + ), + "final_demand": ("float64", "_final_demand_evolution", "industries", np.nan), + "intermediate_demand": ( + "float64", + "_io_demand_evolution", + "industries", + np.nan, + ), + "rebuild_demand": ( + "float64", + "_rebuild_demand_evolution", + "industries", + np.nan, + ), + "overproduction": ( + "float64", + "_overproduction_evolution", "industries", np.nan, ), - "final_demand": ("float64", "final_demand_evolution", "industries", np.nan), - "intermediate_demand": ("float64", "io_demand_evolution", "industries", np.nan), - "rebuild_demand": ("float64", "rebuild_demand_evolution", "industries", np.nan), - "overproduction": ("float64", "overproduction_evolution", "industries", np.nan), "final_demand_unmet": ( "float64", - "final_demand_unmet_evolution", + "_final_demand_unmet_evolution", "industries", np.nan, ), "rebuild_prod": ( "float64", - "rebuild_production_evolution", + "_rebuild_production_evolution", "industries", np.nan, ), - "inputs_stocks": ("float64", "inputs_evolution", "stocks", np.nan), - "limiting_inputs": ("byte", "limiting_inputs_evolution", "stocks", -1), - "kapital_to_recover": ( + "inputs_stocks": ("float64", "_inputs_evolution", "stocks", np.nan), + "limiting_inputs": ("byte", "_limiting_inputs_evolution", "stocks", -1), + "productive_capital_to_recover": ( "float64", - "regional_sectoral_kapital_destroyed_evolution", + "_regional_sectoral_productive_capital_destroyed_evolution", "industries", np.nan, ), @@ -117,7 +139,7 @@ def __init__( save_index: bool = False, save_records: list | str = [], boario_output_dir: str | pathlib.Path = tempfile.mkdtemp(prefix="boario"), - results_dir_name: str = "results", + results_dir_name: Optional[str] = None, ) -> None: """A Simulation instance can be initialized with the following parameters: @@ -135,17 +157,36 @@ def __init__( An optional list of events to run the simulation with [WIP]. separate_sims : bool, default False Whether to run each event separately or during the same simulation [WIP]. + save_events: bool, default False + If True, saves a json file of the list of events to simulate when starting the model loop. + save_params: bool, default False + If True, saves a json file of the parameters used when starting the model loop. + save_index: bool, default False + If True, saves a json file of the list of industries in the model when starting the model loop (convenience). + save_records: list | str, default [] + The list of simulation variable records to save in corresponding files. boario_output_dir : str | pathlib.Path - An optional directory where to save files generated by the simulation. + An optional directory where to save files generated by the simulation. Defaults to a temporary directory prefixed by "boario". results_dir_name : str, default 'results' The name of the folder where simulation results will be stored. - Examples - -------- + """ + self.output_dir = pathlib.Path(boario_output_dir) + """pathlib.Path, optional: Optional path to the directory where output are stored.""" + self.results_storage = ( + self.output_dir.resolve() + if not results_dir_name + else self.output_dir.resolve() / results_dir_name + ) + """str: Name of the folder in `output_dir` where the results will be stored if saved.""" - See #add link to example page. + if not self.results_storage.exists(): + self.results_storage.mkdir(parents=True) - """ + tmp = logging.FileHandler(self.results_storage / "simulation.log") + tmp.setLevel(logging.INFO) + tmp.setFormatter(DEBUGFORMATTER) + logger.addHandler(tmp) if events_list is None: events_list = [] @@ -154,25 +195,27 @@ def __init__( self._save_params = save_params self._save_index = save_index self._register_stocks = register_stocks - self.output_dir = pathlib.Path(boario_output_dir) - """pathlib.Path, optional: Optional path to the directory where output are stored.""" + + # Pre-init record variables + self._production_evolution = np.array([]) + self._production_cap_evolution = np.array([]) + self._final_demand_evolution = np.array([]) + self._io_demand_evolution = np.array([]) + self._rebuild_demand_evolution = np.array([]) + self._overproduction_evolution = np.array([]) + self._final_demand_unmet_evolution = np.array([]) + self._rebuild_production_evolution = np.array([]) + self._inputs_evolution = np.array([]) + self._limiting_inputs_evolution = np.array([]) + self._regional_sectoral_productive_capital_destroyed_evolution = np.array([]) if save_records != [] or save_events or save_params or save_index: self.output_dir.resolve().mkdir(parents=True, exist_ok=True) if save_records != []: - self.results_storage = self.output_dir.resolve() / results_dir_name if not self.results_storage.exists(): self.results_storage.mkdir() - self.results_storage = ( - pathlib.Path(self.output_dir).resolve() / results_dir_name - ) - """str: Name of the folder in `output_dir` where the results will be stored if saved.""" - - if not self.results_storage.exists(): - self.results_storage.mkdir(parents=True) - self.model = model """Union[ARIOBaseModel, ARIOPsiModel] : The model to run the simulation with. See :class:`~boario.model_base.ARIOBaseModel`.""" @@ -184,6 +227,13 @@ def __init__( """list[Event]: A list containing all events that are happening at the current timestep of the simulation.""" self.events_timings = set() + if ( + not isinstance(n_temporal_units_to_sim, (int, np.integer)) + or n_temporal_units_to_sim <= 0 + ): + raise ValueError( + f"n_temporal_units_to_sim should be a positive integer (got {n_temporal_units_to_sim} of type {type(n_temporal_units_to_sim)})" + ) self.n_temporal_units_to_sim = n_temporal_units_to_sim """int: The total number of `temporal_units` to simulate.""" @@ -206,7 +256,7 @@ def __init__( # RECORDS FILES self.records_storage: pathlib.Path = self.results_storage / "records" - """Place where records are stored if stored""" + """pathlib.Path: Place where records are stored if stored""" self._files_to_record = [] @@ -227,15 +277,14 @@ def __init__( f"{impossible_records} are not possible records ({self.__possible_records})" ) logger.info(f"Will save {save_records} records") - logger.info("Records storage is: {}".format(self.records_storage)) + logger.info("Records storage is: {self.records_storage}") self.records_storage.mkdir(parents=True, exist_ok=True) self._save_index = True self._save_events = True self._save_params = True - self.init_records(save_records, register_stocks) + self._init_records(save_records) - Event.temporal_unit_range = self.n_temporal_units_to_sim self.params_dict = { "n_temporal_units_to_sim": self.n_temporal_units_to_sim, "output_dir": str(self.output_dir) @@ -262,8 +311,12 @@ def __init__( """dict: A dictionary saving the parameters the simulation was run with.""" logger.info("Initialized !") + formatted_params_dict = { + key: print_summary(value) if key == "inventory_restoration_tau" else value + for key, value in self.params_dict.items() + } logger.info( - "Simulation parameters:\n{}".format(pformat(self.params_dict, compact=True)) + f"Simulation parameters:\n{pformat(formatted_params_dict, compact=True)}" ) def loop(self, progress: bool = True): @@ -282,21 +335,12 @@ def loop(self, progress: bool = True): If True, shows a progress bar of the loop in the console. """ logger.info( - "Starting model loop for at most {} steps".format( - self.n_temporal_units_to_sim // self.model.n_temporal_units_by_step + 1 - ) + f"Starting model loop for at most {self.n_temporal_units_to_sim // self.model.n_temporal_units_by_step + 1} steps" ) logger.info( - "One step is {}/{} of a year".format( - self.model.n_temporal_units_by_step, - self.model.iotable_year_to_temporal_unit_factor, - ) + f"One step is {self.model.n_temporal_units_by_step}/{self.model.iotable_year_to_temporal_unit_factor} of a year" ) - tmp = logging.FileHandler(self.results_storage / "simulation.log") - tmp.setLevel(logging.INFO) - tmp.setFormatter(DEBUGFORMATTER) - logger.addHandler(tmp) - logger.info("Events : {}".format(self.all_events)) + logger.info("Events : {self.all_events}") run_range = range( 0, @@ -310,9 +354,9 @@ def loop(self, progress: bool = True): ) with ( pathlib.Path(self.results_storage) / "jsons" / "simulated_events.json" - ).open("w") as f: + ).open("w") as ffile: event_dicts = [ev.event_dict for ev in self.all_events] - json.dump(event_dicts, f, indent=4, cls=CustomNumpyEncoder) + json.dump(event_dicts, ffile, indent=4, cls=CustomNumpyEncoder) if progress: widgets = [ "Processed: ", @@ -348,7 +392,7 @@ def loop(self, progress: bool = True): if step_res == 1: self.has_crashed = True logger.warning( - f"""Economy seems to have crashed. + f"""Economy or model seems to have crashed. - At step : {self.current_temporal_unit} """ ) @@ -361,8 +405,8 @@ def loop(self, progress: bool = True): ) break - if self._files_to_record != []: - self.flush_memmaps() + if self._files_to_record: + self._flush_memmaps() if self._save_index: self.model.write_index(self.results_storage / "jsons" / "indexes.json") @@ -373,14 +417,14 @@ def loop(self, progress: bool = True): if self._save_params: with ( pathlib.Path(self.results_storage) / "jsons" / "simulated_params.json" - ).open("w") as f: - json.dump(self.params_dict, f, indent=4, cls=CustomNumpyEncoder) + ).open("w") as ffile: + json.dump(self.params_dict, ffile, indent=4, cls=CustomNumpyEncoder) with ( pathlib.Path(self.results_storage) / "jsons" / "equilibrium_checks.json" - ).open("w") as f: + ).open("w") as ffile: json.dump( {str(k): v for k, v in self.equi.items()}, - f, + ffile, indent=4, cls=CustomNumpyEncoder, ) @@ -399,7 +443,7 @@ def next_step( First it checks if an event is planned to occur at the current step and if so, shocks the model with the - corresponding event. Then it : + corresponding event. Then: 0) If at least one step elapsed, it computes the new overproduction vector for the next step (using :meth:`~boario.model_base.ARIOBaseModel.calc_overproduction`) @@ -407,7 +451,7 @@ def next_step( 2) Distribute the `realised` production towards the different demands (intermediate, final, rebuilding) and compute the changes in the inputs stock matrix (see :meth:`~boario.model_base.ARIOBaseModel.distribute_production`) - Note that it is during this step that the model checks if an event is completely rebuild/recovered. + Note that it is during this step that the model checks if an event is completely rebuild/recovered and removes it from the list the case being. 3) Computes the orders matrix (i.e. the intermediate demand) for the next step (see :meth:`~boario.model_base.ARIOBaseModel.calc_orders`) @@ -423,86 +467,98 @@ def next_step( The minimum number of failing regions required to consider economy has crashed. """ - if min_steps_check is None: - min_steps_check = self.n_temporal_units_to_sim // 5 - if min_failing_regions is None: - min_failing_regions = self.model.n_regions * self.model.n_sectors // 3 - - # Check if there are new events to add, - # if some happening events can start rebuilding (if rebuildable), - # and updates the internal model production_cap decrease and rebuild_demand - self.check_happening_events() - - if "inputs_evolution" in self._files_to_record: - self.write_stocks() - - if self.current_temporal_unit > 1: - self.model.calc_overproduction() - - if "overproduction_evolution" in self._files_to_record: - self.write_overproduction() - if "rebuild_demand_evolution" in self._files_to_record: - self.write_rebuild_demand() - if "final_evolution" in self._files_to_record: - self.write_final_demand() - if "io_demand_evolution" in self._files_to_record: - self.write_io_demand() - - constraints = self.model.calc_production(self.current_temporal_unit) - - if "limiting_inputs_evolution" in self._files_to_record: - self.write_limiting_stocks(constraints) - if "production_evolution" in self._files_to_record: - self.write_production() - if "production_cap_evolution" in self._files_to_record: - self.write_production_max() - if "regional_sectoral_kapital_destroyed_evolution" in self._files_to_record: - self.write_kapital_lost() - try: - rebuildable_events = [ - ev - for ev in self.currently_happening_events - if isinstance(ev, EventKapitalRebuild) and ev.rebuildable - ] - events_to_remove = self.model.distribute_production( - rebuildable_events, self.scheme - ) - if "final_demand_unmet_evolution" in self._files_to_record: - self.write_final_demand_unmet() - if "rebuild_production_evolution" in self._files_to_record: - self.write_rebuild_prod() - except RuntimeError as e: - logger.exception("This exception happened:", e) - return 1 - events_to_remove = events_to_remove + [ - ev for ev in self.currently_happening_events if ev.over - ] - if events_to_remove: - self.currently_happening_events = [ - e for e in self.currently_happening_events if e not in events_to_remove + if min_steps_check is None: + min_steps_check = self.n_temporal_units_to_sim // 5 + if min_failing_regions is None: + min_failing_regions = self.model.n_regions * self.model.n_sectors // 3 + + # Check if there are new events to add, + # if some happening events can start rebuilding (if rebuildable), + # and updates the internal model production_cap decrease and rebuild_demand + self._check_happening_events() + + if "_inputs_evolution" in self._files_to_record: + self._write_stocks() + + # 0) + if self.current_temporal_unit > 1: + self.model.calc_overproduction() + + if "_overproduction_evolution" in self._files_to_record: + self._write_overproduction() + if "_rebuild_demand_evolution" in self._files_to_record: + self._write_rebuild_demand() + if "_final_evolution" in self._files_to_record: + self._write_final_demand() + if "_io_demand_evolution" in self._files_to_record: + self._write_io_demand() + + # 1) + constraints = self.model.calc_production(self.current_temporal_unit) + + if "_limiting_inputs_evolution" in self._files_to_record: + self._write_limiting_stocks(constraints) + if "_production_evolution" in self._files_to_record: + self._write_production() + if "_production_cap_evolution" in self._files_to_record: + self._write_production_max() + if ( + "_regional_sectoral_productive_capital_destroyed_evolution" + in self._files_to_record + ): + self._write_productive_capital_lost() + + # 2) + try: + rebuildable_events = [ + ev + for ev in self.currently_happening_events + if isinstance(ev, EventKapitalRebuild) and ev.rebuildable + ] + events_to_remove = self.model.distribute_production( + rebuildable_events, self.scheme + ) + if "_final_demand_unmet_evolution" in self._files_to_record: + self._write_final_demand_unmet() + if "_rebuild_production_evolution" in self._files_to_record: + self._write_rebuild_prod() + except RuntimeError: + logger.exception("An exception happened:") + self.model.matrix_stock.dump( + self.results_storage / "matrix_stock_dump.pkl" + ) + logger.error( + f"Negative values in the stocks, matrix has been dumped in the results dir : \n {self.results_storage / 'matrix_stock_dump.pkl'}" + ) + return 1 + events_to_remove = events_to_remove + [ + ev for ev in self.currently_happening_events if ev.over ] - for e in events_to_remove: - if isinstance(e, EventKapitalDestroyed): - logger.info( - "Temporal_Unit : {} ~ Event named {} that occured at {} in {} for {} damages is completely rebuilt/recovered".format( - self.current_temporal_unit, - e.name, - e.occurrence, - e.aff_regions, - e.total_kapital_destroyed, + if events_to_remove: + self.currently_happening_events = [ + e + for e in self.currently_happening_events + if e not in events_to_remove + ] + for evnt in events_to_remove: + if isinstance(evnt, EventKapitalDestroyed): + logger.info( + f"""Temporal_Unit : {self.current_temporal_unit} ~ Event named {evnt.name} that occured at {evnt.occurrence} in {evnt.aff_regions.to_list()} for {evnt.total_productive_capital_destroyed} damages is completely rebuilt/recovered""" ) - ) - self.model.calc_orders() + self.model.calc_orders() - n_checks = self.current_temporal_unit // check_period - if n_checks > self._n_checks: - self.check_equilibrium(n_checks) - self._n_checks += 1 + n_checks = self.current_temporal_unit // check_period + if n_checks > self._n_checks: + self.check_equilibrium(n_checks) + self._n_checks += 1 - self.current_temporal_unit += self.model.n_temporal_units_by_step - return 0 + self.current_temporal_unit += self.model.n_temporal_units_by_step + return 0 + except Exception: + logger.exception(f"The following exception happened:") + return 1 def check_equilibrium(self, n_checks: int): """Checks the status of production, stocks and rebuilding demand. @@ -545,43 +601,18 @@ def check_equilibrium(self, n_checks: int): (n_checks, self.current_temporal_unit, "rebuilding") ] = "not finished" - def read_events_from_list(self, events_list: list[dict]): - raise NotImplementedError("I have to redo this one") - - # """Import a list of events (as a list of dictionaries) into the model. - - # Also performs various checks on the events to avoid badly written events. - # See :ref:`How to define Events ` to understand how to write events dictionaries or JSON files. - - # Parameters - # ---------- - # events_list : - # List of events as dictionaries. - - # """ - # logger.info("Reading events from given list and adding them to the model") - # if not isinstance(events_list, list): - # if isinstance(events_list, dict): - # raise TypeError( - # "read_events_from_list() takes a list of event dicts as an argument, not a single event dict, you might want to use read_event(your_event) instead." - # ) - # else: - # raise TypeError( - # "read_events_from_list() takes a list of event dicts as an argument, not a {}".format( - # type(events_list) - # ) - # ) - # for ev_dic in events_list: - # self.read_event(ev_dic) - - def read_event(self, ev_dic: dict): - raise NotImplementedError("I have to redo this one") - - # if ev_dic["aff_sectors"] == "all": - # ev_dic["aff_sectors"] = list(self.model.sectors) - # ev = Event(ev_dic) - # self.all_events.append(ev) - # self.events_timings.add(ev.occurrence) + def add_events(self, events: list[Event]): + """Add a list of events to the simulation. + + Parameters + ---------- + events : list[Event] + The events to add. + """ + if not isinstance(events, list): + raise TypeError(f"list[Event] expected, {type(events)} received.") + for ev in events: + self.add_event(ev) def add_event(self, ev: Event): """Add an event to the simulation. @@ -591,7 +622,8 @@ def add_event(self, ev: Event): ev : Event The event to add. """ - + if not isinstance(ev, Event): + raise ValueError(f"Event expected, {type(ev)} received.") self.all_events.append(ev) self.events_timings.add(ev.occurrence) @@ -603,7 +635,12 @@ def reset_sim_with_same_events(self): self._n_checks = 0 self.n_temporal_units_simulated = 0 self.has_crashed = False - self.reset_records() + self.equi = { + (int(0), int(0), "production"): "equi", + (int(0), int(0), "stocks"): "equi", + (int(0), int(0), "rebuilding"): "equi", + } + self._reset_records() self.model.reset_module() def reset_sim_full(self): @@ -615,29 +652,6 @@ def reset_sim_full(self): self.currently_happening_events = [] self.events_timings = set() - def update_params(self, new_params: dict): - """Update the parameters of the model. - - Replace the ``params`` attribute with ``new_params`` and logs the update. - This method also checks if the directory specified to save the results exists and create it otherwise. - - .. warning:: - Be aware this method calls :meth:`~boario.model_base.ARIOBaseModel.update_params`, which resets the memmap files located in the results directory ! - - Parameters - ---------- - new_params : dict - New dictionnary of parameters to use. - - """ - raise NotImplementedError("To fix") - logger.info("Updating model parameters") - self.params_dict = new_params - results_storage = pathlib.Path(self.results_storage) - if not results_storage.exists(): - results_storage.mkdir(parents=True) - self.model.update_params(self.params_dict) - def write_index(self, index_file: Union[str, pathlib.Path]): """Write the index of the dataframes used in the model in a json file. @@ -651,39 +665,7 @@ def write_index(self, index_file: Union[str, pathlib.Path]): """ self.model.write_index(index_file) - def read_events(self, events_file: Union[str, pathlib.Path]): - """Read events from a json file. - - Parameters - ---------- - events_file : - path to a json file - - Raises - ------ - FileNotFoundError - If file does not exist - - """ - raise NotImplementedError("I have to redo this one") - # logger.info( - # "Reading events from {} and adding them to the model".format(events_file) - # ) - # if isinstance(events_file, str): - # events_file = pathlib.Path(events_file) - # elif not isinstance(events_file, pathlib.Path): - # raise TypeError("Given index file is not an str or a Path") - # if not events_file.exists(): - # raise FileNotFoundError("This file does not exist: ", events_file) - # else: - # with events_file.open("r") as f: - # events = json.load(f) - # if isinstance(events, list): - # self.read_events_from_list(events) - # else: - # self.read_events_from_list([events]) - - def check_happening_events(self) -> None: + def _check_happening_events(self) -> None: """Updates the status of all events. Check the `all_events` attribute and `current_temporal_unit` and @@ -691,112 +673,110 @@ def check_happening_events(self) -> None: rebuild/recover) """ - for ev in self.all_events: - if not ev.happened: + for evnt in self.all_events: + if not evnt.happened: if ( (self.current_temporal_unit - self.model.n_temporal_units_by_step) - <= ev.occurrence + <= evnt.occurrence <= self.current_temporal_unit ): logger.info( - "Temporal_Unit : {} ~ Shocking model with new event".format( - self.current_temporal_unit - ) + f"Temporal_Unit : {self.current_temporal_unit} ~ Shocking model with new event" ) - logger.info("Affected regions are : {}".format(ev.aff_regions)) - ev.happened = True - self.currently_happening_events.append(ev) - for ev in self.currently_happening_events: - if isinstance(ev, EventKapitalRebuild): - ev.rebuildable = self.current_temporal_unit - if isinstance(ev, EventKapitalRecover): - ev.recoverable = self.current_temporal_unit - if ev.recoverable: - ev.recovery(self.current_temporal_unit) + logger.info(f"Affected regions are : {evnt.aff_regions.to_list()}") + evnt.happened = True + self.currently_happening_events.append(evnt) + for evnt in self.currently_happening_events: + if isinstance(evnt, EventKapitalRebuild): + evnt.rebuildable = self.current_temporal_unit + if isinstance(evnt, (EventKapitalRecover, EventArbitraryProd)): + evnt.recoverable = self.current_temporal_unit + if evnt.recoverable: + evnt.recovery(self.current_temporal_unit) self.model.update_system_from_events(self.currently_happening_events) - def write_production(self) -> None: + def _write_production(self) -> None: """Saves the current production vector to the memmap.""" - self.production_evolution[self.current_temporal_unit] = self.model.production + self._production_evolution[self.current_temporal_unit] = self.model.production - def write_rebuild_prod(self) -> None: + def _write_rebuild_prod(self) -> None: """Saves the current rebuilding production vector to the memmap.""" logger.debug( - f"self.rebuild_production_evolution shape : {self.rebuild_production_evolution.shape}, self.model.rebuild_prod shape : {self.model.rebuild_prod.shape}" + f"self._rebuild_production_evolution shape : {self._rebuild_production_evolution.shape}, self.model.rebuild_prod shape : {self.model.rebuild_prod.shape}" ) - self.rebuild_production_evolution[ + self._rebuild_production_evolution[ self.current_temporal_unit ] = self.model.rebuild_prod - def write_kapital_lost(self) -> None: - """Saves the current remaining kapital to rebuild vector to the memmap.""" - self.regional_sectoral_kapital_destroyed_evolution[ + def _write_productive_capital_lost(self) -> None: + """Saves the current remaining productive_capital to rebuild vector to the memmap.""" + self._regional_sectoral_productive_capital_destroyed_evolution[ self.current_temporal_unit - ] = self.model.kapital_lost + ] = self.model.productive_capital_lost - def write_production_max(self) -> None: + def _write_production_max(self) -> None: """Saves the current production capacity vector to the memmap.""" - self.production_cap_evolution[ + self._production_cap_evolution[ self.current_temporal_unit ] = self.model.production_cap - def write_io_demand(self) -> None: + def _write_io_demand(self) -> None: """Saves the current (total per industry) intermediate demand vector to the memmap.""" - self.io_demand_evolution[ + self._io_demand_evolution[ self.current_temporal_unit ] = self.model.matrix_orders.sum(axis=1) - def write_final_demand(self) -> None: + def _write_final_demand(self) -> None: """Saves the current (total per industry) final demand vector to the memmap.""" - self.final_demand_evolution[ + self._final_demand_evolution[ self.current_temporal_unit ] = self.model.final_demand.sum(axis=1) - def write_rebuild_demand(self) -> None: + def _write_rebuild_demand(self) -> None: """Saves the current (total per industry) rebuilding demand vector to the memmap.""" to_write = np.full(self.model.n_regions * self.model.n_sectors, 0.0) - if (r_dem := self.model.tot_rebuild_demand) is not None: - self.rebuild_demand_evolution[self.current_temporal_unit] = r_dem # type: ignore + if len(r_dem := self.model.tot_rebuild_demand) > 0: + self._rebuild_demand_evolution[self.current_temporal_unit] = r_dem # type: ignore else: - self.rebuild_demand_evolution[self.current_temporal_unit] = to_write # type: ignore - - def write_limiting_stocks(self, limiting_stock: np.ndarray) -> None: - """Saves the current limiting inputs matrix to the memmap.""" - self.limiting_inputs_evolution[self.current_temporal_unit] = limiting_stock # type: ignore + self._rebuild_demand_evolution[self.current_temporal_unit] = to_write # type: ignore - def write_overproduction(self) -> None: + def _write_overproduction(self) -> None: """Saves the current overproduction vector to the memmap.""" - self.overproduction_evolution[self.current_temporal_unit] = self.model.overprod + self._overproduction_evolution[self.current_temporal_unit] = self.model.overprod - def write_final_demand_unmet(self) -> None: + def _write_final_demand_unmet(self) -> None: """Saves the unmet final demand (for this step) vector to the memmap.""" - self.final_demand_unmet_evolution[ + self._final_demand_unmet_evolution[ self.current_temporal_unit ] = self.model.final_demand_not_met - def write_stocks(self) -> None: + def _write_stocks(self) -> None: """Saves the current inputs stock matrix to the memmap.""" - self.inputs_evolution[self.current_temporal_unit] = self.model.matrix_stock + self._inputs_evolution[self.current_temporal_unit] = self.model.matrix_stock - def write_limiting_stocks(self, limiting_stock: np.ndarray) -> None: + def _write_limiting_stocks(self, limiting_stock: np.ndarray) -> None: """Saves the current limiting inputs matrix to the memmap.""" - self.limiting_inputs_evolution[self.current_temporal_unit] = limiting_stock # type: ignore + self._limiting_inputs_evolution[self.current_temporal_unit] = limiting_stock # type: ignore - def flush_memmaps(self) -> None: + def _flush_memmaps(self) -> None: """Saves files to record""" - for at in self._files_to_record: - if not hasattr(self, at): + for attr in self._files_to_record: + if not hasattr(self, attr): raise RuntimeError( - f"{at} should be a member yet it isn't. This shouldn't happen." + f"{attr} should be a member yet it isn't. This shouldn't happen." ) else: - getattr(self, at).flush() + getattr(self, attr).flush() - def init_records(self, save_records, register_stocks): + def _init_records(self, save_records): for rec in self.__possible_records: - if rec == "input_stocks" and not register_stocks: - pass + if rec == "inputs_stocks" and not self._register_stocks: + logger.debug("Will not save inputs stocks") else: + if rec == "inputs_stocks": + logger.info( + f"Simulation will save inputs stocks. Estimated size is {sizeof_fmt(self.n_temporal_units_to_sim * self.model.n_sectors * self.model.n_sectors * self.model.n_regions * 64)}" + ) save = rec in save_records filename = rec dtype, attr_name, shapev, fillv = self.__file_save_array_specs[rec] @@ -824,14 +804,181 @@ def init_records(self, save_records, register_stocks): self._files_to_record.append(attr_name) setattr(self, attr_name, memmap_array) - def reset_records( + def _reset_records( self, ): for rec in self.__possible_records: - dtype, attr_name, shapev, fillv = self.__file_save_array_specs[rec] + _, attr_name, _, fillv = self.__file_save_array_specs[rec] if rec == "input_stocks" and not self._register_stocks: pass else: memmap_array = getattr(self, attr_name) memmap_array.fill(fillv) setattr(self, attr_name, memmap_array) + + @property + def production_realised(self) -> pd.DataFrame: + """Return the evolution of the production realised by each industry (region,sector) + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the production realised, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._production_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def production_capacity(self) -> pd.DataFrame: + """Return the evolution of the production capacity of each industry (region,sector) + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the production capacity, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._production_cap_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def final_demand(self) -> pd.DataFrame: + """Return the evolution of final demand asked of each industry. + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the final demand asked, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._final_demand_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def intermediate_demand(self) -> pd.DataFrame: + """Return the evolution of intermediate demand asked of each industry (Total orders). + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the intermediate demand asked, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._io_demand_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def rebuild_demand(self) -> pd.DataFrame: + """Return the evolution of rebuild demand asked of each industry (Total orders). + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the rebuild demand asked, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._rebuild_demand_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def overproduction(self) -> pd.DataFrame: + """Return the evolution of the overproduction factor of each industry (region,sector) + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the overproduction factor, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._overproduction_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def final_demand_unmet(self) -> pd.DataFrame: + """Return the evolution of the final demand that could not be answered by industries. + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the final demand not met, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._final_demand_unmet_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def rebuild_prod(self) -> pd.DataFrame: + """Return the production allocated for the rebuilding demand by each industry (region,sector). + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the production allocated, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._rebuild_production_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") + + @property + def inputs_stocks(self) -> pd.DataFrame: + """Return the evolution of the inventory amount of each input for each industry (region,sector). Not this is not available if "record_stocks" is not set to True, + as the DataFrame can be quite large for "classic" MRIOTs. + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the amount in inventory, the columns are the industries + and the index are the step and input considered (MultiIndex). + """ + return pd.DataFrame( + self._inputs_evolution.reshape( + self.n_temporal_units_to_sim * self.model.n_sectors, -1 + ), + columns=self.model.industries, + copy=True, + index=pd.MultiIndex.from_product( + [list(range(self.n_temporal_units_to_sim)), self.model.sectors], + names=["step", "input"], + ), + ) + + @property + def limiting_inputs(self) -> pd.DataFrame: + """Return the evolution of the inputs lacking for each industry (region,sector) + + Returns: + pd.DataFrame: A pandas DataFrame where the value is a boolean set to 1 if considered input constrains production, the columns are the industries + and the index are the step and input considered (MultiIndex). + """ + return pd.DataFrame( + self._limiting_inputs_evolution.reshape( + self.n_temporal_units_to_sim * self.model.n_sectors, -1 + ), + columns=self.model.industries, + copy=True, + index=pd.MultiIndex.from_product( + [list(range(self.n_temporal_units_to_sim)), self.model.sectors], + names=["step", "input"], + ), + ) + + @property + def productive_capital_to_recover(self) -> pd.DataFrame: + """Return the evolution of remaining capital destroyed/to recover for each industry (region,sector) if it exists. + + Returns: + pd.DataFrame: A pandas DataFrame where the value is the amount of capital (still) destroyed, the columns are the industries + and the index is the step considered. + """ + return pd.DataFrame( + self._regional_sectoral_productive_capital_destroyed_evolution, + columns=self.model.industries, + copy=True, + ).rename_axis("step") diff --git a/boario/utils/misc.py b/boario/utils/misc.py index 90f7f74..c3c77d5 100644 --- a/boario/utils/misc.py +++ b/boario/utils/misc.py @@ -15,6 +15,7 @@ # along with this program. If not, see . from collections.abc import Iterable +import textwrap import pymrio import numpy import json @@ -118,3 +119,42 @@ def lexico_reindex(mrio: pymrio.IOSystem) -> pymrio.IOSystem: mrio.A = mrio.A.reindex(sorted(mrio.A.columns), axis=1) return mrio + + +def sizeof_fmt(num, suffix="B"): + for unit in ["", "Ki", "Mi", "Gi", "Ti", "Pi", "Ei", "Zi"]: + if abs(num) < 1024.0: + return f"{num:3.1f}{unit}{suffix}" + num /= 1024.0 + return f"{num:.1f}Yi{suffix}" + + +def print_summary(my_list): + if my_list: + current_element = None + current_count = 0 + summary = [] + for element in my_list: + if element != current_element: + if current_element is not None: + if current_count == 1: + summary.append(str(current_element)) + else: + summary.append(f"{current_element} (x {current_count})") + current_element = element + current_count = 1 + else: + current_count += 1 + if current_element is not None: + if current_count == 1: + summary.append(str(current_element)) + else: + summary.append(f"{current_element} (x {current_count})") + total_length = len(my_list) + total_sum = sum(my_list) + summary_string = ( + "[" + ", ".join(summary) + f"] (len: {total_length}, sum: {total_sum})" + ) + return textwrap.wrap(summary_string, width=80) + else: + return "" diff --git a/boario/utils/recovery_functions.py b/boario/utils/recovery_functions.py new file mode 100644 index 0000000..d8775f6 --- /dev/null +++ b/boario/utils/recovery_functions.py @@ -0,0 +1,132 @@ +from typing import Optional +import numpy as np + + +def linear_recovery( + elapsed_temporal_unit: int, + init_impact_stock: np.ndarray, + recovery_time: int, +): + r"""Linear Initial impact recovery function + + Initial impact is entirely recovered when `recovery_time` has passed since event + started recovering + + Parameters + ---------- + + init_impact_stock : float + Initial Initial impact destroyed + + elapsed_temporal_unit : int + Elapsed time since event started recovering + + recovery_time : int + Total time it takes the event to fully recover + + """ + + return init_impact_stock * (1 - (elapsed_temporal_unit / recovery_time)) + + +def convexe_recovery( + elapsed_temporal_unit: int, + init_impact_stock: np.ndarray, + recovery_time: int, +): + r"""Convexe Initial impact recovery function + + Initial impact is recovered with characteristic time `recovery_time`. (This doesn't mean Initial impact is fully recovered after this time !) + This function models a recovery similar as the one happening in the rebuilding case, for the same characteristic time. + + Parameters + ---------- + + init_impact_stock : float + Initial Initial impact destroyed + + elapsed_temporal_unit : int + Elapsed time since event started recovering + + recovery_time : int + Total time it takes the event to fully recover + + """ + + return init_impact_stock * (1 - (1 / recovery_time)) ** elapsed_temporal_unit + + +def convexe_recovery_scaled( + elapsed_temporal_unit: int, + init_impact_stock: np.ndarray, + recovery_time: int, + scaling_factor: float = 4, +): + r"""Convexe Initial impact recovery function (scaled to match other recovery duration) + + Initial impact is mostly recovered (>95% by default for most cases) when `recovery_time` has passed since event + started recovering. + + Parameters + ---------- + + init_impact_stock : float + Initial Initial impact destroyed + + elapsed_temporal_unit : int + Elapsed time since event started recovering + + recovery_time : int + Total time it takes the event to fully recover + + scaling_factor: float + Used to scale the exponent in the function so that Initial impact is mostly rebuilt after `recovery_time`. A value of 4 insure >95% of Initial impact is recovered for a reasonable range of `recovery_time` values. + + """ + + return init_impact_stock * (1 - (1 / recovery_time)) ** ( + scaling_factor * elapsed_temporal_unit + ) + + +def concave_recovery( + elapsed_temporal_unit: int, + init_impact_stock: np.ndarray, + recovery_time: int, + steep_factor: float = 0.000001, + half_recovery_time: Optional[int] = None, +): + r"""Concave (s-shaped) Initial impact recovery function + + Initial impact is mostly (>95% in most cases) recovered when `recovery_time` has passed since event started recovering. + + Parameters + ---------- + + init_impact_stock : float + Initial Initial impact destroyed + + elapsed_temporal_unit : int + Elapsed time since event started recovering + + recovery_time : int + Total time it takes the event to fully recover + + steep_factor: float + This coefficient governs the slope of the central part of the s-shape, smaller values lead to a steeper slope. As such it also affect the percentage of Initial impact rebuilt after `recovery_time` has elapsed. A value of 0.000001 should insure 95% of the initial impact is rebuild for a reasonable range of recovery duration. + + half_recovery_time : int + This can by use to change the time the inflexion point of the s-shape curve is attained. By default it is set to half the recovery duration. + + """ + + if half_recovery_time is None: + tau_h = 2 + else: + tau_h = recovery_time / half_recovery_time + exponent = (np.log(recovery_time) - np.log(steep_factor)) / ( + np.log(recovery_time) - np.log(tau_h) + ) + return (init_impact_stock * recovery_time) / ( + recovery_time + steep_factor * (elapsed_temporal_unit**exponent) + ) diff --git a/docs/build/doctrees/autoapi/boario.doctree b/docs/build/doctrees/autoapi/boario.doctree index a2c26a2..ce16be8 100644 Binary files a/docs/build/doctrees/autoapi/boario.doctree and b/docs/build/doctrees/autoapi/boario.doctree differ diff --git a/docs/build/doctrees/autoapi/boario.event.doctree b/docs/build/doctrees/autoapi/boario.event.doctree index 43a95ea..008760b 100644 Binary files a/docs/build/doctrees/autoapi/boario.event.doctree and b/docs/build/doctrees/autoapi/boario.event.doctree differ diff --git a/docs/build/doctrees/autoapi/boario.extended_models.doctree b/docs/build/doctrees/autoapi/boario.extended_models.doctree index a3f6973..dbc1989 100644 Binary files a/docs/build/doctrees/autoapi/boario.extended_models.doctree and b/docs/build/doctrees/autoapi/boario.extended_models.doctree differ diff --git a/docs/build/doctrees/autoapi/boario.indicators.doctree b/docs/build/doctrees/autoapi/boario.indicators.doctree index 2de9efa..942488b 100644 Binary files a/docs/build/doctrees/autoapi/boario.indicators.doctree and b/docs/build/doctrees/autoapi/boario.indicators.doctree differ diff --git a/docs/build/doctrees/autoapi/boario.model_base.doctree b/docs/build/doctrees/autoapi/boario.model_base.doctree index 1243129..2af63f1 100644 Binary files a/docs/build/doctrees/autoapi/boario.model_base.doctree and b/docs/build/doctrees/autoapi/boario.model_base.doctree differ diff --git a/docs/build/doctrees/autoapi/boario.simulation.doctree b/docs/build/doctrees/autoapi/boario.simulation.doctree index 64880dd..747c7a0 100644 Binary files a/docs/build/doctrees/autoapi/boario.simulation.doctree and b/docs/build/doctrees/autoapi/boario.simulation.doctree differ diff --git a/docs/build/doctrees/boario-api-reference.doctree b/docs/build/doctrees/boario-api-reference.doctree index 4065ae0..c965b70 100644 Binary files a/docs/build/doctrees/boario-api-reference.doctree and b/docs/build/doctrees/boario-api-reference.doctree differ diff --git a/docs/build/doctrees/boario-events.doctree b/docs/build/doctrees/boario-events.doctree index 1c35aec..06806f2 100644 Binary files a/docs/build/doctrees/boario-events.doctree and b/docs/build/doctrees/boario-events.doctree differ diff --git a/docs/build/doctrees/boario-installation.doctree b/docs/build/doctrees/boario-installation.doctree index c736dc5..e6df321 100644 Binary files a/docs/build/doctrees/boario-installation.doctree and b/docs/build/doctrees/boario-installation.doctree differ diff --git a/docs/build/doctrees/boario-math.doctree b/docs/build/doctrees/boario-math.doctree index 62bc69c..90e69bb 100644 Binary files a/docs/build/doctrees/boario-math.doctree and b/docs/build/doctrees/boario-math.doctree differ diff --git a/docs/build/doctrees/boario-quickstart.doctree b/docs/build/doctrees/boario-quickstart.doctree index 9bbe1dc..a061220 100644 Binary files a/docs/build/doctrees/boario-quickstart.doctree and b/docs/build/doctrees/boario-quickstart.doctree differ diff --git a/docs/build/doctrees/boario-references.doctree b/docs/build/doctrees/boario-references.doctree index 02f6236..64716a7 100644 Binary files a/docs/build/doctrees/boario-references.doctree and b/docs/build/doctrees/boario-references.doctree differ diff --git a/docs/build/doctrees/boario-terminology.doctree b/docs/build/doctrees/boario-terminology.doctree index 48c6263..93d7485 100644 Binary files a/docs/build/doctrees/boario-terminology.doctree and b/docs/build/doctrees/boario-terminology.doctree differ diff --git a/docs/build/doctrees/boario-what-is.doctree b/docs/build/doctrees/boario-what-is.doctree index 3906bce..3b14bd2 100644 Binary files a/docs/build/doctrees/boario-what-is.doctree and b/docs/build/doctrees/boario-what-is.doctree differ diff --git a/docs/build/doctrees/development.doctree b/docs/build/doctrees/development.doctree index 0876d34..e9c5499 100644 Binary files a/docs/build/doctrees/development.doctree and b/docs/build/doctrees/development.doctree differ diff --git a/docs/build/doctrees/environment.pickle b/docs/build/doctrees/environment.pickle index 97ee411..e0ebe57 100644 Binary files a/docs/build/doctrees/environment.pickle and b/docs/build/doctrees/environment.pickle differ diff --git a/docs/build/doctrees/index.doctree b/docs/build/doctrees/index.doctree index 369755c..c8bc38d 100644 Binary files a/docs/build/doctrees/index.doctree and b/docs/build/doctrees/index.doctree differ diff --git a/docs/build/doctrees/release-note.doctree b/docs/build/doctrees/release-note.doctree index 43e902c..e73f61d 100644 Binary files a/docs/build/doctrees/release-note.doctree and b/docs/build/doctrees/release-note.doctree differ diff --git a/docs/build/doctrees/user-guide.doctree b/docs/build/doctrees/user-guide.doctree index b9926b6..d155b55 100644 Binary files a/docs/build/doctrees/user-guide.doctree and b/docs/build/doctrees/user-guide.doctree differ diff --git a/docs/build/html/.buildinfo b/docs/build/html/.buildinfo index 27d8f68..202edc0 100644 --- a/docs/build/html/.buildinfo +++ b/docs/build/html/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: 451af383ec1f62424fe2b190b4f160ff +config: 9b60168294691cd51a01b397a85f5935 tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/docs/build/html/_modules/boario/event.html b/docs/build/html/_modules/boario/event.html index 6136ee2..0ce6f92 100644 --- a/docs/build/html/_modules/boario/event.html +++ b/docs/build/html/_modules/boario/event.html @@ -1,3 +1,1704 @@ +<<<<<<< HEAD + + + + + + + + + + + boario.event — BoARIO documentation + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+
+ + + +
+
+ +
+ + + + + + + + + + + +
+ +
+ + +
+
+ +
+
+ +
+ +
+ + + + +
+ +
+ + +
+
+ + + + + +
+ +

Source code for boario.event

+# BoARIO : The Adaptative Regional Input Output model in python.
+# Copyright (C) 2022  Samuel Juhel
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+from __future__ import annotations
+import abc
+from typing import Callable, Optional, Union, List, Tuple
+import numpy.typing as npt
+import numpy as np
+import pandas as pd
+from boario import logger
+import math
+import inspect
+from functools import partial
+
+from boario.utils.recovery_functions import (
+    concave_recovery,
+    convexe_recovery,
+    linear_recovery,
+    convexe_recovery_scaled,
+)
+
+__all__ = [
+    "Event",
+    "EventKapitalDestroyed",
+    "EventArbitraryProd",
+    "EventKapitalRecover",
+    "EventKapitalRebuild",
+    "Impact",
+    "IndustriesList",
+    "SectorsList",
+    "RegionsList",
+]
+
+Impact = Union[int, float, list, dict, npt.NDArray, pd.DataFrame, pd.Series]
+IndustriesList = Union[List[Tuple[str, str]], pd.MultiIndex, npt.NDArray]
+SectorsList = Union[List[str], pd.Index, npt.NDArray]
+RegionsList = Union[List[str], pd.Index, npt.NDArray]
+FinalCatList = Union[List[str], pd.Index, npt.NDArray]
+
+rebuilding_finaldemand_cat_regex = r".*[hH]ousehold.*|HFCE"
+
+
+
[docs]class Event(metaclass=abc.ABCMeta): + # Class Attributes + __required_class_attributes = [ + "possible_sectors", + "possible_regions", + "temporal_unit_range", + "z_shape", + "y_shape", + "x_shape", + "regions_idx", + "sectors_idx", + "model_monetary_factor", + "sectors_gva_shares", + "Z_distrib", + "mrio_name", + ] + possible_sectors: pd.Index = pd.Index([]) + r"""List of sectors present in the MRIO used by the model""" + + possible_regions: pd.Index = pd.Index([]) + r"""List of regions present in the MRIO used by the model""" + + possible_final_demand_cat: pd.Index = pd.Index([]) + r"""List of final demand categories present in the MRIO used by the model""" + + temporal_unit_range: int = 0 + r"""Maximum temporal unit simulated""" + + z_shape: tuple[int, int] = (0, 0) + r"""Shape of the Z (intermediate consumption) matrix in the model""" + + y_shape: tuple[int, int] = (0, 0) + r"""Shape of the Y (final demand) matrix in the model""" + + x_shape: tuple[int, int] = (0, 0) + r"""Shape of the x (production) vector in the model""" + + regions_idx: npt.NDArray = np.array([]) + r"""lexicographic region indexes""" + + sectors_idx: npt.NDArray = np.array([]) + r"""lexicographic sector indexes""" + + model_monetary_factor: int = 0 + r"""Amount of unitary currency used in the MRIO (e.g. 1000000 if in € millions)""" + + sectors_gva_shares: npt.NDArray = np.array([]) + r"""Fraction of total (regional) GVA for each sectors""" + + Z_distrib: npt.NDArray = np.array([]) + r"""Normalized intermediate consumption matrix""" + + Y_distrib: npt.NDArray = np.array([]) + r"""Normalized final consumption matrix""" + + mrio_name: str = "" + r"""MRIO identification""" + + def __init__( + self, + *, + impact: Impact, + aff_regions: RegionsList = [], + aff_sectors: SectorsList = [], + aff_industries: IndustriesList = [], + impact_industries_distrib: Optional[npt.ArrayLike] = None, + impact_regional_distrib: Optional[npt.ArrayLike] = None, + impact_sectoral_distrib_type: str = "custom", + impact_sectoral_distrib: Optional[npt.ArrayLike] = None, + name: str = "Unnamed", + occurrence: int = 1, + duration: int = 1, + event_monetary_factor: Optional[int] = None, + ) -> None: + r"""Create an event shocking the model from a dictionary. + + An Event object stores all information about a unique shock during simulation such as + time of occurrence, duration, type of shock, amount of damages. Computation + of recovery or initially requested rebuilding demand is also done in this + class. + + Parameters + ---------- + + event : dict + A dictionary holding the necessary information to define an event. + + Examples + -------- + FIXME: Add docs. + + """ + + self._aff_sectors_idx = np.array([]) + self._aff_sectors = pd.Index([]) + self._aff_regions_idx = np.array([]) + self._aff_regions = pd.Index([]) + impact_regional_distrib = ( + np.array(impact_regional_distrib) + if impact_regional_distrib + else np.array([]) + ) + impact_sectoral_distrib = ( + np.array(impact_sectoral_distrib) + if impact_sectoral_distrib + else np.array([]) + ) + impact_industries_distrib = ( + np.array(impact_industries_distrib) + if impact_sectoral_distrib + else np.array([]) + ) + logger.info("Initializing new Event") + logger.debug("Checking required Class attributes are defined") + + if event_monetary_factor is None: + logger.info( + f"No event monetary factor given. Assuming it is the same as the model ({self.model_monetary_factor})" + ) + self.event_monetary_factor = self.model_monetary_factor + else: + self.event_monetary_factor = event_monetary_factor + if self.event_monetary_factor != self.model_monetary_factor: + logger.warning( + f"Event monetary factor ({self.event_monetary_factor}) differs from model monetary factor ({self.model_monetary_factor}). Be careful to define your impact with the correct unit (ie in event monetary factor)." + ) + + for v in Event.__required_class_attributes: + if Event.__dict__[v] is None: + raise AttributeError( + "Required Event Class attribute {} is not set yet so instantiating an Event isn't possible".format( + v + ) + ) + + self.name: str = name + r"""An identifying name for the event (for convenience mostly)""" + + self.occurrence = occurrence + self.duration = duration + self.impact_df: pd.Series = pd.Series( + 0, + dtype="float64", + index=pd.MultiIndex.from_product( + [self.possible_regions, self.possible_sectors], + names=["region", "sector"], + ), + ) + r"""A pandas Series with all possible industries as index, holding the impact vector of the event. The impact is defined for each sectors in each region.""" + + ################## DATAFRAME INIT ################# + # CASE VECTOR 1 (everything is there and regrouped) (only df creation) + if isinstance(impact, pd.Series): + logger.debug("Given impact is a pandas Series") + self.impact_df.loc[impact.index] = impact + if self.name == "Unnamed" and not impact.name is None: + self.name = str(impact.name) + elif isinstance(impact, dict): + logger.debug("Given impact is a dict, converting it to pandas Series") + impact = pd.Series(impact) + self.impact_df.loc[impact.index] = impact + elif isinstance(impact, pd.DataFrame): + logger.debug("Given impact is a pandas DataFrame, squeezing it to a Series") + impact = impact.squeeze() + if not isinstance(impact, pd.Series): + raise ValueError( + "The given impact DataFrame is not a Series after being squeezed" + ) + self.impact_df.loc[impact.index] = impact + # CASE VECTOR 2 (everything is there but not regrouped) AND CASE SCALAR (Only df creation) + elif ( + isinstance(impact, (int, float, list, np.ndarray)) + and len(aff_industries) > 0 + ): + logger.debug( + f"Given impact is a {type(impact)} and list of impacted industries given. Proceeding." + ) + self.impact_df.loc[aff_industries] = impact + elif ( + isinstance(impact, (int, float, list, np.ndarray)) + and len(aff_regions) > 0 + and len(aff_sectors) > 0 + ): + logger.debug( + f"Given impact is a {type(impact)} and lists of impacted regions and sectors given. Proceeding." + ) + if isinstance(aff_regions, str): + aff_regions = [aff_regions] + if isinstance(aff_sectors, str): + aff_sectors = [aff_sectors] + + self.impact_df.loc[ + pd.MultiIndex.from_product([aff_regions, aff_sectors]) + ] = impact + else: + raise ValueError("Invalid input format. Could not initiate pandas Series.") + + # Check for <0 values and remove 0. + if (self.impact_df < 0).any(): + logger.warning( + "Found negative values in impact vector. This should raise concern" + ) + + # SORT DF + # at this point impact_df is built, and can be sorted. Note that if impact was a scalar, impact_df contains copies of this scalar. + logger.debug("Sorting impact Series") + self.impact_df = self.impact_df.sort_index() + + # Init self.impact_sectoral_distrib_type, + self.impact_sectoral_distrib_type = impact_sectoral_distrib_type + ################################################# + + # SET INDEXES ATTR + # note that the following also sets aff_regions and aff_sectors + assert isinstance(self.impact_df.index, pd.MultiIndex) + # Only look for industries where impact is greater than 0 + self.aff_industries = self.impact_df.loc[self.impact_df > 0].index + + logger.debug( + f"impact df at the moment:\n {self.impact_df.loc[self.aff_industries]}" + ) + + ###### SCALAR DISTRIBUTION ###################### + # if impact_industries_distrib is given, set it. We assume impact is scalar ! + # CASE SCALAR + INDUS DISTRIB + if len(impact_industries_distrib) > 0 and not isinstance( + impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) + ): + logger.debug("impact is Scalar and impact_industries_distrib was given") + self.impact_industries_distrib = np.array(impact_industries_distrib) + self.impact_df.loc[self.aff_industries] = ( + self.impact_df.loc[self.aff_industries] * self.impact_industries_distrib + ) + # if impact_reg_dis and sec_dis are give, deduce the rest. We also assume impact is scalar ! + # CASE SCALAR + REGION and SECTOR DISTRIB + elif ( + len(impact_regional_distrib) > 0 + and len(impact_sectoral_distrib) > 0 + and not isinstance( + impact, + (pd.Series, dict, pd.DataFrame, list, np.ndarray), + ) + ): + logger.debug( + "impact is Scalar and impact_regional_distrib and impact_sectoral_distrib were given" + ) + if len(impact_regional_distrib) != len(self.aff_regions) or len( + impact_sectoral_distrib + ) != len(self.aff_sectors): + raise ValueError( + "Lengths of `impact_regional_distrib` and/or `impact_sectoral_distrib` are incompatible with `aff_regions` and/or `aff_sectors`." + ) + else: + self.impact_regional_distrib = impact_regional_distrib + self.impact_sectoral_distrib = impact_sectoral_distrib + self.impact_industries_distrib = ( + self.impact_regional_distrib[:, np.newaxis] + * self.impact_sectoral_distrib + ).flatten() + self.impact_df.loc[self.aff_industries] = ( + self.impact_df.loc[self.aff_industries] + * self.impact_industries_distrib + ) + # CASE SCALAR + 'gdp' distrib + elif ( + len(impact_regional_distrib) > 0 + and len(impact_sectoral_distrib_type) > 0 + and impact_sectoral_distrib_type == "gdp" + and not isinstance( + impact, + (pd.Series, dict, pd.DataFrame, list, np.ndarray), + ) + ): + logger.debug("impact is Scalar and impact_sectoral_distrib_type is 'gdp'") + + self.impact_regional_distrib = impact_regional_distrib + + shares = self.sectors_gva_shares.reshape( + (len(self.possible_regions), len(self.possible_sectors)) + ) + self.impact_sectoral_distrib = ( + shares[self._aff_regions_idx][:, self._aff_sectors_idx] + / shares[self._aff_regions_idx][:, self._aff_sectors_idx].sum(axis=1)[ + :, np.newaxis + ] + ) + self.impact_industries_distrib = ( + self.impact_regional_distrib[:, np.newaxis] + * self.impact_sectoral_distrib + ).flatten() + self.impact_df.loc[self.aff_industries] = ( + self.impact_df.loc[self.aff_industries] * self.impact_industries_distrib + ) + self.impact_sectoral_distrib_type = "gdp" + # CASE SCALAR + NO DISTRIB + list of industries + # if neither was given, we use default values. Again impact should be scalar here ! + elif isinstance(aff_industries, (list, np.ndarray)) and not isinstance( + impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray) + ): + logger.debug( + "impact is Scalar and no distribution was given but a list of affected industries was given" + ) + self._default_distribute_impact_from_industries_list() + self.impact_sectoral_distrib_type = "default (shared equally between affected regions and then affected sectors)" + # CASE SCALAR + NO DISTRIB + list of region + list of sectors + elif ( + len(aff_regions) > 0 + and len(aff_sectors) > 0 + and not isinstance( + impact, + (pd.Series, dict, pd.DataFrame, list, np.ndarray), + ) + ): + logger.debug( + "impact is Scalar and no distribution was given but lists of regions and sectors affected were given" + ) + self._default_distribute_impact_from_industries_list() + self.impact_sectoral_distrib_type = "default (shared equally between affected regions and then affected sectors)" + elif not isinstance(impact, (pd.Series, dict, pd.DataFrame, list, np.ndarray)): + raise ValueError(f"Invalid input format: Could not compute impact") + + self._finish_init() + ################################################## + + self.happened: bool = False + r"""States if the event happened""" + + self.over: bool = False + r"""States if the event is over""" + + self.event_dict: dict = { + "name": str(self.name), + "occurrence": self.occurrence, + "duration": self.duration, + "aff_regions": list(self.aff_regions), + "aff_sectors": list(self.aff_sectors), + "impact": self.total_impact, + "impact_industries_distrib": list(self.impact_industries_distrib), + "impact_regional_distrib": list(self.impact_regional_distrib), + "impact_sectoral_distrib_type": self.impact_sectoral_distrib_type, + "globals_vars": { + "possible_sectors": list(self.possible_sectors), + "possible_regions": list(self.possible_regions), + "temporal_unit_range": self.temporal_unit_range, + "z_shape": self.z_shape, + "y_shape": self.y_shape, + "x_shape": self.x_shape, + "model_monetary_factor": self.model_monetary_factor, + "event_monetary_factor": self.event_monetary_factor, + "mrio_used": self.mrio_name, + }, + } + r"""Store relevant information about the event""" + + def _default_distribute_impact_from_industries_list(self): + # at this point, impact should still be scalar. + logger.debug("Using default impact distribution to industries") + logger.debug( + f"impact df at the moment:\n {self.impact_df.loc[self.aff_industries]}" + ) + self.impact_regional_distrib = np.full( + len(self.aff_regions), 1 / len(self.aff_regions) + ) + + logger.debug( + f"self.impact_regional_distrib: {list(self.impact_regional_distrib)}" + ) + logger.debug(f"len aff_regions: {len(self.aff_regions)}") + self.impact_df.loc[self.aff_industries] = ( + self.impact_df.loc[self.aff_industries] * 1 / len(self.aff_regions) + ) + impact_sec_vec = np.array( + [ + 1 / len(self.aff_industries.to_series().loc[reg]) + for reg in self.aff_regions + ] + ) + self.impact_df.loc[self.aff_industries] = ( + self.impact_df.loc[self.aff_industries] * impact_sec_vec + ) + logger.debug( + f"impact df after default distrib:\n {self.impact_df.loc[self.aff_industries]}" + ) + + def _finish_init(self): + logger.debug("Finishing Event init") + self.impact_vector = self.impact_df.to_numpy() + self.total_impact = self.impact_vector.sum() + self.impact_industries_distrib = ( + self.impact_vector[self.impact_vector > 0] / self.total_impact + ) + self.impact_regional_distrib = ( + self.impact_df.loc[self.aff_industries].groupby("region").sum().values + / self.total_impact + ) + + @property + def aff_industries(self) -> pd.MultiIndex: + r"""The industries affected by the event. + + Parameters + ---------- + + index : pd.MultiIndex + The affected industries as a pandas MultiIndex + + Returns + ------- + A pandas MultiIndex with the regions affected as first level, and sectors affected as second level + + """ + + return self._aff_industries + + @aff_industries.setter + def aff_industries(self, index: pd.MultiIndex): + if not isinstance(index, pd.MultiIndex): + raise ValueError( + "Given impact vector does not have a (region,sector) MultiIndex" + ) + if index.names != ["region", "sector"]: + raise ValueError("MultiIndex must have levels named 'region' and 'sector'") + + regions = index.get_level_values("region").unique() + sectors = index.get_level_values("sector").unique() + if not set(regions).issubset(self.possible_regions): + raise ValueError( + f"Regions {set(regions) - set(self.possible_regions)} are not in the possible regions list" + ) + if not set(sectors).issubset(self.possible_sectors): + raise ValueError( + f"Sectors {set(sectors) - set(self.possible_sectors)} are not in the possible sectors list" + ) + + self.aff_regions = regions + self.aff_sectors = sectors + logger.debug( + f"Setting _aff_industries. There are {len(index)} affected industries" + ) + self._aff_industries = index + self._aff_industries_idx = np.array( + [ + len(self.possible_sectors) * ri + si + for ri in self._aff_regions_idx + for si in self._aff_sectors_idx + ] + ) + + @property + def occurrence(self) -> int: + r"""The temporal unit of occurrence of the event.""" + + return self._occur + + @occurrence.setter + def occurrence(self, value: int): + if not 0 <= value <= self.temporal_unit_range: + raise ValueError( + "Occurrence of event is not in the range of simulation steps : {} not in [0-{}]".format( + value, self.temporal_unit_range + ) + ) + else: + logger.debug(f"Setting occurence to {value}") + self._occur = value + + @property + def duration(self) -> int: + r"""The duration of the event.""" + + return self._duration + + @duration.setter + def duration(self, value: int): + if not 0 <= self.occurrence + value <= self.temporal_unit_range: + raise ValueError( + "Occurrence + duration of event is not in the range of simulation steps : {} + {} not in [0-{}]".format( + self.occurrence, value, self.temporal_unit_range + ) + ) + else: + logger.debug(f"Setting duration to {value}") + self._duration = value + + @property + def aff_regions(self) -> pd.Index: + r"""The array of regions affected by the event""" + + return self._aff_regions + + @property + def aff_regions_idx(self) -> npt.NDArray: + r"""The array of lexicographically ordered affected region indexes""" + + return self._aff_regions_idx + + @aff_regions.setter + def aff_regions(self, value: npt.ArrayLike | str): + if isinstance(value, str): + value = [value] + value = pd.Index(value) + impossible_regions = np.setdiff1d(value, self.possible_regions) + if impossible_regions.size > 0: + raise ValueError( + "These regions are not in the model : {}".format(impossible_regions) + ) + else: + self._aff_regions = pd.Index(value, name="region") + logger.debug(f"Setting _aff_regions to {self._aff_regions}") + self._aff_regions_idx = np.searchsorted(self.possible_regions, value) + + @property + def aff_sectors(self) -> pd.Index: + r"""The array of affected sectors by the event""" + + return self._aff_sectors + + @property + def aff_sectors_idx(self) -> npt.NDArray: + r"""The array of lexicographically ordered affected sectors indexes""" + + return self._aff_sectors_idx + + @aff_sectors.setter + def aff_sectors(self, value: npt.ArrayLike | str): + if isinstance(value, str): + value = [value] + value = pd.Index(value, name="sector") + impossible_sectors = np.setdiff1d(value, self.possible_sectors) + if impossible_sectors.size > 0: + raise ValueError( + "These sectors are not in the model : {}".format(impossible_sectors) + ) + else: + self._aff_sectors = pd.Index(value, name="sector") + logger.debug(f"Setting _aff_sectors to {self._aff_sectors}") + self._aff_sectors_idx = np.searchsorted(self.possible_sectors, value) + + @property + def impact_regional_distrib(self) -> npt.NDArray: + r"""The array specifying how damages are distributed among affected regions""" + + return self._impact_regional_distrib + + @impact_regional_distrib.setter + def impact_regional_distrib(self, value: npt.ArrayLike): + if self.aff_regions is None: + raise AttributeError("Affected regions attribute isn't set yet") + value = np.array(value) + if value.size != self.aff_regions.size: + raise ValueError( + "There are {} affected regions by the event and length of given damage distribution is {}".format( + self.aff_regions.size, value.size + ) + ) + s = value.sum() + if not math.isclose(s, 1): + raise ValueError( + "Damage distribution doesn't sum up to 1 but to {}, which is not valid".format( + s + ) + ) + self._impact_regional_distrib = value + + @property + def impact_sectoral_distrib_type(self) -> str: + r"""The type of damages distribution among sectors (currently only 'gdp')""" + + return self._impact_sectoral_distrib_type + + @impact_sectoral_distrib_type.setter + def impact_sectoral_distrib_type(self, value: str): + logger.debug(f"Setting _impact_sectoral_distrib_type to {value}") + self._impact_sectoral_distrib_type = value + + def __repr__(self): + # TODO: find ways to represent long lists + return f""" + Event( + name = {self.name}, + occur = {self.occurrence}, + duration = {self.duration} + aff_regions = {self.aff_regions.to_list()}, + aff_sectors = {self.aff_sectors.to_list()}, + ) + """
+ + +
[docs]class EventArbitraryProd(Event): + def __init__( + self, + *, + impact: Impact, + recovery_time: int = 1, + recovery_function: str = "linear", + aff_regions: RegionsList = [], + aff_sectors: SectorsList = [], + aff_industries: IndustriesList = [], + impact_industries_distrib: Optional[npt.ArrayLike] = None, + impact_regional_distrib: Optional[npt.ArrayLike] = None, + impact_sectoral_distrib_type="equally shared", + impact_sectoral_distrib: Optional[npt.ArrayLike] = None, + name: str = "Unnamed", + occurrence: int = 1, + duration: int = 1, + ) -> None: + if not isinstance(impact, pd.Series): + raise NotImplementedError( + "Arbitrary production capacity shock currently require to be setup from a vector indexed by the industries (regions,sectors) affected, and where values are the share of production capacity lost." + ) + + if (impact > 1.0).any(): + raise ValueError( + "Impact is greater than 100% (1.) for at least an industry." + ) + + super().__init__( + impact=impact, + aff_regions=aff_regions, + aff_sectors=aff_sectors, + aff_industries=aff_industries, + impact_industries_distrib=impact_industries_distrib, + impact_regional_distrib=impact_regional_distrib, + impact_sectoral_distrib_type=impact_sectoral_distrib_type, + impact_sectoral_distrib=impact_sectoral_distrib, + name=name, + occurrence=occurrence, + duration=duration, + event_monetary_factor=None, + ) + + self._prod_cap_delta_arbitrary_0 = ( + self.impact_vector.copy() + ) # np.zeros(shape=len(self.possible_sectors)) + self.prod_cap_delta_arbitrary = ( + self.impact_vector.copy() + ) # np.zeros(shape=len(self.possible_sectors)) + self.recovery_time = recovery_time + r"""The characteristic recovery duration after the event is over""" + + self.recovery_function = recovery_function + + logger.info("Initialized") + + @property + def prod_cap_delta_arbitrary(self) -> npt.NDArray: + r"""The optional array storing arbitrary (as in not related to productive_capital destroyed) production capacity loss""" + return self._prod_cap_delta_arbitrary + + @prod_cap_delta_arbitrary.setter + def prod_cap_delta_arbitrary(self, value: npt.NDArray): + self._prod_cap_delta_arbitrary = value + + @property + def recoverable(self) -> bool: + r"""A boolean stating if an event can start recover""" + return self._recoverable + + @recoverable.setter + def recoverable(self, current_temporal_unit: int): + reb = (self.occurrence + self.duration) <= current_temporal_unit + if reb and not self.recoverable: + logger.info( + "Temporal_Unit : {} ~ Event named {} that occured at {} in {} has started recovering (arbitrary production capacity loss)".format( + current_temporal_unit, + self.name, + self.occurrence, + self.aff_regions.to_list(), + ) + ) + self._recoverable = reb + + @property + def recovery_function(self) -> Callable: + r"""The recovery function used for recovery (`Callable`)""" + return self._recovery_fun + + @recovery_function.setter + def recovery_function(self, r_fun: str | Callable | None): + if r_fun is None: + r_fun = "instant" + if self.recovery_time is None: + raise AttributeError( + "Impossible to set recovery function if no recovery time is given." + ) + if isinstance(r_fun, str): + if r_fun == "linear": + fun = linear_recovery + elif r_fun == "convexe": + fun = convexe_recovery_scaled + elif r_fun == "convexe noscale": + fun = convexe_recovery + elif r_fun == "concave": + fun = concave_recovery + else: + raise NotImplementedError( + "No implemented recovery function corresponding to {}".format(r_fun) + ) + elif callable(r_fun): + r_fun_argsspec = inspect.getfullargspec(r_fun) + r_fun_args = r_fun_argsspec.args + r_fun_argsspec.kwonlyargs + if not all( + args in r_fun_args + for args in [ + "init_impact_stock", + "elapsed_temporal_unit", + "recovery_time", + ] + ): + raise ValueError( + "Recovery function has to have at least the following keyword arguments: {}".format( + [ + "init_impact_stock", + "elapsed_temporal_unit", + "recovery_time", + ] + ) + ) + fun = r_fun + + else: + raise ValueError("Given recovery function is not a str or callable") + + r_fun_partial = partial( + fun, + init_impact_stock=self._prod_cap_delta_arbitrary_0, + recovery_time=self.recovery_time, + ) + self._recovery_fun = r_fun_partial + + def recovery(self, current_temporal_unit: int): + elapsed = current_temporal_unit - (self.occurrence + self.duration) + if elapsed < 0: + raise RuntimeError("Trying to recover before event is over") + if self.recovery_function is None: + raise RuntimeError( + "Trying to recover event while recovery function isn't set yet" + ) + res = self.recovery_function(elapsed_temporal_unit=elapsed) + precision = int(math.log10(self.model_monetary_factor)) + 1 + res = np.around(res, precision) + if not np.any(res): + self.over = True + self._prod_cap_delta_arbitrary = res
+ + +
[docs]class EventKapitalDestroyed(Event, metaclass=abc.ABCMeta): + def __init__( + self, + *, + impact: Impact, + households_impact: Optional[Impact] = None, + aff_regions: RegionsList = [], + aff_sectors: SectorsList = [], + aff_industries: IndustriesList = [], + impact_industries_distrib=None, + impact_regional_distrib=None, + impact_sectoral_distrib_type="equally shared", + impact_sectoral_distrib=None, + name="Unnamed", + occurrence=1, + duration=1, + event_monetary_factor=None, + ) -> None: + super().__init__( + impact=impact, + aff_regions=aff_regions, + aff_sectors=aff_sectors, + aff_industries=aff_industries, + impact_industries_distrib=impact_industries_distrib, + impact_regional_distrib=impact_regional_distrib, + impact_sectoral_distrib_type=impact_sectoral_distrib_type, + impact_sectoral_distrib=impact_sectoral_distrib, + name=name, + occurrence=occurrence, + duration=duration, + event_monetary_factor=event_monetary_factor, + ) + # The only thing we have to do is affecting/computing the regional_sectoral_productive_capital_destroyed + self.total_productive_capital_destroyed = self.total_impact + self.total_productive_capital_destroyed *= ( + self.event_monetary_factor / self.model_monetary_factor + ) + logger.info( + f"Total impact on productive capital is {self.total_productive_capital_destroyed} (in model unit)" + ) + self.remaining_productive_capital_destroyed = ( + self.total_productive_capital_destroyed + ) + self._regional_sectoral_productive_capital_destroyed_0 = self.impact_vector * ( + self.event_monetary_factor / self.model_monetary_factor + ) + self.regional_sectoral_productive_capital_destroyed = self.impact_vector * ( + self.event_monetary_factor / self.model_monetary_factor + ) + self.households_impact_df: pd.Series = pd.Series( + 0, + dtype="float64", + index=pd.MultiIndex.from_product( + [self.possible_regions, self.possible_final_demand_cat], + names=["region", "final_demand_cat"], + ), + ) + r"""A pandas Series with all possible (regions, final_demand_cat) as index, holding the households impacts vector of the event. The impact is defined for each region and each final demand category.""" + + if households_impact: + try: + rebuilding_demand_idx = self.possible_final_demand_cat[ + self.possible_final_demand_cat.str.match( + rebuilding_finaldemand_cat_regex + ) + ] # .values[0] + + except IndexError: + rebuilding_demand_idx = self.possible_final_demand_cat[0] + if isinstance(households_impact, pd.Series): + logger.debug("Given household impact is a pandas Series") + self.households_impact_df.loc[ + households_impact.index, rebuilding_demand_idx + ] = households_impact + elif isinstance(households_impact, dict): + logger.debug("Given household impact is a dict") + households_impact = pd.Series(households_impact) + self.households_impact_df.loc[ + households_impact.index, rebuilding_demand_idx + ] = households_impact + elif isinstance(households_impact, pd.DataFrame): + logger.debug("Given household impact is a dataframe") + households_impact = households_impact.squeeze() + if not isinstance(households_impact, pd.Series): + raise ValueError( + "The given households_impact DataFrame is not a Series after being squeezed" + ) + self.households_impact_df.loc[ + households_impact.index, rebuilding_demand_idx + ] = households_impact + elif isinstance(households_impact, (list, np.ndarray)): + if len(households_impact) != len(self.aff_regions): + raise ValueError( + f"Length mismatch between households_impact ({len(households_impact)} and affected regions ({len(self.aff_regions)}))" + ) + else: + self.households_impact_df.loc[ + self.aff_regions, rebuilding_demand_idx + ] = households_impact + elif isinstance(households_impact, (float, int)): + if self.impact_regional_distrib is not None: + logger.warning( + f"households impact given as scalar, distributing among region following `impact_regional_distrib` ({self.impact_regional_distrib}) " + ) + self.households_impact_df.loc[ + self.aff_regions, rebuilding_demand_idx + ] = (households_impact * self.impact_regional_distrib) + self.households_impact_df *= ( + self.event_monetary_factor / self.model_monetary_factor + ) + logger.info( + f"Total impact on households is {self.households_impact_df.sum()} (in model unit)" + ) + + @property + def regional_sectoral_productive_capital_destroyed(self) -> npt.NDArray: + r"""The (optional) array of productive_capital destroyed per industry (ie region x sector)""" + return self._regional_sectoral_productive_capital_destroyed + + @regional_sectoral_productive_capital_destroyed.setter + def regional_sectoral_productive_capital_destroyed(self, value: npt.ArrayLike): + value = np.array(value) + logger.debug( + f"Changing regional_sectoral_productive_capital_destroyed with {value}\n Sum is {value.sum()}" + ) + if value.shape != self.x_shape: + raise ValueError( + "Incorrect shape give for regional_sectoral_productive_capital_destroyed: {} given and {} expected".format( + value.shape, self.x_shape + ) + ) + self._regional_sectoral_productive_capital_destroyed = value
+ + +
[docs]class EventKapitalRebuild(EventKapitalDestroyed): + def __init__( + self, + *, + impact: Impact, + rebuilding_sectors: Union[dict[str, float], pd.Series], + rebuild_tau=60, + households_impact: Optional[Impact] = None, + aff_regions: RegionsList = [], + aff_sectors: SectorsList = [], + aff_industries: IndustriesList = [], + impact_industries_distrib=None, + impact_regional_distrib=None, + impact_sectoral_distrib_type="equally shared", + impact_sectoral_distrib=None, + name="Unnamed", + occurrence=1, + duration=1, + rebuilding_factor: float = 1.0, + event_monetary_factor=None, + ) -> None: + super().__init__( + impact=impact, + households_impact=households_impact, + aff_regions=aff_regions, + aff_sectors=aff_sectors, + aff_industries=aff_industries, + impact_industries_distrib=impact_industries_distrib, + impact_regional_distrib=impact_regional_distrib, + impact_sectoral_distrib_type=impact_sectoral_distrib_type, + impact_sectoral_distrib=impact_sectoral_distrib, + name=name, + occurrence=occurrence, + duration=duration, + event_monetary_factor=event_monetary_factor, + ) + + self._rebuildable = False + self.rebuild_tau = rebuild_tau + self.rebuilding_sectors = rebuilding_sectors + self.rebuilding_factor = rebuilding_factor + + industrial_rebuilding_demand = np.zeros(shape=self.z_shape) + tmp = np.zeros(self.z_shape, dtype="float") + mask = np.ix_( + np.union1d( + self._rebuilding_industries_RoW_idx, self._rebuilding_industries_idx + ), + self._aff_industries_idx, + ) + industrial_rebuilding_demand = np.outer( + self.rebuilding_sectors_shares, + self.regional_sectoral_productive_capital_destroyed, + ) + industrial_rebuilding_demand = industrial_rebuilding_demand * rebuilding_factor + tmp[mask] = self.Z_distrib[mask] * industrial_rebuilding_demand[mask] + + precision = int(math.log10(self.model_monetary_factor)) + 1 + np.around(tmp, precision, tmp) + + self.rebuilding_demand_indus = tmp + + households_rebuilding_demand = np.zeros(shape=self.y_shape) + tmp = np.zeros(self.y_shape, dtype="float") + mask = np.ix_( # type: ignore + np.union1d( # type: ignore + self._rebuilding_industries_RoW_idx, self._rebuilding_industries_idx + ), + list(range(self.y_shape[1])), + ) + + households_rebuilding_demand = np.outer( + self.rebuilding_sectors_shares, self.households_impact_df.to_numpy() + ) + households_rebuilding_demand = households_rebuilding_demand * rebuilding_factor + tmp[mask] = self.Y_distrib[mask] * households_rebuilding_demand[mask] + np.around(tmp, precision, tmp) + self.rebuilding_demand_house = tmp + logger.debug( + f"house rebuilding demand is: {pd.DataFrame(self.rebuilding_demand_house, index=pd.MultiIndex.from_product([self.possible_regions,self.possible_sectors]))}" + ) + + self.event_dict["rebuilding_sectors"] = { + sec: share + for sec, share in zip( + self.rebuilding_sectors, self.rebuilding_sectors_shares + ) + } + + @property + def rebuilding_sectors(self) -> pd.Index: + r"""The (optional) array of rebuilding sectors""" + return self._rebuilding_sectors + + @property + def rebuilding_sectors_idx(self) -> npt.NDArray: + r"""The (optional) array of indexes of the rebuilding sectors (in lexicographic order)""" + return self._rebuilding_sectors_idx + + @property + def rebuilding_sectors_shares(self) -> npt.NDArray: + r"""The array specifying how rebuilding demand is distributed along the rebuilding sectors""" + return self._rebuilding_sectors_shares + + @rebuilding_sectors.setter + def rebuilding_sectors(self, value: dict[str, float] | pd.Series): + if isinstance(value, dict): + reb_sectors = pd.Series(value) + else: + reb_sectors = value + assert np.isclose(reb_sectors.sum(), 1.0) + impossible_sectors = np.setdiff1d(reb_sectors.index, self.possible_sectors) + if impossible_sectors.size > 0: + raise ValueError( + "These sectors are not in the model : {}".format(impossible_sectors) + ) + else: + self._rebuilding_sectors = reb_sectors.index + self._rebuilding_sectors_idx = np.searchsorted( + self.possible_sectors, reb_sectors.index + ) + + self._rebuilding_sectors_shares = np.zeros(self.x_shape) + self._rebuilding_industries_idx = np.array( + [ + len(self.possible_sectors) * ri + si + for ri in self._aff_regions_idx + for si in self._rebuilding_sectors_idx + ] + ) + self._rebuilding_industries_RoW_idx = np.array( + [ + len(self.possible_sectors) * ri + si + for ri in range(len(self.possible_regions)) + if ri not in self._aff_regions_idx + for si in self._rebuilding_sectors_idx + ] + ) + self._rebuilding_sectors_shares[self._rebuilding_industries_idx] = np.tile( + np.array(reb_sectors.values), len(self.aff_regions) + ) + self._rebuilding_sectors_shares[ + self._rebuilding_industries_RoW_idx + ] = np.tile( + np.array(reb_sectors.values), + (len(self.possible_regions) - len(self.aff_regions)), + ) + + @property + def rebuilding_demand_house(self) -> npt.NDArray: + r"""The optional array of rebuilding demand from households""" + return self._rebuilding_demand_house + + @rebuilding_demand_house.setter + def rebuilding_demand_house(self, value: npt.ArrayLike): + value = np.array(value) + if value.shape != self.y_shape: + raise ValueError( + "Incorrect shape give for rebuilding_demand_house: {} given and {} expected".format( + value.shape, self.y_shape + ) + ) + value[value < 10 / self.model_monetary_factor] = 0.0 + self._rebuilding_demand_house = value + + @property + def rebuilding_demand_indus(self) -> npt.NDArray: + r"""The optional array of rebuilding demand from industries (to rebuild their productive_capital)""" + return self._rebuilding_demand_indus + + @rebuilding_demand_indus.setter + def rebuilding_demand_indus(self, value: npt.ArrayLike): + value = np.array(value) + if value.shape != self.z_shape: + raise ValueError( + "Incorrect shape give for rebuilding_demand_indus: {} given and {} expected".format( + value.shape, self.z_shape + ) + ) + value[value < 10 / self.model_monetary_factor] = 0.0 + self._rebuilding_demand_indus = value + # Also update productive_capital destroyed + self.regional_sectoral_productive_capital_destroyed = value.sum(axis=0) * ( + 1 / self.rebuilding_factor + ) + + @property + def rebuildable(self) -> bool: + r"""A boolean stating if an event can start rebuild""" + return self._rebuildable + + @rebuildable.setter + def rebuildable(self, current_temporal_unit: int): + reb = (self.occurrence + self.duration) <= current_temporal_unit + if reb and not self.rebuildable: + logger.info( + "Temporal_Unit : {} ~ Event named {} that occured at {} in {} for {} damages has started rebuilding".format( + current_temporal_unit, + self.name, + self.occurrence, + self.aff_regions.to_list(), + self.total_productive_capital_destroyed, + ) + ) + self._rebuildable = reb
+ + +
[docs]class EventKapitalRecover(EventKapitalDestroyed): + def __init__( + self, + *, + impact: Impact, + recovery_time: int, + recovery_function: str = "linear", + households_impact: Impact | None = None, + aff_regions: RegionsList = [], + aff_sectors: SectorsList = [], + aff_industries: IndustriesList = [], + impact_industries_distrib=None, + impact_regional_distrib=None, + impact_sectoral_distrib_type="equally shared", + impact_sectoral_distrib=None, + name="Unnamed", + occurrence=1, + duration=1, + event_monetary_factor=None, + ) -> None: + super().__init__( + impact=impact, + households_impact=households_impact, + aff_regions=aff_regions, + aff_sectors=aff_sectors, + aff_industries=aff_industries, + impact_industries_distrib=impact_industries_distrib, + impact_regional_distrib=impact_regional_distrib, + impact_sectoral_distrib_type=impact_sectoral_distrib_type, + impact_sectoral_distrib=impact_sectoral_distrib, + name=name, + occurrence=occurrence, + duration=duration, + event_monetary_factor=event_monetary_factor, + ) + self.recovery_time = recovery_time + self.recovery_function = recovery_function + self._recoverable = False + + @property + def recoverable(self) -> bool: + r"""A boolean stating if an event can start recover""" + return self._recoverable + + @recoverable.setter + def recoverable(self, current_temporal_unit: int): + reb = (self.occurrence + self.duration) <= current_temporal_unit + if reb and not self.recoverable: + logger.info( + "Temporal_Unit : {} ~ Event named {} that occured at {} in {} for {} damages has started recovering (no rebuilding demand)".format( + current_temporal_unit, + self.name, + self.occurrence, + self.aff_regions.to_list(), + self.total_productive_capital_destroyed, + ) + ) + self._recoverable = reb + + @property + def recovery_function(self) -> Callable: + r"""The recovery function used for recovery (`Callable`)""" + return self._recovery_fun + + @recovery_function.setter + def recovery_function(self, r_fun: str | Callable | None): + if r_fun is None: + r_fun = "linear" + if self.recovery_time is None: + raise AttributeError( + "Impossible to set recovery function if no recovery time is given." + ) + if isinstance(r_fun, str): + if r_fun == "linear": + fun = linear_recovery + elif r_fun == "convexe": + fun = convexe_recovery_scaled + elif r_fun == "convexe noscale": + fun = convexe_recovery + elif r_fun == "concave": + fun = concave_recovery + else: + raise NotImplementedError( + "No implemented recovery function corresponding to {}".format(r_fun) + ) + elif callable(r_fun): + r_fun_argsspec = inspect.getfullargspec(r_fun) + r_fun_args = r_fun_argsspec.args + r_fun_argsspec.kwonlyargs + if not all( + args in r_fun_args + for args in [ + "init_impact_stock", + "elapsed_temporal_unit", + "recovery_time", + ] + ): + raise ValueError( + "Recovery function has to have at least the following keyword arguments: {}".format( + [ + "init_impact_stock", + "elapsed_temporal_unit", + "recovery_time", + ] + ) + ) + fun = r_fun + + else: + raise ValueError("Given recovery function is not a str or callable") + + r_fun_partial = partial( + fun, + init_impact_stock=self._regional_sectoral_productive_capital_destroyed_0, + recovery_time=self.recovery_time, + ) + self._recovery_fun = r_fun_partial + + def recovery(self, current_temporal_unit: int): + elapsed = current_temporal_unit - (self.occurrence + self.duration) + if elapsed < 0: + raise RuntimeError("Trying to recover before event is over") + if self.recovery_function is None: + raise RuntimeError( + "Trying to recover event while recovery function isn't set yet" + ) + res = self.recovery_function(elapsed_temporal_unit=elapsed) + precision = int(math.log10(self.model_monetary_factor)) + 1 + res = np.around(res, precision) + res[res < 0] = 0.0 + if not np.any(res): + self.over = True + self.regional_sectoral_productive_capital_destroyed = res
+
+ +
+ + + +
+ +
+
+
+ +
+ + + + +
+
+ +
+ +
+
+
+ + + + + +
+ + +
+ + +||||||| 66c6efed +======= @@ -1590,4 +3291,5 @@

Source code for boario.event

 
   
   
-
\ No newline at end of file
+
+>>>>>>> master
diff --git a/docs/build/html/_modules/boario/extended_models.html b/docs/build/html/_modules/boario/extended_models.html
index a1829e3..aa6dd53 100644
--- a/docs/build/html/_modules/boario/extended_models.html
+++ b/docs/build/html/_modules/boario/extended_models.html
@@ -1,5 +1,4 @@
 
-
 
 
 
@@ -18,12 +17,12 @@
   
   
   
-  
-
-
+  
+
+
 
   
-  
+  
   
 
 
@@ -34,14 +33,15 @@
     
   
   
-  
-
+  
+
 
     
     
     
     
     
+    
     
     
     
@@ -157,7 +157,7 @@
 
                     
                 
@@ -182,7 +182,7 @@
       
       
         
       
@@ -269,7 +269,7 @@
 
                     
                 
@@ -283,7 +283,7 @@
       
@@ -622,7 +685,8 @@ 

Defining events from a scalarDifferent types of Event objects
  • Defining Event from an impact vector
  • Defining events from a scalar
  • @@ -657,7 +721,7 @@

    Defining events from a scalar - + @@ -665,8 +729,8 @@

    Defining events from a scalar - + +

    diff --git a/docs/build/html/boario-installation.html b/docs/build/html/boario-installation.html index 9c820f3..de68c46 100644 --- a/docs/build/html/boario-installation.html +++ b/docs/build/html/boario-installation.html @@ -1,5 +1,4 @@ - @@ -7,7 +6,7 @@ - + Installation — BoARIO documentation @@ -19,12 +18,12 @@ - - - + + + - + @@ -35,14 +34,15 @@ - - + + + @@ -160,7 +160,7 @@ @@ -185,7 +185,7 @@ @@ -276,7 +276,7 @@ @@ -290,7 +290,7 @@ @@ -441,7 +457,9 @@

    Requirements @@ -473,7 +491,7 @@

    Requirements - + @@ -481,8 +499,8 @@

    Requirements - + +