Sample code for log-likelihood calibration¶

About this document¶

In [1]:

from swift2.doc_helper import pkg_versions_info

print(pkg_versions_info("This document was generated from a jupyter notebook"))

from swift2.doc_helper import pkg_versions_info print(pkg_versions_info("This document was generated from a jupyter notebook"))

This document was generated from a jupyter notebook on 2025-03-27 17:24:58.815735
    swift2 2.5.1
    uchronia 2.6.2

Setting up a calibration on daily data¶

We will use some sample data from (MMH) included in the package

In [2]:

import numpy as np
import pandas as pd
import xarray as xr

import numpy as np import pandas as pd import xarray as xr

In [3]:

from cinterop.timeseries import as_timestamp
from swift2.doc_helper import get_free_params, sample_series, set_loglik_param_keys
from swift2.parameteriser import (
    concatenate_parameterisers,
    create_parameter_sampler,
    create_parameteriser,
    create_sce_termination_wila,
    extract_optimisation_log,
    get_default_sce_parameters,
    parameteriser_as_dataframe,
    sort_by_score,
)
from swift2.simulation import create_subarea
from swift2.utils import c, mk_full_data_id, paste0
from swift2.vis import OptimisationPlots

s = as_timestamp('1990-01-01')
e = as_timestamp('2005-12-31')

rain = sample_series('MMH', 'rain')[slice(s, e)]
evap = sample_series('MMH', 'evap')[slice(s, e)]
flow = sample_series('MMH', 'flow')[slice(s, e)]

from cinterop.timeseries import as_timestamp from swift2.doc_helper import get_free_params, sample_series, set_loglik_param_keys from swift2.parameteriser import ( concatenate_parameterisers, create_parameter_sampler, create_parameteriser, create_sce_termination_wila, extract_optimisation_log, get_default_sce_parameters, parameteriser_as_dataframe, sort_by_score, ) from swift2.simulation import create_subarea from swift2.utils import c, mk_full_data_id, paste0 from swift2.vis import OptimisationPlots s = as_timestamp('1990-01-01') e = as_timestamp('2005-12-31') rain = sample_series('MMH', 'rain')[slice(s, e)] evap = sample_series('MMH', 'evap')[slice(s, e)] flow = sample_series('MMH', 'flow')[slice(s, e)]

In [4]:

rain.describe()

rain.describe()

Out[4]:

count    5844.000000
mean        3.545405
std         7.737554
min         0.000000
25%         0.000000
50%         0.283600
75%         3.308775
max        97.645500
dtype: float64

In [5]:

flow.describe()

flow.describe()

Out[5]:

count    5844.000000
mean       -1.993059
std        16.361702
min       -99.999000
25%         0.194400
50%         0.438400
75%         0.900200
max        17.221100
dtype: float64

We need to adjust the observed flow, as the SWIFTv1 legacy missing value code is -99.

In [6]:

flow[flow < 0] = np.nan

flow[flow < 0] = np.nan

In [7]:

flow

flow

Out[7]:

1990-01-01    0.2577
1990-01-02    0.2459
1990-01-03    0.2374
1990-01-04    0.2218
1990-01-05    0.2127
               ...  
2005-12-27    0.3477
2005-12-28    0.3314
2005-12-29    0.3333
2005-12-30    0.3066
2005-12-31    0.2896
Length: 5844, dtype: float64

Catchment setup¶

Let's create a single catchment setup, using daily data. We need to specify the simulation time step to be consistent with the daily input data.

In [8]:

ms = create_subarea('GR4J', 1.0)
from cinterop.timeseries import xr_ts_end, xr_ts_start

s = xr_ts_start(rain)
e = xr_ts_end(rain)
ms.set_simulation_span(s, e)
ms.set_simulation_time_step('daily')

ms = create_subarea('GR4J', 1.0) from cinterop.timeseries import xr_ts_end, xr_ts_start s = xr_ts_start(rain) e = xr_ts_end(rain) ms.set_simulation_span(s, e) ms.set_simulation_time_step('daily')

Assign input time series

In [9]:

sa_name = ms.get_subarea_names()[0]
ms.play_subarea_input(rain, sa_name, "P")
ms.play_subarea_input(evap, sa_name, "E")

sa_name = ms.get_subarea_names()[0] ms.play_subarea_input(rain, sa_name, "P") ms.play_subarea_input(evap, sa_name, "E")

Model variables identifiers are hierarchical, with separators '.' and '|' supported. The "dot" notation should now be preferred, as some R functions producing data frames may change the variable names and replace some characters with '.'.

In [10]:

sa_id = paste0("subarea.", sa_name)
root_id = paste0(sa_id, ".")
print(ms.get_variable_ids(sa_id))

sa_id = paste0("subarea.", sa_name) root_id = paste0(sa_id, ".") print(ms.get_variable_ids(sa_id))

['subarea.Subarea.areaKm2', 'subarea.Subarea.P', 'subarea.Subarea.E', 'subarea.Subarea.En', 'subarea.Subarea.LAI', 'subarea.Subarea.runoff', 'subarea.Subarea.S', 'subarea.Subarea.R', 'subarea.Subarea.TWS', 'subarea.Subarea.Eactual', 'subarea.Subarea.Ps', 'subarea.Subarea.Es', 'subarea.Subarea.Pr', 'subarea.Subarea.ech1', 'subarea.Subarea.ech2', 'subarea.Subarea.Perc', 'subarea.Subarea.alpha', 'subarea.Subarea.k', 'subarea.Subarea.x1', 'subarea.Subarea.x2', 'subarea.Subarea.x3', 'subarea.Subarea.x4', 'subarea.Subarea.UHExponent', 'subarea.Subarea.PercFactor', 'subarea.Subarea.OutflowVolume', 'subarea.Subarea.OutflowRate']

In [11]:

gr4_state_names = paste0(root_id, c('runoff', 'S', 'R', 'Perc'))
for name in gr4_state_names: 
    ms.record_state(name)

gr4_state_names = paste0(root_id, c('runoff', 'S', 'R', 'Perc')) for name in gr4_state_names: ms.record_state(name)

Let's check that one simulation runs fine, before we build a calibration definition.

In [12]:

ms.exec_simulation()
sState = ms.get_recorded(gr4_state_names[2])

ms.exec_simulation() sState = ms.get_recorded(gr4_state_names[2])

In [13]:

sState.plot(figsize=(10,4))

sState.plot(figsize=(10,4))

Out[13]:

[<matplotlib.lines.Line2D at 0x7f04f7a3d850>]

Let's build the objective calculator that will guide the calibration process:

In [14]:

w = pd.Timestamp("1992-01-01")

w = pd.Timestamp("1992-01-01")

In [15]:

runoff_depth_varname = 'subarea.Subarea.runoff'
mod_runoff = ms.get_recorded(runoff_depth_varname)
# zoo::index(flow) = zoo::index(mod_runoff)
objective = ms.create_objective(runoff_depth_varname, flow, 'log-likelihood', w, e)

runoff_depth_varname = 'subarea.Subarea.runoff' mod_runoff = ms.get_recorded(runoff_depth_varname) # zoo::index(flow) = zoo::index(mod_runoff) objective = ms.create_objective(runoff_depth_varname, flow, 'log-likelihood', w, e)

In [16]:

mod_runoff.plot()

mod_runoff.plot()

Out[16]:

[<matplotlib.lines.Line2D at 0x7f04f7aed490>]

Parameterisation¶

Define the feasible parameter space, using a generic parameter set for the model parameters. This is 'wrapped' by a log-likelihood parameter set with the extra parameters used in the log likelihood calculation, but which exposes all the parameters as 8 independent degrees of freedom to the optimiser.

In [17]:

pspec_gr4j = get_free_params('GR4J')
pspec_gr4j.Value = c(542.1981111, -0.4127542, 7.7403390, 1.2388548)
pspec_gr4j.Min = c(1,-30, 1,1)
pspec_gr4j.Max = c(3000, 30, 1000, 240)
pspec_gr4j.Name = paste0(root_id, pspec_gr4j.Name)


maxobs = np.max(flow)
p = create_parameteriser(type='Generic', specs=pspec_gr4j)
set_loglik_param_keys(a='a', b='b', m='m', s='s', ct="ct", censopt='censopt')
censor_threshold = maxobs / 100 # TBC
censopt = 0.0

loglik = create_parameteriser(type='no apply')
loglik.add_to_hypercube( 
          pd.DataFrame({ 
          "Name": c('b','m','s','a','maxobs','ct', 'censopt'),
          "Min": c(-30, 0, -10,    -20, maxobs, censor_threshold, censopt),
          "Max":  c(5,   0, 10, 0, maxobs, censor_threshold, censopt),
          "Value": c(-7,  0, 0,  -10, maxobs, censor_threshold, censopt),
          }
          ) )
p = concatenate_parameterisers(p, loglik)
p.as_dataframe()

pspec_gr4j = get_free_params('GR4J') pspec_gr4j.Value = c(542.1981111, -0.4127542, 7.7403390, 1.2388548) pspec_gr4j.Min = c(1,-30, 1,1) pspec_gr4j.Max = c(3000, 30, 1000, 240) pspec_gr4j.Name = paste0(root_id, pspec_gr4j.Name) maxobs = np.max(flow) p = create_parameteriser(type='Generic', specs=pspec_gr4j) set_loglik_param_keys(a='a', b='b', m='m', s='s', ct="ct", censopt='censopt') censor_threshold = maxobs / 100 # TBC censopt = 0.0 loglik = create_parameteriser(type='no apply') loglik.add_to_hypercube( pd.DataFrame({ "Name": c('b','m','s','a','maxobs','ct', 'censopt'), "Min": c(-30, 0, -10, -20, maxobs, censor_threshold, censopt), "Max": c(5, 0, 10, 0, maxobs, censor_threshold, censopt), "Value": c(-7, 0, 0, -10, maxobs, censor_threshold, censopt), } ) ) p = concatenate_parameterisers(p, loglik) p.as_dataframe()

Out[17]:

	Name	Value	Min	Max
0	subarea.Subarea.x1	542.198111	1.000000	3000.000000
1	subarea.Subarea.x2	-0.412754	-30.000000	30.000000
2	subarea.Subarea.x3	7.740339	1.000000	1000.000000
3	subarea.Subarea.x4	1.238855	1.000000	240.000000
4	b	-7.000000	-30.000000	5.000000
5	m	0.000000	0.000000	0.000000
6	s	0.000000	-10.000000	10.000000
7	a	-10.000000	-20.000000	0.000000
8	maxobs	17.221100	17.221100	17.221100
9	ct	0.172211	0.172211	0.172211
10	censopt	0.000000	0.000000	0.000000

Check that the objective calculator works, at least with the default values in the feasible parameter space:

In [18]:

score = objective.get_score(p)
print(score)

score = objective.get_score(p) print(score)

{'scores': {'Log-likelihood': -1e+20}, 'sysconfig':                   Name       Value        Min          Max
0   subarea.Subarea.x1  542.198111   1.000000  3000.000000
1   subarea.Subarea.x2   -0.412754 -30.000000    30.000000
2   subarea.Subarea.x3    7.740339   1.000000  1000.000000
3   subarea.Subarea.x4    1.238855   1.000000   240.000000
4                    b   -7.000000 -30.000000     5.000000
5                    m    0.000000   0.000000     0.000000
6                    s    0.000000 -10.000000    10.000000
7                    a  -10.000000 -20.000000     0.000000
8               maxobs   17.221100  17.221100    17.221100
9                   ct    0.172211   0.172211     0.172211
10             censopt    0.000000   0.000000     0.000000}

In [19]:

mod_runoff = ms.get_recorded(runoff_depth_varname)

mod_runoff = ms.get_recorded(runoff_depth_varname)

In [20]:

from swift2.vis import plot_two_series

from swift2.vis import plot_two_series

In [21]:

plot_two_series(flow, mod_runoff, ylab="obs/mod runoff", start_time = "2000-01-01", end_time = "2002-12-31", names=['observed','modelled'])

plot_two_series(flow, mod_runoff, ylab="obs/mod runoff", start_time = "2000-01-01", end_time = "2002-12-31", names=['observed','modelled'])

Calibration¶

Build the optimiser definition, instrument with a logger.

In [22]:

# term = getMaxRuntimeTermination(max_hours = 0.3/60)  # ~20 second appears enough with SWIFT binaries in Release mode
# term = getMarginalTermination(tolerance = 1e-06, cutoff_no_improvement = 10, max_hours = 0.3/60) 
term = create_sce_termination_wila('relative standard deviation', c('0.005',str(1/60)))

sce_params = get_default_sce_parameters()
urs = create_parameter_sampler(0, p, 'urs')
optimiser = objective.create_sce_optim_swift(term, sce_params, urs)
calib_logger = optimiser.set_calibration_logger('')

# term = getMaxRuntimeTermination(max_hours = 0.3/60) # ~20 second appears enough with SWIFT binaries in Release mode # term = getMarginalTermination(tolerance = 1e-06, cutoff_no_improvement = 10, max_hours = 0.3/60) term = create_sce_termination_wila('relative standard deviation', c('0.005',str(1/60))) sce_params = get_default_sce_parameters() urs = create_parameter_sampler(0, p, 'urs') optimiser = objective.create_sce_optim_swift(term, sce_params, urs) calib_logger = optimiser.set_calibration_logger('')

In [23]:

%%time 
calib_results = optimiser.execute_optimisation()

%%time calib_results = optimiser.execute_optimisation()

CPU times: user 5min 7s, sys: 44.1 ms, total: 5min 7s
Wall time: 1min

In [24]:

opt_log = extract_optimisation_log(optimiser, fitness_name = 'Log-likelihood')
geom_ops = opt_log.subset_by_message(pattern= 'Initial.*|Reflec.*|Contrac.*|Add.*')

opt_log = extract_optimisation_log(optimiser, fitness_name = 'Log-likelihood') geom_ops = opt_log.subset_by_message(pattern= 'Initial.*|Reflec.*|Contrac.*|Add.*')

In [25]:

import matplotlib.pyplot as plt

import matplotlib.pyplot as plt

In [26]:

ll_max = max(geom_ops._data['Log-likelihood'].values)
ll_min = np.median(geom_ops._data['Log-likelihood'].values)

ll_max = max(geom_ops._data['Log-likelihood'].values) ll_min = np.median(geom_ops._data['Log-likelihood'].values)

Parameter plots¶

In [27]:

p_var_ids = p.as_dataframe().Name.values
v = OptimisationPlots(geom_ops)
for pVar in p_var_ids:
    g = v.parameter_evolution(pVar, obj_lims=[ll_min, ll_max])
    plt.gcf().set_size_inches(10,8)

p_var_ids = p.as_dataframe().Name.values v = OptimisationPlots(geom_ops) for pVar in p_var_ids: g = v.parameter_evolution(pVar, obj_lims=[ll_min, ll_max]) plt.gcf().set_size_inches(10,8)

Finally, get a visual of the runoff time series with the best known parameter set (the penultimate entry in the data frame with the log of the calibration process).

In [28]:

sortedResults = sort_by_score(calib_results, 'Log-likelihood')
sortedResults.as_dataframe().head().T

sortedResults = sort_by_score(calib_results, 'Log-likelihood') sortedResults.as_dataframe().head().T

Out[28]:

	0	1	2	3	4
Log-likelihood	3340.915859	3340.771324	3340.698917	3340.632166	3340.481544
subarea.Subarea.x1	142.311803	143.683797	141.209587	138.884465	138.401395
subarea.Subarea.x2	-29.926501	-29.985372	-29.892518	-29.987330	-29.915535
subarea.Subarea.x3	769.304425	771.741061	768.582954	768.664694	768.479620
subarea.Subarea.x4	1.001487	1.000453	1.001558	1.001635	1.002588
b	-1.069403	-1.180434	-1.182628	-0.981440	-1.247725
m	0.000000	0.000000	0.000000	0.000000	0.000000
s	0.394041	0.514402	0.517914	0.318156	0.597338
a	-7.887809	-7.851411	-8.134082	-7.777277	-8.522691
maxobs	17.221100	17.221100	17.221100	17.221100	17.221100
ct	0.172211	0.172211	0.172211	0.172211	0.172211
censopt	0.000000	0.000000	0.000000	0.000000	0.000000

In [29]:

best_pset = calib_results.get_best_score('Log-likelihood').parameteriser
best_pset.apply_sys_config(ms)
ms.exec_simulation()
mod_runoff = ms.get_recorded(runoff_depth_varname)
# joki::plot_two_series(flow, mod_runoff, ylab="obs/mod runoff", startTime = start(flow), endTime = end(flow))

best_pset = calib_results.get_best_score('Log-likelihood').parameteriser best_pset.apply_sys_config(ms) ms.exec_simulation() mod_runoff = ms.get_recorded(runoff_depth_varname) # joki::plot_two_series(flow, mod_runoff, ylab="obs/mod runoff", startTime = start(flow), endTime = end(flow))

In [30]:

mod_runoff

mod_runoff

Out[30]:

<xarray.DataArray (variable_identifiers: 1, ensemble: 1, time: 5844)> Size: 47kB
array([[[0.        , 0.        , 0.        , ..., 0.41868723,
         0.41110903, 0.40375453]]], shape=(1, 1, 5844))
Coordinates:
  * ensemble              (ensemble) int64 8B 0
  * time                  (time) datetime64[ns] 47kB 1990-01-01 ... 2005-12-31
  * variable_identifiers  (variable_identifiers) object 8B 'subarea.Subarea.r...

xarray.DataArray

• variable_identifiers: 1
• ensemble: 1
• time: 5844

• 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.4429 0.4346 0.4265 0.4187 0.4111 0.4038

  array([[[0.        , 0.        , 0.        , ..., 0.41868723,
           0.41110903, 0.40375453]]], shape=(1, 1, 5844))

• Coordinates: (3)
  • ensemble
    (ensemble)
    int64
    0

    array([0])

  • time
    (time)
    datetime64[ns]
    1990-01-01 ... 2005-12-31

    array(['1990-01-01T00:00:00.000000000', '1990-01-02T00:00:00.000000000',
           '1990-01-03T00:00:00.000000000', ..., '2005-12-29T00:00:00.000000000',
           '2005-12-30T00:00:00.000000000', '2005-12-31T00:00:00.000000000'],
          shape=(5844,), dtype='datetime64[ns]')

  • variable_identifiers
    (variable_identifiers)
    object
    'subarea.Subarea.runoff'

    array(['subarea.Subarea.runoff'], dtype=object)

• Indexes: (3)
  • ensemble
    PandasIndex

    PandasIndex(Index([0], dtype='int64', name='ensemble'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1990-01-01', '1990-01-02', '1990-01-03', '1990-01-04',
                   '1990-01-05', '1990-01-06', '1990-01-07', '1990-01-08',
                   '1990-01-09', '1990-01-10',
                   ...
                   '2005-12-22', '2005-12-23', '2005-12-24', '2005-12-25',
                   '2005-12-26', '2005-12-27', '2005-12-28', '2005-12-29',
                   '2005-12-30', '2005-12-31'],
                  dtype='datetime64[ns]', name='time', length=5844, freq='D'))

  • variable_identifiers
    PandasIndex

    PandasIndex(Index(['subarea.Subarea.runoff'], dtype='object', name='variable_identifiers'))

• Attributes: (0)

In [31]:

mod_runoff.squeeze(drop=True).sel(time=slice(e - pd.offsets.DateOffset(years=1), e)).plot(figsize=(16,9))

mod_runoff.squeeze(drop=True).sel(time=slice(e - pd.offsets.DateOffset(years=1), e)).plot(figsize=(16,9))

Out[31]:

[<matplotlib.lines.Line2D at 0x7f04df4f7e50>]

In [32]:

plot_two_series(flow, mod_runoff, ylab="obs/mod runoff", start_time = "2000-01-01", end_time = "2002-12-31", names=['observed','modelled'])

plot_two_series(flow, mod_runoff, ylab="obs/mod runoff", start_time = "2000-01-01", end_time = "2002-12-31", names=['observed','modelled'])