use _lambertw_pvlib etc. (#2732)

* use _lambertw_pvlib etc.

* whatsnew, lint

* use _lambertw_pvlib in ivtools

* Update pvlib/singlediode.py

* avoid warning

* pacify flake8

* another instance replaced

Author: cwhanse

docs/sphinx/source/whatsnew/v0.15.1.rst

@@ -36,6 +36,8 @@ Enhancements
 * Include `ross` and `faiman_rad` in the allowed models within 
   :py:meth:`pvlib.pvsystem.PVSystem.get_cell_temperature` (:issue:`2625`, :pull:`2631`)  
 * Accelerate the internals of :py:func:`~pvlib.solarpostion.ephemeris`. (:pull:`2626`)
+* Accelerate the intervals of :py:func:`~pvlib.pvsystem.singlediode` when
+  `method='lambertw'`. (:pull:`2732`, :pull:`2723`)
 
 Documentation
 ~~~~~~~~~~~~~

pvlib/ivtools/sdm/_fit_desoto_pvsyst_sandia.py

@@ -5,10 +5,9 @@
 import numpy as np
 
 from scipy import optimize
-from scipy.special import lambertw
 
 from pvlib.pvsystem import singlediode, v_from_i
-from pvlib.ivtools.utils import rectify_iv_curve, _numdiff
+from pvlib.ivtools.utils import rectify_iv_curve, _numdiff, _lambertw_pvlib
 from pvlib.pvsystem import _pvsyst_Rsh
 
 
@@ -535,7 +534,7 @@ def _calc_theta_phi_exact(vmp, imp, iph, io, rs, rsh, nnsvth):
             nnsvth == 0,
             np.nan,
             rsh * io / nnsvth * np.exp(rsh * (iph + io - imp) / nnsvth))
-        phi = np.where(argw > 0, lambertw(argw).real, np.nan)
+        phi = np.where(argw > 0, _lambertw_pvlib(argw), np.nan)
 
     # NaN where argw overflows. Switch to log space to evaluate
     u = np.isinf(argw)
@@ -561,7 +560,7 @@ def _calc_theta_phi_exact(vmp, imp, iph, io, rs, rsh, nnsvth):
             np.nan,
             rsh / (rsh + rs) * rs * io / nnsvth * np.exp(
                 rsh / (rsh + rs) * (rs * (iph + io) + vmp) / nnsvth))
-        theta = np.where(argw > 0, lambertw(argw).real, np.nan)
+        theta = np.where(argw > 0, _lambertw_pvlib(argw), np.nan)
 
     # NaN where argw overflows. Switch to log space to evaluate
     u = np.isinf(argw)

pvlib/ivtools/sdm/desoto.py

@@ -2,9 +2,8 @@
 
 from scipy import constants
 from scipy import optimize
-from scipy.special import lambertw
 
-from pvlib.ivtools.utils import rectify_iv_curve
+from pvlib.ivtools.utils import rectify_iv_curve, _lambertw_pvlib
 from pvlib.ivtools.sde import _fit_sandia_cocontent
 
 from pvlib.ivtools.sdm._fit_desoto_pvsyst_sandia import (
@@ -454,7 +453,7 @@ def fit_desoto_batzelis(v_mp, i_mp, v_oc, i_sc, alpha_sc, beta_voc):
     # Equation numbers refer to [1]
     t0 = 298.15  # K
     del0 = (1 - beta_voc * t0) / (50.1 - alpha_sc * t0)  # Eq 9
-    w0 = np.real(lambertw(np.exp(1/del0 + 1)))
+    w0 = _lambertw_pvlib(np.exp(1/del0 + 1))
 
     # Eqs 11-15
     a0 = del0 * v_oc

pvlib/pvarray.py

@@ -11,7 +11,8 @@
 import numpy as np
 import pandas as pd
 from scipy.optimize import curve_fit
-from scipy.special import exp10, lambertw
+from scipy.special import exp10
+from pvlib.ivtools.utils import _lambertw_pvlib
 
 
 def pvefficiency_adr(effective_irradiance, temp_cell,
@@ -483,7 +484,7 @@ def batzelis(effective_irradiance, temp_cell,
 
     # Eq 9-10
     del0 = (1 - beta_voc * t0) / (50.1 - alpha_sc * t0)
-    w0 = np.real(lambertw(np.exp(1/del0 + 1)))
+    w0 = _lambertw_pvlib(np.exp(1/del0 + 1))
 
     # Eqs 27-28
     alpha_imp = alpha_sc + (beta_voc - 1/t0) / (w0 - 1)

pvlib/singlediode.py

@@ -5,9 +5,10 @@
 import numpy as np
 import pandas as pd
 from pvlib.tools import _golden_sect_DataFrame
+from pvlib.ivtools.utils import _lambertw_pvlib, _log_lambertw
 
 from scipy.optimize import brentq, newton
-from scipy.special import lambertw
+
 
 # newton method default parameters for this module
 NEWTON_DEFAULT_PARAMS = {
@@ -801,12 +802,12 @@ def _lambertw_v_from_i(current, photocurrent, saturation_current,
         with np.errstate(over='ignore'):
             argW = I0 / (Gsh * a) * np.exp((-I + IL + I0) / (Gsh * a))
 
-        # lambertw typically returns complex value with zero imaginary part
-        # may overflow to np.inf
-        lambertwterm = lambertw(argW).real
+        lambertwterm = np.zeros_like(argW)
+
+        # Record indices where lambertw input overflowed
+        idx_inf = np.isinf(argW)
 
-        # Record indices where lambertw input overflowed output
-        idx_inf = np.isinf(lambertwterm)
+        lambertwterm[~idx_inf] = _lambertw_pvlib(argW[~idx_inf])
 
         # Only re-compute LambertW if it overflowed
         if np.any(idx_inf):
@@ -815,15 +816,7 @@ def _lambertw_v_from_i(current, photocurrent, saturation_current,
                        np.log(a[idx_inf]) +
                        (-I[idx_inf] + IL[idx_inf] + I0[idx_inf]) /
                        (Gsh[idx_inf] * a[idx_inf]))
-
-            # Three iterations of Newton-Raphson method to solve
-            #  w+log(w)=logargW. The initial guess is w=logargW. Where direct
-            #  evaluation (above) results in NaN from overflow, 3 iterations
-            #  of Newton's method gives approximately 8 digits of precision.
-            w = logargW
-            for _ in range(0, 3):
-                w = w * (1. - np.log(w) + logargW) / (1. + w)
-            lambertwterm[idx_inf] = w
+            lambertwterm[idx_inf] = _log_lambertw(logargW)
 
         # Eqn. 3 in Jain and Kapoor, 2004
         #  V = -I*(Rs + Rsh) + IL*Rsh - a*lambertwterm + I0*Rsh
@@ -866,27 +859,25 @@ def _lambertw_i_from_v(voltage, photocurrent, saturation_current,
     # Explicit solutions where Rs=0
     if np.any(idx_z):
         I[idx_z] = IL[idx_z] - I0[idx_z] * np.expm1(V[idx_z] / a[idx_z]) - \
-                   Gsh[idx_z] * V[idx_z]
+            Gsh[idx_z] * V[idx_z]
 
     # Only compute using LambertW if there are cases with Rs>0
     # Does NOT handle possibility of overflow, github issue 298
     if np.any(idx_p):
         # LambertW argument, cannot be float128, may overflow to np.inf
         argW = Rs[idx_p] * I0[idx_p] / (
-                    a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)) * \
-               np.exp((Rs[idx_p] * (IL[idx_p] + I0[idx_p]) + V[idx_p]) /
-                      (a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)))
+            a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.)) * (
+                np.exp((Rs[idx_p] * (IL[idx_p] + I0[idx_p]) + V[idx_p])
+                       / (a[idx_p] * (Rs[idx_p] * Gsh[idx_p] + 1.))))
 
-        # lambertw typically returns complex value with zero imaginary part
-        # may overflow to np.inf
-        lambertwterm = lambertw(argW).real
+        lambertwterm = _lambertw_pvlib(argW)
 
         # Eqn. 2 in Jain and Kapoor, 2004
         #  I = -V/(Rs + Rsh) - (a/Rs)*lambertwterm + Rsh*(IL + I0)/(Rs + Rsh)
         # Recast in terms of Gsh=1/Rsh for better numerical stability.
         I[idx_p] = (IL[idx_p] + I0[idx_p] - V[idx_p] * Gsh[idx_p]) / \
-                   (Rs[idx_p] * Gsh[idx_p] + 1.) - (
-                               a[idx_p] / Rs[idx_p]) * lambertwterm
+            (Rs[idx_p] * Gsh[idx_p] + 1.) \
+            - (a[idx_p] / Rs[idx_p]) * lambertwterm
 
     if output_is_scalar:
         return I.item()
@@ -1031,7 +1022,7 @@ def batzelis(photocurrent, saturation_current, resistance_series,
         voc = a * np.log(Iph / Is)
 
     # Eqs 5-8
-    w = np.real(lambertw(np.e * Iph / Is))
+    w = _lambertw_pvlib(np.e * Iph / Is)
     # vmp = (1 + Rs/Rsh) * a * (w - 1) - Rs * Iph * (1 - 1/w)  # not needed
     with np.errstate(divide='ignore', invalid='ignore'):  # zero Iph -> zero w
         imp = Iph * (1 - 1/w) - a * (w - 1) / Rsh