Shapes could not be broadcast together with remaped shapes

Hello Numba community.

I’m trying to use the guvectorize decorator to make gufuncs out of some of my functions in NumPy-Financial. Currently, what I have is as follows:

@numba.guvectorize("(),(n)->()")
def _npv_internal(r, values, res):
    acc = 0.0
    for t in range(values.shape[0]):
        acc += values[t] / ((1.0 + r) ** t)
    res[0] = acc

This is called as follows:

r = np.atleast_1d(rate)
v = np.atleast_2d(values)
res = np.empty(shape=(r.shape[0], v.shape[0]))
_npv_internal(r, v, res)

This passes two of my existing unit tests. However, it fails when I try adding a new test:

    def test_npv_broadcast_equals_for_loop(self):
        cashflows = [
            [-15000, 1500, 2500, 3500, 4500, 6000],
            [-25000, 1500, 2500, 3500, 4500, 6000],
            [-35000, 1500, 2500, 3500, 4500, 6000],
            [-45000, 1500, 2500, 3500, 4500, 6000],
        ]
        rates = [-0.05, 0.00, 0.05, 0.10, 0.15]

        res = np.empty((len(rates), len(cashflows)))
        for i, r in enumerate(rates):
            for j, cf in enumerate(cashflows):
                res[i, j] = npf.npv(r, cf)

        assert assert_allclose(npf.npv(rates, cashflows), res)

This gives the error: ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (5,)->(5) (4,6)->(4) (5,4)->(5,4) and requested shape () I’m not sure if there is a neat way of running this to vectorize my operations.

Could someone please help me?

I’m a little confused about what you’re trying to do. Your unit test has a nested loop over each item (scaler) in rates and each row (1D array) in cashflows. And the computation seems to result in a scalar again, yet your results array has a 2D shape, which doesn’t seem to match that pattern/signature.

The shape of cashflows (4,6) also doesn’t match the result shape (5,4).

Perhaps adding some size 1 dimensions can resolve it, like reshaping your inputs make sure the dimensions align correctly:

r = r[:, np.newaxis]     # (5,1)   (l,m)
v = v[np.newaxis, :, :]  # (1,4,6) (l,m,n)

_npv_internal(r, v, res, axis=-1) # res = (5,4) (l,m)

But as said, I’m not entirely sure what output values you expect to get.

Sorry for being unclear. What I’m trying to do is to broadcast (I think I’m using that term correctly) the net present value computation. So for n results and k cashflows (where each cashflow is of a given length) we return an n x k array of results.

I have managed to get the tests to pass by simply using jit, but I would prefer to have it working as a ufunc to understand the concept and to apply it later.

What I have at the moment that passes is:


    r = np.atleast_1d(rate)
    v = np.atleast_2d(values)
    out = np.empty(shape=(r.shape[0], v.shape[0]))

    for i in prange(r.shape[0]):
        for j in prange(v.shape[0]):
            acc = 0.0
            for t in range(v.shape[1]):
                acc += v[j, t] / ((1.0 + r[i]) ** t)
            out[i, j] = acc

Essentially, what I would like to do is tell Numba how to do it in the simplest case (a scalar rate and 1d array of cashflows) and let the decorator figure out how to scale that to the more complicated case (a 1d array of rates and a 2d array of cashflows)

It’s still not clear to me. Changing the decorator from guvectorize to jit of the function in your OP should still result in a similar shape error.

What are rate and values in the code that you share? I assumed it would map to rates and cashflows from your OP. If that’s the case, the reshaping of the inputs I shared would give the exact same result right?

Keeping your original function as it is, but reshaping the inputs:

image

1 Like

Firstly, thank you! Reshaping the input solves the problem.

I’ll try to clarify; what I’m trying to do here is calculate the net present value (NPV) of a cashflow. NPV encapsulates how much an investment is worth in today’s terms by discounting each term in the cashflow by a given rate.

But how does the NPV change when we have different rates? Or what if we have multiple cashflows we want to analyse with the same rate? Or a combination of the two previous cases? I thought that, in this case, I could tell numba how to solve the simplest case, then use the guvectorise decorator to expand to more complicated cases.

1 Like

Hey @Kai ,

Risk disclaimer, I am a Numba user not a developer.
In your described case you could make guvectorize work if it makes sense in your broader context.

You should keep some considerations in mind while adapting guvectorize functions to your specific use case.
If not already in use Numba will be an additional dependency for your package.
Numpy and Numba’s guvectorize do not support keyword arguments in the function signature.
There is currently a limitation where using guvectorize inside a JIT-compiled function is not supported. However, there is a pull request (#8984) that aims to address this issue, and future versions of Numba might provide support for this.
To the code:
The function signature provided to guvectorize should treat the result as a 0D-array for your npv-function, which is not intuitive to me.

["void(float64, float64[:], float64[:])"], "(),(n)->()"

When using the function with guvectorize, it’s important to ensure that the difference in dimensions remains the same. As @Rutger described it is not enough to make sure that the used arrays have a minimum dimension like:

r = np.atleast_1d(rate)
v = np.atleast_2d(values)

This is why the example below works for different scenarios with rates and cashflows, but the function signature and dimensions must be consistent.
To make sure the correct input is provided to the function, you could write an additional wrapper.
Here is an example:

import numpy as np
import numba as nb

@nb.guvectorize(["void(float64, float64[:], float64[:])"], "(),(n)->()")
def _npv_internal(r: float,
                  values: np.ndarray[np.float64],
                  res: float | None = None):
    acc = np.float64(0.0)
    for t in range(len(values)):
        acc += values[t] / ((1.0 + r) ** t)
    res[0] = acc

print(f'Numba version {nb.__version__}\n')
# Numba version 0.58.1

print('base case:')
rates_0D = 0.08
cashflows_1D = np.array([-40_000, 5_000, 8_000, 12_000, 30_000], dtype=np.float64)
npv = _npv_internal(rates_0D, cashflows_1D)
print('rates:\n', rates_0D)
print('cashflows:\n', cashflows_1D)
print('npv:\n', npv)
# base case:
# rates:
#  0.08
# cashflows:
#  [-40000.   5000.   8000.  12000.  30000.]
# npv:
#  3065.2226681795255

print('\none dimension higher than base case (rates over cashflows)')
rates_1D = np.arange(0.08, 0.11, 0.01)
shape = (rates_1D.size, cashflows_1D.size)
cashflows_2D = np.broadcast_to(cashflows_1D, shape)
npv = _npv_internal(rates_1D, cashflows_2D)
print('rates:\n', rates_1D)
print('cashflows:\n', cashflows_2D)
print('npv:\n', npv)
# one dimension higher than base case (rates over cashflows)
# rates:
#  [0.08 0.09 0.1 ]
# cashflows:
#  [[-40000.   5000.   8000.  12000.  30000.]
#  [-40000.   5000.   8000.  12000.  30000.]
#  [-40000.   5000.   8000.  12000.  30000.]]
# npv:
#  [3065.22266818 1839.55400212  663.20606516]

print('\none dimension higher than base case (cashflows over rate)')
rates_1D = np.repeat(rates_0D, 3)
shape = (rates_1D.size, cashflows_1D.size)
cashflows_2D = np.broadcast_to(cashflows_1D, shape) + np.random.rand(*shape) * 10
npv = _npv_internal(rates_1D, cashflows_2D)
print('rates:\n', rates_1D)
print('cashflows:\n', cashflows_2D)
print('npv:\n', npv)
# one dimension higher than base case (cashflows over rate)
# rates:
#  [0.08 0.08 0.08]
# cashflows:
#  [[-39991.38776559   5009.77660016   8000.7422399   12000.2466202
#    30005.12153158]
#  [-39994.35224078   5004.44263713   8005.66121354  12003.29206549
#    30002.4836989 ]
#  [-39992.06003294   5007.14147516   8005.49237424  12001.97664394
#    30007.8053703 ]]
# npv:
#  [3087.4839149  3084.276499   3091.79024176]
2 Likes