Variance in `scipy.special.lambertw` output depending on scipy version and installation source

[kandersolar]

Originally discovered in #1996 (comment)

It seems that scipy.special.lambertw can return slightly different values from 1.12.0, depending on how it is installed (with conda or from PyPI with pip). I observe this difference locally on Windows:

>>> import scipy
>>> argW = 6.588940213030326e-08
>>> scipy.special.lambertw(argW)
(6.588939778889038e-08+0j)  # v1.12.0 (PyPI)
(6.588939778889037e-08+0j)  # v1.12.0 (conda)

As an aside, mpmath says that the one ending in 8 is the more accurate of the two:

>>> import mpmath
>>> mpmath.mp.dps = 25
>>> mpmath.lambertw(argW)
mpf('0.0000000658893977888903808156892075')

The relevance is that the particular argW value above was (implicitly) used to generate negative values of v_oc using _lambertw_v_from_i for testing in test_singlediode.py::test_singlediode_lambert_negative_voc. The change from ...8 to ...7 means that the returned voltage is no longer negative and thus causes the test to fail.

Other information: scipy.special.lambertw was translated from cython to C++ for 1.12.0 in scipy/scipy#19435.

What should we do about it? I think the first step is to ask the scipy folks if the difference is considered a bug in scipy that they will address. I will ask the scipy maintainers about this. If they say it is not a bug, then I suppose we should identify alternative values of argW that result in negative Voc for both versions of the scipy calculation and edit our test.

[adriesse]

Without looking at the test details, my gut feeling would be that the test case is much too close to the edge, so to speak, if the last decimal place matters. Can you move away from the edge?

[kandersolar]

Feedback from scipy maintainers: this level of variance is expected across different scipy versions and distributions and is not something they view as a bug to be fixed. So we should either find a different magic input and/or adjust the test strategy to not rely on the full precision of scipy's output, as @adriesse suggests.

[echedey-ls]

Credits to @cwhanse for the following comment:

I can't reproduce this test error on my Windows system with scipy 1.12 and numpy 1.26.4. There were some changes to scipy's lambertw implementation for scipy 1.12 but I can't see that those changes would affect the values for that failed test. I don't know how to look to see if there's some Windows library that is different for scipy 1.10 vs. 1.12.

Locally it happens the same to me: no failures. However, on a remote Win11 machine, I have installed conda and what I find strange is that, compared to bare python, it is the only one that gives the following failure:

================================== FAILURES ===================================
____________________ test_singlediode_lambert_negative_voc ____________________

    def test_singlediode_lambert_negative_voc():
    
        # Those values result in a negative v_oc out of `_lambertw_v_from_i`
        x = np.array([0., 1.480501e-11, 0.178, 8000., 1.797559])
        outs = pvsystem.singlediode(*x, method='lambertw')
>       assert outs['v_oc'] == 0
E       assert 1.3234889800848443e-23 == 0

pvlib\tests\test_singlediode.py:175: AssertionError

Conda and Bare python installed environment packages versions listed and outputs after running the Conda and then Bare (scipy1.13.0), then restarting the remote and only doing Bare(sicpy1.13.0), then checking again with Bare(scipy1.12.0).

I think we can isolate this problem to only happen on conda, possibly due to some difference in their recipe/runner OS config in contrast to the scipy procedure. Anyways, that does not change the fact that these tests require some rework.

And by the way, there are a few warnings too. IMO we should strive to have clean tests (edit: I'm referring to deprecationwarnings specially).

Edit: output from the stated runs in the remote Win11 machine linked in this Mega folder.

[cwhanse]

Without looking at the test details, my gut feeling would be that the test case is much too close to the edge, so to speak, if the last decimal place matters. Can you move away from the edge?

I don't think we want to do that. The point of that test is to confirm that small negative Voc is converted to 0, without large negative Voc being converted (as those indicate a problem with the inputs). See [these lines](

pvlib-python/pvlib/singlediode.py

Line 786 in ccf74ea

# Set small elements <0 in v_oc to 0

I think we want a test that produces small negative Voc consistently, which is not easy to formulate.

[echedey-ls]

@cwhanse so I should find which range of inputs generate a v_oc in the range(-1e-12, 0) in my PR, right?

[cwhanse]

That test was written to satisfy code coverage requirements after we decided to set small negative Voc to 0, since these Voc values can occur due to round-off error. Options:

1. find a new set of inputs that produces small negative Voc for all environments etc. As @adriesse points out, we're working at the edges of precision, so a new test may fail in the future, and maybe this kind of test wasn't the best choice for testing those lines to begin with.
2. change how we are handling small negative Voc. I don't have any ideas here.
3. change how the single diode equations are solved (default to bishop88?). That would cover up the problem for most users, but not fix it.
4. Change how lambertw_v_from_i method does its calculation. I don't have available time for this, and it is likely to be fiddly work.
5. Live with the test failure (ugh) or delete the test and live with reduced coverage.

Open to other ideas.

[echedey-ls]

I started working on no. 1, but I've realised we can just mock the output of _lambert_v_from_i for the v_oc. What do you think of the latter?

No. 1 script

import numpy as np
from pvlib.singlediode import _lambertw_v_from_i

from functools import partial


test_premise = {
    "photocurrent": 0.0,
    "saturation_current": 1.480501e-11,
    "resistance_series": 0.178,
    "resistance_shunt": 8000.0,
    "nNsVth": 1.797559,
}
variable_sweep = "saturation_current"
percentage_sweep = 0.5
bin_min, bin_max = -1e-12, 0.0
partial_kwargs = {k: v for k, v in test_premise.items() if k != variable_sweep}
voc_lambertw = partial(_lambertw_v_from_i, **partial_kwargs)

param_sweep = np.linspace(
    test_premise[variable_sweep] * (1 - percentage_sweep),
    test_premise[variable_sweep] * (1 + percentage_sweep),
    10000,
)

voc_returned = voc_lambertw(current=0.0, **{variable_sweep: param_sweep})
param_min, param_max = param_sweep[
    (voc_returned < bin_max) & (voc_returned > bin_min)
][(0, -1),]

print(
    f"Input range of {variable_sweep} where _lambertw_v_from_i is between {bin_min} and {bin_max}: {param_min} to {param_max}"
)

fulfills_condition = param_min < test_premise[variable_sweep] < param_max

print(
    f"Condition param_min < {variable_sweep} < param_max: {fulfills_condition}"
)
assert fulfills_condition

I'll be trying this script on a Win11 x64 with conda.

EDIT1: regarding the mocking, I don't expect it to create a problem in the long-run since what I understand is that the v_oc check is a fast path for the calculation. If I'm wrong and it is there to avoid problems, then we should choose another approach.

[echedey-ls]

I ended up doing the sixth option, mocking of _lambert_v_from_i in #2007

Mind giving it a review? Also, thanks for all the insight regarding what the test was about.

[adriesse]

Open to other ideas.

It seems like there's a working test solution already, but could the edge have been softened by using the range range (-1e-6, 0) rather than range (-1e-12, 0) ?

[adriesse]

Open to other ideas.

It seems like there's a working test solution already, but could the edge have been softened by using the range range (-1e-6, 0) rather than range (-1e-12, 0) ?

Still curious, although this is now closed.