Examples

C library

Here we show a selection of examples of using the lintegrate C library functions. A Linux Makefile for the three examples can be found here.

Example 1

In this example the lintegration_qag(), lintegration_qng() and lintegration_cquad() functions are used to evaluate the integral of a Gaussian function in log-space. The outputs are compared to the result of integration the original function using the GSL gsl_integration_qag() function.

/* example using lintegrate functionality */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include <lintegrate.h>

/* create function for integration */
double lintegrand(double x, void *params);

struct intparams {
  double mu;
  double sig;
};

double lintegrand(double x, void *params){
  struct intparams * p = (struct intparams *)params;
  double mu = p->mu;
  double sig = p->sig;

  return -0.5*(mu-x)*(mu-x)/(sig*sig);
}

double integrand(double x, void *params){
  struct intparams * p = (struct intparams *)params;
  double mu = p->mu;
  double sig = p->sig;

  return exp(-0.5*(mu-x)*(mu-x)/(sig*sig));
}

int main( int argv, char **argc ){
  gsl_function F;
  struct intparams params;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc (100);
  gsl_integration_cquad_workspace *cw = gsl_integration_cquad_workspace_alloc(50);
  double qaganswer = 0., qnganswer = 0., cquadanswer = 0., answer = 0.;
  double abserr = 0.;
  size_t neval = 0;

  double minlim = -6.; /* minimum for integration range */
  double maxlim = 6.;  /* maximum for integration range */

  double abstol = 1e-10; /* absolute tolerance */
  double reltol = 1e-10; /* relative tolerance */

  params.mu = 0.;
  params.sig = 1.;

  F.function = &lintegrand;
  F.params = &params;

  /* integrate log of function using QAG */
  lintegration_qag(&F, minlim, maxlim, abstol, reltol, 100, GSL_INTEG_GAUSS31, w, &qaganswer, &abserr);

  /* integrate log of function using QNG */
  lintegration_qng(&F, minlim, maxlim, abstol, reltol, &qnganswer, &abserr, &neval);

  /* integrate log of function using CQUAD */
  lintegration_cquad(&F, minlim, maxlim, abstol, reltol, cw, &cquadanswer, &abserr, &neval);

  /* integrate function using GSL QAG */
  F.function = &integrand;
  gsl_integration_qag(&F, minlim, maxlim, abstol, reltol, 100, GSL_INTEG_GAUSS31, w, &answer, &abserr);

  gsl_integration_workspace_free(w);
  gsl_integration_cquad_workspace_free(cw);

  fprintf(stdout, "Answer \"lintegrate QAG\" = %.8lf\n", qaganswer);
  fprintf(stdout, "Answer \"lintegrate QNG\" = %.8lf\n", qnganswer);
  fprintf(stdout, "Answer \"lintegrate CQUAD\" = %.8lf\n", cquadanswer);
  fprintf(stdout, "Answer \"gsl_integrate_qag\" = %.8lf\n", log(answer));
  fprintf(stdout, "Analytical answer = %.8lf\n", log(sqrt(2.*M_PI)));

  return 0;
}

Example 2

In this example the lintegration_qag(), lintegration_qng() and lintegration_cquad() functions are used to evaluate the integral of flat function in log-space. The outputs are compared to the result of integration the original function using the GSL gsl_integration_qag() function.

/* example using lintegrate functionality */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include <lintegrate.h>

/* create function for integration */
double lintegrand(double x, void *params);

double lintegrand(double x, void *params){
  return 50.;
}

double integrand(double x, void *params){
  return exp(50.);
}

int main( int argv, char **argc ){
  gsl_function F;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc (100);
  gsl_integration_cquad_workspace *cw = gsl_integration_cquad_workspace_alloc(50);
  double qaganswer = 0., qnganswer = 0., cquadanswer = 0., answer = 0.;
  double abserr = 0.;
  size_t neval = 0;

  double minlim = -6.; /* minimum for integration range */
  double maxlim = 6.;  /* maximum for integration range */

  double abstol = 1e-10; /* absolute tolerance */
  double reltol = 1e-10; /* relative tolerance */

  F.function = &lintegrand;

  /* integrate log of function using QAG */
  lintegration_qag(&F, minlim, maxlim, abstol, reltol, 100, GSL_INTEG_GAUSS31, w, &qaganswer, &abserr);

  /* integrate log of function using QNG */
  lintegration_qng(&F, minlim, maxlim, abstol, reltol, &qnganswer, &abserr, &neval);

  /* integrate log of function using CQUAD */
  lintegration_cquad(&F, minlim, maxlim, abstol, reltol, cw, &cquadanswer, &abserr, &neval);

  /* integrate function using GSL QAG */
  F.function = &integrand;
  gsl_integration_qag(&F, minlim, maxlim, abstol, reltol, 100, GSL_INTEG_GAUSS31, w, &answer, &abserr);

  gsl_integration_workspace_free(w);
  gsl_integration_cquad_workspace_free(cw);

  fprintf(stdout, "Answer \"lintegrate QAG\" = %.8lf\n", qaganswer);
  fprintf(stdout, "Answer \"lintegrate QNG\" = %.8lf\n", qnganswer);
  fprintf(stdout, "Answer \"lintegrate CQUAD\" = %.8lf\n", cquadanswer);
  fprintf(stdout, "Answer \"gsl_integrate_qag\" = %.8lf\n", log(answer));
  fprintf(stdout, "Analytic answer = %.8lf\n", log(maxlim-minlim) + 50.);

  return 0;
}

Example 3

In this example the lintegration_qag_split(), lintegration_qng_split() and lintegration_cquad_split() functions are used to evaluate the integral of function that is very tightly peaked log-space. The outputs are compared to an analytical calculation of the integral of original function found using Wolfram Alpha.

/* example using lintegrate_split functionality */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include <time.h>

#include <lintegrate.h>

/* create function for integration */
double lintegrand(double x, void *params);

double lintegrand(double x, void *params){
  return 1e-11/x;
}

int main( int argv, char **argc ){
  gsl_function F;

  double qaganswer = 0., cquadanswer = 0., qnganswer = 0.;
  double abserr = 0.;
  size_t neval = 0;

  double minlim = 1e-13; /* minimum for integration range */
  double maxlim = 1e-5;  /* maximum for integration range */

  /* split the interval logarithmically for integrating */
  size_t nints = 50, i = 0;
  double splits[nints+1];
  for ( i=0; i<nints+1; i++ ){
    splits[i] = pow(10., (log10(minlim) + i*(log10(maxlim)-log10(minlim))/(double)nints));
  }

  double abstol = 1e-10; /* absolute tolerance */
  double reltol = 1e-10; /* relative tolerance */

  F.function = &lintegrand;

  clock_t t1, t2;

  /* integrate log of function using QAG */
  t1 = clock();
  lintegration_qag_split(&F, splits, nints+1, abstol, reltol, 100, GSL_INTEG_GAUSS31, &qaganswer, &abserr);
  t2 = clock();
  fprintf(stderr, "lintegration_qag_split: run time = %ld mus\n", (t2-t1)*1000000/CLOCKS_PER_SEC);

  /* integrate log of function using QNG */
  t1 = clock();
  lintegration_qng_split(&F, splits, nints+1, abstol, reltol, &qnganswer, &abserr, &neval);
  t2 = clock();
  fprintf(stderr, "lintegration_qng_split: run time = %ld mus\n", (t2-t1)*1000000/CLOCKS_PER_SEC);

  /* integrate log of function using CQUAD */
  t1 = clock();
  lintegration_cquad_split(&F, splits, nints+1, abstol, reltol, 50, &cquadanswer, &abserr, &neval);
  t2 = clock();
  fprintf(stderr, "lintegration_cquad_split: run time = %ld mus\n", (t2-t1)*1000000/CLOCKS_PER_SEC);

  fprintf(stdout, "Answer \"lintegrate QAG\" = %.8lf\n", qaganswer);
  fprintf(stdout, "Answer \"lintegrate QNG\" = %.8lf\n", qnganswer);
  fprintf(stdout, "Answer \"lintegrate CQUAD\" = %.8lf\n", cquadanswer);
  fprintf(stdout, "(Approximate) Analytical answer = %.8lf\n", log(2.74356e28)); // via Wolfram Alpha

  return 0;
}

Python interface

Here we show and example of using the Python interface for lintegrate.

from lintegrate import lqag, lqng, lcquad, logtrapz
import numpy as np

# define the log of the function to be integrated
def integrand(x, args):
    mu, sig = args # unpack extra arguments
    return -0.5*((x-mu)/sig)**2

# set integration limits
xmin = -6.
xmax = 6.

# set additional arguments
mu = 0.
sig = 1.

resqag = lqag(integrand, xmin, xmax, args=(mu, sig))
resqng = lqng(integrand, xmin, xmax, args=(mu, sig))
rescquad = lcquad(integrand, xmin, xmax, args=(mu, sig))
restrapz = logtrapz(integrand, np.linspace(xmin, xmax, 100), args=(mu, sig))

Using R

The lintegrate Python interface can be accessed using R through the reticulate package. The above example would be:

library(reticulate)
py_install("lintegrate", pip = TRUE) ## run once to make sure lintegrate is installed and visible to reticulate.
lint <- import("lintegrate", convert =FALSE)
integrand <- function(x, args){
  mu = args[1]
  sig = args[2]
  return(-.5 * ((x-mu)/sig)^2 )
}
integrand <- Vectorize(integrand)
mu <- 0
sig <- 1
mmin <- -10
mmax <- 10
lint$lqag(py_func(integrand), r_to_py(mmin), r_to_py(mmax), c(mu, sig))