stfmath.h

Go to the documentation of this file.

// Header file for the stimfit namespace
// General-purpose routines
// last revision: 08-08-2006
// C. Schmidt-Hieber, christsc@gmx.de

// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// 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 General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

#ifndef _CORE_H
#define _CORE_H

#ifdef _WINDOWS
#pragma warning( disable : 4251 )  // Disable warning messages
#endif

#include <vector>
#include <complex>
#include <deque>
#ifdef _OPENMP
#include <omp.h>
#endif

#ifdef _MSC_VER
#define INFINITY (DBL_MAX+DBL_MAX)
#define NAN (INFINITY-INFINITY)
#endif

#include "./../stf.h"

// header for the fourier transform:
#include "fftw3.h"
#include "./spline.h"

namespace stf {



typedef boost::function<double(double, const Vector_double&)> Func;

typedef boost::function<Vector_double(double, const Vector_double&)> Jac;

typedef boost::function<double(double, double, double, double, double)> Scale;

Vector_double nojac( double x, const Vector_double& p);

double noscale(double param, double xscale, double xoff, double yscale, double yoff);
 

struct parInfo {
    parInfo()
    : desc(""),toFit(true), constrained(false), constr_lb(0), constr_ub(0), scale(noscale), unscale(noscale) {}


  parInfo( const std::string& desc_, bool toFit_, bool constrained_ = false, 
             double constr_lb_ = 0, double constr_ub_ = 0, Scale scale_ = noscale, Scale unscale_ = noscale)
    : desc(desc_),toFit(toFit_),
        constrained(false), constr_lb(constr_lb_), constr_ub(constr_ub_),
        scale(scale_), unscale(unscale_)
    {}

    std::string desc; 
    bool toFit;    
    bool constrained; 
    double constr_lb; 
    double constr_ub; 
    Scale scale; 
    Scale unscale; 
};

typedef boost::function<Table(const Vector_double&,const std::vector<stf::parInfo>,double)> Output;
 
Table defaultOutput(const Vector_double& pars, 
                    const std::vector<parInfo>& parsInfo,
                    double chisqr);

typedef boost::function<void(const Vector_double&, double, double, double, double, double, Vector_double&)> Init;


struct StfDll storedFunc {


    storedFunc( const std::string& name_, const std::vector<parInfo>& pInfo_,
            const Func& func_, const Init& init_, const Jac& jac_, bool hasJac_ = true,
            const Output& output_ = defaultOutput /*,
            bool hasId_ = true*/
    ) : name(name_),pInfo(pInfo_),func(func_),init(init_),jac(jac_),hasJac(hasJac_),output(output_) /*, hasId(hasId_)*/
    {
/*        if (hasId) {
            id = NextId();
            std::string new_name;
            new_name << id << ": " << name;
            name = new_name;
        } else
            id = 0;
*/    }
     
    ~storedFunc() { }

//    static int n_funcs;          /*!< Static function counter */
//    int id;                      /*!< Function id; set automatically upon construction, so don't touch. */
    std::string name;            
    std::vector<parInfo> pInfo;  
    Func func;                   
    Init init;                   
    Jac jac;                     
    bool hasJac;                 
    Output output;               
//    bool hasId;                  /*!< Determines whether a function should have an id. */

};


template <typename T>
T SQR (T a);


Vector_double
filter(
        const Vector_double& toFilter,
        std::size_t filter_start,
        std::size_t filter_end,  
        const Vector_double &a,
        int SR,
        stf::Func func,
        bool inverse = false
);


std::map<double, int>
histogram(const Vector_double& data, int nbins=-1);


StfDll Vector_double
deconvolve(const Vector_double& data, const Vector_double& templ,
           int SR, double hipass, double lopass, stfio::ProgressInfo& progDlg);


template <class T>
std::vector<T>
cubicSpline(
        const std::vector<T>& y,
        T oldF,
        T newF
);

/* \param input The valarray to be differentiated.
 * \param x_scale The sampling interval.
 * \return The result of the differentiation.
 */
template <class T>
std::vector<T> diff(const std::vector<T>& input, T x_scale);


double integrate_simpson(
        const Vector_double& input,
        std::size_t a,
        std::size_t b,
        double x_scale
);


double integrate_trapezium(
        const Vector_double& input,
        std::size_t a,
        std::size_t b,
        double x_scale
);


int
linsolv(
        int m,
        int n,
        int nrhs,
        Vector_double& A,
        Vector_double& B
);


Vector_double
quad(const Vector_double& data, std::size_t begin, std::size_t end);
 


StfDll Vector_double
detectionCriterion(
        const Vector_double& data,
        const Vector_double& templ,
        stfio::ProgressInfo& progDlg
);

// TODO: Add negative-going peaks.

StfDll std::vector<int> peakIndices(const Vector_double& data, double threshold, int minDistance);


StfDll Vector_double linCorr(const Vector_double& va1, const Vector_double& va2, stfio::ProgressInfo& progDlg); 


double fgaussColqu(double x, const Vector_double& p);


double fboltz(double x, const Vector_double& p);


double fbessel(double x, int n);


double fbessel4(double x, const Vector_double& p);


int fac(int arg);


int pow2(int arg);

}

typedef std::vector< stf::storedFunc >::const_iterator c_stfunc_it; 
inline int stf::pow2(int arg) {return 1<<arg;}


template <typename T>
void SWAP(T s1, T s2) {
    T aux=s1;
    s1=s2;
    s2=aux;
}

template <class T>
std::vector<T>
stf::cubicSpline(const std::vector<T>& y,
        T oldF,
        T newF)
{
    double factor_i=newF/oldF;
    int size=(int)y.size();
    // size of interpolated data:
    int size_i=(int)(size*factor_i);
    Vector_double x(size);
    Vector_double y_d(size);
    for (int n_p=0; n_p < size; ++n_p) {
        x[n_p]=n_p;
        y_d[n_p]=y[n_p];
    }
    Vector_double y_i(stf::spline_cubic_set(x,y_d,0,0,0,0));

    std::vector<T> y_if(size_i);
    Vector_double x_i(size_i);

    //Cubic spline interpolation:
    for (int n_i=0; n_i < size_i; ++n_i) {
        x_i[n_i]=(double)n_i * (double)size/(double)size_i;
        double yp, ypp;
        y_if[n_i]=(T)stf::spline_cubic_val(x,x_i[n_i],y_d,y_i,yp,ypp);
    }
    return y_if;
}

template <class T>
std::vector<T> stf::diff(const std::vector<T>& input, T x_scale) {
    std::vector<T> diffVA(input.size()-1);
    for (unsigned n=0;n<diffVA.size();++n) {
        diffVA[n]=(input[n+1]-input[n])/x_scale;
    }
    return diffVA;
}

template <typename T>
inline T stf::SQR(T a) {return a*a;}

#endif