PyCuba

Cuba is a library for multidimensional numerical integration. The manual in Cuba is helpful as a reference. Parallelisation is currently not supported.

The algorithms available are

  • Cuhre

  • Divonne

  • Suave

  • Vegas

Make sure you installed everything correctly.

PyCuba Demo Program

Take a look at the demo

PyCuba API

The return value is always a dictionary:

{ 'neval': number of evaluations,
  'fail': 0 or 1,
  'comp': number of components, usually 1,
  'nregions': number of regions used,
  'results': [ # a list of results for each component
     {
       'integral': value of the integral,
       'error':  uncertainty,
       'prob': probability (see manual),
     },
     ...
  ]
}

Define your integrand, like so:

def Integrand(ndim, xx, ncomp, ff, userdata):
        # access the current parameters
        x,y,z = [xx[i] for i in range(ndim.contents.value)]
        # compute the result
        result = math.sin(x)*math.cos(y)*math.exp(z)
        # store the result (here only one component)
        ff[0] = result
        return 0

It will be called with xx in the interval from 0 to 1 (so scale to your borders in this function).

The pycuba.demo() in the pycuba source shows how to run all of the algorithms and access the results:

  import math

  def Integrand(ndim, xx, ncomp, ff, userdata):
    x,y,z = [xx[i] for i in range(ndim.contents.value)]
    result = math.sin(x)*math.cos(y)*math.exp(z)
    ff[0] = result
    return 0
    

  NDIM = 3
  NCOMP = 1

  NNEW = 1000
  NMIN = 2
  FLATNESS = 50.

  KEY1 = 47
  KEY2 = 1
  KEY3 = 1
  MAXPASS = 5
  BORDER = 0.
  MAXCHISQ = 10.
  MINDEVIATION = .25
  NGIVEN = 0
  LDXGIVEN = NDIM
  NEXTRA = 0
  MINEVAL = 0
  MAXEVAL = 50000


  KEY = 0
  
  def print_header(name):
    print('-------------------- %s test -------------------' % name)
  def print_results(name, results):
    keys = ['nregions', 'neval', 'fail']
    keys = list(filter(results.has_key, keys))
    text = ["%s %d" % (k, results[k]) for k in keys]
    print("%s RESULT:\t" % name.upper() + "\t".join(text))
    for comp in results['results']:
      print("%s RESULT:\t" % name.upper() + \
	"%(integral).8f +- %(error).8f\tp = %(prob).3f\n" % comp)
  
  from os import environ as env
  verbose = 2
  if 'CUBAVERBOSE' in env:
    verbose = int(env['CUBAVERBOSE'])
  
  print_header('Vegas')
  print_results('Vegas', Vegas(Integrand, NDIM, verbose=2))
  
  print_header('Suave')
  print_results('Suave', Suave(Integrand, NDIM, NNEW, FLATNESS, verbose=2 | 4))
  
  print_header('Divonne')
  print_results('Divonne', Divonne(Integrand, NDIM, 
    mineval=MINEVAL, maxeval=MAXEVAL,
    key1=KEY1, key2=KEY2, key3=KEY3, maxpass=MAXPASS,
    border=BORDER, maxchisq=MAXCHISQ, mindeviation=MINDEVIATION,
    ldxgiven=LDXGIVEN, verbose=2))
  
  print_header('Cuhre')
  print_results('Cuhre', Cuhre(Integrand, NDIM, key=KEY, verbose=2 | 4))

API documentation