Issue: Specifying Multiple Arrays As Limits

[jackrolph]

Issue:

When using variable limits in more than one dimension (i.e. for binning) for vectorised integration, the integral is incorrect compared to evaluating in a loop.

Example: Code to integrate a 2D gaussian in bins, using vectorised and loop based methods.

import numpy as np
from scipy.stats import norm
from itertools import product
import cubepy as cp

def binstobcs(bins):
    return (bins[:-1] + bins[1:]) / 2.

def Integrated2DNormal(bins_x, bins_y):
    # Two uncorrelated normal distributions
    norm_x = norm(loc=0, scale=1)
    norm_y = norm(loc=0, scale=2)
    dx = bins_x[1] - bins_x[0]
    dy = bins_y[1] - bins_y[0]
    N_bin_x = len(bins_x)
    N_bin_y = len(bins_y)

    # Bin centers
    x = binstobcs(bins_x)
    y = binstobcs(bins_y)

    # PDF function
    pdf = lambda x, y: np.exp(norm_x.logpdf(x) + norm_y.logpdf(y))

    # Collect all bounds for vectorized integration
    x_lows, x_his, y_lows, y_his = [], [], [], []
    for i_x, i_y in product(range(N_bin_x-1), range(N_bin_y-1)):
        bin_x_low, bin_x_hi = bins_x[i_x], bins_x[i_x+1]
        bin_y_low, bin_y_hi = bins_y[i_y], bins_y[i_y+1]
        
        x_lows.append(bin_x_low)
        x_his.append(bin_x_hi)
        y_lows.append(bin_y_low)
        y_his.append(bin_y_hi)

    # Convert bounds to arrays
    low = [x_lows, y_lows]
    hi = [x_his, y_his]

    # Vectorized integration
    vs_test, err = cp.integrate(pdf, low, hi)
    vs_test = np.array(vs_test).reshape(N_bin_x-1, N_bin_y-1) / (dx * dy)

    # Loop-based integration
    vs = []
    for lx, ly, hx, hy in zip(x_lows, y_lows, x_his, y_his):
        v, error = cp.integrate(pdf, [lx, ly], [hx, hy])
        vs.append(v)

    vs = np.array(vs).reshape(N_bin_x-1, N_bin_y-1) / (dx * dy)

    # Compare the two methods
    print("vs_test:")
    print(vs_test)
    print("vs:")
    print(vs)
    print("Ratio (vs_test / vs):")
    print(vs_test / vs)

    # Integrate over x and y using trapezoidal rule
    I = np.trapz(vs, x, axis=0)
    I = np.trapz(I, y)
    print("Total Integrated Value:", I)

# Define bins
bins_x = np.linspace(-10, 10, 11)
bins_y = np.linspace(-10, 10, 10)

# Call the function
Integrated2DNormal(bins_x, bins_y)

OUTPUT:

Ratio (vs_test / vs):
[[1.         1.         1.         1.         1.         1.
  1.         1.         1.        ]
 [1.         1.         1.         1.         1.         1.
  1.         1.         1.        ]
 [1.         1.         1.0090048  1.         1.00902335 1.
  1.0090048  1.         1.        ]
 [1.05784913 1.         1.21738983 1.03748608 0.99999999 1.03748608
  1.21738983 1.         1.05784913]
 [1.7152618  1.75537945 1.         1.8160293  1.71523276 1.8160293
  1.         1.75537945 1.7152618 ]
 [1.7152618  1.75537945 1.         1.8160293  1.71523276 1.8160293
  1.         1.75537945 1.7152618 ]
 [1.05784913 1.         1.21738983 1.03748608 0.99999999 1.03748608
  1.21738983 1.         1.05784913]
 [1.         1.         1.0090048  1.         1.00902335 1.
  1.0090048  1.         1.        ]
 [1.         1.         1.         1.         1.         1.
  1.         1.         1.        ]
 [1.         1.         1.         1.         1.         1.
  1.         1.         1.        ]]
Total Integrated Value: 0.9999493567999094

[Areustle]

Will try to refamiliarize myself with the codebase and investigate.