Response function plugin for Pastas

This notebook shows usage of a custom response function from the responses plugin for Pastas.

import pandas as pd
import pastas as ps

import pastas_plugins as pp
from pastas_plugins import responses

pp.show_plugin_versions()

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[2], line 1
----> 1 pp.show_plugin_versions()

File ~/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas_plugins/__init__.py:27, in show_plugin_versions()
     25 for plugin in plugins:
     26     try:
---> 27         module = import_module(f"pastas_plugins.{plugin}.version")
     28         version = module.__version__
     29     except ModuleNotFoundError:

File ~/.asdf/installs/python/3.11.12/lib/python3.11/importlib/__init__.py:126, in import_module(name, package)
    124             break
    125         level += 1
--> 126 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1204, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1176, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1126, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:241, in _call_with_frames_removed(f, *args, **kwds)

File <frozen importlib._bootstrap>:1204, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1176, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1147, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:690, in _load_unlocked(spec)

File <frozen importlib._bootstrap_external>:940, in exec_module(self, module)

File <frozen importlib._bootstrap>:241, in _call_with_frames_removed(f, *args, **kwds)

File ~/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas_plugins/modflow/__init__.py:9
      1 # ruff: noqa: F401
      2 from pastas_plugins.modflow.modflow import (
      3     ModflowDrn,
      4     ModflowDrnSto,
   (...)      7     ModflowUzf,
      8 )
----> 9 from pastas_plugins.modflow.stressmodels import ModflowModel
     10 from pastas_plugins.modflow.version import __version__

File ~/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas_plugins/modflow/stressmodels.py:7
      5 from pastas.stressmodels import StressModelBase
      6 from pastas.timeseries import TimeSeries
----> 7 from pastas.typing import ArrayLike, TimestampType
      9 from .modflow import Modflow
     11 logger = getLogger(__name__)

ImportError: cannot import name 'TimestampType' from 'pastas.typing' (/home/docs/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas/typing/__init__.py)

index = pd.date_range("2014-01-01", "2024-01-01", freq="D")
df = pd.read_csv("data/aftopping.csv", index_col=0, parse_dates=True)
stress = (df["Precipitation"] - df["Evaporation"]).multiply(1e3)


def get_model(rfunc: ps.rfunc.RfuncBase, oseries: pd.Series | None = None) -> ps.Model:
    # Create a Pastas model
    oseries = pd.Series(0.0, stress.index) if oseries is None else oseries

    ml = ps.Model(oseries)

    # Create a stress model with recharge response
    sm = ps.StressModel(stress, rfunc=rfunc, name="sm", settings="evap")
    ml.add_stressmodel(sm)

    return ml


ml_synth = get_model(ps.Exponential())
ml_synth.set_parameter("sm_A", initial=0.1)
ml_synth.set_parameter("sm_a", initial=100.0)
ml_synth.set_parameter("constant_d", initial=5.0)
head = ml_synth.simulate()

Model is not optimized yet, initial parameters are used.

The Theis response function

The Theis response function is based on the analytical solution for pumping between two infinitely long ditches.

ml = get_model(responses.Theis(), head)
ml.solve()

Fit report Simulation              Fit Statistics
=================================================
nfev    58                     EVP          96.02
nobs    5351                   R2            0.96
noise   False                  RMSE          0.01
tmin    2002-06-10 00:00:00    AICc     -46532.24
tmax    2017-02-01 00:00:00    BIC      -46505.91
freq    D                      Obj           0.45
warmup  3650 days 00:00:00     ___               
solver  LeastSquares           Interp.         No

Parameters (4 optimized)
=================================================
               optimal     initial  vary
sm_A          0.038215    0.258227  True
sm_a        198.906792  100.000000  True
sm_b          0.162112    0.001000  True
constant_d    4.998249    5.057874  True

axes = ml.plots.results()

The Edelman response function

ml_em = get_model(responses.Edelman(), head)
ml_em.solve()

Fit report Simulation              Fit Statistics
=================================================
nfev    32                     EVP           0.56
nobs    5351                   R2            0.01
noise   False                  RMSE          0.06
tmin    2002-06-10 00:00:00    AICc     -29314.06
tmax    2017-02-01 00:00:00    BIC      -29300.90
freq    D                      Obj          11.17
warmup  3650 days 00:00:00     ___               
solver  LeastSquares           Interp.         No

Parameters (2 optimized)
=================================================
             optimal   initial  vary
sm_beta     0.006867  1.000000  True
constant_d  5.051175  5.057874  True

Warnings! (2)
=================================================
Response tmax for 'sm' > than calibration period.
Response tmax for 'sm' > than warmup period.

ml_em.plot()

<Axes: ylabel='Head'>