LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Private Member Functions | List of all members
gshf::MarqFitAlg Class Reference

#include "MarqFitAlg.h"

Public Member Functions

 MarqFitAlg ()
 
virtual ~MarqFitAlg ()
 
int cal_perr (float p[], float y[], const int nParam, const int nData, float perr[])
 
int mrqdtfit (float &lambda, float p[], float y[], const int nParam, const int nData, float &chiSqr, float &dchiSqr)
 
int mrqdtfit (float &lambda, float p[], float plimmin[], float plimmax[], float y[], const int nParam, const int nData, float &chiSqr, float &dchiSqr)
 

Private Member Functions

void fgauss (const float yd[], const float p[], const int npar, const int ndat, std::vector< float > &res)
 
void dgauss (const float p[], const int npar, const int ndat, std::vector< float > &dydp)
 
float cal_xi2 (const std::vector< float > &res, const int ndat)
 
void setup_matrix (const std::vector< float > &res, const std::vector< float > &dydp, const int npar, const int ndat, std::vector< float > &beta, std::vector< float > &alpha)
 
void solve_matrix (const std::vector< float > &beta, const std::vector< float > &alpha, const int npar, std::vector< float > &dp)
 
float invrt_matrix (std::vector< float > &alphaf, const int npar)
 

Detailed Description

Definition at line 21 of file MarqFitAlg.h.

Constructor & Destructor Documentation

gshf::MarqFitAlg::MarqFitAlg ( )
explicit
virtual gshf::MarqFitAlg::~MarqFitAlg ( )
inlinevirtual

Definition at line 24 of file MarqFitAlg.h.

References cal_perr(), cal_xi2(), dgauss(), fgauss(), invrt_matrix(), mrqdtfit(), setup_matrix(), solve_matrix(), and y.

24 {}

Member Function Documentation

int gshf::MarqFitAlg::cal_perr ( float  p[],
float  y[],
const int  nParam,
const int  nData,
float  perr[] 
)

Definition at line 262 of file MarqFitAlg.cxx.

References dgauss(), fgauss(), invrt_matrix(), and setup_matrix().

Referenced by ~MarqFitAlg().

263  {
264  int i, j;
265  float det;
266 
267  std::vector<float> res(nData);
268  std::vector<float> dydp(nData * nParam);
269  std::vector<float> beta(nParam);
270  std::vector<float> alpha(nParam * nParam);
271  std::vector<std::vector<float>> alpsav(nParam, std::vector<float>(nParam));
272 
273  fgauss(y, p, nParam, nData, res);
274  dgauss(p, nParam, nData, dydp);
275  setup_matrix(res, dydp, nParam, nData, beta, alpha);
276  for (i = 0; i < nParam; i++) {
277  for (j = 0; j < nParam; j++) {
278  alpsav[i][j] = alpha[i * nParam + j];
279  }
280  }
281  det = invrt_matrix(alpha, nParam);
282 
283  if (det == 0) return 1;
284  for (i = 0; i < nParam; i++) {
285  if (alpha[i * nParam + i] >= 0.) { perr[i] = sqrt(alpha[i * nParam + i]); }
286  else {
287  perr[i] = alpha[i * nParam + i];
288  }
289  }
290 
291  return 0;
292  }
y
Float_t y
Definition: compare.C:6
gshf::MarqFitAlg::fgauss
void fgauss(const float yd[], const float p[], const int npar, const int ndat, std::vector< float > &res)
Definition: MarqFitAlg.cxx:7
gshf::MarqFitAlg::dgauss
void dgauss(const float p[], const int npar, const int ndat, std::vector< float > &dydp)
Definition: MarqFitAlg.cxx:26
gshf::MarqFitAlg::invrt_matrix
float invrt_matrix(std::vector< float > &alphaf, const int npar)
Definition: MarqFitAlg.cxx:149
gshf::MarqFitAlg::setup_matrix
void setup_matrix(const std::vector< float > &res, const std::vector< float > &dydp, const int npar, const int ndat, std::vector< float > &beta, std::vector< float > &alpha)
Definition: MarqFitAlg.cxx:59
float gshf::MarqFitAlg::cal_xi2 ( const std::vector< float > &  res,
const int  ndat 
)
private

Definition at line 47 of file MarqFitAlg.cxx.

Referenced by mrqdtfit(), and ~MarqFitAlg().

48  {
49  int i;
50  float xi2;
51  xi2 = 0.;
52  for (i = 0; i < ndat; i++) {
53  xi2 += res[i] * res[i];
54  }
55  return xi2;
56  }
void gshf::MarqFitAlg::dgauss ( const float  p[],
const int  npar,
const int  ndat,
std::vector< float > &  dydp 
)
private

Definition at line 26 of file MarqFitAlg.cxx.

Referenced by cal_perr(), mrqdtfit(), and ~MarqFitAlg().

27  {
28 #if defined WITH_OPENMP
29 #pragma omp simd
30 #ifdef __INTEL_COMPILER
31 #pragma ivdep
32 #endif
33 #endif
34  for (int i = 0; i < ndat; i++) {
35  for (int j = 0; j < npar; j += 3) {
36  const float xmu = float(i) - p[j + 1];
37  const float xmu_sg = xmu / p[j + 2];
38  const float xmu_sg2 = xmu_sg * xmu_sg;
39  dydp[i * npar + j] = std::exp(-0.5 * xmu_sg2);
40  dydp[i * npar + j + 1] = p[j] * dydp[i * npar + j] * xmu_sg / p[j + 2];
41  dydp[i * npar + j + 2] = dydp[i * npar + j + 1] * xmu_sg;
42  }
43  }
44  }
void gshf::MarqFitAlg::fgauss ( const float  yd[],
const float  p[],
const int  npar,
const int  ndat,
std::vector< float > &  res 
)
private

Definition at line 7 of file MarqFitAlg.cxx.

Referenced by cal_perr(), mrqdtfit(), and ~MarqFitAlg().

12  {
13 #if defined WITH_OPENMP
14 #pragma omp simd
15 #endif
16  for (int i = 0; i < ndat; i++) {
17  float yf = 0.;
18  for (int j = 0; j < npar; j += 3) {
19  yf = yf + p[j] * std::exp(-0.5 * std::pow((float(i) - p[j + 1]) / p[j + 2], 2));
20  }
21  res[i] = yd[i] - yf;
22  }
23  }
float gshf::MarqFitAlg::invrt_matrix ( std::vector< float > &  alphaf,
const int  npar 
)
private

Definition at line 149 of file MarqFitAlg.cxx.

Referenced by cal_perr(), and ~MarqFitAlg().

150  {
151  /*
152  Inverts the curvature matrix alpha using Gauss-Jordan elimination and
153  returns the determinant. This is based on the implementation in "Data
154  Reduction and Error Analysis for the Physical Sciences" by P. R. Bevington.
155  That implementation, in turn, uses the algorithm of the subroutine MINV,
156  described on page 44 of the IBM System/360 Scientific Subroutine Package
157  Program Description Manual (GH20-0586-0). This only needs to be called
158  once after the fit has converged if the parameter errors are desired.
159  */
160 
161  //turn input alphas into doubles
162  std::vector<double> alpha(npar * npar);
163 
164  int i, j, k;
165  std::vector<int> ik(npar);
166  std::vector<int> jk(npar);
167  double aMax, save, det;
168  float detf;
169 
170  for (i = 0; i < npar * npar; i++) {
171  alpha[i] = alphaf[i];
172  }
173 
174  det = 0;
175  /* ... search for the largest element which we will then put in the diagonal */
176  for (k = 0; k < npar; k++) {
177  aMax = 0;
178  for (i = k; i < npar; i++) {
179  for (j = k; j < npar; j++) {
180 
181  // alpha[i*npar+j]=alphaf[i*npar+j];
182 
183  if (fabs(alpha[i * npar + j]) > fabs(aMax)) {
184  aMax = alpha[i * npar + j];
185  ik[k] = i;
186  jk[k] = j;
187  }
188  }
189  }
190  if (aMax == 0) return (det); /* return 0 determinant to signal problem */
191  det = 1;
192  /* ... interchange rows if necessary to put aMax in diag */
193  i = ik[k];
194  if (i > k) {
195  for (j = 0; j < npar; j++) {
196  save = alpha[k * npar + j];
197  alpha[k * npar + j] = alpha[i * npar + j];
198  alpha[i * npar + j] = -save;
199  }
200  }
201  /* ... interchange columns if necessary to put aMax in diag */
202  j = jk[k];
203  if (j > k) {
204  for (i = 0; i < npar; i++) {
205  save = alpha[i * npar + k];
206  alpha[i * npar + k] = alpha[i * npar + j];
207  alpha[i * npar + j] = -save;
208  }
209  }
210  /* ... accumulate elements of inverse matrix */
211  for (i = 0; i < npar; i++) {
212  if (i != k) alpha[i * npar + k] = -alpha[i * npar + k] / aMax;
213  }
214 #if defined WITH_OPENMP
215 #pragma omp simd
216 #ifdef __INTEL_COMPILER
217 #pragma ivdep
218 #endif
219 #endif
220  for (i = 0; i < npar; i++) {
221  for (j = 0; j < npar; j++) {
222  if ((i != k) && (j != k))
223  alpha[i * npar + j] = alpha[i * npar + j] + alpha[i * npar + k] * alpha[k * npar + j];
224  }
225  }
226  for (j = 0; j < npar; j++) {
227  if (j != k) alpha[k * npar + j] = alpha[k * npar + j] / aMax;
228  }
229  alpha[k * npar + k] = 1 / aMax;
230  det = det * aMax;
231  }
232 
233  /* ... restore ordering of matrix */
234  for (k = npar - 1; k >= 0; k--) {
235  j = ik[k];
236  if (j > k) {
237  for (i = 0; i < npar; i++) {
238  save = alpha[i * npar + k];
239  alpha[i * npar + k] = -alpha[i * npar + j];
240  alpha[i * npar + j] = save;
241  }
242  }
243  i = jk[k];
244  if (i > k) {
245  for (j = 0; j < npar; j++) {
246  save = alpha[k * npar + j];
247  alpha[k * npar + j] = -alpha[i * npar + j];
248  alpha[i * npar + j] = save;
249  }
250  }
251  }
252 
253  for (i = 0; i < npar * npar; i++) {
254  alphaf[i] = alpha[i];
255  }
256 
257  detf = det;
258  return (detf);
259  }
int gshf::MarqFitAlg::mrqdtfit ( float &  lambda,
float  p[],
float  y[],
const int  nParam,
const int  nData,
float &  chiSqr,
float &  dchiSqr 
)

Definition at line 294 of file MarqFitAlg.cxx.

Referenced by util::LArFFTW::PeakCorrelation(), and ~MarqFitAlg().

301  {
302  std::vector<float> plimmin(nParam, std::numeric_limits<float>::lowest());
303  std::vector<float> plimmax(nParam, std::numeric_limits<float>::max());
304  return mrqdtfit(lambda, p, &plimmin[0], &plimmax[0], y, nParam, nData, chiSqr, dchiSqr);
305  }
y
Float_t y
Definition: compare.C:6
gshf::MarqFitAlg::mrqdtfit
int mrqdtfit(float &lambda, float p[], float y[], const int nParam, const int nData, float &chiSqr, float &dchiSqr)
Definition: MarqFitAlg.cxx:294
int gshf::MarqFitAlg::mrqdtfit ( float &  lambda,
float  p[],
float  plimmin[],
float  plimmax[],
float  y[],
const int  nParam,
const int  nData,
float &  chiSqr,
float &  dchiSqr 
)

Definition at line 307 of file MarqFitAlg.cxx.

References cal_xi2(), dgauss(), fgauss(), setup_matrix(), and solve_matrix().

316  {
317  int j;
318  float nu, rho, lzmlh, amax, chiSq0;
319 
320  std::vector<float> res(nData);
321  std::vector<float> beta(nParam);
322  std::vector<float> dp(nParam);
323  std::vector<float> alpsav(nParam);
324  std::vector<float> psav(nParam);
325  std::vector<float> dydp(nData * nParam);
326  std::vector<float> alpha(nParam * nParam);
327 
328  bool haslimits = false;
329  for (j = 0; j < nParam; j++) {
330  if (plimmin[j] > std::numeric_limits<float>::lowest() ||
331  plimmax[j] < std::numeric_limits<float>::max()) {
332  haslimits = true;
333  break;
334  }
335  }
336 
337  fgauss(y, p, nParam, nData, res);
338  chiSq0 = cal_xi2(res, nData);
339  dgauss(p, nParam, nData, dydp);
340  setup_matrix(res, dydp, nParam, nData, beta, alpha);
341  if (lambda < 0.) {
342  amax = -999.;
343  for (j = 0; j < nParam; j++) {
344  if (alpha[j * nParam + j] > amax) amax = alpha[j * nParam + j];
345  }
346  lambda = 0.001 * amax;
347  }
348  for (j = 0; j < nParam; j++) {
349  alpsav[j] = alpha[j * nParam + j];
350  alpha[j * nParam + j] = alpsav[j] + lambda;
351  }
352  solve_matrix(beta, alpha, nParam, dp);
353 
354  nu = 2.;
355  rho = -1.;
356 
357  do {
358  for (j = 0; j < nParam; j++) {
359  psav[j] = p[j];
360  p[j] = p[j] + dp[j];
361  }
362  fgauss(y, p, nParam, nData, res);
363  chiSqr = cal_xi2(res, nData);
364  if (haslimits) {
365  for (j = 0; j < nParam; j++) {
366  if (p[j] <= plimmin[j] || p[j] >= plimmax[j])
367  chiSqr *= 10000; //penalty for going out of limits!
368  }
369  }
370 
371  lzmlh = 0.;
372  for (j = 0; j < nParam; j++) {
373  lzmlh += dp[j] * (lambda * dp[j] + beta[j]);
374  }
375  rho = 2. * (chiSq0 - chiSqr) / lzmlh;
376  if (rho < 0.) {
377  for (j = 0; j < nParam; j++)
378  p[j] = psav[j];
379  chiSqr = chiSq0;
380  lambda = nu * lambda;
381  nu = 2. * nu;
382  for (j = 0; j < nParam; j++) {
383  alpha[j * nParam + j] = alpsav[j] + lambda;
384  }
385  solve_matrix(beta, alpha, nParam, dp);
386  }
387  } while (rho < 0.);
388  lambda = lambda * fmax(0.333333, 1. - pow(2. * rho - 1., 3));
389  dchiSqr = chiSqr - chiSq0;
390  return 0;
391  }
gshf::MarqFitAlg::solve_matrix
void solve_matrix(const std::vector< float > &beta, const std::vector< float > &alpha, const int npar, std::vector< float > &dp)
Definition: MarqFitAlg.cxx:101
y
Float_t y
Definition: compare.C:6
gshf::MarqFitAlg::fgauss
void fgauss(const float yd[], const float p[], const int npar, const int ndat, std::vector< float > &res)
Definition: MarqFitAlg.cxx:7
gshf::MarqFitAlg::cal_xi2
float cal_xi2(const std::vector< float > &res, const int ndat)
Definition: MarqFitAlg.cxx:47
gshf::MarqFitAlg::dgauss
void dgauss(const float p[], const int npar, const int ndat, std::vector< float > &dydp)
Definition: MarqFitAlg.cxx:26
gshf::MarqFitAlg::setup_matrix
void setup_matrix(const std::vector< float > &res, const std::vector< float > &dydp, const int npar, const int ndat, std::vector< float > &beta, std::vector< float > &alpha)
Definition: MarqFitAlg.cxx:59
void gshf::MarqFitAlg::setup_matrix ( const std::vector< float > &  res,
const std::vector< float > &  dydp,
const int  npar,
const int  ndat,
std::vector< float > &  beta,
std::vector< float > &  alpha 
)
private

Definition at line 59 of file MarqFitAlg.cxx.

Referenced by cal_perr(), mrqdtfit(), and ~MarqFitAlg().

65  {
66  int i, j, k;
67 
68 /* ... Calculate beta */
69 #if defined WITH_OPENMP
70 #pragma omp simd
71 #ifdef __INTEL_COMPILER
72 #pragma ivdep
73 #endif
74 #endif
75  for (j = 0; j < npar; j++) {
76  beta[j] = 0.0;
77  for (i = 0; i < ndat; i++) {
78  beta[j] += res[i] * dydp[i * npar + j];
79  }
80  }
81 
82 /* ... Calculate alpha */
83 #if defined WITH_OPENMP
84 #pragma omp simd
85 #ifdef __INTEL_COMPILER
86 #pragma ivdep
87 #endif
88 #endif
89  for (j = 0; j < npar; j++) {
90  for (k = j; k < npar; k++) {
91  alpha[j * npar + k] = 0.0;
92  for (i = 0; i < ndat; i++) {
93  alpha[j * npar + k] += dydp[i * npar + j] * dydp[i * npar + k];
94  }
95  if (k != j) alpha[k * npar + j] = alpha[j * npar + k];
96  }
97  }
98  }
void gshf::MarqFitAlg::solve_matrix ( const std::vector< float > &  beta,
const std::vector< float > &  alpha,
const int  npar,
std::vector< float > &  dp 
)
private

Definition at line 101 of file MarqFitAlg.cxx.

Referenced by mrqdtfit(), and ~MarqFitAlg().

105  {
106  int i, j, k, imax;
107  float hmax, hsav;
108 
109  std::vector<std::vector<float>> h(npar, std::vector<float>(npar + 1, 0));
110 
111  /* ... set up augmented N x N+1 matrix */
112  for (i = 0; i < npar; i++) {
113  h[i][npar] = beta[i];
114  for (j = 0; j < npar; j++) {
115  h[i][j] = alpha[i * npar + j];
116  }
117  }
118 
119  /* ... diagonalize N x N matrix but do only terms required for solution */
120  for (i = 0; i < npar; i++) {
121  hmax = h[i][i];
122  imax = i;
123  for (j = i + 1; j < npar; j++) {
124  if (h[j][i] > hmax) {
125  hmax = h[j][i];
126  imax = j;
127  }
128  }
129  if (imax != i) {
130  for (k = 0; k <= npar; k++) {
131  hsav = h[i][k];
132  h[i][k] = h[imax][k];
133  h[imax][k] = hsav;
134  }
135  }
136  for (j = 0; j < npar; j++) {
137  if (j == i) continue;
138  for (k = i; k < npar; k++) {
139  h[j][k + 1] -= h[i][k + 1] * h[j][i] / h[i][i];
140  }
141  }
142  }
143  /* ... scale (N+1)'th column with factor which normalizes the diagonal */
144  for (i = 0; i < npar; i++) {
145  dp[i] = h[i][npar] / h[i][i];
146  }
147  }
The documentation for this class was generated from the following files: