Stimfit  0.13.15
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines
stfmath.h
Go to the documentation of this file.
00001 // Header file for the stimfit namespace
00002 // General-purpose routines
00003 // last revision: 08-08-2006
00004 // C. Schmidt-Hieber, christsc@gmx.de
00005 
00006 // This program is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU General Public License
00008 // as published by the Free Software Foundation; either version 2
00009 // of the License, or (at your option) any later version.
00010 
00011 // This program is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014 // GNU General Public License for more details.
00015 
00016 // You should have received a copy of the GNU General Public License
00017 // along with this program; if not, write to the Free Software
00018 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
00019 
00026 #ifndef _CORE_H
00027 #define _CORE_H
00028 
00029 #ifdef _WINDOWS
00030 #pragma warning( disable : 4251 )  // Disable warning messages
00031 #endif
00032 
00033 #include <vector>
00034 #include <complex>
00035 #include <deque>
00036 #ifdef _OPENMP
00037 #include <omp.h>
00038 #endif
00039 
00040 #ifdef _MSC_VER
00041 #define INFINITY (DBL_MAX+DBL_MAX)
00042 #define NAN (INFINITY-INFINITY)
00043 #endif
00044 
00045 #include "./../stf.h"
00046 
00047 // header for the fourier transform:
00048 #include "fftw3.h"
00049 #include "./spline.h"
00050 
00051 namespace stf {
00052 
00057 
00058 
00062 typedef boost::function<double(double, const Vector_double&)> Func;
00063 
00065 typedef boost::function<Vector_double(double, const Vector_double&)> Jac;
00066 
00068 typedef boost::function<double(double, double, double, double, double)> Scale;
00069 
00071 Vector_double nojac( double x, const Vector_double& p);
00072 
00074 double noscale(double param, double xscale, double xoff, double yscale, double yoff);
00075  
00077 
00082 struct parInfo {
00084     parInfo()
00085     : desc(""),toFit(true), constrained(false), constr_lb(0), constr_ub(0), scale(noscale), unscale(noscale) {}
00086 
00088 
00097   parInfo( const std::string& desc_, bool toFit_, bool constrained_ = false, 
00098              double constr_lb_ = 0, double constr_ub_ = 0, Scale scale_ = noscale, Scale unscale_ = noscale)
00099     : desc(desc_),toFit(toFit_),
00100         constrained(false), constr_lb(constr_lb_), constr_ub(constr_ub_),
00101         scale(scale_), unscale(unscale_)
00102     {}
00103 
00104     std::string desc; 
00105     bool toFit;    
00106     bool constrained; 
00107     double constr_lb; 
00108     double constr_ub; 
00109     Scale scale; 
00110     Scale unscale; 
00111 };
00112 
00114 typedef boost::function<Table(const Vector_double&,const std::vector<stf::parInfo>,double)> Output;
00115  
00117 Table defaultOutput(const Vector_double& pars, 
00118                     const std::vector<parInfo>& parsInfo,
00119                     double chisqr);
00120 
00122 typedef boost::function<void(const Vector_double&, double, double, double, double, double, Vector_double&)> Init;
00123 
00125 
00130 struct StfDll storedFunc {
00131 
00133 
00141     storedFunc( const std::string& name_, const std::vector<parInfo>& pInfo_,
00142             const Func& func_, const Init& init_, const Jac& jac_, bool hasJac_ = true,
00143             const Output& output_ = defaultOutput /*,
00144             bool hasId_ = true*/
00145     ) : name(name_),pInfo(pInfo_),func(func_),init(init_),jac(jac_),hasJac(hasJac_),output(output_) /*, hasId(hasId_)*/
00146     {
00147 /*        if (hasId) {
00148             id = NextId();
00149             std::string new_name;
00150             new_name << id << ": " << name;
00151             name = new_name;
00152         } else
00153             id = 0;
00154 */    }
00155      
00157     ~storedFunc() { }
00158 
00159 //    static int n_funcs;          /*!< Static function counter */
00160 //    int id;                      /*!< Function id; set automatically upon construction, so don't touch. */
00161     std::string name;            
00162     std::vector<parInfo> pInfo;  
00163     Func func;                   
00164     Init init;                   
00165     Jac jac;                     
00166     bool hasJac;                 
00167     Output output;               
00168 //    bool hasId;                  /*!< Determines whether a function should have an id. */
00169 
00170 };
00171 
00173 
00176 template <typename T>
00177 T SQR (T a);
00178 
00180 
00189 Vector_double
00190 filter(
00191         const Vector_double& toFilter,
00192         std::size_t filter_start,
00193         std::size_t filter_end,  
00194         const Vector_double &a,
00195         int SR,
00196         stf::Func func,
00197         bool inverse = false
00198 );
00199 
00201 
00205 std::map<double, int>
00206 histogram(const Vector_double& data, int nbins=-1);
00207 
00209 
00216 StfDll Vector_double
00217 deconvolve(const Vector_double& data, const Vector_double& templ,
00218            int SR, double hipass, double lopass, stfio::ProgressInfo& progDlg);
00219 
00221 
00226 template <class T>
00227 std::vector<T>
00228 cubicSpline(
00229         const std::vector<T>& y,
00230         T oldF,
00231         T newF
00232 );
00233 
00235 /* \param input The valarray to be differentiated.
00236  * \param x_scale The sampling interval.
00237  * \return The result of the differentiation.
00238  */
00239 template <class T>
00240 std::vector<T> diff(const std::vector<T>& input, T x_scale);
00241 
00243 
00249 double integrate_simpson(
00250         const Vector_double& input,
00251         std::size_t a,
00252         std::size_t b,
00253         double x_scale
00254 );
00255 
00257 
00263 double integrate_trapezium(
00264         const Vector_double& input,
00265         std::size_t a,
00266         std::size_t b,
00267         double x_scale
00268 );
00269 
00271 
00283 int
00284 linsolv(
00285         int m,
00286         int n,
00287         int nrhs,
00288         Vector_double& A,
00289         Vector_double& B
00290 );
00291 
00293 
00298 Vector_double
00299 quad(const Vector_double& data, std::size_t begin, std::size_t end);
00300  
00301 
00303 
00307 StfDll Vector_double
00308 detectionCriterion(
00309         const Vector_double& data,
00310         const Vector_double& templ,
00311         stfio::ProgressInfo& progDlg
00312 );
00313 
00314 // TODO: Add negative-going peaks.
00316 
00321 StfDll std::vector<int> peakIndices(const Vector_double& data, double threshold, int minDistance);
00322 
00324 
00328 StfDll Vector_double linCorr(const Vector_double& va1, const Vector_double& va2, stfio::ProgressInfo& progDlg); 
00329 
00331 
00339 double fgaussColqu(double x, const Vector_double& p);
00340 
00342 
00349 double fboltz(double x, const Vector_double& p);
00350 
00352 
00359 double fbessel(double x, int n);
00360 
00362 
00371 double fbessel4(double x, const Vector_double& p);
00372 
00374 
00377 int fac(int arg);
00378 
00380 
00383 int pow2(int arg);
00384 
00387 }
00388 
00389 typedef std::vector< stf::storedFunc >::const_iterator c_stfunc_it; 
00391 inline int stf::pow2(int arg) {return 1<<arg;}
00392 
00394 
00397 template <typename T>
00398 void SWAP(T s1, T s2) {
00399     T aux=s1;
00400     s1=s2;
00401     s2=aux;
00402 }
00403 
00404 template <class T>
00405 std::vector<T>
00406 stf::cubicSpline(const std::vector<T>& y,
00407         T oldF,
00408         T newF)
00409 {
00410     double factor_i=newF/oldF;
00411     int size=(int)y.size();
00412     // size of interpolated data:
00413     int size_i=(int)(size*factor_i);
00414     Vector_double x(size);
00415     Vector_double y_d(size);
00416     for (int n_p=0; n_p < size; ++n_p) {
00417         x[n_p]=n_p;
00418         y_d[n_p]=y[n_p];
00419     }
00420     Vector_double y_i(stf::spline_cubic_set(x,y_d,0,0,0,0));
00421 
00422     std::vector<T> y_if(size_i);
00423     Vector_double x_i(size_i);
00424 
00425     //Cubic spline interpolation:
00426     for (int n_i=0; n_i < size_i; ++n_i) {
00427         x_i[n_i]=(double)n_i * (double)size/(double)size_i;
00428         double yp, ypp;
00429         y_if[n_i]=(T)stf::spline_cubic_val(x,x_i[n_i],y_d,y_i,yp,ypp);
00430     }
00431     return y_if;
00432 }
00433 
00434 template <class T>
00435 std::vector<T> stf::diff(const std::vector<T>& input, T x_scale) {
00436     std::vector<T> diffVA(input.size()-1);
00437     for (unsigned n=0;n<diffVA.size();++n) {
00438         diffVA[n]=(input[n+1]-input[n])/x_scale;
00439     }
00440     return diffVA;
00441 }
00442 
00443 template <typename T>
00444 inline T stf::SQR(T a) {return a*a;}
00445 
00446 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines