28 #ifndef _math_optimize_gaussianfittimpl_h
29 #define _math_optimize_gaussianfittimpl_h
33 #include <math/optimize/levmar/lm.h>
34 #include <math/optimize/gaussianfit.h>
35 #include <util/misc/scexception.h>
38 void __eval_slater(
double* params,
double* f,
int nparam,
int np,
40 void __eval_slater_dfdp(
double* params,
double* dfdp,
int nparam,
int np,
42 void __eval_slater_pgauss(
double* params,
double* f,
int nparam,
int np,
44 void __eval_slater_dfdp_pgauss(
double* params,
double* dfdp,
int nparam,
45 int np,
void *extraparams);
52 (*eval_f_ptr)(
double *p,
double *hx,
int m,
int n,
void *adata);
53 typedef void (*eval_dfdp_ptr)(
double *p,
double *hx,
int m,
int n,
67 weight_(W), gaussians_(), left_(left), right_(right), npts_(NP), p_(new double[2*N]), scratch_(new double[NP])
69 if (left < 0.0 || right < 0.0 || left >= right || N < 1 || NP < 2)
70 throw sc::ProgrammingError(
"GaussianFit::GaussianFit() -- invalid parameters",__FILE__,__LINE__);
78 for (
int i=N/2; i>=0; --i, gamma/=3.0)
79 gaussians_[i] = std::make_pair(gamma, 1.0);
81 for (
int i=N/2+1; i<N; ++i, gamma*=3.0)
82 gaussians_[i] = std::make_pair(gamma, 1.0);
85 double gamma = 3.0/sqrt(3.0);
86 for (
int i=N/2-1; i>=0; --i, gamma/=3.0)
87 gaussians_[i] = std::make_pair(gamma, 1.0);
88 gamma = 3.0*sqrt(3.0);
89 for (
int i=N/2; i<N; ++i, gamma*=3.0)
90 gaussians_[i] = std::make_pair(gamma, 1.0);
94 for (
unsigned int p=0; p<npts_; ++p)
105 template <
class Function,
class Weight>
struct ExtraParams {
114 template <
class Function,
class Weight>
void eval_f(
double* params,
115 double* f,
int nparam,
118 typedef ExtraParams<Function,Weight> XPS;
119 typedef GaussianFit<Function,Weight> GaussFit;
120 XPS* XParams = static_cast<XPS*>(extraparams);
121 const double lb = XParams->lb;
122 const double ub = XParams->ub;
123 const double dx = (ub - lb)/(np - 1);
124 const double npts = XParams->npts;
125 const Function& F = *(XParams->f);
126 const Weight& W = *(XParams->w);
127 const int k = XParams->k;
130 for (
unsigned int p=0; p<npts; ++p, x+=dx) {
131 double fitvalue = 0.0;
132 for (
unsigned int i=0; i<nparam;) {
133 const double gamma = params[i++];
134 const double coef = params[i++];
135 fitvalue += coef * std::pow(x, k) * std::exp(-gamma*x*x);
138 if (GaussFit::weigh_F) {
139 const double df = fitvalue - W(x) * F(x);
142 const double df = fitvalue - F(x);
145 f[p] = df * std::sqrt(W(x));
149 template <
class Function,
class Weight>
void eval_dfdp(
double* params,
153 typedef ExtraParams<Function,Weight> XPS;
154 typedef GaussianFit<Function,Weight> GaussFit;
155 XPS* XParams = static_cast<XPS*>(extraparams);
156 const double lb = XParams->lb;
157 const double ub = XParams->ub;
158 const double dx = (ub - lb)/(np - 1);
159 const double npts = XParams->npts;
160 const Function& F = *(XParams->f);
161 const Weight& W = *(XParams->w);
162 const int k = XParams->k;
166 for (
unsigned int p=0, j=0; p<npts; ++p, x+=dx) {
167 const double sqrt_W = std::sqrt(W(x));
168 for (
unsigned int i=0; i<nparam;) {
169 const unsigned int p1 = i++;
170 const unsigned int p2 = i++;
171 const double gamma = params[p1];
172 const double coef = params[p2];
173 const double expgxx = std::pow(x, k) * std::exp(-gamma*x*x);
174 if (GaussFit::weigh_F) {
175 dfdp[j++] = - (x * x * coef * expgxx);
178 dfdp[j++] = -sqrt_W * (x * x * coef * expgxx);
179 dfdp[j++] = sqrt_W * expgxx;
188 static eval_f_ptr f_ptr;
189 static eval_dfdp_ptr dfdp_ptr;
194 template <
typename Function,
typename Weight>
const typename GaussianFit<Function,Weight>::Gaussians&
GaussianFit<
197 ExtraParams<Function,Weight> xp;
204 const int ngaussians = gaussians_.size();
210 p_, scratch_, 2*ngaussians, npts_, 100000, NULL, NULL, NULL, NULL, static_cast<void*>(&xp));
213 std::cout <<
"GaussianFit() -- converged in " << niter <<
" iterations" << std::endl;
214 for(
unsigned int g=0, p=0; g<ngaussians; ++g) {
215 std::cout <<
" gamma = " << p_[p++];
216 std::cout <<
" coef = " << p_[p++] << std::endl;
234 const unsigned int ngaussians = gaussians_.size();
236 for(
unsigned int p=0, pp=0; p<ngaussians; ++p) {
237 p_[pp++] = gaussians_[p].first;
238 p_[pp++] = gaussians_[p].second;
242 template <
typename Function,
245 GaussianFit<Function,Weight>::assign_params()
const
247 const unsigned int ngaussians = gaussians_.size();
249 for(
unsigned int p=0, pp=0; p<ngaussians; ++p) {
250 gaussians_[p].first = p_[pp++];
251 gaussians_[p].second = p_[pp++];
256 #endif // include guards