Update to Hector 3

[rgieseke]

• fix TODOs
• update and check tests/comparisons with vanilla Hector
• finish conversion of variable renamings, see http://jgcri.github.io/hector/articles/v3_output_changes.html
• fix versioneer config

[rgieseke]

Test script

import pyhector
from copy import deepcopy
from pyhector.default_config import _default_config

print(pyhector.__hector_version__)

h = pyhector.Hector()

h.set_emissions(pyhector.ssp126)

parameters = deepcopy(_default_config)
h.config(parameters)

h.run(until=2100)

Prints the expected Hector version and runs.

However, adding an observable fails:

import pyhector
from copy import deepcopy
from pyhector.default_config import _default_config

print(pyhector.__hector_version__)

h = pyhector.Hector()

h.set_emissions(pyhector.ssp126)

parameters = deepcopy(_default_config)
h.config(parameters)

h.add_observable(
    "simpleNbox",
    "atmos_co2",
    False
)

h.run(until=2100)

Error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/Work/pyhector/session.py:20
     12 h.config(parameters)
     14 h.add_observable(
     15     "simpleNbox",
     16     "atmos_co2",
     17     False
     18 )
---> 20 h.run(until=2100)

ValueError: vector::_M_fill_insert

[swillner]

(hint for enabling debug compilation: add to Pybind11Extension in setup.py: extra_compile_args=["-g3"])

[rgieseke]

With the now compiling and running version only a single value is returned in the array.

import pyhector
from copy import deepcopy
from pyhector.default_config import _default_config

print(pyhector.__hector_version__)

h = pyhector.Hector()

h.set_emissions(pyhector.ssp126)

parameters = deepcopy(_default_config)
h.config(parameters)

h.add_observable(
	"simpleNbox",
	"CO2_concentration",
	True
)

h.run()

out = h.get_observable("simpleNbox", "CO2_concentration")

print(out)

[rgieseke]

Pyhector v2.x results and output streams from the C++ command line version are in good agreement, but the tests are not catching larger differences, like volcanic forcing changes, missing albedo forcing, or parameter changes introduced with the Hector v3 update.

With Pyhector 2.5.0.2

import pyhector

df = pyhector.run(pyhector.rcp26)

outstream = pyhector.read_hector_output(
    "https://github.com/openclimatedata/pyhector/raw/main/tests/data/outputstream_rcp26.csv"
)

df["temperature.Tgav"].plot(color="orange")
outstream["Tgav"].plot(linestyle=":", color="blue")

I think the Pandas testing function uses numpy's isclose which actually does

absolute(a - b) <= (atol + rtol * absolute(b))

The test only changes the rtol default.

assert_series_equal(
    output["temperature.global_tas"],
    original.global_tas,
    check_names=False,
    check_exact=False,
    rtol=3,
)

[rgieseke]

I'm not entirely sure how command line Hector sets its output precision for the CSV file, but it appears to have four significant figures.
One way to make a comparison between Pyhector output and the CSV streams could be to compare the difference between these significant figures, ignoring the exponent and checking that there are only floating point rounding differences.

With 88dcebf

import os
import pyhector
from pyhector import (
    ssp119,
    ssp126,
    ssp245,
    ssp370,
    ssp434,
    ssp460,
    ssp534_over,
    ssp585,
    read_hector_output
)

import matplotlib.pyplot as plt


ssps = {
    "ssp119": ssp119,
    "ssp126": ssp126,
    "ssp245": ssp245,
    "ssp370": ssp370,
    "ssp434": ssp434,
    "ssp460": ssp460,
    "ssp534-over": ssp534_over,
    "ssp585": ssp585,
}


for name, scenario in list(ssps.items())[-2:-1]:
    df = pyhector.run(scenario)
    out = read_hector_output(os.path.join("tests/data/outputstream_{}.csv".format(name)))
    
    fig, ax = plt.subplots(1,3, figsize=(14,6))
    out["global_tas"].plot(label="cli outputstream", ax=ax[0])
    df["temperature.global_tas"].plot(label="pyhector", ax=ax[0])
    ax[0].legend()
    ax[0].set_title("Temperature")
    
    out["CO2_concentration"].plot(label="cli outputstream", ax=ax[1])
    df["simpleNbox.CO2_concentration"].plot(label="pyhector", ax=ax[1])
    ax[1].legend()
    ax[1].set_title("Concentration")

    (out["global_tas"].map(lambda x: float(f'{x:.5e}'[:7])) - 
     df["temperature.global_tas"].map(lambda x: float(f'{x:.5e}'[:7]))).plot(ax=ax[2])
    ax[2].set_title('Subtracting significant digits (Temp.)')
    
    fig.suptitle(name)

    assert (out["global_tas"].map(lambda x: float(f'{x:.5e}'[:7])) - 
        df["temperature.global_tas"].map(lambda x: float(f'{x:.5e}'[:7])).abs() < 6e-4).all()

This should catch small changes, e.g. like

df = pyhector.run(scenario, {"temperature": {"S": 3.01}})

[swillner]

awesome!