MPQC  3.0.0-alpha
gaussianfit.timpl.h
1 //
2 // gaussianfit.timpl.h
3 //
4 // Copyright (C) 2007 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _math_optimize_gaussianfittimpl_h
29 #define _math_optimize_gaussianfittimpl_h
30 
31 #include <cmath>
32 #include <algorithm>
33 #include <math/optimize/levmar/lm.h>
34 #include <math/optimize/gaussianfit.h>
35 #include <util/misc/scexception.h>
36 
37 extern "C" {
38  void __eval_slater(double* params, double* f, int nparam, int np,
39  void *extraparams);
40  void __eval_slater_dfdp(double* params, double* dfdp, int nparam, int np,
41  void *extraparams);
42  void __eval_slater_pgauss(double* params, double* f, int nparam, int np,
43  void *extraparams);
44  void __eval_slater_dfdp_pgauss(double* params, double* dfdp, int nparam,
45  int np, void *extraparams);
46 }
47 ;
48 
49 namespace {
50  extern "C" {
51  typedef void
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,
54  void *adata);
55  }
56  ;
57 }
58 
59 namespace sc {
60 
61  template <typename Function, typename Weight> GaussianFit<Function, Weight>::GaussianFit(
62  unsigned int N,
63  const Weight& W,
64  double left,
65  double right,
66  unsigned int NP) :
67  weight_(W), gaussians_(), left_(left), right_(right), npts_(NP), p_(new double[2*N]), scratch_(new double[NP])
68  {
69  if (left < 0.0 || right < 0.0 || left >= right || N < 1 || NP < 2)
70  throw sc::ProgrammingError("GaussianFit::GaussianFit() -- invalid parameters",__FILE__,__LINE__);
71 
72  // Set initial parameters
73  gaussians_.resize(N);
74  // center exponents at 3 and use exp. ratio of 3
75  if (N%2) {
76  // Odd N
77  double gamma = 3.0;
78  for (int i=N/2; i>=0; --i, gamma/=3.0)
79  gaussians_[i] = std::make_pair(gamma, 1.0);
80  gamma = 9.0;
81  for (int i=N/2+1; i<N; ++i, gamma*=3.0)
82  gaussians_[i] = std::make_pair(gamma, 1.0);
83  } else {
84  // Even N
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);
91  }
92 
93  // zero out scratch
94  for (unsigned int p=0; p<npts_; ++p)
95  scratch_[p] = 0.0;
96  }
97 
98  template <typename Function, typename Weight> GaussianFit<Function, Weight>::~GaussianFit() {
99  delete[] p_;
100  delete[] scratch_;
101  }
102 
103  namespace {
104 
105  template <class Function, class Weight> struct ExtraParams {
106  const Function* f; // functor
107  const Weight* w; // weight
108  unsigned int npts; // number of points
109  double lb; // lower bound
110  double ub; // upper bound
111  int k; // fit to x^k exp(-a*x^2) functions
112  };
113 
114  template <class Function, class Weight> void eval_f(double* params,
115  double* f, int nparam,
116  int np,
117  void *extraparams) {
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;
128 
129  double x = 0.0;
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);
136  }
137 
138  if (GaussFit::weigh_F) {
139  const double df = fitvalue - W(x) * F(x);
140  f[p] = df;
141  } else {
142  const double df = fitvalue - F(x);
143  // lm will minimize the difference between the target function (identically 0.0) and the fitting function (Function - fit)
144  // to weigh the sum of squares by W need to multiply this by the square root of W
145  f[p] = df * std::sqrt(W(x));
146  }
147  }
148  }
149  template <class Function, class Weight> void eval_dfdp(double* params,
150  double* dfdp,
151  int nparam, int np,
152  void *extraparams) {
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;
163 
164  double x = 0.0;
165  unsigned int j = 0;
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);
176  dfdp[j++] = expgxx;
177  } else {
178  dfdp[j++] = -sqrt_W * (x * x * coef * expgxx);
179  dfdp[j++] = sqrt_W * expgxx;
180  }
181  }
182  }
183  }
184  } // namespace
185 
186  namespace detail {
187  template <class Function, class Weight> struct __to_extern_C_eval {
188  static eval_f_ptr f_ptr;
189  static eval_dfdp_ptr dfdp_ptr;
190  };
191 
192  }
193 
194  template <typename Function, typename Weight> const typename GaussianFit<Function,Weight>::Gaussians& GaussianFit<
195  Function, Weight>::operator()(const Function& F) const
196  {
197  ExtraParams<Function,Weight> xp;
198  xp.f = &F;
199  xp.w = &weight_;
200  xp.npts = npts_;
201  xp.lb = left_;
202  xp.ub = right_;
203  xp.k = k_;
204  const int ngaussians = gaussians_.size();
205 
206  extract_params();
207 
208  const int niter = dlevmar_der(detail::__to_extern_C_eval<Function,Weight>::f_ptr,
210  p_, scratch_, 2*ngaussians, npts_, 100000, NULL, NULL, NULL, NULL, static_cast<void*>(&xp));
211  if (niter> 0) {
212  if (classdebug_) {
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;
217  }
218  }
219 
220  assign_params();
221  }
222  else {
223  throw AlgorithmException("GaussianFit() -- fit failed",__FILE__,__LINE__);
224  }
225 
226  return gaussians_;
227  }
228 
229  template <typename Function,
230  typename Weight>
231  void
233  {
234  const unsigned int ngaussians = gaussians_.size();
235  unsigned int pp;
236  for(unsigned int p=0, pp=0; p<ngaussians; ++p) {
237  p_[pp++] = gaussians_[p].first;
238  p_[pp++] = gaussians_[p].second;
239  }
240  }
241 
242  template <typename Function,
243  typename Weight>
244  void
245  GaussianFit<Function,Weight>::assign_params() const
246  {
247  const unsigned int ngaussians = gaussians_.size();
248  unsigned int pp;
249  for(unsigned int p=0, pp=0; p<ngaussians; ++p) {
250  gaussians_[p].first = p_[pp++];
251  gaussians_[p].second = p_[pp++];
252  }
253  }
254 };
255 
256 #endif // include guards
257 
259 // Local Variables:
260 // mode: c++
261 // c-file-style: "CLJ-CONDENSED"
262 // End:
sc::GaussianFit::GaussianFit
GaussianFit(unsigned int N, const Weight &W, double left, double right, unsigned int NP)
Will fit a function F with weight W evenly sampled on NP points in an interval [left,...
Definition: gaussianfit.timpl.h:61
sc::AlgorithmException
This exception is thrown whenever a problem with an algorithm is encountered.
Definition: scexception.h:342
sc::detail::__to_extern_C_eval
Definition: gaussianfit.timpl.h:187
sc::Function
The Function class is an abstract base class that, given a set of coordinates, will compute a value a...
Definition: function.h:44
sc::GaussianFit
GaussianFit<Function> is a fit of Function(x)*Weight(x) to N Gaussians on range [left,...
Definition: gaussianfit.h:41
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14

Generated at Sun Jan 26 2020 23:24:00 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.