stfnum.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 _STFNUM_H
#define _STFNUM_H
 
#ifdef _WINDOWS
#pragma warning( disable : 4251 )  // Disable warning messages
#endif

#if defined(_WINDOWS) && !defined(__MINGW32__)
    #ifdef STFNUMDLL
        #define StfnumDll __declspec( dllexport )
    #else
        #define StfnumDll __declspec( dllimport )
    #endif
#else
    #define StfnumDll
#endif
 
#include <vector>
#include <complex>
#include <deque>
 
#if (__cplusplus < 201103)
#  include <boost/function.hpp>
#else
#  include <algorithm>
#  include <cassert>
#  include <functional>
#endif
 
#ifdef _OPENMP
#include <omp.h>
#endif
#include <fftw3.h>
 
#ifdef _MSC_VER
#ifndef INFINITY
#define INFINITY (DBL_MAX+DBL_MAX)
#endif
#ifndef NAN
        static const unsigned long __nan[2] = {0xffffffff, 0x7fffffff};
        #define NAN (*(const float *) __nan)
#endif
#endif
 
#include "../libstfio/stfio.h"
#include "./spline.h"
 
namespace stfnum {



#if (__cplusplus < 201103)
 
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;
 
#else
 
typedef std::function<double(double, const Vector_double&)> Func;

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

typedef std::function<double(double, double, double, double, double)> Scale;
 
#endif
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(constrained_), 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; 
};


class StfnumDll Table {
public:

    Table(std::size_t nRows,std::size_t nCols);


    Table(const std::map< std::string, double >& map);


    double at(std::size_t row,std::size_t col) const;


    double& at(std::size_t row,std::size_t col);
    

    bool IsEmpty(std::size_t row,std::size_t col) const;


    void SetEmpty(std::size_t row,std::size_t col,bool value=true);


    void SetRowLabel(std::size_t row,const std::string& label);


    void SetColLabel(std::size_t col,const std::string& label);


    const std::string& GetRowLabel(std::size_t row) const;


    const std::string& GetColLabel(std::size_t col) const;


    std::size_t nRows() const { return rowLabels.size(); }


    std::size_t nCols() const { return colLabels.size(); }
    

    void AppendRows(std::size_t nRows);
 
private:
    // row major order:
    std::vector< std::vector<double> > values;
    std::vector< std::deque< bool > > empty;
    std::vector< std::string > rowLabels;
    std::vector< std::string > colLabels;
};
 
#if (__cplusplus < 201103)
typedef boost::function<Table(const Vector_double&,const std::vector<stfnum::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;
#else
typedef std::function<Table(const Vector_double&,const std::vector<stfnum::parInfo>,double)> Output;

Table defaultOutput(const Vector_double& pars,
                    const std::vector<parInfo>& parsInfo,
                    double chisqr);

typedef std::function<void(const Vector_double&, double, double, double, double, double, Vector_double&)> Init;
#endif


struct StfnumDll 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);


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


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


StfnumDll 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);


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


StfnumDll
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
);


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


StfnumDll Vector_double
detectionCriterion(
        const Vector_double& data,
        const Vector_double& templ,
        stfio::ProgressInfo& progDlg
);
 
// TODO: Add negative-going peaks.

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


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


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


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


double fbessel(double x, int n);


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


int fac(int arg);


int pow2(int arg);
 
enum direction {
    up,                 
    down,               
    both,               
    undefined_direction 
};

enum baseline_method {
    mean_sd   = 0,  
    median_iqr = 1   
};
 
 
}
 
typedef std::vector< stfnum::storedFunc >::const_iterator c_stfunc_it; 
 
inline int stfnum::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>
stfnum::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(stfnum::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)stfnum::spline_cubic_val(x,x_i[n_i],y_d,y_i,yp,ypp);
    }
    return y_if;
}
 
template <class T>
std::vector<T> stfnum::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 stfnum::SQR(T a) {return a*a;}
 
#endif