Added rounding of aa and b in _schumaker_qspline.py, see Issue #1311

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

@@ -0,0 +1,32 @@
+.. _whatsnew_0910:
+
+v0.9.1 (December 1, 2021)
+--------------------------
+
+Breaking changes
+~~~~~~~~~~~~~~~~
+
+Deprecations
+~~~~~~~~~~~~
+
+Enhancements
+~~~~~~~~~~~~
+
+Bug fixes
+~~~~~~~~~
+* Address round-off effects in :py:func:`pvlib.ivtools.utils._schumaker_qspline`
+  (:issue:`1311`, :pull:`1315`)
+
+Testing
+~~~~~~~
+
+Documentation
+~~~~~~~~~~~~~
+
+Requirements
+~~~~~~~~~~~~
+
+Contributors
+~~~~~~~~~~~~
+* Cliff Hansen (:ghuser:`cwhanse`)
+* :ghuser:`Antoine-0`

pvlib/ivtools/utils.py

@@ -9,7 +9,9 @@
 
 
 # A small number used to decide when a slope is equivalent to zero
-EPS = np.finfo('float').eps**(1/3)
+EPS_slope = np.finfo('float').eps**(1/3)
+# A small number used to decide when a slope is equivalent to zero
+EPS_val = np.finfo('float').eps
 
 
 def _numdiff(x, f):
@@ -263,7 +265,7 @@ def _schumaker_qspline(x, y):
     # [2], Algorithm 4.1 first 'if' condition of step 5 defines intervals
     # which won't get internal knots
     tests = s[:-1] + s[1:]
-    u = np.isclose(tests, 2.0 * delta, atol=EPS)
+    u = np.isclose(tests, 2.0 * delta, atol=EPS_slope)
     # u = true for an interval which will not get an internal knot
 
     k = n + sum(~u)  # total number of knots = original data + inserted knots
@@ -313,6 +315,11 @@ def _schumaker_qspline(x, y):
     aa = s[:-1] - delta
     b = s[1:] - delta
 
+    # Since the above two lines can lead to numerical errors, aa and b
+    # are rounded to 0.0 is their absolute value is small enough.
+    aa[np.isclose(aa, 0., atol=EPS_val)] = 0.
+    b[np.isclose(b, 0., atol=EPS_val)] = 0.
+
     sbar = np.zeros(k)
     eta = np.zeros(k)
     # will contain mapping from the left points of intervals containing an

pvlib/tests/ivtools/test_sde.py

@@ -53,14 +53,6 @@ def test_fit_sandia_simple_bad_iv(get_bad_iv_curves):
 
 
 @pytest.mark.parametrize('i,v,nsvth,expected',  [
-    (np.array(
-        [4., 3.95, 3.92, 3.9, 3.89, 3.88, 3.82, 3.8, 3.75, 3.7, 3.68, 3.66,
-         3.65, 3.5, 3.2, 2.7, 2.2, 1.3, .6, 0.]),
-     np.array(
-        [0., .2, .4, .6, .8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.2, 2.4, 2.6, 2.7,
-         2.76, 2.78, 2.81, 2.85, 2.88]),
-     2.,
-     (-96695.792, 96699.876, 7.4791, .0288, -.1413)),
     (np.array([3., 2.9, 2.8, 2.7, 2.6, 2.5, 2.4, 1.7, 0.8, 0.]),
      np.array([0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.45, 1.5]),
      10.,