Stimfit
0.13.15
|
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