ENH: NPV: Move npv to Cython

benchmarks/benchmarks.py

@@ -1,40 +1,46 @@
-import numpy as np
+from decimal import Decimal
 
+import numpy as np
 import numpy_financial as npf
 
 
-class Npv1DCashflow:
+class Npv2D:
 
-    param_names = ["cashflow_length"]
+    param_names = ["n_cashflows", "cashflow_lengths", "rates_lengths"]
     params = [
         (1, 10, 100, 1000),
+        (1, 10, 100, 1000),
+        (1, 10, 100, 1000),
     ]
 
     def __init__(self):
+        self.rates_decimal = None
+        self.rates = None
+        self.cashflows_decimal = None
         self.cashflows = None
 
-    def setup(self, cashflow_length):
+    def setup(self, n_cashflows, cashflow_lengths, rates_lengths):
         rng = np.random.default_rng(0)
-        self.cashflows = rng.standard_normal(cashflow_length)
-
-    def time_1d_cashflow(self, cashflow_length):
-        npf.npv(0.08, self.cashflows)
+        cf_shape = (n_cashflows, cashflow_lengths)
+        self.cashflows = rng.standard_normal(cf_shape)
+        self.rates = rng.standard_normal(rates_lengths)
+        self.cashflows_decimal = rng.standard_normal(cf_shape, dtype=Decimal)
+        self.rates_decimal = rng.standard_normal(rates_lengths, dtype=Decimal)
 
+    def time_broadcast(self, n_cashflows, cashflow_lengths, rates_lengths):
+        npf.npv(self.rates, self.cashflows)
 
-class Npv2DCashflows:
+    def time_for_loop(self, n_cashflows, cashflow_lengths, rates_lengths):
+        for i, rate in enumerate(self.rates):
+            for j, cashflow in enumerate(self.cashflows):
+                npf.npv(rate, cashflow)
 
-    param_names = ["n_cashflows", "cashflow_lengths"]
-    params = [
-        (1, 10, 100, 1000),
-        (1, 10, 100, 1000),
-    ]
+    def time_broadcast_decimal(self, n_cashflows, cashflow_lengths, rates_lengths):
+        npf.npv(self.rates_decimal, self.cashflows_decimal)
 
-    def __init__(self):
-        self.cashflows = None
+    def time_for_loop_decimal(self, n_cashflows, cashflow_lengths, rates_lengths):
+        for i, rate in enumerate(self.rates_decimal):
+            for j, cashflow in enumerate(self.cashflows_decimal):
+                npf.npv(rate, cashflow)
 
-    def setup(self, n_cashflows, cashflow_lengths):
-        rng = np.random.default_rng(0)
-        self.cashflows = rng.standard_normal((n_cashflows, cashflow_lengths))
 
-    def time_2d_cashflow(self, n_cashflows, cashflow_lengths):
-        npf.npv(0.08, self.cashflows)

numpy_financial/_cy_financial.pyx

@@ -0,0 +1,35 @@
+cimport cython
+from cython.parallel cimport  prange
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+@cython.cpow(True)
+cdef double npv_inner_loop(const double rate, const double[::1] cashflow) noexcept nogil:
+    cdef:
+        long cashflow_len = cashflow.shape[0]
+        long t
+        double acc
+
+    acc = 0.0
+    for t in range(cashflow_len):
+        acc += cashflow[t] / ((1.0 + rate) ** t)
+    return acc
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cpdef void cy_npv(
+        const double[::1] rates,
+        const double[:, ::1] cashflows,
+        double[:, ::1] out
+) noexcept nogil:
+    cdef:
+        long rate_len = rates.shape[0]
+        long no_of_cashflows = cashflows.shape[0]
+        long i, j
+
+    for i in prange(rate_len):
+        for j in prange(no_of_cashflows):
+            out[i, j] = npv_inner_loop(rates[i], cashflows[j])

numpy_financial/_financial.py

@@ -15,6 +15,9 @@
 
 import numpy as np
 
+from numpy_financial._cy_financial import cy_npv
+
+
 __all__ = ['fv', 'pmt', 'nper', 'ipmt', 'ppmt', 'pv', 'rate',
            'irr', 'npv', 'mirr',
            'NoRealSolutionError', 'IterationsExceededError']
@@ -46,6 +49,30 @@ def _convert_when(when):
         return [_when_to_num[x] for x in when]
 
 
+def _return_ufunc_like(array):
+    """Follow the ufunc convention of returning scalars for size 1 arrays"""
+    if array.size == 1:
+        return array.item()
+    return array
+
+
+def _make_out_array(*arrays):
+    """Make an ``out`` array
+
+    Output arrays have the following properties:
+
+    * Are of type decimal if any of the input arrays are object arrays
+    * Have shape of the first dimension of each input array
+    """
+    def _is_object_dtype(array):
+        return array.dtype == np.dtype("O")
+
+    shape = tuple(array.shape[0] for array in arrays)
+    if any(_is_object_dtype(array) for array in arrays):
+        return np.empty(shape, dtype=Decimal)
+    return np.empty(shape)
+
+
 def fv(rate, nper, pmt, pv, when='end'):
     """Compute the future value.
 
@@ -825,14 +852,29 @@ def irr(values, *, guess=None, tol=1e-12, maxiter=100, raise_exceptions=False):
     return np.nan
 
 
+def _npv_decimal(rates, cashflows, result):
+    r"""Version of the ``npv`` function supporting ``decimal.Decimal`` types
+
+    Warnings
+    --------
+    For internal use only, note that this function performs no error checking.
+    """
+    for i in range(rates.shape[0]):
+        for j in range(cashflows.shape[0]):
+            acc = Decimal("0.0")
+            for t in range(cashflows.shape[1]):
+                acc += cashflows[j, t] / ((Decimal("1.0") + rates[i]) ** t)
+            result[i, j] = acc
+
+
 def npv(rate, values):
     r"""Return the NPV (Net Present Value) of a cash flow series.
 
     Parameters
     ----------
-    rate : scalar
+    rate : scalar or array_like, shape(K, )
         The discount rate.
-    values : array_like, shape(M, )
+    values : array_like, shape(M, ) or shape(M, N)
         The values of the time series of cash flows.  The (fixed) time
         interval between cash flow "events" must be the same as that for
         which `rate` is given (i.e., if `rate` is per year, then precisely
@@ -843,7 +885,7 @@ def npv(rate, values):
 
     Returns
     -------
-    out : float
+    out : scalar or array_like, shape(K, M)
         The NPV of the input cash flow series `values` at the discount
         `rate`.
 
@@ -891,16 +933,34 @@ def npv(rate, values):
     >>> np.round(npf.npv(rate, cashflows) + initial_cashflow, 5)
     3065.22267
 
+    The NPV calculation may be applied to several ``rates`` and ``cashflows``
+    simulatneously. This produces an array of shape
+    ``(len(rates), len(cashflows))``.
+
+    >>> rates = np.array([0.00, 0.05, 0.10])
+    >>> cashflows = np.array([[-4_000, 500, 800], [-5_000, 600, 900]])
+    >>> npf.npv(rates, cashflows).round(2)
+    array([[-2700.  , -3500.  ],
+           [-2798.19, -3612.24],
+           [-2884.3 , -3710.74]])
+
     """
+    rates = np.atleast_1d(rate)
     values = np.atleast_2d(values)
-    timestep_array = np.arange(0, values.shape[1])
-    npv = (values / (1 + rate) ** timestep_array).sum(axis=1)
-    try:
-        # If size of array is one, return scalar
-        return npv.item()
-    except ValueError:
-        # Otherwise, return entire array
-        return npv
+
+    if rates.ndim != 1:
+        raise ValueError("invalid shape for rates. Rate must be either a scalar or 1d array")
+
+    if values.ndim != 2:
+        raise ValueError("invalid shape for values. Values must be either a 1d or 2d array")
+
+    out = _make_out_array(rates, values)
+    if out.dtype == np.dtype("O"):
+        _npv_decimal(rates, values, out)
+    else:
+        cy_npv(rates.astype(np.float64), values.astype(np.float64), out)
+
+    return _return_ufunc_like(out)
 
 
 def mirr(values, finance_rate, reinvest_rate, *, raise_exceptions=False):

tests/test_financial.py

@@ -164,7 +164,7 @@ def test_rate_maximum_iterations_exception_array(self):
 class TestNpv:
     def test_npv(self):
         assert_almost_equal(
-            npf.npv(0.05, [-15000, 1500, 2500, 3500, 4500, 6000]),
+            npf.npv(0.05, [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0]),
             122.89, 2)
 
     def test_npv_decimal(self):
@@ -174,17 +174,50 @@ def test_npv_decimal(self):
 
     def test_npv_broadcast(self):
         cashflows = [
-            [-15000, 1500, 2500, 3500, 4500, 6000],
-            [-15000, 1500, 2500, 3500, 4500, 6000],
-            [-15000, 1500, 2500, 3500, 4500, 6000],
-            [-15000, 1500, 2500, 3500, 4500, 6000],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
+            [-15000.0, 1500.0, 2500.0, 3500.0, 4500.0, 6000.0],
         ]
         expected_npvs = [
-            122.8948549, 122.8948549, 122.8948549, 122.8948549
+            [122.8948549, 122.8948549, 122.8948549, 122.8948549]
         ]
         actual_npvs = npf.npv(0.05, cashflows)
         assert_allclose(actual_npvs, expected_npvs)
 
+    @pytest.mark.parametrize("dtype", [Decimal, float])
+    def test_npv_broadcast_equals_for_loop(self, dtype):
+        cashflows_str = [
+            ["-15000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+            ["-25000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+            ["-35000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+            ["-45000.0", "1500.0", "2500.0", "3500.0", "4500.0", "6000.0"],
+        ]
+        rates_str = ["-0.05", "0.00", "0.05", "0.10", "0.15"]
+
+        cashflows = numpy.array([[dtype(x) for x in cf] for cf in cashflows_str])
+        rates = numpy.array([dtype(x) for x in rates_str])
+
+        expected = numpy.empty((len(rates), len(cashflows)), dtype=dtype)
+        for i, r in enumerate(rates):
+            for j, cf in enumerate(cashflows):
+                expected[i, j] = npf.npv(r, cf)
+
+        actual = npf.npv(rates, cashflows)
+        assert_equal(actual, expected)
+
+    @pytest.mark.parametrize("rates", ([[1, 2, 3]], numpy.empty(shape=(1,1,1))))
+    def test_invalid_rates_shape(self, rates):
+        cashflows = [1, 2, 3]
+        with pytest.raises(ValueError):
+            npf.npv(rates, cashflows)
+
+    @pytest.mark.parametrize("cf", ([[[1, 2, 3]]], numpy.empty(shape=(1, 1, 1))))
+    def test_invalid_cashflows_shape(self, cf):
+        rates = [1, 2, 3]
+        with pytest.raises(ValueError):
+            npf.npv(rates, cf)
+
 
 class TestPmt:
     def test_pmt_simple(self):
@@ -632,7 +665,7 @@ def test_npv_irr_congruence(self):
         # a series of cashflows to be zero, so we should have
         #
         # NPV(IRR(x), x) = 0.
-        cashflows = numpy.array([-40000, 5000, 8000, 12000, 30000])
+        cashflows = numpy.array([-40000.0, 5000.0, 8000.0, 12000.0, 30000.0])
         assert_allclose(
             npf.npv(npf.irr(cashflows), cashflows),
             0,