Nicholaswogan / NumbaMinpack

Python wrapper of Minpack (root finding) which can be called from within numba functions.
MIT License
31 stars 2 forks source link

Can NumbaMinpack be used with Cuda? #3

Open jimkring opened 2 years ago

jimkring commented 2 years ago

Hi @Nicholaswogan.

Thanks for the great tools! I was looking at this stack overflow discussion where you posted an example of using NumbaMinpack for curve fitting.

I was wondering if NumbaMinpack (and this general approach you showed in that thread -- and copied below) can be made to work with Numba's cuda functionality. Do you know if that's possible/straightforward?

from NumbaQuadpack import quadpack_sig, dqags
from NumbaMinpack import minpack_sig, lmdif
import numpy as np
import numba as nb
import timeit
np.random.seed(0)

x = np.linspace(0,2*np.pi,100)
y = np.sin(x)+ np.random.rand(100)

@nb.jit
def fitfunction(x, A, B):
    return A*np.sin(B*x)

@nb.cfunc(minpack_sig)
def fitfunction_optimize(u_, fvec, args_):
    u = nb.carray(u_,(2,))
    args = nb.carray(args_,(200,))
    A, B = u
    x = args[:100]
    y = args[100:]
    for i in range(100):
        fvec[i] = fitfunction(x[i], A, B) - y[i] 
optimize_ptr = fitfunction_optimize.address

@nb.cfunc(quadpack_sig)
def fitfunction_integrate(x, data):
    A = data[0]
    B = data[1]
    return fitfunction(x, A, B)
integrate_ptr = fitfunction_integrate.address

@nb.njit
def fast_function():  
    try:
        neqs = 100
        u_init = np.array([2.0,.8],np.float64)
        args = np.append(x,y)
        fitparam, fvec, success, info = lmdif(optimize_ptr , u_init, neqs, args)
        if not success:
            raise Exception

        lower = 0.0
        uppers = np.linspace(np.pi,np.pi*2.0,200)
        solutions = np.empty(len(uppers))
        for i in range(len(uppers)):
            solutions[i], abserr, success = dqags(integrate_ptr, lower, uppers[i], data = fitparam)
            if not success:
                raise Exception
    except:
        print('doing something else')

fast_function()
iters = 1000
t_nb = timeit.Timer(fast_function).timeit(number=iters)/iters
print(t_nb)

Thanks.

Nicholaswogan commented 2 years ago

Thanks! Unfortunately, I do not think you can’t use numba minpack with cuda. These extensions are compiled for the CPU, so I don’t think there is any way it can run on the GPU. There would need to be a cuda GPU implementation of minpack, which I do not think exists.