extern crate libc;
use libc::{c_char, c_double, c_int, c_longlong, c_void};
use std::ffi::CString;
use std::ptr;
use std::slice;
pub enum CubaVerbosity {
Silent = 0,
Progress = 1,
Input = 2,
Subregions = 3,
}
macro_rules! gen_setter {
($setr:ident, $r:ident, $t: ty) => {
pub fn $setr(&mut self, $r: $t) -> &mut Self {
self.$r = $r;
self
}
};
}
#[link(name = "cuba")]
extern "C" {
fn cubacores(n: c_int, p: c_int);
fn llVegas(
ndim: c_int,
ncomp: c_int,
integrand: Option<IntegrandC>,
userdata: *mut c_void,
nvec: c_longlong,
epsrel: c_double,
epsabs: c_double,
flags: c_int,
seed: c_int,
mineval: c_longlong,
maxeval: c_longlong,
nstart: c_longlong,
nincrease: c_longlong,
nbatch: c_longlong,
gridno: c_int,
statefile: *const c_char,
spin: *mut c_void,
neval: *mut c_longlong,
fail: *mut c_int,
integral: *mut c_double,
error: *mut c_double,
prob: *mut c_double,
);
fn llSuave(
ndim: c_int,
ncomp: c_int,
integrand: Option<IntegrandC>,
userdata: *mut c_void,
nvec: c_longlong,
epsrel: c_double,
epsabs: c_double,
flags: c_int,
seed: c_int,
mineval: c_longlong,
maxeval: c_longlong,
nnew: c_longlong,
nmin: c_longlong,
flatness: c_double,
statefile: *const c_char,
spin: *mut c_void,
nregions: *mut c_int,
neval: *mut c_longlong,
fail: *mut c_int,
integral: *mut c_double,
error: *mut c_double,
prob: *mut c_double,
);
fn llDivonne(
ndim: c_int,
ncomp: c_int,
integrand: Option<IntegrandC>,
userdata: *mut c_void,
nvec: c_longlong,
epsrel: c_double,
epsabs: c_double,
flags: c_int,
seed: c_int,
mineval: c_longlong,
maxeval: c_longlong,
key1: c_int,
key2: c_int,
key3: c_int,
maxpass: c_int,
border: c_double,
maxchisq: c_double,
mindeviation: c_double,
ngiven: c_longlong,
lxdgiven: c_int,
xgiven: *const c_double,
nextra: c_longlong,
peakfinder: Option<PeakfinderC>,
statefile: *const c_char,
spin: *mut c_void,
nregions: *mut c_int,
neval: *mut c_longlong,
fail: *mut c_int,
integral: *mut c_double,
error: *mut c_double,
prob: *mut c_double,
);
fn llCuhre(
ndim: c_int,
ncomp: c_int,
integrand: Option<IntegrandC>,
userdata: *mut c_void,
nvec: c_longlong,
epsrel: c_double,
epsabs: c_double,
flags: c_int,
mineval: c_longlong,
maxeval: c_longlong,
key: c_int,
statefile: *const c_char,
spin: *mut c_void,
nregions: *mut c_int,
neval: *mut c_longlong,
fail: *mut c_int,
integral: *mut c_double,
error: *mut c_double,
prob: *mut c_double,
);
}
type IntegrandC = extern "C" fn(
ndim: *const c_int,
x: *const c_double,
ncomp: *const c_int,
f: *mut c_double,
userdata: *mut c_void,
nvec: *const c_int,
core: *const c_int,
) -> c_int;
type PeakfinderC = extern "C" fn(
ndim: *const c_int,
b: *const c_double,
n: *mut c_int,
x: *mut c_double,
userdata: *mut c_void,
);
pub type Integrand<T> = fn(
x: &[f64],
f: &mut [f64],
user_data: &mut T,
nvec: usize,
core: i32,
) -> Result<(), &'static str>;
#[repr(C)]
struct CubaUserData<T> {
integrand: Integrand<T>,
user_data: T,
}
#[derive(Debug)]
pub struct CubaResult {
pub neval: i64,
pub fail: i32,
pub result: Vec<f64>,
pub error: Vec<f64>,
pub prob: Vec<f64>,
}
pub struct CubaIntegrator<T> {
integrand: Integrand<T>,
mineval: i64,
maxeval: i64,
nstart: i64,
nincrease: i64,
epsrel: f64,
epsabs: f64,
batch: i64,
cores: usize,
max_points_per_core: usize,
seed: i32,
use_only_last_sample: bool,
save_state_file: String,
keep_state_file: bool,
reset_vegas_integrator: bool,
key: i32,
key1: i32,
key2: i32,
key3: i32,
maxpass: i32,
border: f64,
maxchisq: f64,
mindeviation: f64,
}
impl<T> CubaIntegrator<T> {
pub fn new(integrand: Integrand<T>) -> CubaIntegrator<T> {
CubaIntegrator {
integrand,
mineval: 0,
maxeval: 50000,
nstart: 1000,
nincrease: 500,
epsrel: 0.001,
epsabs: 1e-12,
batch: 1000,
cores: 1,
max_points_per_core: 1000,
seed: 0,
use_only_last_sample: false,
save_state_file: String::new(),
keep_state_file: false,
reset_vegas_integrator: false,
key: 0,
key1: 47,
key2: 1,
key3: 1,
maxpass: 5,
border: 0.,
maxchisq: 10.,
mindeviation: 0.25,
}
}
pub fn set_cores(&mut self, cores: usize, max_points_per_core: usize) -> &mut Self {
self.cores = cores;
self.max_points_per_core = max_points_per_core;
unsafe {
cubacores(cores as c_int, max_points_per_core as c_int);
}
self
}
gen_setter!(set_mineval, mineval, i64);
gen_setter!(set_maxeval, maxeval, i64);
gen_setter!(set_nstart, nstart, i64);
gen_setter!(set_nincrease, nincrease, i64);
gen_setter!(set_epsrel, epsrel, f64);
gen_setter!(set_epsabs, epsabs, f64);
gen_setter!(set_batch, batch, i64);
gen_setter!(set_seed, seed, i32);
gen_setter!(set_use_only_last_sample, use_only_last_sample, bool);
gen_setter!(set_save_state_file, save_state_file, String);
gen_setter!(set_keep_state_file, keep_state_file, bool);
gen_setter!(set_reset_vegas_integrator, reset_vegas_integrator, bool);
gen_setter!(set_key, key, i32);
gen_setter!(set_key1, key1, i32);
gen_setter!(set_key2, key2, i32);
gen_setter!(set_key3, key3, i32);
gen_setter!(set_maxpass, maxpass, i32);
gen_setter!(set_border, border, f64);
gen_setter!(set_maxchisq, maxchisq, f64);
gen_setter!(set_mindeviation, mindeviation, f64);
extern "C" fn c_integrand(
ndim: *const c_int,
x: *const c_double,
ncomp: *const c_int,
f: *mut c_double,
userdata: *mut c_void,
nvec: *const c_int,
core: *const c_int,
) -> c_int {
unsafe {
let k: &mut CubaUserData<T> = &mut *(userdata as *mut _);
match (k.integrand)(
&slice::from_raw_parts(x, *ndim as usize * *nvec as usize),
&mut slice::from_raw_parts_mut(f, *ncomp as usize * *nvec as usize),
&mut k.user_data,
*nvec as usize,
*core as i32,
) {
Ok(_) => 0,
Err(e) => {
println!("Error during integration: {}. Aborting.", e);
-999
}
}
}
}
extern "C" fn c_peakfinder(
ndim: *const c_int,
b: *const c_double,
n: *mut c_int,
x: *mut c_double,
userdata: *mut c_void,
) {
}
pub fn vegas(
&mut self,
ndim: usize,
ncomp: usize,
nvec: usize,
verbosity: CubaVerbosity,
gridno: i32,
user_data: T,
) -> CubaResult {
let mut out = CubaResult {
neval: 0,
fail: 0,
result: vec![0.; ncomp],
error: vec![0.; ncomp],
prob: vec![0.; ncomp],
};
assert!(
gridno >= -10 && gridno <= 10,
"Grid number needs to be between -10 and 10."
);
assert!(
nvec <= self.batch as usize && nvec <= self.max_points_per_core,
"nvec needs to be at most the vegas batch size or the max points per core"
);
let mut x = CubaUserData {
integrand: self.integrand,
user_data: user_data,
};
let user_data_ptr = &mut x as *mut _ as *mut c_void;
let mut cubaflags = verbosity as i32;
if self.use_only_last_sample {
cubaflags |= 0b100;
}
if self.keep_state_file {
cubaflags |= 0b10000;
}
if self.reset_vegas_integrator {
cubaflags |= 0b100000;
}
let c_str = CString::new(self.save_state_file.as_str()).expect("CString::new failed");
unsafe {
llVegas(
ndim as c_int,
ncomp as c_int,
Some(CubaIntegrator::<T>::c_integrand),
user_data_ptr,
nvec as c_longlong,
self.epsrel,
self.epsabs,
cubaflags as c_int,
self.seed,
self.mineval,
self.maxeval,
self.nstart,
self.nincrease,
self.batch,
gridno,
c_str.as_ptr(),
ptr::null_mut(),
&mut out.neval,
&mut out.fail,
&mut out.result[0],
&mut out.error[0],
&mut out.prob[0],
);
}
out
}
pub fn suave(
&mut self,
ndim: usize,
ncomp: usize,
nvec: usize,
nnew: usize,
nmin: usize,
flatness: f64,
verbosity: CubaVerbosity,
user_data: T,
) -> CubaResult {
let mut out = CubaResult {
neval: 0,
fail: 0,
result: vec![0.; ncomp],
error: vec![0.; ncomp],
prob: vec![0.; ncomp],
};
assert!(
nvec <= self.max_points_per_core && nvec <= nnew,
"nvec needs to be at most the max points per core and nnew"
);
let mut x = CubaUserData {
integrand: self.integrand,
user_data: user_data,
};
let user_data_ptr = &mut x as *mut _ as *mut c_void;
let mut cubaflags = verbosity as i32;
if self.use_only_last_sample {
cubaflags |= 0b100;
}
if self.keep_state_file {
cubaflags |= 0b10000;
}
if self.reset_vegas_integrator {
cubaflags |= 0b100000;
}
let mut nregions = 0;
let c_str = CString::new(self.save_state_file.as_str()).expect("CString::new failed");
unsafe {
llSuave(
ndim as c_int,
ncomp as c_int,
Some(CubaIntegrator::<T>::c_integrand),
user_data_ptr,
nvec as c_longlong,
self.epsrel,
self.epsabs,
cubaflags as c_int,
self.seed,
self.mineval,
self.maxeval,
nnew as c_longlong,
nmin as c_longlong,
flatness,
c_str.as_ptr(),
ptr::null_mut(),
&mut nregions,
&mut out.neval,
&mut out.fail,
&mut out.result[0],
&mut out.error[0],
&mut out.prob[0],
);
}
out
}
pub fn divonne(
&mut self,
ndim: usize,
ncomp: usize,
nvec: usize,
xgiven: &[f64],
verbosity: CubaVerbosity,
user_data: T,
) -> CubaResult {
let mut out = CubaResult {
neval: 0,
fail: 0,
result: vec![0.; ncomp],
error: vec![0.; ncomp],
prob: vec![0.; ncomp],
};
assert!(
nvec <= self.max_points_per_core,
"nvec needs to be at most the max points per core"
);
let mut x = CubaUserData {
integrand: self.integrand,
user_data: user_data,
};
let user_data_ptr = &mut x as *mut _ as *mut c_void;
let mut cubaflags = verbosity as i32;
if self.use_only_last_sample {
cubaflags |= 0b100;
}
if self.keep_state_file {
cubaflags |= 0b10000;
}
if self.reset_vegas_integrator {
cubaflags |= 0b100000;
}
let mut nregions = 0;
let c_str = CString::new(self.save_state_file.as_str()).expect("CString::new failed");
unsafe {
llDivonne(
ndim as c_int,
ncomp as c_int,
Some(CubaIntegrator::<T>::c_integrand),
user_data_ptr,
nvec as c_longlong,
self.epsrel,
self.epsabs,
cubaflags as c_int,
self.seed,
self.mineval,
self.maxeval,
self.key1,
self.key2,
self.key3,
self.maxpass,
self.border,
self.maxchisq,
self.mindeviation,
(xgiven.len() / ndim) as c_longlong,
ndim as c_int,
if xgiven.len() == 0 {
ptr::null_mut()
} else {
&xgiven[0]
},
0,
None,
c_str.as_ptr(),
ptr::null_mut(),
&mut nregions,
&mut out.neval,
&mut out.fail,
&mut out.result[0],
&mut out.error[0],
&mut out.prob[0],
);
}
out
}
pub fn cuhre(
&mut self,
ndim: usize,
ncomp: usize,
nvec: usize,
verbosity: CubaVerbosity,
user_data: T,
) -> CubaResult {
let mut out = CubaResult {
neval: 0,
fail: 0,
result: vec![0.; ncomp],
error: vec![0.; ncomp],
prob: vec![0.; ncomp],
};
assert!(nvec <= 32, "nvec needs to be at most 32");
let mut x = CubaUserData {
integrand: self.integrand,
user_data: user_data,
};
let user_data_ptr = &mut x as *mut _ as *mut c_void;
let mut cubaflags = verbosity as i32;
if self.use_only_last_sample {
cubaflags |= 0b100;
}
if self.keep_state_file {
cubaflags |= 0b10000;
}
let mut nregions = 0;
let c_str = CString::new(self.save_state_file.as_str()).expect("CString::new failed");
unsafe {
llCuhre(
ndim as c_int,
ncomp as c_int,
Some(CubaIntegrator::<T>::c_integrand),
user_data_ptr,
nvec as c_longlong,
self.epsrel,
self.epsabs,
cubaflags as c_int,
self.mineval,
self.maxeval,
self.key,
c_str.as_ptr(),
ptr::null_mut(),
&mut nregions,
&mut out.neval,
&mut out.fail,
&mut out.result[0],
&mut out.error[0],
&mut out.prob[0],
);
}
out
}
}