Merge pull request #118 from Kai-Striega/enh/nper/broadcast

ENH: nper: broadcast rework with Cython

Author: Kai-Striega

numpy_financial/_cfinancial.pyx

@@ -1,6 +1,59 @@
-from libc.math cimport NAN
+from libc.math cimport NAN, INFINITY, log
 cimport cython
 
+
+cdef double nper_inner_loop(
+    const double rate_,
+    const double pmt_,
+    const double pv_,
+    const double fv_,
+    const double when_
+) nogil:
+    if rate_ == 0.0 and pmt_ == 0.0:
+        return INFINITY
+
+    if rate_ == 0.0:
+        return -(fv_ + pv_) / pmt_
+
+    if rate_ <= -1.0:
+        return NAN
+
+    z = pmt_ * (1.0 + rate_ * when_) / rate_
+    return log((-fv_ + z) / (pv_ + z)) / log(1.0 + rate_)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def nper(
+    const double[::1] rates,
+    const double[::1] pmts,
+    const double[::1] pvs,
+    const double[::1] fvs,
+    const double[::1] whens,
+    double[:, :, :, :, ::1] out):
+
+    cdef:
+        Py_ssize_t rate_, pmt_, pv_, fv_, when_
+
+    for rate_ in range(rates.shape[0]):
+        for pmt_ in range(pmts.shape[0]):
+            for pv_ in range(pvs.shape[0]):
+                for fv_ in range(fvs.shape[0]):
+                    for when_ in range(whens.shape[0]):
+                        # We can have several ``ZeroDivisionErrors``s here
+                        # At the moment we want to replicate the existing function as
+                        # closely as possible however we should return financially
+                        # sensible results here.
+                        try:
+                            res = nper_inner_loop(
+                                rates[rate_], pmts[pmt_], pvs[pv_], fvs[fv_], whens[when_]
+                            )
+                        except ZeroDivisionError:
+                            res = NAN
+
+                        out[rate_, pmt_, pv_, fv_, when_] = res
+
+
 @cython.boundscheck(False)
 @cython.cdivision(True)
 def npv(const double[::1] rates, const double[:, ::1] values, double[:, ::1] out):

numpy_financial/_financial.py

@@ -306,35 +306,33 @@ def nper(rate, pmt, pv, fv=0, when='end'):
     The same analysis could be done with several different interest rates
     and/or payments and/or total amounts to produce an entire table.
 
-    >>> npf.nper(*(np.ogrid[0.07/12: 0.08/12: 0.01/12,
-    ...                     -150   : -99    : 50    ,
-    ...                     8000   : 9001   : 1000]))
-    array([[[ 64.07334877,  74.06368256],
-            [108.07548412, 127.99022654]],
+    >>> rates = [0.05, 0.06, 0.07]
+    >>> payments = [100, 200, 300]
+    >>> amounts = [7_000, 8_000, 9_000]
+    >>> npf.nper(rates, payments, amounts).round(3)
+    array([[[-30.827, -32.987, -34.94 ],
+            [-20.734, -22.517, -24.158],
+            [-15.847, -17.366, -18.78 ]],
     <BLANKLINE>
-           [[ 66.12443902,  76.87897353],
-            [114.70165583, 137.90124779]]])
+           [[-28.294, -30.168, -31.857],
+            [-19.417, -21.002, -22.453],
+            [-15.025, -16.398, -17.67 ]],
+    <BLANKLINE>
+           [[-26.234, -27.891, -29.381],
+            [-18.303, -19.731, -21.034],
+            [-14.311, -15.566, -16.722]]])
     """
     when = _convert_when(when)
-    rate, pmt, pv, fv, when = np.broadcast_arrays(rate, pmt, pv, fv, when)
-    nper_array = np.empty_like(rate, dtype=np.float64)
-
-    zero = rate == 0
-    nonzero = ~zero
-
-    with np.errstate(divide='ignore'):
-        # Infinite numbers of payments are okay, so ignore the
-        # potential divide by zero.
-        nper_array[zero] = -(fv[zero] + pv[zero]) / pmt[zero]
-
-    nonzero_rate = rate[nonzero]
-    z = pmt[nonzero] * (1 + nonzero_rate * when[nonzero]) / nonzero_rate
-    nper_array[nonzero] = (
-            np.log((-fv[nonzero] + z) / (pv[nonzero] + z))
-            / np.log(1 + nonzero_rate)
-    )
-
-    return nper_array
+    rates = np.atleast_1d(rate).astype(np.float64)
+    pmts = np.atleast_1d(pmt).astype(np.float64)
+    pvs = np.atleast_1d(pv).astype(np.float64)
+    fvs = np.atleast_1d(fv).astype(np.float64)
+    whens = np.atleast_1d(when).astype(np.float64)
+
+    out_shape = _get_output_array_shape(rates, pmts, pvs, fvs, whens)
+    out = np.empty(out_shape)
+    _cfinancial.nper(rates, pmts, pvs, fvs, whens, out)
+    return _ufunc_like(out)
 
 
 def _value_like(arr, value):

numpy_financial/tests/test_financial.py

@@ -45,12 +45,17 @@ def uint_dtype():
     cashflow_list_strategy,
 )
 
-short_scalar_array = npst.arrays(
+short_scalar_array_strategy = npst.arrays(
     dtype=real_scalar_dtypes,
     shape=npst.array_shapes(min_dims=0, max_dims=1, min_side=0, max_side=5),
 )
 
 
+when_strategy = st.sampled_from(
+    ['end', 'begin', 'e', 'b', 0, 1, 'beginning', 'start', 'finish']
+)
+
+
 def assert_decimal_close(actual, expected, tol=Decimal("1e-7")):
     # Check if both actual and expected are iterable (like arrays)
     if hasattr(actual, "__iter__") and hasattr(expected, "__iter__"):
@@ -280,7 +285,7 @@ def test_npv(self):
             rtol=1e-2,
         )
 
-    @given(rates=short_scalar_array, values=cashflow_array_strategy)
+    @given(rates=short_scalar_array_strategy, values=cashflow_array_strategy)
     @settings(deadline=None)
     def test_fuzz(self, rates, values):
         npf.npv(rates, values)
@@ -426,6 +431,16 @@ def test_broadcast(self):
             npf.nper(0.075, -2000, 0, 100000.0, [0, 1]), [21.5449442, 20.76156441], 4
         )
 
+    @given(
+        rates=short_scalar_array_strategy,
+        payments=short_scalar_array_strategy,
+        present_values=short_scalar_array_strategy,
+        future_values=short_scalar_array_strategy,
+        whens=when_strategy,
+    )
+    def test_fuzz(self, rates, payments, present_values, future_values, whens):
+        npf.nper(rates, payments, present_values, future_values, whens)
+
 
 class TestPpmt:
     def test_float(self):