###################################################################################
# MIT License
# Copyright (c) 2015-2025 Gregory Paradis
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
###################################################################################
"""
This module implements functions for building and running the wood supply simulation
models.
The :py:class:`ws3.forest.ForestModel` and :py:class:`ws3.forest.DevelopmentType` classes
constitute the core functional units of this module, and of the :py:mod:`ws3` package in general.
"""
from __future__ import annotations
import copy
import itertools
import operator
import re
import sys
from collections import defaultdict as dd # noqa: E402
from collections.abc import Callable
from concurrent.futures import ProcessPoolExecutor, as_completed # noqa: E402
from functools import reduce
from itertools import chain
from multiprocessing import get_context # noqa: E402
from typing import (
Any,
)
import dill # noqa: E402
import pandas as pd # noqa: E402
try: # noqa: E402
from ws3 import common, core, opt
except ImportError: # "__main__" case # noqa: E402
from ws3 import common, core, opt
from ws3.forest_helper import ( # noqa: E402
MP_CONTEXT,
PersistentWorkerPool,
auto_batch,
init_worker_gen_vars,
sanitize_func,
worker_cmp_cflw_batch,
worker_cmp_cflw_phase3,
worker_cmp_cflw_phase3_batch,
worker_cmp_cgen_phase3,
worker_cmp_cgen_phase3_batch,
worker_gen_vars,
worker_summarize_tree_batch,
)
def _search(pattern: str, string: str, what: str, flags: int = 0) -> re.Match[str]:
"""
``re.search`` that fails with a useful message instead of ``None``.
The Woodstock parsers call ``re.search(...).group(...)`` in many places. When
the input is malformed the search returns ``None`` and the chained ``.group()``
raises ``AttributeError: 'NoneType' object has no attribute 'group'`` -- which
says nothing about which file, which construct, or what was expected.
:param pattern: Regular expression to apply.
:param string: Text to search.
:param what: Human-readable description of the construct being parsed, used in
the error message.
:param flags: Optional ``re`` flags.
:return: The match object, guaranteed non-None.
:raises ValueError: If the pattern does not match.
"""
match = re.search(pattern, string, flags)
if match is None:
raise ValueError(
f"Could not parse {what}.\n"
f" expected pattern: {pattern}\n"
f" actual input : {string.strip()!r}"
)
return match
[docs]
class GreedyAreaSelector:
"""
Default AreaSelector implementation. Selects areas for treatment from oldest age classes.
"""
parent: ForestModel
def __init__(self, parent: ForestModel):
self.parent = parent
[docs]
def operate(self, period: int, acode: str, target_area: float, mask: tuple[Any, ...] | None = None,
commit_actions: bool = True, verbose: bool = False) -> float:
"""
Greedily operate on oldest operable age classes.
Returns missing area (i.e., difference between target and operated areas).
:param int period: The time period for the operation.
:param str acode: The action code to specify the action.
:param float target_area: The desired area to be achieved through operation.
:param tuple mask: Tuple of values for development types.
:param bool commit_actions: Flag indicating whether to commit actions. Defaults to True.
:param bool verbose: Verbosity flag. Defaults to False.
"""
def key(item):
return max(item[1])
odt = sorted(self.parent.operable_dtypes(acode, period, mask).items(), key=key)
if verbose:
print(' entering selector.operate()', len(odt), 'operable dtypes')
while target_area > 0 and odt:
while target_area > 0 and odt:
popped = odt.pop()
try:
dtk, ages = popped #odt.pop()
except Exception:
print(odt)
print(popped)
raise
age = sorted(ages)[-1]
oa = self.parent.dtypes[dtk].operable_area(acode, period, age)
if not oa:
continue # nothing to operate
area = min(oa, target_area)
target_area -= area
if area < 0:
print('negative area', area, oa, target_area, acode, period, age)
raise AssertionError()
if verbose:
print(' selector found area', [' '.join(dtk)], acode, period, age, area)
self.parent.apply_action(dtk, acode, period, age, area, compile_c_ycomps=True,
fuzzy_age=False, recourse_enabled=False, verbose=verbose)
odt = sorted(self.parent.operable_dtypes(acode, period, mask).items(), key=key)
self.parent.commit_actions(period, repair_future_actions=True)
if verbose:
print(f'GreedyAreaSelector.operate done (remaining target_area: {target_area:0.1f})')
return target_area
[docs]
class Action:
"""
Encapsulates data for an action.
"""
code: str
targetage: int | None
descr: str
lockexempt: bool
oper_a: Any | None
oper_p: Any | None
components: list[str]
partial: list[str]
is_compiled: bool
is_harvest: int
is_sticky: int
treatment_type: Any | None
def __init__(self,
code: str,
targetage: int | None = None,
descr: str = '',
lockexempt: bool = False,
components: list[str] | None = None,
partial: list[str] | None = None,
is_harvest: int = 0,
is_sticky: int = 0):
self.code = code
self.targetage = targetage
self.descr = descr
self.lockexempt = lockexempt
self.oper_a = None
self.oper_p = None
self.components = components or []
self.partial = partial or []
self.is_compiled = False
self.is_harvest = is_harvest
self.is_sticky = is_sticky
self.treatment_type = None
[docs]
class DevelopmentType:
"""
Encapsulates development type data (curves, age, area), and provides methods
to operate on the data. This is the core class in this module, with respect
to tracking forest inventory and simulating growth and actions.
"""
_bo: dict[str, Callable[[Any, Any], Any]] = {'AND':operator.and_, '&':operator.and_, 'OR':operator.or_, '|':operator.or_}
key: tuple[str, ...]
parent: ForestModel
_rc: Callable[[core.Curve], core.Curve]
_max_age: int
_ycomps: dict[str, core.Curve | None]
_complex_ycomps: dict[str, str]
_zero_curve: core.Curve
_unit_curve: core.Curve
_ages_curve: core.Curve
_resolvers: dict[str, Callable[[str, str], tuple[str, core.Curve]]]
transitions: dict[tuple[str, int], list[Any]]
_areas: dict[int, dd[float]] # type: ignore[type-arg]
oper_expr: dd[list] # type: ignore[type-arg]
operability: dict[str, dict[int, tuple[int, int] | None]]
def __init__(self,
key: tuple[str, ...],
parent: ForestModel):
"""
:param tuple key: Development type key (a unique combination of theme values).
Tuple length must match the number of themes in the parent forest model.
:param parent: Parent forest model.
:type parent: :py:class:`ws3.forest.ForestModel`
"""
self.key = key
self.parent = parent
self._rc = parent.register_curve # shorthand
self._max_age = parent.max_age
self._ycomps = {}
self._complex_ycomps = {}
self._zero_curve = parent.common_curves['zero']
self._unit_curve = parent.common_curves['unit']
self._ages_curve = parent.common_curves['ages']
self._resolvers = {'MULTIPLY':self._resolver_multiply,
'DIVIDE':self._resolver_divide,
'SUM':self._resolver_sum,
'CAI':self._resolver_cai,
'MAI':self._resolver_mai,
'YTP':self._resolver_ytp,
'RANGE':self._resolver_range}
self.transitions = {} # keys are (acode, age) tuples
#######################################################################
# Use period 0 slot to store starting inventory.
self._areas = {p:dd(float) for p in range(0, self.parent.horizon+1)}
#######################################################################
self.oper_expr = dd(list)
self.operability = {}
[docs]
def operable_ages(self, acode: str, period: int) -> list[int] | None:
"""
Finds list of ages at which ``self`` is operable, given an action code and period index.
Takes into account both action operability age range and current inventory at the specified period.
:param str acode: Action code for which to compile operable ages.
:param int period: Period at which to compile operable ages.
:return list: List of ages at which the specified action is operable in the specified period.
"""
if acode not in self.oper_expr: # action not defined for this development type
return None
if acode not in self.operability: # action not compiled yet...
if self.compile_action(acode) == -1:
return None # never operable
if period not in self.operability[acode]:
return None
else:
period_oper = self.operability[acode][period]
if period_oper is None:
return None
lo, hi = period_oper
return list(set(range(lo, hi+1)).intersection(list(self._areas[period].keys())))
[docs]
def is_operable(self, acode: str, period: int, age: int | None = None, verbose: bool = False) -> bool | tuple[int, int]:
"""
Test hypothetical operability, given an action code, a period, and optional age.
Does not imply that there is any operable area in current inventory at the specified period.
:param str acode: The action code for which to test operability.
:param int period: The period in which to test operability.
:param int age: The age at which to test operability. If ``None``, only checks operability for the period.
:param bool verbose: Verbosity flag.
"""
if acode not in self.oper_expr: # action not defined for this development type
if verbose:
print('acode operability undefined', acode, self.oper_expr)
return False
if acode not in self.operability: # action not compiled yet...
if self.compile_action(acode) == -1:
if verbose:
print('never operable', acode)
return False # never operable
if period not in self.operability[acode]:
return False
else:
period_oper = self.operability[acode][period]
if period_oper is None:
return False
lo, hi = period_oper
if age is not None:
return age >= lo and age <= hi
else:
return lo, hi
[docs]
def operable_area(self, acode: str, period: int, age: int | None = None, cleanup: bool = True) -> float:
"""
Compiles operable area, given an action code, a period, and optional age.
:param str acode: The action code to determine operability.
:param int period: The period to determine operability for.
:param int age: The age to determine operability for. If None, only checks operability for the period.
:param bool cleanup: If ``True`` (default), removes the age class from the inventory dict if operable
area is less than ``self.parent.area_epsilon``.
:return float: Operable area. Returns 0 if inoperable or no current inventory, and operable area otherwise.
"""
if acode not in self.oper_expr: # action not defined for this development type
return 0.
if acode not in self.operability: # action not xf yet...
if self.compile_action(acode) == -1:
return 0. # never operable
if age is None: # return total operable area
return sum(self.operable_area(acode, period, a) for a in list(self._areas[period].keys()))
if age not in self._areas[period]: # age class not in inventory
return 0.
elif abs(self._areas[period][age]) < self.parent.area_epsilon: # negligible area
if cleanup: # remove ageclass from dict (frees up memory)
del self._areas[period][age]
return 0.0 # type: ignore[no-any-return]
elif self.is_operable(acode, period, age):
return self._areas[period][age] # type: ignore[no-any-return]
else:
return 0.0 # type: ignore[no-any-return]
raise AssertionError() # type: ignore[unreachable]
[docs]
def area(self, period: int, age: int | None = None, area: float | None = None, delta: bool = True) -> float | None:
"""
If area not specified, returns area inventory for period (with optional age filter), else sets area for period and age.
If delta switch active (default ``True``), area value is interpreted as an increment on current inventory
(otherwise will clobber current inventory).
:param int period: The period for which the area is being retrieved or set.
:param int age: The age for which the area is being retrieved or set. If None, returns total area.
:param float area: The area value to set. If None, returns the area inventory.
:param bool delta: If True (default), interprets the area value as an increment on the current inventory.
If False, sets the area value directly.
"""
if area is None: # return area for period and age
if age is not None:
try:
return float(self._areas[period][age]) # type: ignore[no-any-return]
except Exception as e:
print(e)
return 0.0 # type: ignore[no-any-return]
else: # return total area
return float(sum(self._areas[period][a] for a in self._areas[period])) # type: ignore[no-any-return]
else:
if delta:
self._areas[period][age] += area
else:
self._areas[period][age] = area
return None
return None
[docs]
def resolve_condition(self, yname: str, lo: float, hi: float) -> list[int]:
"""
Compile list of ages corresponding to lower- and and upper-bound values of specified yield component
Yield bounds are interpreted as first occurence of lower-bound value (reading curve from left to right)
and first occurrence of upper-bound value (reading curve from right to left).
:param str yname: Yield component name to use for age lookup
:param float lo: Yield lower-bound value to use for age lookup
:param float hi: Yield upper-bound value to use for age lookup
:return list: List of ages corresponding to specified yield bounds on specified yield curve.
"""
return [x for x, y in enumerate(self.ycomp(yname)) if y >= lo and y <= hi] # type: ignore[arg-type]
[docs]
def reset_areas(self, period: int | None = None) -> None:
"""
Reset areas dictionary. By default will reset all periods (except for period 0),
unless period is specified.
:param int period: Period for which to reset areas dictionary.
"""
periods = self.parent.periods if period is None else [period]
for period in periods:
self._areas[period] = dd(float)
[docs]
def ycomps(self) -> list[str]:
"""
:return: List of yield component names.
"""
return list(self._ycomps.keys())
[docs]
def ycomp(self, yname: str, silent_fail: bool = True) -> core.Curve | None:
"""
Returns the yield components associated with the given yield name.
:param str yname: The name of the yield to retrieve components for.
:param bool silent_fail: If ``True`` (default), returns ``None`` if the yield
name is not found. If ``False``, raises a ``KeyError`` that yield name is not found.
:return: Returns ``None`` if the yield name is not found and ``silent_fail`` is ``True``,
otherwise returns the requested yield component.
:rtype: :py:class:`ws3.core.Curve`
"""
if yname in self._ycomps:
if not self._ycomps[yname]: # complex ycomp not compiled yet
self._compile_complex_ycomp(yname)
return self._ycomps[yname]
else: # not a valid yname
if silent_fail:
return None
else:
raise KeyError("ycomp '{}' not in development type '{}'".format(yname, ' '.join(self.key)))
def _o(self, s: str, default_ycomp: Any = None) -> Any: # resolve string operands
if not default_ycomp:
default_ycomp = self._zero_curve
if common.is_num(s):
return float(s)
elif s.startswith('#'):
return self.parent.constants[s[1:]]
else:
s = s.lower() # just to be safe
ycomp = self.ycomp(s)
return ycomp if ycomp else default_ycomp
def _resolver_multiply(self, yname: str, d: str) -> Any:
args = [self._o(s.lower()) for s in re.split(r'\s?,\s?', _search(r'(?<=\().*(?=\))', d, f'MULTIPLY arguments for yield {yname!r}').group(0))]
##################################################################################################
# NOTE: Not consistent with Remsoft documentation on 'complex-compound yields' (fix me)...
ytype_set = {a.type for a in args if isinstance(a, core.Curve)}
return ytype_set.pop() if len(ytype_set) == 1 else 'c', self._rc(reduce(lambda x, y: x*y, args))
##################################################################################################
def _resolver_divide(self, yname: str, d: str) -> Any:
_tmp = list(zip(re.split(r'\s?,\s?', _search(r'(?<=\().*(?=\))', d, f'DIVIDE arguments for yield {yname!r}').group(0)),
(self._zero_curve, self._unit_curve), strict=False))
args = [self._o(s, default_ycomp) for s, default_ycomp in _tmp]
return args[0].type if not args[0].is_special else args[1].type, self._rc(args[0] / args[1])
def _resolver_sum(self, yname, d):
args = [self._o(s.lower()) for s in re.split(r'\s?,\s?', _search(r'(?<=\().*(?=\))', d, f'SUM arguments for yield {yname!r}').group(0))]
ytype_set = {a.type for a in args if isinstance(a, core.Curve)}
return ytype_set.pop() if len(ytype_set) == 1 else 'c', self._rc(reduce(lambda x, y: x+y, list(args)))
def _resolver_cai(self, yname, d):
arg = self._o(re.split(r'\s?,\s?', _search(r'(?<=\().*(?=\))', d, f'CAI arguments for yield {yname!r}').group(0))[0])
return arg.type, self._rc(arg.mai())
def _resolver_mai(self, yname, d):
arg = self._o(re.split(r'\s?,\s?', _search(r'(?<=\().*(?=\))', d, f'MAI arguments for yield {yname!r}').group(0))[0])
return arg.type, self._rc(arg.mai())
def _resolver_ytp(self, yname, d):
arg = self._o(_search(r'(?<=\().*(?=\))', d, f'YTP argument for yield {yname!r}').group(0).lower())
return arg.type, self._rc(arg.ytp())
def _resolver_range(self, yname, d):
args = [self._o(s.lower()) for s in re.split(r'\s?,\s?', _search(r'(?<=\().*(?=\))', d, f'RANGE arguments for yield {yname!r}').group(0))]
arg_triplets = [args[i:i+3] for i in range(0, len(args), 3)]
return args[0].type, self._rc(reduce(lambda x, y: x*y, [t[0].range(t[1], t[2]) for t in arg_triplets]))
def _compile_complex_ycomp(self, yname: str) -> None:
expression = self._complex_ycomps[yname]
keyword = _search(r'(?<=_)[A-Z]+(?=\()', expression,
f'complex yield keyword in {expression.strip()!r}').group(0)
try:
ytype, ycomp = self._resolvers[keyword](yname, expression)
ycomp.label = yname
ycomp.type = ytype
self._ycomps[yname] = ycomp
except KeyError:
raise ValueError(f'Problem compiling complex yield: {yname}, {expression}') from None
[docs]
def compile_actions(self, verbose: bool = False) -> None:
"""
Compile all actions.
:param bool verbose: Verbosity flag. Defaults to False.
"""
for acode in self.oper_expr:
self.compile_action(acode, verbose)
[docs]
def compile_action(self, acode: str, verbose: bool = False) -> int | None:
"""
Compile action, given action code.
This mostly involves resolving operability expression strings into
lower and upper operability limits, defined as ``(alo, ahi)`` age pair for each period.
Deletes action from ``self`` if never operable.
:param str acode: Action code.
:param bool verbose: Verbosity flag. Defaults to ``False``.
"""
self.operability[acode] = {}
for expr in self.oper_expr[acode]:
self._compile_oper_expr(acode, expr, verbose)
is_operable = False
for p in self.operability[acode]:
if self.operability[acode][p] is not None:
is_operable = True
if not is_operable:
if verbose:
print('not operable (deleting):', acode)
del self.operability[acode]
del self.oper_expr[acode]
return -1
else:
if verbose:
print('operable:', acode)
return 0
def _compile_oper_expr(self, acode: str, expr: str, verbose: bool = False) -> None:
expr = expr.replace('&', 'and').replace('|', 'or')
oper = None
plo, phi = 1, self.parent.horizon # count periods from 1, as in Forest...
alo, ahi = 0, self._max_age
if 'and' in expr:
oper = 'and'
elif 'or' in expr:
oper = 'or'
alo, ahi = self._max_age+1, -1
cond_comps = expr.split(f' {oper} ')
lhs, rel_operators, rhs = list(zip(*[cc.split(' ') for cc in cond_comps], strict=False)) # type: ignore[assignment]
rhs = list(map(float, rhs)) # type: ignore[assignment]
_plo, _phi, _alo, _ahi = None, None, None, None # type: ignore[assignment] # type: ignore[assignment] # type: ignore[assignment]
for i, o in enumerate(lhs):
if o == '_cp':
period = int(rhs[i])
assert period <= self.parent.horizon # sanity check
#################################################################
# Nonsense to relate time-based and age-based conditions with OR?
# Recondider if this actually ever comes up...
assert oper != 'or'
#################################################################
if rel_operators[i] == '=':
_plo, _phi = period, period
elif rel_operators[i] == '>=':
_plo = period
elif rel_operators[i] == '<=':
_phi = period
else:
raise ValueError('Bad relational operator.')
# Guard for None, matching the _age handling below: a one-sided
# comparison (>= or <=) leaves the opposite bound unset.
if _plo is not None:
plo = max(_plo, plo)
if _phi is not None:
phi = min(_phi, phi)
elif o == '_age':
age = int(rhs[i])
if rel_operators[i] == '=':
_alo, _ahi = age, age
elif rel_operators[i] == '>=':
_alo = age
elif rel_operators[i] == '<=':
_ahi = age
else:
raise ValueError('Bad relational operator.')
else: # must be yname
ycomp = self.ycomp(o)
if ycomp is None:
raise ValueError(
f"Operability expression for action {acode!r} references yield "
f"component {o!r}, which is not defined for development type "
f"{' '.join(self.key)}."
)
if rel_operators[i] == '=':
_alo = _ahi = ycomp.lookup(rhs[i])
elif rel_operators[i] == '>=':
_alo = ycomp.lookup(rhs[i])
elif rel_operators[i] == '<=':
_ahi = ycomp.lookup(rhs[i])
else:
raise ValueError('Bad relational operator.')
if oper == 'and' or not oper:
if _alo is not None:
alo = max(_alo, alo)
if _ahi is not None:
ahi = min(_ahi, ahi)
else: # or
if _alo is not None:
alo = min(_alo, alo)
if _ahi is not None:
ahi = max(_ahi, ahi)
assert plo <= phi # should never explicitly declare infeasible period range...
for p in range(plo, phi+1):
self.operability[acode][p] = (alo, ahi) if alo <= ahi else None
[docs]
def add_ycomp(self, ytype: str, yname: str, ycomp: Any, first_match: bool = True) -> None:
"""
Adds a yield component.
:param str ytype: Type of yield component to add (``'c'`` for complex).
:param str yname: Name of the yield component.
:param str ycomp: Yield component to add.
:param bool first_match: Flag indicating whether to only add the component
if it does not already exist. Defaults to ``True``.
"""
if first_match and yname in self._ycomps:
return # already exists (reject)
if ytype == 'c':
self._complex_ycomps[yname] = ycomp
self._ycomps[yname] = None
if isinstance(ycomp, core.Curve):
self._ycomps[yname] = ycomp
[docs]
def grow(self, start_period: int = 1, cascade: bool = True) -> None:
"""
Grow self (default starting period 1, and cascading to end of planning horizon).
Growing basically just increments age and bumps inventory area to the next period.
:param int start_period: The starting period for growth (default is 1).
:param bool cascade: If ``True``, growth cascades to the end of the planning horizon, otherwise
only grows the specified period. Default is ``True``.
"""
end_period = start_period + 1 if not cascade else self.parent.horizon
for p in range(start_period, end_period):
self.reset_areas(p+1) #, self._areas[p], self._areas[p+1] # WTF?
for age, area in list(self._areas[p].items()):
self._areas[p+1][age+self.parent.period_length] = area
[docs]
def overwrite_initial_areas(self, period):
"""
Overwrites the initial areas with area from a specified period.
Basically rolls the planning horizon forward to the specified period.
:param int period: Source period from which to copy initial areas.
"""
self._areas[0] = copy.copy(self._areas[period])
self.initialize_areas()
[docs]
def initialize_areas(self) -> None:
"""
Copy initial inventory to period-1 inventory.
"""
self._areas[1] = copy.copy(self._areas[0])
[docs]
class Output:
"""
Encapsulates data and methods to operate on aggregate outputs from the model.
Emulates behaviour of Forest outputs.
.. warning::
Behaviour of Forest outputs is quite complex. This class needs more work
before it is used in a production setting (i.e., resolution of some complex
output cases is buggy).
"""
def __init__(self,
parent,
code=None,
expression=None,
factor=(1., 1),
description='',
theme_index=-1,
is_basic=False,
is_level=False):
self.parent = parent
self.code = code
self.expression = expression
self._factor = factor
self.description = description
self.theme_index = theme_index
self.is_themed = True if (theme_index is not None and theme_index > -1) else False
self.is_basic = is_basic
if is_basic:
self._compile_basic(expression) # shortcut
elif not is_level:
self._compile(expression) # will detect is_basic
self.is_level = is_level
def _lval(self, s: str) -> Any:
"""
Resolve left operand in sub-expression.
"""
if s.lower() in self.parent.outputs:
return self.parent.outputs[s.lower()]
else: # expression
return s.lower()
def _rval(self, s: str) -> Any:
"""
Resolve right operand in sub-expression.
"""
if common.is_num(s):
return float(s)
elif s.startswith('#'):
return self.parent.constants[s[1:].lower()]
else: # time-based ycomp code
return s.lower()
def _compile(self, expression: str) -> None:
"""
Resolve operands in expression to the extent possible.
Can be basic or summary.
Assuming operand pattern:
lval_1 [*|/ rval_1] +|- .. +|- lval_n [*|/ rval_n]
where
lval := ocode or expression
rval := number or #constant or ycomp
"""
t = re.split(r'\s+(\+|-)\s+', expression)
ocomps = t[::2] # output component sub-expressions
signs = [1.] # implied + in front of expression
signs.extend(1. if s == '+' else -1 for s in t[1::2])
factors = [(1., 1) for i in ocomps]
for i, s in enumerate(ocomps):
tt = re.split(r'\s+(\*|/)\s+', s) # split on */ operator
lval = self._lval(tt[0])
if len(tt) > 1:
factors[i] = self._rval(tt[2]), 1 if tt[1] == '*' else -1
if not isinstance(lval, Output):
if len(ocomps) == 1: # simple basic output (special case)
self.is_basic = True
self._factor = factors[0]
self._compile_basic(s)
return
else: # compound basic output
ocomps[i] = Output(parent=self.parent, # type: ignore[no-untyped-call]
expression=lval,
factor=factors[i],
is_basic=True)
else: # summary output
ocomps[i] = lval #self.parent.outputs[lval]
self._ocomps = ocomps
self._signs = signs
self._factors = factors
def _compile_basic(self, expression: str) -> None:
# clean up (makes parsing easier)
s = re.sub(r'\s+', ' ', expression) # separate tokens by single space
s = s.replace(' (', '(') # remove space to left of left parentheses
t = s.lower().split(' ')
# filter dtypes, if starts with mask
if not (t[0] == '@' or t[0] == '_' or t[0] in self.parent.actions):
t = t[self.parent.nthemes():] # pop
# extract @AGE or @YLD condition, if present
self._ages = None
self._condition = None
if t[0].startswith('@age'):
lo, hi = [int(a)+i for i, a in enumerate(t[0][5:-1].split('..'))]
hi = min(hi, self.parent.max_age+1) # they get carried away with range bounds...
self._ages = list(range(lo, hi))
t = t[1:] # pop
elif t[0].startswith('@yld'):
ycomp, args = t[0][5:-1].split(',')
self._condition = tuple([ycomp] + [float(a) for a in args.split('..')])
self._ages = None
t = t[1:] # pop
if not self._ages and not self._condition:
self._ages = self.parent.ages
# extract _INVENT or acode
if t[0].startswith('_'): # _INVENT
self._is_invent = True
self._invent_acodes = t[0][8:-1].split(',') if len(t[0]) > 7 else None
self._acode = None
else: # acode
self._is_invent = False
self._invent_acodes = None
self._acode = t[0]
t = t[1:] # pop
# extract _AREA or ycomp
if t[0].startswith('_'): # _AREA
self._is_area = True
self._ycomp = None
else: # acode
self._is_area = False
self._ycomp = t[0]
t = t[1:] # pop
def _evaluate_basic(self, period: int, factors: list[tuple[Any, int]], verbose: int = 0, cut_corners: bool = True) -> Any:
result = 0.
if self._invent_acodes:
acodes = [acode for acode in self._invent_acodes
if self.parent.applied_actions[period][acode]]
if cut_corners and not acodes:
return 0. # area will be 0...
for k in list(self.parent.dtypes.keys()):
dt = self.parent.dtypes[k]
if cut_corners and not self._is_invent and k not in self.parent.applied_actions[period][self._acode]:
if verbose:
print('bailing on', period, self._acode, ' '.join(k))
continue # area will be 0...
if isinstance(self._factor[0], float):
f = pow(*self._factor)
else:
f = pow(dt.ycomp(self._factor[0])[period], self._factor[1])
for factor in factors:
if isinstance(factor[0], float):
f *= pow(*factor)
else:
f *= pow(dt.ycomp(factor[0])[period], factor[0])
if cut_corners and not f:
if verbose:
print('f is null', f)
continue # one of the factors is 0, no point calculating area...
ages = self._ages if not self._condition else dt.resolve_condition(*self._condition)
if ages is None:
continue # no ages resolved for this condition, nothing to accumulate
for age in ages:
area = 0.
if self._is_invent:
if cut_corners and not dt.area(period, age):
continue
if self._invent_acodes:
any_operable = False
for acode in acodes:
if acode not in dt.operability:
continue
if dt.is_operable(acode, period, age):
any_operable = True
if any_operable:
area += dt.area(period, age)
else:
area += dt.area(period, age)
else:
raise AssertionError() # not implemented yet...
y = 1. if self._is_area else dt.ycomp(self._ycomp)[age]
result += y * area * f
return result
def _evaluate_summary(self, period: int, factors: list[tuple[Any, int]]) -> Any:
result = 0.
for i, ocomp in enumerate(self._ocomps):
result += ocomp(period, [self._factors[i]] + factors) # type: ignore[operator]
return result
def _evaluate_basic_themed(self, period):
pass
def _evaluate_summed_themed(self, period):
pass
def __call__(self, period, factors=None):
if factors is None:
factors = [(1.0, 1)]
if self.is_basic:
return self._evaluate_basic(period, factors)
else:
return self._evaluate_summary(period, factors)
def __add__(self, other):
# assume Output + Output
if self.is_themed:
return [i + j for i, j in zip(self(), other(), strict=False)] # type: ignore[call-arg]
else:
return self() + other() # type: ignore[call-arg]
def __sub__(self, other):
# assume Output - Output
if self.is_themed:
return [i - j for i, j in zip(self(), other(), strict=False)] # type: ignore[call-arg]
else:
return self() - other() # type: ignore[call-arg]
def _normalize_cflw_e(e: Any, periods: list[int]) -> dict[int, tuple[float | None, float | None, int | None]]:
"""Normalize a ``cflw_e`` value into ``{period: (alpha, beta, ref_period)}``.
Two forms are supported for each even-flow output value ``e``:
- Legacy tuple ``(eps_dict, ref_period)``: a symmetric +/-``eps`` band tying
each period's output to the single anchor period ``ref_period`` (int).
``eps_dict`` maps period -> epsilon. This is the historical behaviour and
is unchanged.
- Extended dict ``{"decrease": d, "increase": i, "ref": r}``: separate
period-keyed tolerances. ``decrease`` (``alpha``) is the maximum allowed
fractional period-over-period *decrease* (``H_t - (1-alpha) H_ref >= 0``);
``increase`` (``beta``) is the maximum allowed fractional *increase*
(``H_t - (1+beta) H_ref <= 0``). Either may be ``None`` (constraint not
added). ``ref`` is an int anchor period or ``"consecutive"`` (each period
anchored to the previous period), enabling the classic FORPLAN
sequential-flow forms (non-declining yield, bounded decline,
bounded deviation) used e.g. in Daugherty (1991, eq. 3-4/3-5, Table 5.6).
Returns a dict keyed on period of ``(alpha, beta, ref_period)`` tuples;
periods with no constraint (both tolerances ``None``) are omitted.
"""
if isinstance(e, dict):
dec = e.get("decrease")
inc = e.get("increase")
ref = e.get("ref", "consecutive")
else: # legacy tuple (eps_dict, ref_period)
eps_dict, ref = e
dec = eps_dict
inc = eps_dict
spec: dict[int, tuple[float | None, float | None, int | None]] = {}
for k, t in enumerate(periods):
if ref == "consecutive":
ref_t: int | None = periods[k - 1] if k > 0 else None
else:
ref_t = ref
alpha = None if dec is None else dec.get(t)
beta = None if inc is None else inc.get(t)
if alpha is None and beta is None:
continue
spec[t] = (alpha, beta, ref_t)
return spec
[docs]
class ForestModel:
"""
This is the core class of the ws3 package.
Includes methods import data from various sources, simulate growth and apply actions.
The model can be used in either a (prescriptive) simulation-based approach or a (descriptive)
optimization-based approach.
This class encapsulates all the information used to simulate scenarios from a given dataset (i.e.,
stratified intial inventory, growth and yield functions, action eligibility, transition matrix, action
schedule, etc.), as well as a large collection of functions to import and export data, generate activity
schedules, and simulate application of these schedules (i.e., run scenarios).
At the heart of the ``ForestModel`` class is a list of ``DevelopentType`` instances. Each ``DevelopmentType``
instance encapsulates information about one development type (i.e., a forest stratum, which is an aggregate of
smaller *stands* that make up the raw forest inventory input data). The ``DevelopmentType`` class also
stores a list of operable *actions*, maps *state variable transitions* to these actions, stores growth and
yield functions, and knows how to *grow itself* when time is incremented during a simulation.
A typical use case starts with creating an instance of the ``ForestModel`` class. Then, we need to load data
into this instance, define one or more scenarios (using a mix of heuristic and optimization approaches),
run the scenarios, and export output data to a format suitable for analysis (or link to the next model in
a larger modelling pipeline).
"""
_ytypes = {'*Y':'a', '*YT':'t', '*YC':'c'}
tree: Any = (lambda f: f(f))(lambda a: (lambda: dd(a(a))))
def __init__(self,
model_name,
model_path,
base_year,
horizon=common.HORIZON_DEFAULT,
period_length=common.PERIOD_LENGTH_DEFAULT,
max_age=common.MAX_AGE_DEFAULT,
area_epsilon=common.AREA_EPSILON_DEFAULT,
curve_epsilon=common.CURVE_EPSILON_DEFAULT):
"""
Initializes the ``ForestModel`` with the provided parameters.
:param str model_name: The name of model.
:param str model_path: The path to input data of model.
:param int base_year: The base year of teh model.
:param int horizon: The length (in number of periods) of the simulation horizon.
:param int max_age: The maximum age considered in the model.
:param int area_epsilon:
:param int curve_epsilon:
"""
self.model_name = model_name
self.model_path = model_path
self.base_year = base_year
self.set_horizon(horizon)
self.period_length = period_length
self.max_age = max_age
self.ages = list(range(max_age+1))
self._period_to_years_factor = None
self.yields = []
self.ynames = set()
self.actions: dict[str, Any] = {} # type: ignore[var-annotated]
self.transitions: dict[Any, Any] = {} # type: ignore[var-annotated]
self.oper_expr = {}
self._themes = []
self._theme_basecodes = []
self.dtypes = {}
self.constants = {}
self.output_groups: dict[str, Any] = {} # type: ignore[var-annotated]
self.outputs: dict[str, Any] = {} # type: ignore[var-annotated]
self.applied_actions: dict[int, dict[str, Any]] = {p:{acode:{} for acode in list(self.actions.keys())} for p in self.periods} # type: ignore[var-annotated]
self.reset_actions()
self.curves = {}
self.problems = {}
c_zero = self.register_curve(core.Curve('zero',
is_special=True,
type=''))
c_unit = self.register_curve(core.Curve('unit',
points=[(0, 1)],
is_special=True,
type=''))
c_ages = self.register_curve(core.Curve('ages',
points=[(0, 0), (max_age, max_age)],
is_special=True,
type=''))
self.common_curves = {'zero':c_zero,
'unit':c_unit,
'ages':c_ages}
self.area_epsilon = area_epsilon
self.curve_epsilon = curve_epsilon
self.areaselector = GreedyAreaSelector(self)
self.inoperable_dtypes = []
[docs]
def nthemes(self) -> int:
"""
:return: Number of themes
"""
return len(self._themes)
[docs]
def reset(self) -> None:
"""
Resets the forest model by clearing applied actions and reinitializing areas.
"""
self.reset_actions()
self.initialize_areas()
[docs]
def set_horizon(self, horizon: int) -> None:
"""
Sets the horizon of the model.
This method updates the horizon of the model to the specified value and adjusts the list of periods accordingly.
"""
self.horizon = int(horizon)
self.periods = list(range(1, horizon+1))
def _resolve_period_multiplier(self, convert_periods_to_years: Any) -> int:
"""Resolve the Woodstock period-to-year multiplier for import helpers."""
if convert_periods_to_years is None:
return self._period_to_years_factor or 1
try:
multiplier = int(convert_periods_to_years)
except Exception as exc:
raise ValueError("convert_periods_to_years must be an integer number of years per period") from exc
if multiplier <= 0:
raise ValueError("convert_periods_to_years must be positive")
if self._period_to_years_factor is None:
self._period_to_years_factor = multiplier
elif self._period_to_years_factor != multiplier:
raise ValueError(
f"convert_periods_to_years={multiplier} conflicts with existing multiplier {self._period_to_years_factor}"
)
return multiplier
[docs]
def compile_actions(self, mask: Any = None, verbose: bool = False) -> None:
"""
Compile actions for the development types filtered by mask.
"""
dtype_keys = self.unmask(mask) if mask else list(self.dtypes.keys())
for dtk in dtype_keys:
dt = self.dtypes[dtk]
dt.compile_actions(verbose=verbose)
def _compile_schedule_from_problem(self, problem: Any, formulation: int = 1, skip_null: str = 'null') -> Any:
"""
Compiles a ``ws3``-compatible schedule data object from a solved ``ws3.opt.Problem`` instance.
This is just a dispatcher function---the actual compilation is done by a formulation-specific function
(assumes *Model I* formulation if not specified).
"""
cmp_sch_dsp = {1:self._cmp_sch_m1, 2:self._cmp_sch_m2}
return cmp_sch_dsp[formulation](problem, skip_null) # type: ignore[operator]
[docs]
def add_problem(self, name: str, coeff_funcs: Any, cflw_e: Any = None, cgen_data: Any = None,
solver: Any = opt.SOLVER_HIGHS, formulation: int = 1, z_coeff_key: str = 'z', acodes: Any = None,
sense: int = opt.SENSE_MAXIMIZE, mask: Any = None, workers: int = 1, verbose: bool = False) -> Any:
"""
Add an optimization problem to the model.
:param str name: Used as key to store :py:class:`ws3.opt.Problem` instances
in a dict in the :py:class:`ws3.forest.ForestModel` instanace, so make
sure it is unique within a given model or you will overwrite dict values
(assuming you want to stuff multiple problems, and their solutions, into
your model at the same time).
:param dict coeff_funcs: Dict of function references, keyed on row name strings. These are the functions that generate
the LP optimization problem matrix coefficients (for the objective function and constraint rows). This one gets
complicated, and is a likely source of bugs. Make sure the row name key strings are all unique or you will make
a mess. You can name the constraint rows anything you want, but the objective function row has to be named 'z'.
All coefficient functions must accept exactly two args, in this order: a ``ws3.forest.ForestModel`` instance and a
path (a tuple of ``ws3.core.Node`` object instances). The 'z' coefficient function is special in that it must return a single float value.
All other (i.e., constraint) coefficient functions just return a dict of floats, keyed on period ints (can be
sparse, i.e., not necessary to include key:value pairs in output dict if value is 0.0). It is useful (but not
necessary) to use ``functools.partial`` to specialize a smaller number of more general function definitions (with
more args, that get "locked down" and hidden by ``partial``) as we have done in the example in this notebook.
:param dict cflw_e: Even-flow (flow-constraint) specification, keyed on row name strings (must match row name key
values used to define coefficient functions for flow constraints in coeff_func dict). Two value forms are
supported:
- Legacy ``(eps_dict, ref_period)`` tuple: a symmetric +/-``eps`` band tying each period's output to the
single anchor period ``ref_period`` (int). ``eps_dict`` maps period -> epsilon (must include all periods).
``{'foo':({1:0.01, ..., 10:0.01}, 1), 'bar':({1:0.05, ..., 10:0.05}, 1)}``
- Extended ``{"decrease": d, "increase": i, "ref": r}`` dict: separate period-keyed tolerances.
``decrease`` (alpha) bounds the fractional period-over-period *decrease* (``H_t - (1-alpha) H_ref >= 0``);
``increase`` (beta) bounds the fractional *increase* (``H_t - (1+beta) H_ref <= 0``); either may be
``None`` to omit that bound. ``ref`` is an int anchor period or ``"consecutive"`` (each period anchored
to the previous period). This enables the classic FORPLAN sequential-flow policies — e.g. non-declining
yield ``{"decrease": {t: 0.0}, "increase": None, "ref": "consecutive"}`` and bounded deviation
``{"decrease": {t: eps}, "increase": {t: eps}, "ref": "consecutive"}`` (cf. Daugherty 1991, Table 5.6).
:param dict cgen_data: Dict of dict of dicts. The outer-level dict is keyed on row name strings (must match row names used
in coeff_funcs. The middle second level of dicts always has keys 'lb' and 'ub', and the inner level of dicts
specifies lower- and upper-bound general constraint RHS (float) values, keyed on period (int). See example below.
``{'foo':{'lb':{1:1., ..., 10:1.}, 'ub':{1:2., ..., 10:2.}}, 'bar':{{'lb':{1:1., ..., 10:1.}, 'ub':{1:2., ..., 10:4.}}}}``
:param int acodes: List of strings. Action codes to be included in optimization
problem formulation (actions must defined in the :py:class:`ws3.forest.ForestModel`
instance, but can be only a subset).
:param int sense: Must be one of :py:attr:`ws3.opt.SENSE_MAXIMIZE` or
:py:attr:`ws3.opt.SENSE_MINIMIZE`, or equivalent int values
(use the constants to keep code more legible).
:param tuple mask: Tuple of strings constituting a valid mask for your
:py:class:`ws3.forest.ForestModel` instance. Can be ``None`` if you do not want
to filter :py:class:`ws3.forest.DevelopmentType` instances.
:param int workers: Number of worker threads to use for parallel processing.
:return: ws3.opt.Problem. Reference to a new Problem instance that was created. Also stored in the ForestModel instance (problems attribute,
keyed on problem name).
"""
# --- Prepare serialization for parallel execution ---
if workers > 1:
problems_backup = self.problems
self.problems = None # type: ignore[assignment]
blob_bytes = dill.dumps(self) # Serialize model
self.problems = problems_backup
rebased_funcs = {k: sanitize_func(f) for k, f in coeff_funcs.items()}
serialized_funcs = {k: dill.dumps(f) for k, f in rebased_funcs.items()}
else:
blob_bytes = None
serialized_funcs = None
# --- Reset model state for problem creation ---
self.reset()
# --- Dispatch maps for formulation type ---
bld_p_dsp = {1: self._bld_p_m1, 2: self._bld_p_m2}
cmp_cflw_dsp = {1: self._cmp_cflw_m1, 2: self._cmp_cflw_m2}
cmp_cgen_dsp = {1: self._cmp_cgen_m1, 2: self._cmp_cgen_m2}
assert formulation == 1, "Only Model I supported for now"
# --- Persistent worker pool (None if workers=1) ---
with PersistentWorkerPool(workers, blob_bytes, serialized_funcs) as executor:
if verbose:
print('add_problem: build problem')
p = bld_p_dsp[formulation]( # type: ignore[operator]
name,
coeff_funcs,
solver,
z_coeff_key,
acodes,
sense,
mask,
workers,
executor, # None if serial
verbose
)
if verbose:
print('add_problem: compile flow constraints')
cmp_cflw_dsp[formulation](p, cflw_e, workers=workers, executor=executor, verbose=verbose) # type: ignore[operator]
if verbose:
print('add_problem: compile general constraints')
cmp_cgen_dsp[formulation](p, cgen_data, workers=workers, executor=executor, verbose=verbose) # type: ignore[operator]
# --- Save and return the problem ---
self.problems[name] = p
return p
def _bld_p_m1(
self,
name,
coeff_funcs,
solver,
z_coeff_key='z',
acodes=None,
sense=opt.SENSE_MAXIMIZE,
mask=None,
workers=1,
executor=None,
verbose=False
):
"""
Build a Model I optimization problem with batched, parallel tree processing.
"""
p = opt.Problem(name, sense=sense, solver=solver)
# Step 1: Generate trees and variables
if verbose:
print('generate trees using', workers, 'workers')
p.trees, p._vars, p._leaf_ids = self._gen_vars_m1( # type: ignore[assignment]
# type: ignore[assignment]
coeff_funcs,
acodes=acodes,
mask=mask,
workers=workers,
executor=executor,
verbose=verbose
)
# Step 2: Process trees into coverage constraints
if verbose:
print('process trees')
tree_items = list(p.trees.items()) # type: ignore[attr-defined]
batches = list(auto_batch(tree_items, workers, max_batch_factor=4))
tasks = [(batch, z_coeff_key) for batch in batches]
results = []
if workers == 1:
for task in tasks:
results.extend(worker_summarize_tree_batch(task)) # type: ignore[arg-type]
else:
exec_ctx = executor or ProcessPoolExecutor(max_workers=workers, mp_context=get_context(MP_CONTEXT))
futures = [exec_ctx.submit(worker_summarize_tree_batch, task) for task in tasks] # type: ignore[arg-type]
for f in as_completed(futures):
results.extend(f.result())
if executor is None:
exec_ctx.shutdown()
# Step 3: Apply results to the Problem object
if verbose:
print('_bld_p_m1: build problem')
for cname, coeffs, z_coeffs in results:
p.add_constraint(name=cname, coeffs=coeffs, sense=opt.SENSE_EQ, rhs=1.0)
p._z.update(z_coeffs)
p.coeff_funcs = coeff_funcs
p.formulation = 1
if verbose:
print('_bld_p_m1: done building problem')
return p
def _bld_p_m2(self, problem):
"""
Builds optimization problem, using Model II (m2) formulation.
.. warning::
Not implemented yet.
"""
pass # not implemented
def _gen_vars_m1(self, coeff_funcs: Any, acodes: Any = None, mask: Any = None, workers: int = 1, executor: Any = None, verbose: bool = False) -> list[tuple[Any, ...]]:
"""
Generate trees, variables, and leaf IDs for Model I problems.
Parallelized with model and coeff_funcs preloaded per worker.
"""
dtype_keys = self.dtypes.keys() if not mask else self.unmask(mask)
# --- Step 1: Build (dtk, age) task list ---
self.reset()
tract_tasks = [
(dtk, age)
for dtk in dtype_keys
for age in self.dtypes[dtk]._areas[1].keys()
if self.dtypes[dtk].area(1, age)
]
# --- Step 3: Serial or Parallel execution ---
results = []
if workers == 1:
# Serial mode (do not use worker functions or it breaks _cbm_sit_yield add_problem hack)
for (dtk, age) in tract_tasks:
self.reset()
area = self.dtypes[dtk].area(1, age)
if not area:
continue
tree = self._bld_tree_m1(
area, dtk, age, coeff_funcs,
tree=None, period=1,
acodes=acodes, compile_c_ycomps=True,
verbose=verbose
)
results.append((dtk, age, tree))
else:
task_batches = auto_batch(tract_tasks, workers, max_batch_factor=4)
if not executor:
# prepare serialized model and coeff_funcs
problems_backup = self.problems
self.problems = None # type: ignore[assignment]
blob_bytes = dill.dumps(self)
self.problems = problems_backup
rebased_funcs = {k: sanitize_func(f) for k, f in coeff_funcs.items()}
serialized_funcs = {k: dill.dumps(f) for k, f in rebased_funcs.items()}
with ProcessPoolExecutor(
max_workers=workers,
mp_context=get_context(MP_CONTEXT),
initializer=init_worker_gen_vars,
initargs=(blob_bytes, serialized_funcs, workers),
) as executor:
futures = [executor.submit(worker_gen_vars, batch, acodes) for batch in task_batches]
for f in as_completed(futures):
res = f.result()
for item in res:
if isinstance(item, Exception):
raise item
if item is not None:
results.append(item)
else: # use executor that was passed in as arg
futures = [executor.submit(worker_gen_vars, batch, acodes) for batch in task_batches]
for f in as_completed(futures):
res = f.result()
for item in res:
if isinstance(item, Exception):
raise item
if item is not None:
results.append(item)
# --- Step 4: Restore problems and rebuild trees/vars ---
trees, vars, leaf_ids = {}, {}, {}
for dtk, age, tree in results:
i = (dtk, age)
trees[i] = tree
for path in tree.paths():
j = tuple(n.data('acode') for n in path)
leaf_id = path[-1].data('leaf_id')
vname = f"x_{leaf_id}"
leaf_ids[(i, j)] = leaf_id
vars[vname] = opt.Variable(vname, opt.VTYPE_CONTINUOUS, 0.0, 1.0)
return trees, vars, leaf_ids # type: ignore[return-value]
def _gen_vars_m2(self):
pass
def _bld_tree_m1(
self,
area: float,
dtk: tuple[str, ...],
age: int,
coeff_funcs: Any,
tree: Any = None,
period: int = 1,
acodes: Any = None,
compile_c_ycomps: bool = True,
verbose: bool = False) -> Any:
"""
Build a tree of feasible action sequences (full-length paths = |periods|).
"""
# --- Step 0: Initialize tree if needed ---
if tree is None:
dt = self.dtypes[dtk]
dt.reset_areas()
self.dtypes[dtk]._areas[1][age] = area
self.reset_actions()
tree = core.Tree() # type: ignore[no-untyped-call]
acodes = list(self.actions.keys()) if not acodes else acodes
# --- Step 1: Depth-First Search (DFS) to build the tree structure ---
for acode in acodes:
if self.dt(dtk).is_operable(acode, period, age):
# Reset actions for this period and grow stand if needed
self.reset_actions(period)
if period > 1:
self.dt(dtk).grow(period - 1, False)
# Apply action to get next state
errorcode, missingarea, tstate = self.apply_action(
dtk, acode, period, age, area,
compile_c_ycomps=compile_c_ycomps,
override_operability=False,
fuzzy_age=False,
recourse_enabled=False
)
if errorcode and verbose:
print(
'apply_action error', dtk, acode, period, age, area,
errorcode, missingarea, tstate
)
_dtk, tprop, _age = tstate[0]
assert tprop == 1. # no split handling yet
# Push node to the tree
tree.grow({
'dtk': dtk, '_dtk': _dtk, 'acode': acode, 'period': period,
'age': age, '_age': _age, 'products': None, 'area': area
})
# Recurse deeper or compute leaf coefficients
if period < self.periods[-1]: # not at last period, continue DFS
self.dt(_dtk).grow(period, False)
self._bld_tree_m1(
area, _dtk, _age + self.period_length,
coeff_funcs, tree, period + 1, acodes,
compile_c_ycomps=compile_c_ycomps)
elif period == self.periods[-1]: # reached a leaf
path = tree.path()
leaf = path[-1]
assert leaf.is_leaf()
# Serial leaf coefficient computation
leaf._data.update({k: coeff_funcs[k](self, path) for k in coeff_funcs})
i = (path[0].data('dtk'), path[0].data('age'))
j = tuple(node.data('acode') for node in path)
leaf._data['leaf_id'] = common.hex_id((i, j))
# Pop node from the tree (DFS backtrack)
tree.ungrow()
return tree
def _cmp_cflw_m1(self, problem, cflw_e, workers=1, executor=None, verbose=False):
"""
Compile flow (even-flow) constraints in parallel using batched workers.
Optimized for less overhead while respecting the original (i, j) API.
"""
if not cflw_e:
return
periods = self.periods
cflw_keys = list(cflw_e.keys())
if verbose:
print("_cmp_cflw_m1: phase 1")
# Phase 1: Compute mu values in parallel
tree_items = list(problem.trees.items())
batches = auto_batch(
tree_items,
workers,
size_fn=lambda x: len(x[1].nodes()),
max_batch_factor=1
)
tasks = [(batch, cflw_keys, periods) for batch in batches] # type: ignore[arg-type]
if workers == 1:
results = []
for task in tasks:
results.extend(worker_cmp_cflw_batch(task)) # type: ignore[arg-type]
else:
# Use existing executor if passed
exec_ctx = executor or ProcessPoolExecutor(max_workers=workers, mp_context=get_context(MP_CONTEXT))
futures = [exec_ctx.submit(worker_cmp_cflw_batch, task) for task in tasks] # type: ignore[arg-type]
# Collect results without repeated extend() overhead
results_nested = [f.result() for f in as_completed(futures)]
results = [item for batch in results_nested for item in batch]
if executor is None:
exec_ctx.shutdown()
# Phase 2: Merge results into mu dict
if verbose:
print("_cmp_cflw_m1: phase 2")
mu: dict[int, dict[str, dict[tuple[int, int], float]]] = {t: {o: {} for o in cflw_keys} for t in periods} # type: ignore[var-annotated]
for t, o, i, j, val in results: # type: ignore[index]
mu[t][o][(i, j)] = val # type: ignore[index]
# Phase 3: Build constraints (parallel with batching)
if verbose:
print("_cmp_cflw_m1: phase 3")
leaf_ids = problem._leaf_ids
xnames = {k: f"x_{v}" for k, v in leaf_ids.items()}
add_constraint = problem.add_constraint
# Build Phase 3 tasks
tasks = []
for o, e in cflw_e.items():
spec = _normalize_cflw_e(e, periods) # {t: (alpha, beta, ref_t)}
for t in periods:
if t not in spec:
continue
alpha, beta, ref_t = spec[t]
if ref_t is None or ref_t not in mu:
continue
mu_t_o = mu[t][o]
mu_ref_o = mu[ref_t][o]
tasks.append((t, o, mu_t_o, mu_ref_o, alpha, beta, xnames)) # type: ignore[arg-type]
results = []
if workers == 1:
# Serial processing
for task in tasks:
results.extend(worker_cmp_cflw_phase3(task)) # type: ignore[arg-type]
else:
# Create batches of tasks for more efficient multiprocessing
batches = auto_batch(tasks, workers, max_batch_factor=2)
exec_ctx = executor or ProcessPoolExecutor(max_workers=workers, mp_context=get_context(MP_CONTEXT))
futures = [exec_ctx.submit(worker_cmp_cflw_phase3_batch, batch) for batch in batches]
for f in as_completed(futures):
results.extend(f.result())
if executor is None:
exec_ctx.shutdown()
# Add constraints sequentially
for name, coeffs, sense, rhs in results: # type: ignore[misc]
add_constraint(name=name, coeffs=coeffs, sense=sense, rhs=rhs)
def _cmp_cflw_m2(self):
"""
Compiles flow constraints (lb and ub, per targeted output, per targeted period) for a Model I problem.
.. warning:
Not implemented yet.
"""
pass # not implemented
# def _cmp_cgen_m1(self, problem, cgen_data, workers=1, executor=None, verbose=False):
# print('foo')
# if not cgen_data: return
# mu = {t:{o:{} for o in list(cgen_data.keys())} for t in self.periods}
# for i, tree in list(problem.trees.items()):
# for path in tree.paths():
# j = tuple(n.data('acode') for n in path)
# for o in list(cgen_data.keys()):
# _mu = path[-1].data(o)
# for t in self.periods:
# mu[t][o][i, j] = _mu[t] if t in _mu else 0.
# for o, b in list(cgen_data.items()):
# for t in self.periods:
# _mu = {'x_%s' % common.hex_id((i, j)):mu[t][o][i, j] for i, j in mu[t][o]}
# if b['lb'] is not None and t in b['lb']:
# problem.add_constraint(name='gen-lb_%03d_%s' % (t, o), coeffs=_mu, sense=opt.SENSE_GEQ, rhs=b['lb'][t])
# if b['ub'] is not None and t in b['ub']:
# problem.add_constraint(name='gen-ub_%03d_%s' % (t, o), coeffs=_mu, sense=opt.SENSE_LEQ, rhs=b['ub'][t])
def _cmp_cgen_m1(self, problem, cgen_data, workers=1, executor=None, verbose=False):
"""
Compile general (CGEN) constraints.
Phase 1+2 are kept IDENTICAL to the known-good serial logic to avoid coefficient drift.
Only Phase 3 (row emission) is parallelized/batched.
"""
if not cgen_data:
return
periods = self.periods
cgen_keys = list(cgen_data.keys())
# --- Phase 1+2: build mu exactly like the reference implementation ---
if verbose:
print("_cmp_cgen_m1: phase 1 and 2")
mu: dict[int, dict[str, dict[Any, Any]]] = {t: {o: {} for o in cgen_keys} for t in periods} # type: ignore[var-annotated]
for i, tree in problem.trees.items():
for path in tree.paths():
j = tuple(n.data("acode") for n in path)
leaf = path[-1]
for o in cgen_keys:
_mu = leaf.data(o) # dict(period -> value)
for t in periods:
mu[t][o][(i, j)] = _mu[t] if t in _mu else 0.0
# --- Phase 3: build rows (can parallelize safely) ---
# Build tasks
if verbose:
print("_cmp_cgen_m1: phase 1")
tasks = []
for o, bounds in cgen_data.items():
lb, ub = bounds.get('lb'), bounds.get('ub')
for t in periods:
tasks.append((t, o, mu[t][o], lb, ub))
results = []
if workers == 1:
for task in tasks:
results.extend(worker_cmp_cgen_phase3(task))
else:
batches = auto_batch(tasks, workers) # use your tuned auto_batch; no forced factor
exec_ctx = executor or ProcessPoolExecutor(max_workers=workers, mp_context=get_context(MP_CONTEXT))
futures = [exec_ctx.submit(worker_cmp_cgen_phase3_batch, batch) for batch in batches]
for f in as_completed(futures):
results.extend(f.result())
if executor is None:
exec_ctx.shutdown()
# Add rows
add_constraint = problem.add_constraint
for name, coeffs, sense, rhs in results:
add_constraint(name=name, coeffs=coeffs, sense=sense, rhs=rhs)
def _cmp_cgen_m2(self):
pass # not implemented
def _cmp_sch_m1(self, problem, skip_null):
_sch: list[list[Any]] = [[] for t in self.periods] # type: ignore[var-annotated]
sln = problem.solution()
if not sln:
return None
for i, tree in list(problem.trees.items()):
for path in tree.paths():
x = 'x_{}'.format(common.hex_id((i, tuple(n.data('acode') for n in path))))
if not sln[x]:
continue
for t, n in enumerate(path):
d = n.data()
if skip_null and d['acode'] == skip_null:
continue
etype = '_existing' if self.dt(i[0]).area(0) else '_future'
_sch[t].append((d['dtk'], d['age'], d['area'] * sln[x], d['acode'], d['period'], etype))
return list(itertools.chain.from_iterable(_sch))
def _cmp_sch_m2(self, problem):
pass
[docs]
def add_null_action(self, acode='null', maxage=None):
"""
Adds a null action with the specified action code, minimum age (default is None),
and maximum age (default is None).
:param str acode: Action code for the new null action. Defaults to ``'null'``.
:param int maxage: Maximum age at which the new null action is operable.
"""
mask = tuple(['?' for _ in range(self.nthemes())])
oe = f'_age >= 0 and _age <= {self.max_age}'
target = [(mask, 1.0, None, None, None, None, None)]
self.actions[acode] = Action(acode)
self.oper_expr[acode] = {mask:oe}
self.transitions[acode] = {mask:{'':target}}
for dtk in self.dtypes:
self.dtypes[dtk].oper_expr[acode] = [oe]
self.dtypes[dtk].transitions[acode, -1] = target
for p in self.applied_actions:
self.applied_actions[p][acode] = {}
[docs]
def is_harvest(self, acode):
"""
Returns the value of :py:attr:`ws3.forest.Action.is_harvest` for a given action code.
:param str acode: The action code for which to look up the ``is_harvest`` attribute value.
Should return ``True`` or ``False`` (not guaranteed---must have been correctly set
when action defined).
"""
return self.actions[acode].is_harvest
[docs]
def dt(self, dtype_key: tuple[str, ...]) -> Any:
"""
Returns development type, given key (returns None on invalid key).
:param tuple dtype_key: Development type key.
:return: Development type.
:rtype: :py:class:`ws3.forest.DevelopmentType`
"""
try:
return self.dtypes[dtype_key]
except Exception:
return None
[docs]
def age_class_distribution(self, period, mask=None, omit_null=False):
"""
Returns age class distribution (dict of areas, keys on age).
:param int period: The period for which to retrieve the age class distribution.
:param tuple mask: A mask to filter development types. Default is ``None``.
:param bool omit_null: If ``True``, omits null areas from the distribution. Default is ``False``.
:return dict: A dictionary where keys are ages and values are the corresponding area distributions.
"""
result = dict.fromkeys(self.ages, 0.0)
dtype_keys = self.unmask(mask) if mask else list(self.dtypes.keys())
for dtk in dtype_keys:
dt = self.dtypes[dtk]
for age in dt._areas[period]:
result[age] += dt._areas[period][age]
if omit_null:
result = {k:v for k, v in result.items() if v}
return result
[docs]
def operable_dtypes(self, acode: str, period: int, mask: Any = None) -> dict[tuple[str, ...], list[int]]:
"""
Looks up operable development types for a given action code and period (and optional mask).
:param acode: Action code for which to look up operable development types.
:return: Dictionary keyed on development type key, values are lists of operable ages.
"""
result = {}
dtype_keys = self.unmask(mask) if mask else list(self.dtypes.keys())
for dtk in dtype_keys:
dt = self.dtypes[dtk]
operable_ages = dt.operable_ages(acode, period)
if operable_ages:
result[dt.key] = operable_ages
return result
[docs]
def inventory(self, period, yname=None, age=None, mask=None, dtype_keys=None, verbose=0):
"""
Flexible method that compiles inventory at given period. Unit of return data defaults to
area if ``yname`` not given, but takes on unit of specificed yield component otherwise.
Can optionally be constrained by age and development type mask.
:param int period: Period for which to compile inventory.
:param str yname: Name of yield component to use when compiling inventory.
:param int age: Optional age filter.
:param tuple mask: Optional development type mask filter (``dtype_keys`` must be ``None`` if this is used).
:param list dtype_keys: Optional development type key filter (``mask`` must be ``None`` if this is used).
:param int verbose: Optional verbosity setting (passed to call to :py:meth:`ws3.forest.ForestModel.unmask`).
:return float: Result of compiling inventory query.
"""
result = 0.
assert not (mask and dtype_keys) # too confusing to allow both to be specified...
if mask:
_dtype_keys = self.unmask(mask, verbose=verbose)
elif dtype_keys:
_dtype_keys = dtype_keys
else:
_dtype_keys = list(self.dtypes.keys())
shift = self.period_length
for dtk in _dtype_keys:
dt = self.dtypes[dtk]
if period == 0:
inventory_map = dt._areas[0]
else:
inventory_map = dd(float)
for src_age, src_area in dt._areas[period].items():
aged_age = src_age + shift
inventory_map[aged_age] += src_area
if yname:
ycomp = dt.ycomp(yname)
if not ycomp:
continue
ymax = getattr(ycomp, "xmax", None)
else:
ycomp = None
ymax = None
if age is not None:
value = 0.0
if age in inventory_map:
if ycomp:
lookup_age = min(age, ymax) if ymax is not None else age
factor = ycomp[lookup_age]
else:
factor = 1.0
value = inventory_map[age] * factor
result += value
else:
if ycomp:
if ymax is not None:
result += sum(area * ycomp[min(a, ymax)] for a, area in inventory_map.items())
else:
result += sum(area * ycomp[a] for a, area in inventory_map.items())
else:
result += sum(inventory_map.values())
return result
[docs]
def operable_area(self, acode, period, age=None, mask=None):
"""
Returns total operable area, given action code and period (and optionally age).
:param str acode: Action code for which to compile operable area.
:param int period: Period for which to compile operable area.
:param int age: Optional age filter.
:param tuple mask: Optional development type mask.
:return float: Result of operable area query.
"""
dtype_keys = list(self.dtypes.keys()) if not mask else self.unmask(mask)
return sum(self.dtypes[dtk].operable_area(acode, period, age) for dtk in dtype_keys)
[docs]
def overwrite_initial_areas(self, period):
"""
Overwrites the initial areas for all development types, for the specified period.
:param int period: Period for which to overwrite initial areas.
"""
for dt in list(self.dtypes.values()):
dt.overwrite_initial_areas(period)
[docs]
def initialize_areas(self, reset_areas: bool = True) -> None:
"""
Copies areas from period 0 to period 1.
:param bool reset_areas: Optionally calls ``self.reset_areas()`` if ``True``.
"""
if reset_areas:
self.reset_areas()
for dtk in self.dtypes:
self.dtypes[dtk].initialize_areas()
[docs]
def reset_areas(self, period: int | None = None) -> None:
"""
Reset areas for all development types.
:param int period: Optional period for which to reset areas.
Defaults to ``None`` (in which case resets all periods).
"""
for dtk in self.dtypes:
self.dtypes[dtk].reset_areas(period)
[docs]
def register_curve(self, curve: Any) -> Any:
"""
Add curve to global curve dictionary (uses result of ``curve.points()`` to construct key).
"""
key = tuple(curve.points()) # TO DO: use builtin common.hex_id() function to convert curves to hashed valued?
if key not in self.curves:
# new curve (lock and register)
curve.is_locked = True # points list must not change, else not valid key
self.curves[key] = curve
return self.curves[key]
[docs]
def reset_actions(self, period: int | None = None, acode: str | None = None, override_sticky: bool = False) -> None:
"""
Resets actions. By default resets, all actions in all periods (except for
sticky actions, unless overridden), unless ``period`` or ``acode`` specified.
:param int period: Optional period for which to reset actions.
:param str acode: Optional action code for which to reset actions.
:param bool override_sticky: Will override sticky actions if ``True``.
"""
periods = [period] if period else self.periods
acodes = [acode] if acode else list(self.actions.keys())
for p in periods:
if p not in self.applied_actions:
self.applied_actions[p] = {}
for a in acodes:
if a in self.actions and self.actions[a].is_sticky and not override_sticky:
continue
self.applied_actions[p][a] = {}
[docs]
def compile_product(self,
period,
expr,
acode=None,
dtype_keys=None,
age=None,
coeff=False,
verbose=False):
"""
Compiles products from applied actions in given period. Parses string expression,
which resolves to a single coefficient. Operated area can be filtered on action code,
development type key list, and age. Result is product of sum of filtered area and coefficient.
:param int period: Period for which to compile product query.
:param str expr: String expression to use when compiling product query. Must be a valid
expression string (see documentation for more details on what types of expressions
:py:mod:`ws3` can parse). Yield component names used in expressions will be automatically
resolved to the corresponding float yield values (by development type, with age lookup
corresponding to the age at which an action was applied in the current schedule).
:param str acode: Optional action code filter.
:param list dtype_keys: Optional list of development type keys on which to filter the query.
:param int age: Optional age filter.
:param bool coeff: Will force areas used to compile product query to 1 if ``True``,
else will use the actionned areas in the current schedule. Mostly only for use
when validating or debugging a model (and does not really make to use unless a single
development type key is specified in ``dtype_keys`` filter).
:param bool verbose: Verbosity flag.
:return float: Result of product query.
"""
aa = self.applied_actions
if acode is None:
acodes = list(self.actions.keys())
else:# elif type(acode) == str:
acodes = [acode] if not self.actions[acode].components else self.actions[acode].components
tokens = expr.split(' ')
result = 0.
for _acode in acodes:
if _acode not in list(aa[period].keys()):
continue # acode not in solution
_dtype_keys = list(aa[period][_acode].keys()) if dtype_keys is None else dtype_keys
for dtk in _dtype_keys:
if dtk not in list(aa[period][_acode].keys()):
continue
ages = list(aa[period][_acode][dtk].keys()) if age is None else [age]
for _age in ages:
aaa = aa[period][_acode][dtk][_age]
_tokens = []
for token in tokens:
if token in self.ynames: # found reference to ycomp
if token in aaa[1]: # token is yname in products (replace with value)
_tokens.append(str(aaa[1][token]))
else: # assume null value if ycomp exists but not stored in solution
_tokens.append('0.')
else:
_tokens.append(token)
_expr = ' '.join(_tokens)
area = aaa[0] if not coeff else 1.
try:
result += eval(_expr) * area
except ZeroDivisionError:
pass # let this one go...
except Exception:
print(("Unexpected error:", sys.exc_info()[0]))
print("evaluating expression '{}' for case:".format(' '.join(_tokens)), period, [' '.join(dtk)], _acode, _age)
raise
return result
[docs]
def operated_area(self, acode, period, dtype_key=None, age=None):
"""
Compiles operated area, given action code and period (and optionally list
of development type keys or age).
:param str acode: Action code for which to compile operated area.
:param int period: Period for which to compile operated area.
:param list dtype_keys: Optional list of development type keys to use as filter for operated area query.
:param int age: Optional age filter.
:return float: Result of operated area query.
"""
aa = self.applied_actions
acodes = [acode] if not self.actions[acode].components else self.actions[acode].components
result = 0.
for _acode in acodes:
if not aa[period][_acode]:
continue # acode not in solution
dtype_keys = list(aa[period][_acode].keys()) if dtype_key is None else [dtype_key]
for _dtype_key in dtype_keys:
ages = list(aa[period][_acode][_dtype_key].keys()) if age is None else [age]
for _age in ages:
result += aa[period][_acode][_dtype_key][_age][0]
return result
[docs]
def repair_actions(self, period: int, areaselector: Any = None, verbose: bool = False) -> None:
"""
Attempts to repair the action schedule for given period, using :py:class:`ws3.forest.AreaSelector`
object (defaults to class-default greed oldest-first area selector).
:param int period: Period for which to attempt to repair the action schedule.
:param areaselector: Area selector to use when attempting to repair action schedule.
:type areaselector: :py:class:`ws3.forest.AreaSelector`
:param book verbose: Verbosity flag.
"""
if areaselector is None: # use default (greedy) selector
areaselector = self.areaselector
aa = copy.copy(self.applied_actions[period])
self.reset_actions(period)
for acode in aa:
if not aa[acode]:
continue # null solution, move along...
old_area = 0.
new_area = 0.
# start by re-applying as much of the old solution as possible
for dtype_key in aa[acode]:
for age in aa[acode][dtype_key]:
aaa = aa[acode][dtype_key][age][0]
old_area += aaa
oa = self.dtypes[dtype_key].operable_area(acode, period, age)
if not oa:
continue
applied_area = min(aaa, oa)
#print ' applying old area', applied_area
new_area += applied_area
self.apply_action(dtype_key, acode, period, age, applied_area)
# try to make up for missing area...
target_area = old_area - new_area
if verbose:
print(f' patched {int(new_area)} of {int(old_area)} solution hectares, missing', target_area)
if areaselector is None: # use default area selector
areaselector = self.areaselector
areaselector.operate(period, acode, target_area)
[docs]
def commit_actions(self, period: int = 1, repair_future_actions: bool = False, verbose: bool = False) -> None:
"""
Commits applied actions (i.e., apply transitions and grow, default starting at period 1).
By default, will attempt to repair broken (infeasible) future actions, attempting to
replace infeasiblea operated area using default area selector.
:param int period: Period at which to start committing actions. Defaults to 1.
:param bool repair_future_actions: If ``True`` will attempt to repair future actions
(i.e., actions currently scheduled for periods after ``period``), else resets actions
in future periods. Defaults to ``False``.
"""
while period < self.horizon:
if verbose:
print('growing period', period)
self.grow(period, cascade=False) # type: ignore[no-untyped-call]
period += 1
if repair_future_actions:
if verbose:
print('repairing actions in period', period)
self.repair_actions(period)
else:
self.reset_actions(period)
[docs]
def resolve_replace(self, dtk: tuple[str, ...], expr: str) -> str:
"""
Enables the creation of new development types by replacing an existing
attribute code with a new value for a specific theme, instead of directly
coding the attribute change in transition.
:param tuple dtk: Source development type key.
:param str expr: Woodstock REPLACE expression parsed from an TRANSITIONS
model section using the :py:meth:`ws3.forest.ForestModel.import_transitions_section`
method.
:return: Updated development type key tuple.
:rtype: tuple
"""
# HACK ####################################################################
# Too lazy to implement all the use cases.
# This should work OK for BFEC models (TO DO: confirm).
tokens = re.split(r'\s+', expr)
i = int(tokens[0][3]) - 1
try:
return str(eval(expr.replace(tokens[0], dtk[i])))
except Exception:
print('dtk', ' '.join(dtk))
print('expr', expr)
raise
[docs]
def resolve_append(self, dtk: tuple[str, ...], expr: str) -> Any:
"""
Not been implemented yet.
"""
raise AssertionError() # brick wall (deal with this case later, as needed)
[docs]
def resolve_targetage(self, dtk: tuple[str, ...], tyield: Any, sage: int, tage: int | None, acode: str, verbose: bool = False) -> int:
"""
Determines the target age for a transition.
:param tuple dtk: Development type key tuple
:param str tyield: Target yield component name
:param int sage: Source age
:param int tage: Target age
:param str acode: Action code
:param bool verbose: Verbosity flag
:return int: Target age
"""
action = self.actions[acode]
if tyield is not None: # yield-based age definition
if verbose:
print('yield-based age definition', tyield, self.dt(dtk).ycomp(tyield[0]).lookup(tyield[1], roundx=True))
try:
targetage = self.dt(dtk).ycomp(tyield[0]).lookup(tyield[1], roundx=True)
except Exception:
print(' '.join(dtk), tyield[0], self.dt(dtk).ycomps())
raise AssertionError() from None
elif tage is not None: # target age override specifed in transition
if verbose:
print('_AGE override', tage)
targetage = tage
elif action.targetage is None: # use source age
if verbose:
print('source age', sage)
targetage = sage
else: # default: age reset to 0
if verbose:
print('default age reset to 0')
targetage = 0
return targetage # type: ignore[no-any-return]
[docs]
def apply_action(self,
dtype_key: tuple[str, ...],
acode: str,
period: int,
age: int,
area: float,
override_operability: bool = False,
fuzzy_age: bool = True,
recourse_enabled: bool = True,
areaselector: Any = None,
compile_t_ycomps: bool = False,
compile_c_ycomps: bool = False,
verbose: bool = False) -> tuple[int, float, list[tuple[tuple[str, ...], Any, int]]]:
"""
Applies action, given action code, development type, period, age, area.
Can optionally override operability limits, optionally use fuzzy age (i.e., attempt
to apply action to proximal age class if specified age is not operable), optionally use
default AreaSelector to patch missing area (if recourse enabled). Applying an action is
a rather complex process, involving testing for operability (JIT-compiling operability
expression as required), checking that valid transitions are defined, checking that area
is available (possibly using fuzzy age and area selector functions to find missing area),
generate list of target development types (from source development type and transition
expressions [which may need to be JIT-compiled]), creating new development types
(as needed), doing the area accounting correctly (without creating or destroying any area)
and compiling the products from the action (which gets a bit complicated in the case of
partial cuts...).
:param tuple dtype_key: The key identifying the development type.
:param str acode: The action code to apply.
:param int period: The period in which to apply the action.
:param int age: The age at which to apply the action.
:param float area: The area to apply the action on.
:param bool override_operability: If True, overrides operability limits. Default is False.
:param bool fuzzy_age: If True, attempts to apply action to proximal age class if specified age is not operable.
:param bool recourse_enabled: If True, uses default AreaSelector to patch missing area. Default is True.
:param bool areaselector: The AreaSelector object to use for patching missing area. Default is None.
:param bool compile_t_ycomps: If True, compiles time-indexed yield components. Default is False.
:param bool compile_c_ycomps: If True, compiles complex yield components. Default is False.
:param bool verbose: If True, prints additional information for debugging purposes. Default is False.
:return tuple: Returns ``(errorcode, missing_area, target_dt)`` triplet, where ``errorcode`` is an error code,
``missing_area`` is the missing area, and ``target_dt`` is a list of ``(dtk, tprop, targetage)``
triplets (one triplet per target development type).
**Error codes:**
1. Invalid area argument
2. Requested action not defined for development type
3. Requested action defined, but never operable
4. Action not operable
5. Transitions not defined for action
"""
if area <= 0. and not override_operability:
return 1, None, None # type: ignore[return-value]
if verbose > 1:
print('applying action', [' '.join(dtype_key)], acode, period, age, area)
dt = self.dtypes[dtype_key]
if acode not in dt.oper_expr:
print('requested action not defined for development type...')
print(' ', [' '.join(dtype_key)], acode, period, age, area)
return 2, None, None # type: ignore[return-value]
if acode not in dt.operability: # action not compiled yet...
if dt.compile_action(acode) == -1:
print('requested action is defined, but never not operable...')
print(' ', [' '.join(dtype_key)], acode, period, age, area)
return 3, None, None # type: ignore[return-value]
if not dt.is_operable(acode, period, age) and not override_operability:
print('not operable')
print(' '.join(dt.key), acode, period, age)
print(dt.operability[acode][period])
return 4, None, None # type: ignore[return-value]
if not any((acode, __age) in dt.transitions for __age in (age, -1)): # sanity check...
print('transitions not defined...')
print(' ', [' '.join(dtype_key)], acode, period, age, area)
print(dt.oper_expr)
print(dt.operability)
return 5, None, None # type: ignore[return-value]
if dt.area(period, age) - area < self.area_epsilon:
# tweak area if slightly over or under, so we don't get any accounting drift...
area = dt.area(period, age)
missing_area = 0.
if dt.area(period, age) < area:
# insufficient area in dt to operate (infeasible)
# apply action to operable area, then look for missing area in adjacent ageclasses
if dt.area(period, age) > 0: # operate available area before applying recourse
print('insufficient area in dt to operate (infeasible)', dtype_key, period, age)
self.apply_action(dtype_key, acode, period, age, dt.area(period, age),
False, False, False, None, True)
missing_area = area - dt.area(period, age)
if fuzzy_age and missing_area:
for age_delta in [+1, -1, +2, -2]:
_age = age + age_delta
if dt.area(period, _age) > 0 and any((acode, __age) in dt.transitions for __age in (_age, -1)):
_area = min(missing_area, dt.area(period, _age))
self.apply_action(dtype_key, acode, period, _age, _area,
False, False, False, None, True)
missing_area -= _area
if missing_area < self.area_epsilon:
missing_area = 0.
break
if recourse_enabled and missing_area:
areaselector = self.areaselector if areaselector is None else areaselector
missing_area = areaselector.operate(period, acode, missing_area)
if missing_area < self.area_epsilon:
missing_area = 0.
action = self.actions[acode]
###########################################################################
dt.area(period, age, -area)
target_dt = []
__age = age if (acode, age) in dt.transitions else -1
for target in dt.transitions[acode, __age]:
tmask, tprop, tyield, tage, tlock, treplace, tappend = target # unpack tuple
dtk = list(dtype_key) # start with source key
###########################################################################
# DO TO: Confirm correct order for evaluating mask, _APPEND and _REPLACE...
dtk = [t if tmask[i] == '?' else tmask[i] for i, t in enumerate(dtk)]
if treplace:
dtk[treplace[0]] = self.resolve_replace(dtk, treplace[1]) # type: ignore[arg-type]
if tappend:
dtk[tappend[0]] = self.resolve_append(dtk, tappend[1]) # type: ignore[arg-type, assignment]
dtk = tuple(dtk) # type: ignore[assignment]
###########################################################################
# import pdb; pdb.set_trace()
_dt = self.create_dtype_fromkey(dtk) if dtk not in self.dtypes else self.dtypes[dtk] # type: ignore[arg-type]
targetage = self.resolve_targetage(dtk, tyield, age, tage, acode) # type: ignore[arg-type]
_dt.area(period, targetage, area*tprop)
target_dt.append([dtk, tprop, targetage])
aa = self.applied_actions[period][acode]
if dtype_key not in aa:
aa[dtype_key] = {}
if age not in aa[dtype_key]:
aa[dtype_key][age] = [0., {}]
aa[dtype_key][age][0] += area
for yname in dt.ycomps():
ycomp = dt.ycomp(yname)
if ycomp.type == 't' and not compile_t_ycomps:
continue # skip time-indexed ycomps
if ycomp.type == 'c' and not compile_c_ycomps:
continue # skip complex ycomps
if yname in action.partial:
value = 0.
for dtk, tprop, targetage in target_dt:
_dt = self.dtypes[dtk]
_value = 0.
if yname in dt.ycomps():
if yname in _dt.ycomps():
_value = (dt.ycomp(yname)[age] - _dt.ycomp(yname)[targetage])
else:
_value = dt.ycomp(yname)[age]
if _value > 0.:
value += _value * tprop
else:
if verbose:
if _value < 0:
print('negative partial value', acode, yname, tprop, _value)
print(' ', ''.join(dtype_key), age)
print(' ', ''.join(dtk), targetage)
print()
else: # not partial action
value = dt.ycomp(yname)[age]
if value != 0.:
aa[dtype_key][age][1][yname] = value
return 0, missing_area, target_dt # type: ignore[return-value]
[docs]
def create_dtype_fromkey(self, key: tuple[str, ...]) -> Any:
"""
Creates a new development type, given a development type key (checks for existing,
auto-assigns yield compompontents, auto-assign actions and transitions, checks for
operability (filed under inoperable if applicable).
:param tuple key: Development type key
:return: New development type
:rtype: :py:class:`ws3.forest.DevelopmentType`
"""
assert key not in self.dtypes # should not be creating new dtypes from existing key
dt = DevelopmentType(key, self)
self.dtypes[key] = dt
for mask, t, ycomps in self.yields:
if self.match_mask(mask, key):
for yname, ycomp in ycomps:
dt.add_ycomp(t, yname, ycomp)
# assign actions and transitions
for acode in self.oper_expr:
for mask in self.oper_expr[acode]:
if self.match_mask(mask, key):
dt.oper_expr[acode].append(self.oper_expr[acode][mask])
for mask in self.transitions[acode]:
if self.match_mask(mask, key):
for scond in self.transitions[acode][mask]:
for x in self.resolve_condition(scond, key):
dt.transitions[acode, x] = self.transitions[acode][mask][scond]
if not dt.transitions:
self.inoperable_dtypes.append(key)
return dt
def _resolve_outputs_buffer(self, s: str, for_flag: Any = None) -> Any:
"""
Resolve outputs buffer.
:param str s: String to resolve
:param int for_flag: Flag indicating for loop iteration (this is a recursive function)
"""
group = 'no_group' # outputs declared at top of file assigned to 'no_group'
self.output_groups[group] = set()
ocode = ''
buffering_for = False
# Loop-carried state: bound in earlier iterations before first use, but
# initialized here so the bindings are explicit and statically provable.
for_var, for_lo, for_hi = None, 0, 0
for_buffer: list[Any] = [] # type: ignore[var-annotated]
expression, description, theme_index = '', '', None
s = re.sub(r'\{.*?\}', '', s, flags=re.M|re.S) # remove curly-bracket comments
for line_ in re.split(r'[\r\n]+', s, flags=re.M|re.S):
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
matches = re.findall(r'#[A-Za-z0-9_]*', line_)
for m in matches: # replace CONSTANTS variables with value
try:
line_ = line_.replace(m, str(self.constants[m[1:].lower()]))
except Exception:
import sys
print(sys.exc_info()[0])
print(line_)
print(matches, m)
raise AssertionError() from None
if buffering_for:
if line_.strip().startswith('ENDFOR'):
for i in range(for_lo, for_hi+1):
ss = '\n'.join(for_buffer).replace(for_var, str(i)) # type: ignore[arg-type]
self._resolve_outputs_buffer(ss, for_flag=i)
buffering_for = False
continue
else:
for_buffer.append(line_)
continue
line_ = re.sub(r'\s+', ' ', line_) # separate tokens by single space
line_ = line_.strip().partition(';')[0].strip()
line_ = line_.replace(' (', '(') # remove space to left of left parentheses
##################################################
# HACK ###########################################
# substitute ugly symbols have in ocodes...
line_ = line_.replace(r'%', 'p')
line_ = line_.replace(r'$', 's')
##################################################
tokens = line_.lower().split(' ')
if line_.startswith('*GROUP'):
keyword = 'group'
group = tokens[1].lower()
self.output_groups[group] = set()
# Parse output codes from the rest of the line (comma-separated)
if len(tokens) > 2:
codes_str = ' '.join(tokens[2:]) # rejoin in case of spaces
for code in codes_str.split(','):
code = code.strip().lower()
if code:
self.output_groups[group].add(code)
elif line_.startswith('FOR'):
# pattern matching may not be very robust, but works for now with:
# 'FOR XX := 1 to 99'
# TO DO: implement DOWNTO, etc.
for_var = _search(r'(?<=FOR\s).+(?=:=)', line_, 'FOR loop variable').group(0).strip()
for_lo = int(_search(r'(?<=:=).+(?=to)', line_, 'FOR loop lower bound').group(0))
for_hi = int(_search(r'(?<=to).+', line_, 'FOR loop upper bound').group(0))
for_buffer = []
buffering_for = True
continue
if line_.startswith('*OUTPUT') or line_.startswith('*LEVEL'):
keyword = 'output' if line_.startswith('*OUTPUT') else 'level'
if ocode: # flush data collected from previous lines
self.outputs[ocode] = Output(parent=self, # type: ignore[no-untyped-call]
code=ocode,
expression=expression,
description=description,
theme_index=theme_index)
tt = tokens[1].split('(')
ocode = tt[0]
theme_index = int(tt[1][:-1]) - 1 if len(tt) > 1 else None
description = ' '.join(tokens[2:])
expression = ''
self.output_groups[group].add(ocode)
if keyword == 'level':
self.outputs[ocode] = Output(parent=self, # type: ignore[no-untyped-call]
code=ocode,
expression=expression,
description=description,
theme_index=theme_index,
is_level=True)
ocode = ''
elif line_.startswith('*SOURCE'):
keyword = 'source'
expression += line_[8:]
elif keyword == 'source': # continuation line of SOURCE expression
expression += ' '
expression += line_
# Flush any remaining output after the loop ends.
if ocode:
self.outputs[ocode] = Output(parent=self, # type: ignore[no-untyped-call]
code=ocode,
expression=expression,
description=description,
theme_index=theme_index)
[docs]
def import_outputs_section(self, filename_suffix='out'):
"""
Imports OUTPUTS section from a Woodstock-formatted model input dataset.
Model path and file base name assumed from ``self.model_path`` and ``self.model_path``.
:param str filename_suffix: Filename suffix in which to look for OUTPUTS section data.
"""
with open(f'{self.model_path}/{self.model_name}.{filename_suffix}') as f:
s = f.read()
self._resolve_outputs_buffer(s)
[docs]
def add_theme(self, name, basecodes=None, aggs=None, description=''):
"""
Adds a theme to the model.
:param str name: The name of theme.
:param list basecodes: List of base codes for the theme.
:param dict aggs: Dictionary containing aggregate values for the theme.
:param str description: Description of the theme.
"""
if aggs is None:
aggs = {}
if basecodes is None:
basecodes = []
self._themes.append({})
self._themes[-1]['__name__'] = name
self._themes[-1]['__description__'] = description
if basecodes:
self._theme_basecodes.append([])
for c in basecodes:
self._themes[-1][c] = c
self._theme_basecodes[-1].append(c)
for c in aggs:
self._themes[-1][c] = aggs[c]
[docs]
def import_landscape_section(self, filename_suffix='lan', ti_offset=0):
"""
Imports LANDSCAPE section from a Woodstock-formatted model input dataset.
Model path and file base name assumed from ``self.model_path`` and ``self.model_path``.
:param str filename_suffix: Filename suffix in which to look for LANDSCAPE section data
:param str ti_offset: Theme index offset. Will bump theme index by this value. Defaults to 0.
"""
with open(f'{self.model_path}/{self.model_name}.{filename_suffix}') as f:
data = f.read()
_data = _search(r'\*THEME.*', data, 'landscape section (no *THEME declaration found)',
re.M | re.S).group(0) # strip leading junk
# Capture the text trailing each *THEME declaration rather than discarding it.
# That text is the modeller's own description of the theme (e.g. "Analysis Unit
# (AU)"), and it is the only thing in the dataset that says what a theme position
# means -- theme order and meaning are entirely user-defined in this format.
_chunks = re.split(r'\*THEME(.*)\n', _data)[1:] # [desc, body, desc, body, ...]
t_data = list(zip(_chunks[::2], _chunks[1::2], strict=False))
for ti, (t_description, t) in enumerate(t_data, start=ti_offset):
self._themes.append({})
self._themes[-1]['__name__'] = f'theme{ti}'
self._themes[-1]['__description__'] = t_description.strip().lstrip(';').strip()
self._theme_basecodes.append([])
defining_aggregates = False
for line_ in [line_ for line_ in t.split('\n') if not re.match(r'^\s*(;|{|$)', line_)]:
if re.match(r'^\s*\*AGGREGATE', line_): # aggregate theme attribute code
tac = re.split(r'\s+', line_.strip())[1].lower()
self._themes[ti][tac] = []
defining_aggregates = True
continue
if not defining_aggregates: # line defines basic theme attribute code
tac = _search(r'\S+', line_.strip(), 'theme attribute code').group(0).lower()
self._themes[ti][tac] = tac
self._theme_basecodes[ti].append(tac)
else: # line defines aggregate values (parse out multiple values before comment)
_tacs = [_tac.lower() for _tac in re.split(r'\s+', line_.strip().partition(';')[0].strip())]
self._themes[ti][tac].extend(_tacs)
[docs]
def theme_basecodes(self, theme_index: int) -> list[str]:
"""
Return list of base codes, given theme index.
:param int theme_index: Theme index for which to return basecodes.
:return list: List of theme basecodes.
"""
return self._theme_basecodes[theme_index] # type: ignore[no-any-return]
[docs]
def import_areas_section(self,
model_path=None,
model_name=None,
filename_suffix='are',
import_empty=False,
convert_periods_to_years=None):
"""
Imports AREAS section from a Woodstock-formatted model input dataset.
Each line in the section represents an area for a unique combination of
development type key and age class. Empty areas (with values less than
``area_epsilon``) will be skipped if ``import_empty`` is ``False``.
:param str filename_suffix: Suffix for AREAS section file name.
:param bool import_empty: Whether or not to import empty areas (with values less than ``area_epsilon``).
:return int: 0 if succcess, 1 otherwise.
"""
n = self.nthemes()
model_path = self.model_path if not model_path else model_path
model_name = self.model_name if not model_name else model_name
multiplier = self._resolve_period_multiplier(convert_periods_to_years)
with open(f'{model_path}/{model_name}.{filename_suffix}') as f:
for line_ in f:
try:
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
line_ = line_.lower().strip().partition(';')[0] # strip leading whitespace and trailing comments
t = re.split(r'\s+', line_)
key = tuple(_t for _t in t[1:n+1])
age = int(t[n+1]) * multiplier
area = float(t[n+2].replace(',', ''))
if area < self.area_epsilon and not import_empty:
continue
if key not in self.dtypes:
self.dtypes[key] = DevelopmentType(key, self)
self.dtypes[key].area(0, age, area)
except Exception:
print(f'Failed AREAS import on line: \n{line_}')
return 1
return 0
def _expand_theme(self, t: Any, c: str, verbose: int = 0) -> list[str]:
"""
Depth-first search recursive aggregate theme code expansion.
:param t: Target theme.
:param c: Theme code to expand.
:param verbose: Verbosity level (0=silent, 1=normal, 2=verbose).
:return list: List of expanded theme code values.
"""
if verbose > 1:
print('ws3.forest.ForestModel._expand_theme', t, c)
print(c)
return [c] if t[c] == c else list(chain.from_iterable(self._expand_theme(t, c) for c in t[c]))
[docs]
def match_mask(self, mask: tuple[str, ...], key: tuple[str, ...]) -> bool:
"""
Checks if a development type key matches a development type mask.
:param tuple mask: Development type mask
:param tuple key: Development type key
:return bool: Returns ``True`` if key matches mask, ``False`` otherwise.
"""
for ti, tac in enumerate(mask):
if tac == '?':
continue # wildcard matches all keys
tacs = self._expand_theme(self._themes[ti], tac)
if key[ti] not in tacs:
return False # reject key
return True # key matches
[docs]
def unmask(self, mask: Any, verbose: int = 0) -> list[tuple[str, ...]]:
"""
Iteratively filter list of development type keys using mask values.
Accepts Woodstock-style string masks to facilitate cut-and-paste testing.
:param tuple or str mask: Development type mask (tuple or Woodstock-style string format)
:param int verbose: Verbosity level (passed to ``self._expand_theme``).
:return list: List of development type keys that match the mask.
"""
if isinstance(mask, str): # Woodstock-style string mask format
mask = tuple(re.sub(r'\s+', ' ', mask).lower().split(' '))
assert len(mask) == self.nthemes() # must be bad mask if wrong theme count
else:
try:
assert isinstance(mask, tuple) and len(mask) == self.nthemes()
except Exception:
print(len(mask), type(mask), mask)
raise AssertionError() from None
dtype_keys = copy.copy(list(self.dtypes.keys())) # filter this
for ti, tac in enumerate(mask):
if tac == '?':
continue # wildcard matches all
tacs = self._expand_theme(self._themes[ti], tac, verbose=verbose) if tac in self._themes[ti] else []
dtype_keys = [dtk for dtk in dtype_keys if dtk[ti] in tacs] # exclude bad matches
return dtype_keys
[docs]
def import_constants_section(self, filename_suffix='con'):
"""
Imports CONSTANTS section from a Woodstock-formatted input dataset.
Each line in the section represents a constant and its value.
Constants are stored in a dictionary where the keys are the constant names and the values are their respective values.
:param str filename_suffix: Suffix for CONSTANTS section file name.
"""
with open(f'{self.model_path}/{self.model_name}.{filename_suffix}') as f:
for _lnum, line_ in enumerate(f):
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
line_ = line_.strip().partition(';')[0].strip() # strip leading whitespace, trailing comments
t = re.split(r'\s+', line_)
self.constants[t[0].lower()] = float(t[1])
[docs]
def import_yields_section(self,
filename_suffix='yld',
mask_func=None,
verbose=False,
convert_periods_to_years=None):
"""
Imports YIELDS section from a Forest model.
:param str filename_suffix: Suffix for CONSTANTS section file name.
:param function mask_func: Custom mask function (deprecate?)
:param bool verbose: Verbosity flag
"""
def flush_ycomps(t, m, n, c):
if t == 'a': # age-based ycomps
def _c(y):
return self.register_curve(core.Curve(y,
points=c[y],
type='a',
period_length=self.period_length))
ycomps = [(y, _c(y)) for y in n] # type: ignore[no-untyped-call]
elif t == 't': # time-based ycomps (skimp on x range)
def _c(y):
return self.register_curve(core.Curve(y,
points=c[y],
type='t',
xmax=self.horizon,
period_length=self.period_length))
ycomps = [(y, _c(y)) for y in n] # type: ignore[no-untyped-call]
else: # complex ycomps
ycomps = [(y, c[y]) for y in n]
self.yields.append((m, t, ycomps)) # stash for creating new dtypes at runtime...
self.ynames.update(n)
for k in self.unmask(m):
for yname, ycomp in ycomps:
self.dtypes[k].add_ycomp(t, yname, ycomp)
self._resolve_period_multiplier(convert_periods_to_years)
age_multiplier = self._period_to_years_factor or 1
period_step = self._period_to_years_factor or self.period_length
ytype = ''
mask = ('?',) * self.nthemes()
ynames: list[str] = [] # type: ignore[var-annotated]
data: dict[str, list[Any]] = {} # type: ignore[assignment]
with open(f'{self.model_path}/{self.model_name}.{filename_suffix}') as f:
for lnum, line_ in enumerate(f):
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
line_ = line_.strip().partition(';')[0].strip() # strip leading whitespace and trailing comments
t = re.split(r'\s+', line_)
if t[0].startswith('*Y'): # new yield definition
newyield = True
flush_ycomps(ytype, mask, ynames, data) # type: ignore[no-untyped-call]
ytype = self._ytypes[t[0]]
mask = tuple(_t.lower() for _t in t[1:])
mask = mask_func(mask) if mask_func else mask
if verbose:
print(lnum, ' '.join(mask))
continue
if newyield:
if t[0] == '_AGE':
is_tabular = True
ynames = [_t.lower() for _t in t[1:]]
data = {yname:[] for yname in ynames}
newyield = False
continue
else:
is_tabular = False
ynames = []
data = {}
newyield = False
else:
if t[0] == '_AGE': # same yield block, new table
flush_ycomps(ytype, mask, ynames, data) # type: ignore[no-untyped-call]
is_tabular = True
ynames = [_t.lower() for _t in t[1:]]
data = {yname:[] for yname in ynames}
newyield = False
continue
if is_tabular:
try:
x = int(t[0]) * age_multiplier
except Exception:
print(lnum, line_)
for i, yname in enumerate(ynames):
data[yname].append((x, float(t[i+1])))
else:
if ytype in 'at': # standard or time-based yield (extract xy values)
if not common.is_num(t[0]): # first line of row-based yield component
yname = t[0].lower()
ynames.append(yname)
start = int(t[1]) * period_step
data[yname] = [(start + (i * period_step), float(t[i+2]))
for i in range(len(t) - 2)]
else: # continuation of row-based yield compontent
x_last = data[yname][-1][0]
if len(data[yname]) >= 2:
step = data[yname][1][0] - data[yname][0][0]
else:
step = period_step
if not step:
step = period_step
data[yname].extend([(x_last + step * (i + 1), float(t[i]))
for i in range(len(t))])
else:
yname = t[0].lower()
ynames.append(yname)
data[yname] = ' '.join(t[1:]) # type: ignore[assignment]
flush_ycomps(ytype, mask, ynames, data) # type: ignore[no-untyped-call]
[docs]
def import_actions_section(self,
filename_suffix='act',
mask_func=None,
nthemes=None,
convert_periods_to_years=None):
"""
Imports ACTIONS section from a Woodstock-formatted model input dataset.
:param str filename_suffix: Suffix for CONSTANTS section file name.
:param function mask_func: Custom mask function (deprecate?)
:param int nthemes: Number of themes
"""
nthemes = nthemes if nthemes else self.nthemes()
multiplier = self._resolve_period_multiplier(convert_periods_to_years)
def scale_expr(expr):
if multiplier == 1:
return expr
pattern = re.compile(r'(_age\s*(?:>=|<=|=|<|>)\s*)(\d+)')
def repl(match):
return f"{match.group(1)}{int(int(match.group(2)) * multiplier)}"
return pattern.sub(repl, expr)
partials: dict[str, Any] = {} # type: ignore[var-annotated]
keyword = ''
with open(f'{self.model_path}/{self.model_name}.{filename_suffix}') as f:
s = f.read().lower()
s = re.sub(r'\{.*?\}', '', s, flags=re.M|re.S) # remove curly-bracket comments
for line_ in re.split(r'[\r\n]+', s, flags=re.M|re.S):
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
line_ = line_.strip().partition(';')[0].strip() # strip leading whitespace and trailing comments
line_ = re.sub(r'r\s+', ' ', line_) # separate tokens by single space
tokens = line_.split(' ')
if line_.startswith('*action'):
keyword = 'action'
acode = tokens[1]
targetage = 0 if tokens[2] == 'y' else None
descr = ' '.join(tokens[3:])
lockexempt = '_lockexempt' in tokens
self.actions[acode] = Action(acode, targetage, descr, lockexempt)
self.oper_expr[acode] = {}
elif line_.startswith('*operable'):
keyword = 'operable'
acode = tokens[1]
elif line_.startswith('*aggregate'):
keyword = 'aggregate'
acode = tokens[1]
self.actions[acode] = Action(acode)
elif line_.startswith('*partial'):
keyword = 'partial'
acode = tokens[1]
partials[acode] = []
else: # continuation of OPERABLE, AGGREGATE, or PARTIAL block
if keyword == 'operable':
mask = tuple(tokens[:nthemes])
mask = mask_func(mask) if mask_func else mask
expression = ' '.join(tokens[nthemes:])
self.oper_expr[acode][mask] = scale_expr(expression) # type: ignore[no-untyped-call]
elif keyword == 'aggregate':
self.actions[acode].components.extend(tokens)
elif keyword == 'partial':
self.actions[acode].partial.extend(tokens)
for acode, a in list(self.actions.items()):
if a.components:
continue # aggregate action, skip
for mask, expression in list(self.oper_expr[acode].items()):
for k in self.unmask(mask):
self.dtypes[k].oper_expr[acode].append(expression)
[docs]
def resolve_treplace(self, dt: Any, treplace: str) -> str:
"""
Resolves a theme replace expression, in the context of defining a new development type
when implementing a transition (following application of an action to a source development type).
:param :py:class:`ws3.forest.DevelopmentType` dt: Source development type
:param str treplace: Theme replace expression to apply
:return str: New theme value string
"""
if '_TH' in treplace: # assume incrementing integer theme value
i = int(_search(r'(?<=_TH)\w+', treplace, '_TH theme index in _REPLACE expression').group(0))
return eval(re.sub(f'_TH{i}', str(dt.key[i-1]), treplace)) # type: ignore[no-any-return]
else:
raise AssertionError() # many other possible arguments (see Woodstock documentation for details)
[docs]
def resolve_tappend(self, dt, tappend):
"""
Resolves a theme append expression, in the context of defining a new development type
when implementing a transition (following application of an action to a source development type).
.. warning::
Not implemented yet.
"""
raise AssertionError() # brick wall (not implemented yet)
[docs]
def resolve_tmask(self, dt, tmask, treplace, tappend):
"""
Returns new developement type key (tuple of values, one per theme),
given developement type, theme mask, theme replace expression, and theme append expression.
:param :py:class:`ws3.forest.DevelopmentType` dt: Source development type
:param tuple tmask: Theme mask to apply to source development type key
:param str treplace: Theme replace expresion to apply source development type key
:param str tappend: Theme append expression to apply to soruce development type key
:return tuple: New development type key
"""
key = list(dt.key)
if treplace:
key[treplace[0]] = self.resolve_treplace(dt, treplace[1])
if tappend:
key[tappend[0]] = self.resolve_tappend(dt, tappend[1]) # type: ignore[no-untyped-call]
for i, val in enumerate(tmask):
if val == '?':
continue # wildcard (skip it)
key[i] = val
return tuple(key)
[docs]
def resolve_condition(self, condition: Any, dtype_key: tuple[str, ...] | None = None) -> list[int]: # type: ignore[return]
"""
Expands ``@AGE`` or ``@YLD`` conditions to list of age values.
``@AGE`` condition specifies lower- and upper-bound ages in a range,
so just expands to age values that fall within that range.
``@YLD`` condition specifies lower- and upper-bound yield values for a
given yield component name, so needs to do a reverse age lookup on the
specified yield component (so ``dtype_key`` must be specified or in this
case or the method will crash).
:param str condition: Condition expression string
:param tuple dtype_key: Development type key
:return list: List of age values
"""
if not condition:
return [-1]
elif condition.startswith('@AGE'):
lo, hi = [int(a) for a in condition[5:-1].split('..')]
multiplier = self._period_to_years_factor or 1
if multiplier != 1:
lo *= multiplier
hi *= multiplier
step = multiplier
else:
step = 1
return list(range(lo, hi + step, step))
elif condition.startswith('@YLD'):
args = re.split(r'\s?,\s?', condition[5:-1])
yname = args[0].lower() # type: ignore[assignment]
lo, hi = [float(y) for y in args[1].split('..')] # type: ignore[assignment]
if not dtype_key:
raise AssertionError() # to do: add better error handling
dt = self.dtypes[dtype_key]
lo_age, hi_age = dt.ycomp(yname).range(lo, hi, as_bounds=True) # type: ignore[assignment]
return list(range(int(lo_age), int(hi_age)+1))
[docs]
def import_transitions_section(self,
filename_suffix='trn',
mask_func=None,
nthemes=None,
convert_periods_to_years=None):
"""
Imports TRANSITIONS section from a Woodstock-formatted model input dataset.
:param str filename_suffix: Suffix for CONSTANTS section file name.
:param function mask_func: Custom mask function
:param int nthemes: Number of themes
"""
nthemes = nthemes if nthemes else self.nthemes()
self._resolve_period_multiplier(convert_periods_to_years)
age_multiplier = self._period_to_years_factor or 1
def flush_transitions(acode, sources):
if not acode:
return # nothing to flush on first loop
self.transitions[acode] = {}
for smask, scond in sources:
# store transition data for future dtypes creation
if smask not in self.transitions[acode]:
self.transitions[acode][smask] = {}
self.transitions[acode][smask][scond] = sources[smask, scond]
# assign transitions to existing dtypes
for k in self.unmask(smask):
dt = self.dtypes[k]
for x in self.resolve_condition(scond, k): # store targets
dt.transitions[acode, x] = sources[smask, scond]
acode = None
sources: dict[Any, Any] = {} # type: ignore[var-annotated]
with open(f'{self.model_path}/{self.model_name}.{filename_suffix}') as f:
s = f.read()
s = re.sub(r'\{.*?\}', '', s, flags=re.M|re.S) # remove curly-bracket comments
for line_ in re.split(r'[\r\n]+', s, flags=re.M|re.S):
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
line_ = line_.strip().partition(';')[0].strip() # strip leading whitespace, trailing comments
tokens = re.split(r'\s+', line_)
if line_.startswith('*CASE'):
if acode:
flush_transitions(acode, sources) # type: ignore[no-untyped-call]
acode = tokens[1].lower()
sources = {}
elif line_.startswith('*SOURCE'):
smask = tuple(t.lower() for t in tokens[1:nthemes+1])
smask = mask_func(smask) if mask_func else smask
match = re.search(r'@.+\)', line_)
scond = match.group(0) if match else ''
sources[(smask, scond)] = []
elif line_.startswith('*TARGET'):
tmask = tuple(t.lower() for t in tokens[1:nthemes+1])
tmask = mask_func(tmask) if mask_func else tmask
tprop = float(tokens[nthemes+1]) * 0.01
tyield = None
if len(tokens) > nthemes+2 and tokens[nthemes+2].lower() in self.ynames:
tyield = (tokens[nthemes+2].lower(), float(tokens[nthemes+3]))
try: # _AGE keyword
tage = int(tokens[tokens.index('_AGE')+1]) * age_multiplier
except Exception:
tage = None
try: # _LOCK keyword
tlock = int(tokens[tokens.index('_LOCK')+1]) * age_multiplier
except Exception:
tlock = None
try: # _REPLACE keyword (TO DO: implement other cases)
args = re.split(r'\s?,\s?', _search(r'(?<=_REPLACE\().*(?=\))', line_, '_REPLACE arguments').group(0))
theme_index = int(args[0][3]) - 1
treplace = theme_index, args[1]
except Exception:
treplace = None
try: # _APPEND keyword (TO DO: implement other cases)
args = re.split(r'\s?,\s?', _search(r'(?<=_APPEND\().*(?=\))', line_, '_APPEND arguments').group(0))
theme_index = int(args[0][3]) - 1
tappend = theme_index, args[1]
except Exception:
tappend = None
sources[(smask, scond)].append((tmask, tprop, tyield, tage, tlock, treplace, tappend))
flush_transitions(acode, sources) # type: ignore[no-untyped-call]
[docs]
def import_optimize_section(self, filename_suffix='opt'):
"""
Imports OPTIMIZE section from a Woodstock-formatted model input dataset.
.. warning::
Not implemented yet.
"""
pass
[docs]
def import_graphics_section(self, filename_suffix='gra'):
"""
Imports GRAPHICS section from a Woodstock-formatted model input dataset.
.. warning::
Not implemented yet.
"""
pass
[docs]
def import_lifespan_section(self, filename_suffix='lif'):
"""
Imports LIFESPAN section from a Woodstock-formatted model input dataset.
.. warning::
Not implemented yet.
"""
pass
[docs]
def import_schedule_section(self,
filename_suffix='seq',
replace_commas=True,
filename_prefix=None,
convert_periods_to_years=None):
"""
Imports SCHEDULE section from a Woodstock-formatted model input dataset.
:param str filename_suffix: Suffix for SCHEDULE section file name.
:param bool replace_commas: Remove commas from area value string tokens
:param str filename_prefix: Prefix for SCEDULE section file name.
"""
filename_prefix = self.model_name if filename_prefix is None else filename_prefix
age_multiplier = self._resolve_period_multiplier(convert_periods_to_years)
schedule = []
n = self.nthemes()
with open(f'{self.model_path}/{filename_prefix}.{filename_suffix}') as f:
for _lnum, line_ in enumerate(f):
if re.match(r'^\s*(;|$)', line_):
continue # skip comments and blank lines
line_ = line_.lower().strip().partition(';')[0].strip() # strip leading whitespace and trailing comments
t = re.split(r'\s+', line_)
if len(t) not in (n + 4, n + 5):
break
dtype_key = tuple(t[:n])
age = int(t[n]) * age_multiplier
area = float(t[n+1].replace(',', '')) if replace_commas else float(t[n+1])
acode = t[n+2]
period = int(t[n+3])
etype = t[n+4] if len(t) == n + 5 else ''
schedule.append((dtype_key, age, area, acode, period, etype))
if area <= 0:
print('area <= 0', line_)
return schedule
[docs]
def compile_schedule(self, problem: Any = None) -> list[tuple[Any, ...]]:
"""
Compiles an action schedule.
If a :py:class;`ws3.opt.Problem` instance is specified compiles the schedule from the
optimal solution of the problem instance, otherwise compiles the schedule from the current
solution (i.e., ``self.applied_actions``).
:param :py:class:`ws3.opt.Problem` problem: Optimization problem from which to extract an action schedule
:return list: Action schedule as list of ``(dtk, age, area, acode, period, etype)`` tuples.
"""
if problem is not None:
return self._compile_schedule_from_problem(problem) # type: ignore[no-any-return]
else: # use data in self.applied_actions
return self._compile_schedule_from_actions()
def _compile_schedule_from_actions(self) -> list[tuple[Any, ...]]:
"""
Compiles an action schedule from the current solution.
:return list: Action schedule as list of ``(dtype_key, age, area, acode, period, etype)`` tuples.
"""
result = []
for period in self.periods:
aa = self.applied_actions[period]
for acode in aa.keys():
for dtk in aa[acode].keys():
etype = '_existing' if self.dt(dtk).area(0) else '_future'
for age in aa[acode][dtk].keys():
area = aa[acode][dtk][age][0]
result.append((dtk, age, area, acode, period, etype))
return result
[docs]
def apply_schedule(self, schedule: Any, max_period: int | None = None, verbose: bool = False,
fail_on_missingarea: bool = False, force_integral_area: bool = False,
override_operability: bool = False, fuzzy_age: bool = True,
recourse_enabled: bool = True, areaselector: Any = None,
compile_t_ycomps: bool = False, compile_c_ycomps: bool = False,
rounding_bias: float = 0.15, scale_area: Any = None, reset: bool = True,
crash_on_action_error: bool = False) -> None:
"""
Assumes schedule in format returned by import_schedule_section().
That is: list of ``(dtype_key, age, area, acode, period, etype)`` tuples.
Also assumes that actions in list are sorted by applied period.
:param list schedule: The schedule of actions to apply.
:param int max_period: The maximum period to apply actions for.
If ``None``, defaults to the ``self.horizon``.
:param bool verbose: If True, prints additional information for debugging
purposes. Default is ``False``.
:param bool fail_on_missingarea: If True, raises an exception if missing
area is encountered. Default is ``False``.
:param bool force_integral_area: If True, forces the area to be integral.
Default is ``False``.
:param bool override_operability: If True, overrides operability limits.
Default is ``False``.
:param bool fuzzy_age: If True, attempts to apply action to proximal age
class if specified age is not operable. Default is ``True``.
:param bool recourse_enabled: If True, uses default AreaSelector to patch
missing area. Default is ``True``.
:param object areaselector: The AreaSelector object to use for patching
missing area. Default is ``None``.
:param bool compile_t_ycomps: If True, compiles time-indexed yield components.
Default is ``False``.
:param bool compile_c_ycomps: If True, compiles complex yield components.
Default is ``False``.
:param float rounding_bias: The rounding bias to use when forcing integral area.
Default is 0.15.
:param scale_area: The scaling factor to apply to the area. Default is ``None``.
:param bool reset: If True, resets the model before applying the schedule.
Default is ``True``.
:param bool crash_on_action_error: Crash on error applying action. Default is
``False``.
:return: The missing area (float) after applying the schedule.
"""
if max_period is None:
max_period = self.horizon
if reset:
self.reset()
missing_area = 0.
for _period in self.periods:
for dtype_key, age, area, acode, period, _etype in schedule:
if period != _period:
continue
if scale_area:
area = area * scale_area
if period > _period:
if verbose:
print('apply_schedule: committing actions for period', _period, f'(missing area {missing_area:0.1f})')
self.commit_actions(_period)
if period > max_period:
return
if force_integral_area:
area = round(area+rounding_bias)
if not area:
continue
assert not area % 1.
e, _aa, _ = self.apply_action(dtype_key,
acode,
period,
age,
area,
override_operability=override_operability,
fuzzy_age=fuzzy_age,
recourse_enabled=recourse_enabled,
areaselector=areaselector,
compile_t_ycomps=compile_t_ycomps,
compile_c_ycomps=compile_c_ycomps,
verbose=verbose)
if crash_on_action_error:
assert not e # crash on error (TO DO: better error handling)
else:
if e:
print('apply action error', e, dtype_key, acode, period, age, area)
if isinstance(_aa, float):
if fail_on_missingarea and missing_area:
raise
else:
missing_area += _aa
if verbose:
print(f'missing area {_aa:0.1f} ({_aa/area:0.2f})')
_period = period
self.commit_actions(_period)
return missing_area # type: ignore[return-value]
[docs]
def import_control_section(self, filename_suffix='run'):
"""
Imports CONTROL section from a Woodstock-formatted model input dataset.
.. warning::
Not implemented yet.
"""
pass
[docs]
def grow(self, start_period=1, cascade=True):
"""
Simulates growth (default start at period 1 and cascading to the end of the planning horizon).
Basically just calls :py:meth:`ws3.forest.DevelopmentType.grow` on all development types.
Growth in `ws3` just increments age---any other consequences of aging (e.g., tree height growth,
diameter growth, volume growth) is all in implicitly embedded in the yield curves.
:param int start_period: Period at which to start aging inventory
:param bool cascade: Will cascade growth to all future periods if ``True``, else only grows the specified period.
"""
for dt in list(self.dtypes.values()):
dt.grow(start_period, cascade)
def _cbm_sit_classifiers(self) -> Any:
"""
Compile ``sit_classifiers`` dataframe (part of CBM SIT input dataset).
:return: Dataframe containing CBM SIT classifiers table
:rtype: :py:class:`pandas.DataFrame`
"""
data: dict[str, list[Any]] = {'classifier_id':[], 'name':[], 'description':[]} # type: ignore[var-annotated]
for i, theme in enumerate(self._themes):
data['classifier_id'].append(i+1)
data['name'].append('_CLASSIFIER')
data['description'].append(theme['__name__'])
for v in self.theme_basecodes(i):
data['classifier_id'].append(i+1)
data['name'].append(v)
data['description'].append(v) # not great descriptions, but no effect on CBM output
data['classifier_id'].append(self.nthemes()+1)
data['name'].append('_CLASSIFIER')
data['description'].append('species')
data['classifier_id'].append(self.nthemes()+1)
data['name'].append('softwood')
data['description'].append('softwood')
data['classifier_id'].append(self.nthemes()+1)
data['name'].append('hardwood')
data['description'].append('hardwood')
result = pd.DataFrame(data)
return result
def _cbm_sit_disturbance_types(self) -> Any:
"""
Compile ``sit_disturbance_types`` dataframe (part of CBM SIT input dataset).
:return: Dataframe containing CBM SIT disturbance types table
:rtype: :py:class:`pandas.DataFrame`
"""
acodes = [acode for acode in self.actions.keys() if acode != 'null'] + ['fire']
data = {'id':acodes, 'name':acodes}
result = pd.DataFrame(data)
return result
def _cbm_sit_age_classes(self) -> Any:
"""
Compile ``sit_age_classes`` dataframe (part of CBM SIT input dataset).
:return: Dataframe containing CBM SIT age class table
:rtype: :py:class:`pandas.DataFrame`
"""
data: dict[str, list[Any]] = {'name': ['age_0'],
'class_size': [0],
'start_year': [0],
'end_year': [0]}
for i, ac in enumerate(range(self.period_length,
self.max_age+self.period_length,
self.period_length)):
data['name'].append(f'age_{i+1}')
data['class_size'].append(self.period_length)
data['start_year'].append(ac - self.period_length + 1)
data['end_year'].append(ac)
result = pd.DataFrame(data)
return result
def _cbm_sit_inventory(self, softwood_volume_yname, hardwood_volume_yname, default_last_pass_disturbance='fire',
include_empty_dtypes=False):
"""
Compile sit_inventory dataframe (part of CBM SIT input dataset).
:param str softwood_volume_yname: Softwood yield component name
:param str hardwood_volume_yname: Hardwood yield component name
:param str default_last_pass_disturbance: Default name of last pass disturbance.
Defaults to ``'fire'``.
:param int default_landclass: Default land classification code. The landclass column
in the CBM SIT inventory table should contain integers in the range [0, 22], which CBM
maps to one of 23 UNFCCC land classes (see Table 3-1 in the the CBM-CFS3 user guide for details).
Uses the value of the ``landclass`` attribute of the corresponding development type (if defined),
otherwise default to 0 (i.e., 'forest land remaining forest land').
:return: Dataframe containing CBM SIT inventory table
:rtype: :py:class:`pandas.DataFrame`
"""
def leading_species(r):
"""
Determine if softwood or hardwood leading species by comparing softwood and hardwood
volume at peak mean annual increment (MAI) age (largest volume wins).
"""
dt = self.dtypes[tuple(r.tolist())]
svol_curve, hvol_curve = dt.ycomp(softwood_volume_yname), dt.ycomp(hardwood_volume_yname)
tvol_curve = svol_curve + hvol_curve
x_cmai = tvol_curve.mai().ytp().lookup(0)
return 'softwood' if svol_curve[x_cmai] > hvol_curve[x_cmai] else 'hardwood'
def landclass(r):
"""
Use the value of the ``landclass`` attribute of the corresponding development type (if defined),
otherwise default to ``default_landclass``.
"""
dt = self.dtypes[tuple(r.tolist())]
if hasattr(dt, 'landclass'):
return dt.landclass
else:
return '0'
def last_pass_disturbance(r):
"""
We use the value of the ``last_pass_disturbance`` attribute of the corresponding development
type (if defined), otherwise default to ``default_last_pass_disturbance``.
"""
dt = self.dtypes[tuple(r.tolist())]
if hasattr(dt, 'last_pass_disturbance'):
return dt.last_pass_disturbance
else:
return default_last_pass_disturbance
theme_cols = [theme['__name__'] for theme in self._themes]
data: dict[str, list[Any]] = {**{c:[] for c in theme_cols}, **{c:[] for c in ['age', 'area']}} # type: ignore[var-annotated]
for dtype_key in self.dtypes:
dt = self.dtypes[dtype_key]
if include_empty_dtypes:
if not dt._areas[0]:
dt._areas[0] = {0:0.} # force a non-null initial inventory
else:
if not dt._areas[0]:
continue # developement type not in initial inventory
for age, area in dt._areas[0].items():
for i, c in enumerate(theme_cols):
data[c].append(dtype_key[i])
data['age'].append(age)
data['area'].append(area)
result = pd.DataFrame(data)
age_series, area_series = result['age'], result['area']
result.drop(['age', 'area'], axis=1, inplace=True) # wrong column order
result['species'] = result[theme_cols].apply(leading_species, axis=1)
result['using_age_class'] = 'FALSE'
result['age'] = age_series
result['area'] = area_series
result['delay'] = 0
result['landclass'] = result[theme_cols].apply(landclass, axis=1)
result['historic_disturbance'] = 'fire'
result['last_pass_disturbance'] = result[theme_cols].apply(last_pass_disturbance, axis=1)
return result
def _cbm_sit_yield(self, softwood_volume_yname, hardwood_volume_yname, n_yield_vals):
"""
Compile ``sit_yields`` dataframe (part of CBM SIT input dataset).
:param str softwood_volume_yname: Softwood yield component name
:param str hardwood_volume_yname: Hardwood yield component name
:param int n_yield_vals: Number of yield values to compile
:return: Dataframe containing CBM SIT yield table
:rtype: :py:class:`pandas.DataFrame`
.. note:
This is where the ``ws3`` and CBM data models diverge and things get a bit messy.
The hack is as follows. First, define a bogus Model I optimization problem
in the current ``ForestModel`` instance. A side effect of generating
the Model I problem matrix is to dynamically create all possible
``DevelopmentType`` cases and add them to the current ``ForestModel`` instance.
Then we can can use ``ForestModel.unmask`` to get at least one valid
``DevelopmentType`` instance for each yield mask in ``ForestModel.yields``
(there could be multiple, but we can safely just grab the first one in the list).
From there, use the same species-grokking logic we used to compile
the species column in ``self.sit_inventory``. Ugly, but should work
as long as ``ForestModel.yields`` includes full-wildcard softwood and
hardwood complex yield curves defined with the ``_SUM(...)`` function (because
this way all ``DevelopmentType`` instances, initially present or dynamically
created, will automatically have softwood and hardwood yield curves defined).
The bogus Model I problem can just use a bogus z coefficient
function (that always returns 0) and no flow or general constraints.
Maybe there is a better way?
"""
def leading_species(dt):
"""
Determine if softwood or hardwood leading species by comparing softwood and hardwood
volume at peak mean annual increment (MAI) age.
"""
svol_curve, hvol_curve = dt.ycomp(softwood_volume_yname), dt.ycomp(hardwood_volume_yname)
tvol_curve = svol_curve + hvol_curve
x_cmai = tvol_curve.mai().ytp().lookup(0)
return 'softwood' if svol_curve[x_cmai] > hvol_curve[x_cmai] else 'hardwood'
schedule = self.compile_schedule()
self.add_problem('__cbm_sit_bogus', {'z':(lambda forestmodel, path: 0.)})
theme_cols = [theme['__name__'] for theme in self._themes]
data: dict[str, list[Any]] = {**{c:[] for c in theme_cols}, # type: ignore[var-annotated]
'species':[], 'leading_species':[],
**{f'v{i}':[] for i in range(n_yield_vals + 1)}}
for dtype_key in self.dtypes:
dt = self.dt(dtype_key)
dt.leading_species = leading_species(dt) # type: ignore[no-untyped-call]
for species, yname in zip(('softwood', 'hardwood'), (softwood_volume_yname, hardwood_volume_yname), strict=False):
for i, c in enumerate(theme_cols):
data[c].append(dtype_key[i])
data['species'].append(species)
data['leading_species'].append(dt.leading_species)
for i in range(n_yield_vals + 1):
data[f'v{i}'].append(dt.ycomp(yname)[i * self.period_length])
result = pd.DataFrame(data)
self.apply_schedule(schedule) # running add_problem above broke the schedule so restore from backup we stashed
return result
def _cbm_sit_events(self) -> Any:
"""
Compile ``sit_events`` dataframe (part of CBM SIT input dataset).
:return: Dataframe containing CBM SIT events table
:rtype: :py:class:`pandas.DataFrame`
"""
theme_cols = [theme['__name__'] for theme in self._themes]
columns = theme_cols.copy()
columns += ['species',
'using_age_class',
'min_softwood_age',
'max_softwood_age',
'min_hardwood_age',
'max_hardwood_age',
'MinYearsSinceDist',
'MaxYearsSinceDist',
'LastDistTypeID',
'MinTotBiomassC',
'MaxTotBiomassC',
'MinSWMerchBiomassC',
'MaxSWMerchBiomassC',
'MinHWMerchBiomassC',
'MaxHWMerchBiomassC',
'MinTotalStemSnagC',
'MaxTotalStemSnagC',
'MinSWStemSnagC',
'MaxSWStemSnagC',
'MinHWStemSnagC',
'MaxHWStemSnagC',
'MinTotalStemSnagMerchC',
'MaxTotalStemSnagMerchC',
'MinSWMerchStemSnagC',
'MaxSWMerchStemSnagC',
'MinHWMerchStemSnagC',
'MaxHWMerchStemSnagC',
'efficiency',
'sort_type',
'target_type',
'target',
'disturbance_type',
'disturbance_year']
data: dict[str, list[Any]] = {c:[] for c in columns} # type: ignore[var-annotated]
for dtype_key, _age, area, acode, period, _ in self.compile_schedule():
#set_trace()
for i, c in enumerate(theme_cols):
data[c].append(dtype_key[i])
data['species'].append(self.dt(dtype_key).leading_species)
data['using_age_class'].append('FALSE')
#############################################################################
# might need to be more flexible with age range, to avoid OBO bugs and such?
data['min_softwood_age'].append(-1)
data['max_softwood_age'].append(-1)
data['min_hardwood_age'].append(-1)
data['max_hardwood_age'].append(-1)
#############################################################################
for c in columns[len(theme_cols)+6:-6]:
data[c].append(-1)
data['efficiency'].append(1)
data['sort_type'].append(3) # oldest first (see Table 3-3 in the CBM-CFS3 user guide)
data['target_type'].append('A') # area target
data['target'].append(area)
data['disturbance_type'].append(acode)
data['disturbance_year'].append(period*self.period_length)
result = pd.DataFrame(data)
return result
def _cbm_sit_transitions(self, null_acode: str = 'null') -> Any:
"""
Compile ``sit_transitions`` dataframe (part of CBM SIT input dataset).
:param str null_acode: Null action acode. Defaults to ``'null'``.
:return: Dataframe containing CBM SIT transitions table
:rtype: :py:class:`pandas.DataFrame`
"""
def resolve_target(dtype_key, target, sage):
tmask, tprop, tyield, tage, tlock, treplace, tappend = target # unpack tuple
dtk = list(dtype_key) # start with source key
dtk = [t if tmask[i] == '?' else tmask[i] for i, t in enumerate(dtk)]
if treplace:
dtk[treplace[0]] = self.resolve_replace(dtk, treplace[1]) # type: ignore[arg-type]
if tappend:
dtk[tappend[0]] = self.resolve_append(dtk, tappend[1]) # type: ignore[arg-type, assignment] # type: ignore[assignment]
dtk[tappend[0]] = self.resolve_append(dtk, tappend[1]) # type: ignore[arg-type, assignment] # type: ignore[assignment]
dtk = tuple(dtk) # type: ignore[assignment]
self.resolve_targetage(dtk, tyield, sage, tage, acode) # type: ignore[arg-type]
theme_cols = [theme['__name__'] for theme in self._themes]
columns = theme_cols.copy()
columns += ['species',
'using_age_class',
'min_softwood_age',
'max_softwood_age',
'min_hardwood_age',
'max_hardwood_age',
'disturbance_type']
columns += [f'to_{c}' for c in theme_cols]
columns += ['to_species',
'regen_delay',
'reset_age',
'percent']
data: dict[str, list[Any]] = {c:[] for c in columns} # type: ignore[var-annotated]
for dtype_key in self.dtypes:
dt = self.dt(dtype_key)
for acode, sage in dt.transitions:
if acode == null_acode:
continue
for target in dt.transitions[acode, sage]:
for i, c in enumerate(theme_cols):
data[c].append(dtype_key[i])
data['species'].append(dt.leading_species)
data['using_age_class'].append('FALSE')
data['min_softwood_age'].append(sage)
data['max_softwood_age'].append(sage)
data['min_hardwood_age'].append(sage)
data['max_hardwood_age'].append(sage)
data['disturbance_type'].append(acode)
to_dtype_key, to_age = resolve_target(dtype_key, target, sage) # type: ignore[no-untyped-call]
for i in range(len(theme_cols)):
data[f'to_theme{i}'].append(to_dtype_key[i]) # Monkey patch
target_dt = self.dt(to_dtype_key)
target_species = target_dt.leading_species if target_dt else dt.leading_species
data['to_species'].append(target_species)
data['regen_delay'].append(0)
data['reset_age'].append(to_age)
data['percent'].append(int(target[1] * 100))
result = pd.DataFrame(data)
return result
# --- orphaneed code? delete? ---
# for acode in fm.transitions:
# if acode == null_acode: continue
# for smask in self.transitions[acode]:
# tmask, tprop, _, _, _, _, _ = self.transitions[acode][smask][''][0]
# for i, c in enumerate(theme_cols): data[c].append(smask[i])
# data['species'].append('softwood' if au_table1.loc[int(smask[2])].canfi_species < 1200 else 'hardwood')
# data['using_age_class'].append('FALSE')
# data['min_softwood_age'].append(1)
# data['max_softwood_age'].append(999)
# data['min_hardwood_age'].append(1)
# data['max_hardwood_age'].append(999)
# data['disturbance_type'].append('harvest')
# for i in range(len(theme_cols)): data['to_theme%i' % i].append(tmask[i])
# data['to_%s' % species_classifier_colname].append('softwood' if au_table2.loc[int(tmask[4])].canfi_species < 1200 else 'hardwood')
# data['regen_delay'].append(0)
# data['reset_age'].append(0)
# data['percent'].append(100)
# result = pd.DataFrame(data)
# return result
[docs]
def to_cbm_sit(self, softwood_volume_yname, hardwood_volume_yname, admin_boundary, eco_boundary,
disturbance_type_mapping,
export_csv=False, sit_data_path='', default_last_pass_disturbance='fire', n_yield_vals=100,
include_empty_dtypes=False):
"""
Exports model data in a CBM standard import tool (SIT) data exchange format.
Calls several private methods to compile individual CBM SIT tables for the current
:py:class:`ws3.forest.ForestModel` instance.
:param str softwood_volume_yname: The yield component name for softwood volume.
:param str hardwood_volume_yname: The yield component name for hardwood volume.
:param str admin_boundary: The administrative boundary for spatial units mapping.
:param str eco_boundary: The ecological boundary for spatial units mapping.
:param dict disturbance_type_mapping: A dictionary containing disturbance type mapping information.
:param bool export_csv: Flag indicating whether to export data to CSV files. Default is False.
:param str sit_data_path: The path to export CSV files. Default is empty string.
:param str default_last_pass_disturbance: The default last pass disturbance type. Default is 'fire'.
:param int n_yield_vals: The number of yield values. Default is 100.
:return tuple: Tuple of ``sit_config`` (JSON-like dict namespace) and ``sit_tables`` (dict of
:py:class:`pandas.DataFrame` objects).
"""
sit_config = {
'mapping_config': {
'nonforest': None,
'species': {
'species_classifier': 'species',
'species_mapping': [
{'user_species': 'softwood', 'default_species': 'Softwood forest type'},
{'user_species': 'hardwood', 'default_species': 'Hardwood forest type'}
]
},
'spatial_units': {
'mapping_mode': 'SingleDefaultSpatialUnit',
'admin_boundary': admin_boundary,
'eco_boundary': eco_boundary
},
'disturbance_types': {
'disturbance_type_mapping': disturbance_type_mapping
}
}
}
sit_yield = self._cbm_sit_yield(softwood_volume_yname='swdvol', # type: ignore[no-untyped-call]
hardwood_volume_yname='hwdvol', n_yield_vals=100)
sit_tables = {'sit_classifiers':self._cbm_sit_classifiers(),
'sit_disturbance_types':self._cbm_sit_disturbance_types(),
'sit_age_classes':self._cbm_sit_age_classes(),
'sit_yield':sit_yield,
'sit_inventory':self._cbm_sit_inventory(softwood_volume_yname='swdvol', # type: ignore[no-untyped-call]
hardwood_volume_yname='hwdvol',
include_empty_dtypes=include_empty_dtypes),
'sit_events':self._cbm_sit_events(),
'sit_transitions':self._cbm_sit_transitions()
}
return sit_config, sit_tables
if __name__ == '__main__':
pass