Speed up Brent how-methods in single-diode using Scipy Cython Optimize

[mikofski]

Problem

As described in #968 singlediode.bishop88_i_from_v() and singlediode.bishop88_v_from_i() can be very slow using how method=brentq

Solution

Scipy made changes to access Brent using Cython which can be compiled without the GIL and with prange, a Cython command that parallelizes long loops. Compiled code is already typically very fast, so this could result in 10-100x speed improvements. Here's a gist with a Jupyter notebook benchmark that demonstrates about 30x reduction in runtime, but YMMV.

Alternatives

1. An alternative to Brent is Newton which is already vectorized in Scipy using NumPy arrays. A Newton how-method is already available in singlediode but Newton is unbounded, and so it's not guaranteed to converge like Brent, which is a bounded search method.
2. write a custom vectorized bisection method - actually there already is one, I think Rob Andrews may have authored it in pvlib.tools: _golden_sect_DataFrame, I think it's a golden search algorithm, not as fast as Brent, but still very good. Should we test it? Might need modification to work because I think it was specifically written for solving LambertW problems.
3. is there a numba way to solve this problem?
4. is there a pythran solution to this?
5. wait til next summer and make this a GSoC project 😉

Additional context
There might be some unpleasant hurdles:

• how do you maintain cython code in your repo?
• this should work fine for conda-forge & pvlib Anaconda channels, but what about PyPI users? We don't have a binary distribution on PyPI,
• The closest experience we've had is with the NREL SPA algorithm, which we've sort of abandoned because the NumPy & numba versions are so much faster and easier
• would we consider making binaries for deployment on multiple platforms like Windows (win32) & Mac OS X (darwin)
• would we make Linux users build their own wheels or figure out how to get involved with manylinux?

[kandersolar]

I am involved with another project that has a cython module. I am setting up cibuildwheel with GH Actions for it and am happy with it so far -- it seems to make the PyPI binary distribution process pretty painless, at least for the three major OSs, though I only have a few hours of experience with it at this point. If we start using cython in pvlib-python at some point, cibuildwheel seems like a good way to keep PyPI users from having to compile things themselves, and I'll volunteer to do the CI fiddling to get it set up.

[kandersolar]

SciPy is adding a vectorized bounded root finder that is accessible directly from python: https://github.com/scipy/scipy/blob/0cec3111f825019598a28b9a03858866d6908d39/scipy/optimize/_zeros_py.py#L1399

I installed a nightly scipy wheel (1.12.0.dev0+1691.26d4738) to compare it with brentq and newton in bishop88_i_from_v and found that it is faster than brentq by 1-2 orders of magnitude, and less than an order of magnitude away from newton:

Benchmark script

import pvlib
import numpy as np
import pandas as pd
import time
import matplotlib.pyplot as plt

params = {
    'photocurrent': 9,
    'saturation_current': 1e-5,
    'resistance_series': 0.1,
    'resistance_shunt': 1000,
    'nNsVth': 3,
}

results = []

for n in np.logspace(2, 5, 10).astype(int):
    v = np.linspace(0.1, 40, n)
    
    for method in ['newton', 'brentq', 'chandrupatla']:
        st = time.perf_counter()
        i = pvlib.singlediode.bishop88_i_from_v(v, method=method, **params)
        ed = time.perf_counter()
        results.append({
            'n': n,
            'method': method,
            'elapsed': ed - st
        })

out = pd.DataFrame(results)
out.pivot(index='n', columns='method', values='elapsed').plot(logy=True, logx=True)
plt.ylabel('Elapsed time [s]')
plt.xlabel('Number of curve points')
plt.grid()
plt.show()

bishop88_i_from_v patch

$ git diff
diff --git a/pvlib/singlediode.py b/pvlib/singlediode.py
index 203b20cd..cffc34ff 100644
--- a/pvlib/singlediode.py
+++ b/pvlib/singlediode.py
@@ -324,6 +324,17 @@ def bishop88_i_from_v(voltage, photocurrent, saturation_current,
         vd = newton(func=lambda x, *a: fv(x, voltage, *a), x0=x0,
                     fprime=lambda x, *a: bishop88(x, *a, gradients=True)[4],
                     args=args, **method_kwargs)
+    elif method == 'chandrupatla':
+        from scipy.optimize._zeros_py import _chandrupatla
+
+        voc_est = estimate_voc(photocurrent, saturation_current, nNsVth)
+        vlo = np.zeros_like(voltage)
+        vhi = np.full_like(voltage, voc_est)
+
+        result = _chandrupatla(func=lambda x, *a: fv(x, *a, *args),
+                               a=vlo, b=vhi, args=(voltage,),
+                               **method_kwargs)
+        vd = result.x
     else:
         raise NotImplementedError("Method '%s' isn't implemented" % method)

Not sure when a public function for this root finding method will appear in scipy's API, but anyway it is something to keep an eye on.

[cwhanse]

I would be quite pleased to drop our custom golden mean search function in favor of a vectorized scipy function that guarantees convergence (given that a root is in the bounded interval). Would solve issues in #1757

[adriesse]

This will be useful for reverse transposition as well.