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 = ¶ms;
/* 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))