giordano / Cuba.jl

Library for multidimensional numerical integration with four independent algorithms: Vegas, Suave, Divonne, and Cuhre.
https://giordano.github.io/Cuba.jl/stable
MIT License
75 stars 9 forks source link

Controlling Size of X Matrix using nvec #37

Closed CharlesRSmith44 closed 3 years ago

CharlesRSmith44 commented 3 years ago

Hello,

I'm trying to use the nvec option to choose a large matrix for sampling. However, when I set nvec = 10000 the integration function seems to generate a matrix that's only Mx 500 when it should be M x nvec. Is there anyway to ensure that the nvec option is implemented exactly? Why does it deviate from the selected nvec?

Code below:

import Pkg
Pkg.activate("joint_timing")
Pkg.instantiate()
using Cuba, Distributions
using BenchmarkTools, Test, CUDA
using FLoops, FoldsCUDA
using SpecialFunctions
using Suppressor

# User Inputs
M= 5 # number of independent uniform random variables
atol=1e-10
rtol=1e-10
nvec=100000
maxevals=1000000

## Setting up functions 
function beta_pdf_cpu(x,a,b,dim)
    prod(x.^(a-1.0f0) .* (1.0f0 .- x).^(b-1.0f0)./(gamma(a)*gamma(b)/gamma(a+b)),dims=dim)
end

function int_cpu(x, f)
display(size(x)) 
  f[1] = beta_pdf_cpu(x,1.0,2.0,1)[1]
end

function int_thread_col_cpu(x, f)
    display(size(x))
    @floop for i in 1:size(x,2)
         f[i] = beta_pdf_cpu(x[:,i],1.0,2.0,1)[1]
    end
end

#Integrating   

cuhre(int_cpu, M, 1, atol=atol, rtol=rtol,nvec = nvec)
suave(int_thread_col_cpu, M, 1, atol=atol, rtol=rtol, nvec = nvec, maxevals = maxevals)

## You can see here that the size of X is only 5x500 when it should be 5xNvec = 5 x 100000
giordano commented 3 years ago

Quoting from the upstream library documentation

Note that nvec indicates the actual number of points passed to the integrand here and may be smaller than the nvec given to the integrator.

There is nothing we can do from this wrapper.