Integrate over multiple limits at once using Quadpy

Question

My goal is to integrate over multiple limits at once using quadpy. The documentation says the following:

quadpy is fully vectorized, so if you like to compute the integral of a function on many domains at once, you can provide them all in one integrate() call, e.g.,

However, its still unclear to me which quadpy-function I should use for this objective. My script looks as follows:

import quadpy
import numpy as np


class Kelly:
    def __init__(self):
        self.odds = 1.952
        self.kelly = 0.08961344537815132
        self.i = 0.001
        self.f = np.arange(0, 1 + self.i, self.i).flatten()
        self.c1 = 1
        self.c2 = 2
        self.k = 1.5

    def loss_function(self, p):
        p = p[:, 0]
        loss_function = np.where(p[:, None] < (self.f * self.odds - 1) + 1 / self.odds,
                                 (self.c1 + self.c2) * abs(self.f - self.kelly) ** self.k, 0)
        return loss_function

    def integrate(self):
        xmin = np.zeros(len(self.f))
        xmax = np.array([self.f * (self.odds - 1) + 1 / self.odds]).flatten()

        # vals, errors = quadpy.quad(self.loss_function, xmin, xmax)
        
        return vals, errors


kelly = Kelly()
vals, errors = kelly.integrate()
print(vals, errors)

This results in the following error:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Please advice

Answer 1

The function signature should be adjusted properly, and different libraries have different notion of "vectorization". What looks problematic if we consider quadpy:

• First evidence is the error trace, which shows confusion with many array coordinates

• Second evidence is your function definition which works only with rank-2 input arrays. This shall not be the case, from the documentation, the inputs should be of rank 1:

  Quadpy provides integration schemes for many different 1D, 2D, even nD domains.

In any case, it would be helpful to see what is the precise function definition and input/output shapes. Try to provide an example how to call your function with known inputs/outputs, e.g.:

kelly = Kelly()
test_input = np.array([0.1,0.1,0.1]) # rank 1
kelly.loss_function(test_input.reshape(-1,1)) # passes, but rank 2
kelly.loss_function(test_input) # error

---

I will demonstrate how to bridge the two versions. The use-case here is to integrate vectorized output, namely function values for different hyperparams (here: f). To accomplish this in quadpy, we need to pass a list of individual components.

The full code:

class Kelly:
    def __init__(self):
        self.odds = 1.952
        self.kelly = 0.08961344537815132
        self.i = 0.01
        #self.f = np.arange(0.01, 1 + self.i, self.i).flatten()
        self.f = np.array([0.1,0.2,0.3,0.4,0.5,0.9])
        self.c1 = 1
        self.c2 = 2
        self.k = 1.5

    def loss_function_cub(self, p):
        p = p[:, 0]
        loss_function = np.where(np.less(p[:, None], (self.f * self.odds - 1) + 1 / self.odds),
                                (self.c1 + self.c2) * abs(self.f - self.kelly) ** self.k, 
                                 0)
        return loss_function

    def loss_function(self, f):
        def loss_fn(p):
          loss_function = np.where(np.less(p, (f * self.odds - 1) + 1 / self.odds),
                                  (self.c1 + self.c2) * abs(f - self.kelly) ** self.k, 
                                    0)
          return loss_function
        return loss_fn

    def loss_function_quad(self, p):
        """ outputs a list of functions parametrized by self.f """
        loss_functions = [self.loss_function(f)(p) for f in self.f]
        return loss_functions

    def integrate_cub(self):
        xmin = np.zeros(len(self.f))
        #xmax = np.array([self.f * (self.odds - 1) + 1 / self.odds]).flatten()
        xmax = np.ones(len(self.f))
        vals, errors = cubature(func=self.loss_function_cub,
                                ndim=1,
                                fdim=len(self.f),
                                xmin=xmin,
                                xmax=xmax,
                                vectorized=True)
        return vals, errors

    def integrate_quad(self):
        xmin = 0.0
        xmax = 1.0
        vals, errors = quadpy.quad(
            self.loss_function_quad,
            xmin,
            xmax,
            limit=1000,
        )
        return vals, errors


kelly = Kelly()

vals, errors = kelly.integrate_cub()
print(vals)

vals, errors = kelly.integrate_quad()
print(vals)

This gives

[0.         0.         0.0283406  0.1520492  0.38510662 2.18856573]
[0.         0.         0.0283406  0.1520492  0.38511753 2.18856573]

NOTE: If you like to integrate a bunch of functions over multiple domains (not clear from the question), this is another use-case and possibly harder to implement in quadpy.

Answer 2

The error may be occurring because you are using the comparison operator < between two arrays:

p[:, None] < (self.f * self.odds - 1) + 1 / self.odds.

This comparison generates a boolean array with multiple elements, resulting in an ambiguous truth value.

To fix this issue, adjust the comparison in the loss_function method to ensure that it operates correctly on the arrays. You can use the np.less() function to perform element-wise comparisons

# UNTESTED
def loss_function(self, p):
        p = p[:, 0]
        loss_function = np.where(np.less(p[:, None], (self.f * self.odds - 1) + 1 / self.odds),
                                 (self.c1 + self.c2) * abs(self.f - self.kelly) ** self.k, 0)
        return loss_function