one model of mulitmodel ensemble behaving strangely

[zentouro]

Hi, i'm not sure if this is an issue with my underlying data or input error or something else, but I am attempting to calculate SPI for 19 models and for every model except for CESM2 (and its sibling CESM2-WACCM) the outputs look reasonable. For the CESM2 models, the variance increases dramatically (and exactly at the switch from historical to future period/SSP).

At first I assumed this was an issue with the underlying data, but I also tried the xclim package, which gives a reasonable result.

I'd much rather use the climate_indices package than xclimate as it doesn't have as many dropped values and runs much quicker.
Happy to share the data and/or code I'm using, but was wondering if this is an issue anyone else has run into and/or if I've possibly made a silly error somewhere.
Thanks!

[monocongo]

Thanks for your issue submission, Miriam (@zentouro). It's always gratifying to see this code being used in the wild. I'm afraid I have little time to properly maintain this package so I can't promise a quick resolution. It's interesting you mention that this package is faster and produces fewer missing values than the implementation offered by @Ouranosinc/xclim. I have recently considered mothballing this repository and directing the users to xclim since it seems to be a much better codebase (it's actively developed and maintained, and it's just more robust). Still, I assumed that their SPI implementation worked at least as well as ours, so this is an interesting observation. They're currently working on an implementation of PET (see Ouranosinc/xclim#1306) and if that one's off as well then I'm not so sure...

[monocongo]

I've also just been reminded that the Logger class isn't supported by numba either, so we need to better isolate the computation-heavy code from any logging -- not impossible but not trivial either since I've included a good bit of logging throughout.

[monocongo]

So I have now gone down this rabbit hole a bit and preliminary results look like numba isn't helping since it was never being used in pure "no Python" mode. There are several reasons for this, but mostly because numba doesn't fully support some of the numpy functions being used (hstack, convolve), making some refactorization necessary if we want to get any real benefits from numba. Or at least that's my take so far. One simple test I've run is to use numba in object mode (the fallback that we ended up using by default in our previous version) with the loop-lifting feature enabled, to see how much if anything we gain with that. It looks like that's slower than using no numba at all. Here are the results (testing code listed below):

without numba decoorators:
Timings: [1.5961637909931596, 1.5906572910025716, 1.582102083004429, 1.581416042012279, 1.5790853749931557]

with numba.jit(forceobj=True, looplift=True):
Timings: [2.721167041992885, 1.7240668750018813, 1.7313793329958571, 1.7260585420008283, 1.7234851669927593]

with plain numba.jit
Timings: [4.165839624998625, 1.7441993749962421, 1.733445917008794, 1.7340991250093793, 1.7374107089999598]

This is the benchmarking code:

import os
import timeit
import numpy as np

from climate_indices import compute, indices

_CALIBRATION_YEAR_START_DAILY = 1998
_CALIBRATION_YEAR_END_DAILY = 2016
_DATA_YEAR_START_DAILY = 1998
_DATA_YEAR_START_MONTHLY = 1895
_DATA_YEAR_END_MONTHLY = 2017


def get_precips_mm_daily():
    return np.load(
        os.path.join(os.path.split(__file__)[0],
                     "fixture",
                     "precips_mm_daily.npy")
    )


def get_precips_mm_monthly():
    return np.load(
        os.path.join(os.path.split(__file__)[0],
                     "fixture",
                     "precips_mm_monthly.npy")
    )


def run_several_invocations(
    precips_mm_monthly: np.ndarray,
    precips_mm_daily: np.ndarray,
):
    # compute SPI/gamma 1-month
    indices.spi(
        precips_mm_monthly,
        1,
        indices.Distribution.gamma,
        _DATA_YEAR_START_MONTHLY,
        _DATA_YEAR_START_MONTHLY,
        _DATA_YEAR_END_MONTHLY,
        compute.Periodicity.monthly,
    )

    # compute SPI/gamma 6-month
    indices.spi(
        precips_mm_monthly.flatten(),
        6,
        indices.Distribution.gamma,
        _DATA_YEAR_START_MONTHLY,
        _DATA_YEAR_START_MONTHLY,
        _DATA_YEAR_END_MONTHLY,
        compute.Periodicity.monthly,
    )

    # compute SPI/gamma 30-day
    indices.spi(
        precips_mm_daily,
        30,
        indices.Distribution.gamma,
        _DATA_YEAR_START_DAILY,
        _CALIBRATION_YEAR_START_DAILY,
        _CALIBRATION_YEAR_END_DAILY,
        compute.Periodicity.daily,
    )

precips_mm_daily = get_precips_mm_daily()
precips_mm_monthly = get_precips_mm_monthly()

time_measurement = timeit.repeat(
    lambda: run_several_invocations(
        precips_mm_monthly=precips_mm_monthly,
        precips_mm_daily=precips_mm_daily,
    ),
    repeat=5,
    number=1000,
)

print(f"Timings: {time_measurement}")

So for now I'm reasonably convinced that we were getting no benefit from using numba before. But I'm certainly no expert at this, please advise if anything's off here, as I could easily be missing something obvious.

[monocongo]

@coxipi It seems that you've written your code so that it can be optimized with numba whereas mine has various gotchas that keep that from happening. So your implementation is likely to be much faster on big datasets. My concern however is that the results don't line up, i.e. the xclim SPI results are somehow different from the climate_indices SPI results, or at least that's what I have been led to believe. This goes against our initial motivation for providing a reference implementation, i.e. to facilitate standardization (otherwise every scientist can cook up some version of SPI and call it that, and then nobody's comparing apples to apples).

[coxipi]

Yes, I was mostly focused on performance, but getting the right results is important. We do have some tests in our testing suite comparing our values to yours, but there are only a few of them. I just took a 150 year, 20x20 grid for some simulation I had lying around and checked how the results of a monthly SPI-6 of both libraries compare (so in the end, 20x20x150x12 = 720 000 values to compare). I find that 97.6% of values have a relative difference below 0.001%, and 99.5% of values below 0.01% relative error, etc. So it seems fairly good. Also, I have the same number of NaNs in both cases (2000 = 20x20x5, i.e. the first 5 months in the first year), no more NaNs in xclim.

I could go deeper into this, but before, it would be good to know: @zentouro were you using a specific method ("ML" or "APP")? If you didn't specify it in your input, then the default should be "APP", but the example in the documentation uses method="ML". I would not be surprised if you got more NaNs with "ML" which relies on a scipy routine and I have seen it fail in some instances. In comparison, the "APP" method re-implements the nice and simple well-behaved approximation used in monocongo library.

By the way, @monocongo, I saw your comment that the routine to get gamma parameters should be better documented. A reference that finds the same formula is Thom, 1958 in eq.22, and it's a saddle point approximation of the Maximum likelihood. Not sure if that's how they do it or call it, but that's what I've done myself.

[monocongo]

@coxipi I'm curious as to why you're re-writing a new SPI for xclim -- what's wrong with the implementation provided by this package? This is disheartening because the entire reason for this package is to provide standardized implementations of commonly used indices so folks don't keep reinventing these wheels.