0

I am trying to replicate the technique descrived in the paragraph 1. of the selected answer to another question: How to pass additional parameters to numba cfunc passed as LowLevelCallable to scipy.integrate.quad.

However, I don't know how to modify the implementation so that xx[1] is an array of float and not a unique float.

PabloRdrRbl
  • 247
  • 2
  • 10
  • 1
    I'll try to have a look at this but in essence, xx[1] should be a unique float followed by the rest of the array in xx[2], xx[3], xx[4], etc. – Jacques Gaudin Jun 07 '21 at 20:52
  • I have posted the way I got this working. If you don't mind, you can test this and add it to your answer so we can close/delete my question. – PabloRdrRbl Jun 08 '21 at 07:46
  • 1
    I would keep this question/answer as it is. There's quite a bit to take in in the other answer so I'd rather keep it separate. – Jacques Gaudin Jun 08 '21 at 08:27
  • 1
    Or maybe take this answer and add it to the other question as a generalisation? – Jacques Gaudin Jun 08 '21 at 08:33

1 Answers1

2

I solved the issue by modifiying the code by Jacques Gaudin in https://stackoverflow.com/a/49732825/3925704 to:

import numpy as np
import scipy.integrate as si
import numba
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable


def jit_integrand_function(integrand_function):
    jitted_function = numba.jit(integrand_function, nopython=True)
    
    @cfunc(float64(intc, CPointer(float64)))
    def wrapped(n, xx):
        values = carray(xx, n)
        return jitted_function(values)
    return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def integrand(args):
    t = args[0]
    a = args[1]
    return np.exp(-t/a) / t**2

def do_integrate(func, a):
    """
    Integrate the given function from 1.0 to +inf with additional argument a.
    """
    return si.quad(func, 1, np.inf, args=(a,))
PabloRdrRbl
  • 247
  • 2
  • 10