1 /***************************************************************************
2 File : fitclasses.cpp
3 Project : SciDAVis
4 --------------------------------------------------------------------
5 Copyright : (C) 2006 by Ion Vasilief, Tilman Benkert
6 Email (use @ for *) : ion_vasilief*yahoo.fr, thzs*gmx.net
7 Description : Exponential fit classes
8
9 ***************************************************************************/
10
11 /***************************************************************************
12 * *
13 * This program is free software; you can redistribute it and/or modify *
14 * it under the terms of the GNU General Public License as published by *
15 * the Free Software Foundation; either version 2 of the License, or *
16 * (at your option) any later version. *
17 * *
18 * This program is distributed in the hope that it will be useful, *
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
21 * GNU General Public License for more details. *
22 * *
23 * You should have received a copy of the GNU General Public License *
24 * along with this program; if not, write to the Free Software *
25 * Foundation, Inc., 51 Franklin Street, Fifth Floor, *
26 * Boston, MA 02110-1301 USA *
27 * *
28 ***************************************************************************/
29 #include "ExponentialFit.h"
30 #include "fit_gsl.h"
31 using namespace std;
32 #include <assert.h>
33
34 /*****************************************************************************
35 *
36 * Class ExponentialFit
37 *
38 *****************************************************************************/
39
ExponentialFit(ApplicationWindow * parent,Graph * g,bool expGrowth)40 ExponentialFit::ExponentialFit(ApplicationWindow *parent, Graph *g, bool expGrowth)
41 : Fit(parent, g), is_exp_growth(expGrowth)
42 {
43 init();
44 }
45
ExponentialFit(ApplicationWindow * parent,Graph * g,const QString & curveTitle,bool expGrowth)46 ExponentialFit::ExponentialFit(ApplicationWindow *parent, Graph *g, const QString &curveTitle,
47 bool expGrowth)
48 : Fit(parent, g), is_exp_growth(expGrowth)
49 {
50 init();
51 setDataFromCurve(curveTitle);
52 }
53
ExponentialFit(ApplicationWindow * parent,Graph * g,const QString & curveTitle,double start,double end,bool expGrowth)54 ExponentialFit::ExponentialFit(ApplicationWindow *parent, Graph *g, const QString &curveTitle,
55 double start, double end, bool expGrowth)
56 : Fit(parent, g), is_exp_growth(expGrowth)
57 {
58 init();
59 setDataFromCurve(curveTitle, start, end);
60 }
61
init()62 void ExponentialFit::init()
63 {
64 d_f = exp_f;
65 d_df = exp_df;
66 d_fdf = exp_fdf;
67 d_fsimplex = exp_d;
68 d_p = 3;
69 d_min_points = d_p;
70 d_param_init = gsl_vector_alloc(d_p);
71 gsl_vector_set_all(d_param_init, 1.0);
72
73 covar = gsl_matrix_alloc(d_p, d_p);
74 d_results.resize(d_p);
75 d_param_names << "A"
76 << "t"
77 << "y0";
78
79 if (is_exp_growth) {
80 setObjectName("ExpGrowth");
81 d_explanation = tr("Exponential growth");
82 d_formula = "y0+A*exp(x/t)";
83 d_param_explain << tr("(amplitude)") << tr("(lifetime)") << tr("(offset)");
84 } else {
85 setObjectName("ExpDecay");
86 d_explanation = tr("Exponential decay");
87 d_formula = "y0+A*exp(-x/t)";
88 d_param_explain << tr("(amplitude)") << tr("(e-folding time)") << tr("(offset)");
89 }
90 }
91
storeCustomFitResults(const vector<double> & par)92 void ExponentialFit::storeCustomFitResults(const vector<double> &par)
93 {
94 d_results = par;
95 assert(d_results.size() >= 2);
96
97 if (is_exp_growth)
98 d_results[1] = -1.0 / d_results[1];
99 else
100 d_results[1] = 1.0 / d_results[1];
101 }
102
calculateFitCurveData(const vector<double> & par,std::vector<double> & X,std::vector<double> & Y)103 bool ExponentialFit::calculateFitCurveData(const vector<double> &par, std::vector<double> &X,
104 std::vector<double> &Y)
105 {
106 generateX(X);
107 for (int i = 0; i < d_points; i++)
108 Y[i] = par[0] * exp(-par[1] * X[i]) + par[2];
109 return true;
110 }
111
112 /*****************************************************************************
113 *
114 * Class TwoExpFit
115 *
116 *****************************************************************************/
117
TwoExpFit(ApplicationWindow * parent,Graph * g)118 TwoExpFit::TwoExpFit(ApplicationWindow *parent, Graph *g) : Fit(parent, g)
119 {
120 init();
121 }
122
TwoExpFit(ApplicationWindow * parent,Graph * g,const QString & curveTitle)123 TwoExpFit::TwoExpFit(ApplicationWindow *parent, Graph *g, const QString &curveTitle)
124 : Fit(parent, g)
125 {
126 init();
127 setDataFromCurve(curveTitle);
128 }
129
TwoExpFit(ApplicationWindow * parent,Graph * g,const QString & curveTitle,double start,double end)130 TwoExpFit::TwoExpFit(ApplicationWindow *parent, Graph *g, const QString &curveTitle, double start,
131 double end)
132 : Fit(parent, g)
133 {
134 init();
135 setDataFromCurve(curveTitle, start, end);
136 }
137
init()138 void TwoExpFit::init()
139 {
140 setObjectName("ExpDecay");
141 d_f = expd2_f;
142 d_df = expd2_df;
143 d_fdf = expd2_fdf;
144 d_fsimplex = expd2_d;
145 d_p = 5;
146 d_min_points = d_p;
147 d_param_init = gsl_vector_alloc(d_p);
148 gsl_vector_set_all(d_param_init, 1.0);
149 covar = gsl_matrix_alloc(d_p, d_p);
150 d_results.resize(d_p);
151 d_param_names << "A1"
152 << "t1"
153 << "A2"
154 << "t2"
155 << "y0";
156 d_explanation = tr("Exponential decay");
157 d_formula = "A1*exp(-x/t1)+A2*exp(-x/t2)+y0";
158 d_param_explain << tr("(first amplitude)") << tr("(first lifetime)") << tr("(second amplitude)")
159 << tr("(second lifetime)") << tr("(offset)");
160 }
161
storeCustomFitResults(const vector<double> & par)162 void TwoExpFit::storeCustomFitResults(const vector<double> &par)
163 {
164 d_results = par;
165 assert(d_results.size() > 3);
166 d_results[1] = 1.0 / d_results[1];
167 d_results[3] = 1.0 / d_results[3];
168 }
169
calculateFitCurveData(const vector<double> & par,std::vector<double> & X,std::vector<double> & Y)170 bool TwoExpFit::calculateFitCurveData(const vector<double> &par, std::vector<double> &X,
171 std::vector<double> &Y)
172 {
173 generateX(X);
174 for (int i = 0; i < d_points; i++)
175 Y[i] = par[0] * exp(-par[1] * X[i]) + par[2] * exp(-par[3] * X[i]) + par[4];
176 return true;
177 }
178
179 /*****************************************************************************
180 *
181 * Class ThreeExpFit
182 *
183 *****************************************************************************/
184
ThreeExpFit(ApplicationWindow * parent,Graph * g)185 ThreeExpFit::ThreeExpFit(ApplicationWindow *parent, Graph *g) : Fit(parent, g)
186 {
187 init();
188 }
189
ThreeExpFit(ApplicationWindow * parent,Graph * g,const QString & curveTitle)190 ThreeExpFit::ThreeExpFit(ApplicationWindow *parent, Graph *g, const QString &curveTitle)
191 : Fit(parent, g)
192 {
193 init();
194 setDataFromCurve(curveTitle);
195 }
196
ThreeExpFit(ApplicationWindow * parent,Graph * g,const QString & curveTitle,double start,double end)197 ThreeExpFit::ThreeExpFit(ApplicationWindow *parent, Graph *g, const QString &curveTitle,
198 double start, double end)
199 : Fit(parent, g)
200 {
201 init();
202 setDataFromCurve(curveTitle, start, end);
203 }
204
init()205 void ThreeExpFit::init()
206 {
207 setObjectName("ExpDecay");
208 d_f = expd3_f;
209 d_df = expd3_df;
210 d_fdf = expd3_fdf;
211 d_fsimplex = expd3_d;
212 d_p = 7;
213 d_min_points = d_p;
214 d_param_init = gsl_vector_alloc(d_p);
215 gsl_vector_set_all(d_param_init, 1.0);
216 covar = gsl_matrix_alloc(d_p, d_p);
217 d_results.resize(d_p);
218 d_param_names << "A1"
219 << "t1"
220 << "A2"
221 << "t2"
222 << "A3"
223 << "t3"
224 << "y0";
225 d_explanation = tr("Exponential decay");
226 d_formula = "A1*exp(-x/t1)+A2*exp(-x/t2)+A3*exp(-x/t3)+y0";
227 d_param_explain << tr("(first amplitude)") << tr("(first lifetime)") << tr("(second amplitude)")
228 << tr("(second lifetime)") << tr("(third amplitude)") << tr("(third lifetime)")
229 << tr("(offset)");
230 }
231
storeCustomFitResults(const vector<double> & par)232 void ThreeExpFit::storeCustomFitResults(const vector<double> &par)
233 {
234 d_results = par;
235 assert(d_results.size() > 5);
236 d_results[1] = 1.0 / d_results[1];
237 d_results[3] = 1.0 / d_results[3];
238 d_results[5] = 1.0 / d_results[5];
239 }
240
calculateFitCurveData(const vector<double> & par,std::vector<double> & X,std::vector<double> & Y)241 bool ThreeExpFit::calculateFitCurveData(const vector<double> &par, std::vector<double> &X,
242 std::vector<double> &Y)
243 {
244 generateX(X);
245 for (int i = 0; i < d_points; i++)
246 Y[i] = par[0] * exp(-X[i] * par[1]) + par[2] * exp(-X[i] * par[3])
247 + par[4] * exp(-X[i] * par[5]) + par[6];
248 return true;
249 }
250