catima/integrator.h

/*
 *  Author: Andrej Prochazka
 *  Copyright(C) 2017
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Affero General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Affero General Public License for more details.

 *  You should have received a copy of the GNU Affero General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef INTEGRATOR_H
#define INTEGRATOR_H
#include "catima/build_config.h"
#include "gsl/gsl_integration.h"
#include <functional>
#include <array>
#ifdef USE_THREADS
#include <mutex>
#endif

namespace catima{
    
/// helper class to integrate functions using the GSL library
class IntegratorGSL{
	public:
	IntegratorGSL(bool adapt=true);
	~IntegratorGSL();

	double integrate(std::function<double(double)> f, double min, double  max, double precision=0.001);

	private:
    gsl_integration_workspace *w;
    bool adaptive;
	double error;
	double result;
	double min;
	double max;
	#ifdef USE_THREADS
	std::mutex integration_mutex;
	#endif
    };

// built in integrator 
template<int order>
class GaussLegendreIntegration{
public:
    template<typename F>
    double integrate(F f, double a, double b) const;
    template<typename F>
    double operator()(F f, double a, double b) const {return integrate(f, a, b);}
    double get_w(int i) const {return w[i];}
    double get_x(int i) const {return x[i];}
    int n() const {return order;}
    std::array<double,order> get_points(double a = -1.0, double b = 1.0)const;

public:
    static std::array<double,order/2> w;
    static std::array<double,order/2> x;
};

template<int order>
template<typename F>
double GaussLegendreIntegration<order>::integrate(F f, double a, double b) const{
    double res=0.0;
    double p = 0.5*(b-a);
    double q = 0.5*(b+a);

    for(int i=0;i<order/2;i++){
        res += w[i] * (f(p*x[i] + q) + f(-p*x[i] + q));
    }
    return p*res;
}

template<int order>
std::array<double,order> GaussLegendreIntegration<order>::get_points(double a,  double b)const{
    std::array<double,order> points;
    double p = 0.5*(b-a);
    double q = 0.5*(b+a);
    
    int num = (order/2);
    for(int i=0;i< num;i++){
        points[num-i-1] = -p*x[i] + q;
        points[num+i] = p*x[i] + q;
    }
    return points;
}

#ifdef GSL_INTEGRATION
using integrator_type = IntegratorGSL;
#else
using integrator_type = GaussLegendreIntegration<8>;
#endif

extern integrator_type integrator;
extern IntegratorGSL integratorGSL;
}

#endif