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.

24 {}

Member Function Documentation

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

Definition at line 241 of file MarqFitAlg.cxx.

242  {
243  int i,j;
244  float det;
245 
246  std::vector<float> res(nData);
247  std::vector<float> dydp(nData*nParam);
248  std::vector<float> beta(nParam);
249  std::vector<float> alpha(nParam*nParam);
250  std::vector<std::vector<float>> alpsav(nParam,std::vector<float>(nParam));
251 
252  fgauss(y, p, nParam, nData, res);
253  dgauss(p, nParam, nData, dydp);
254  setup_matrix(res, dydp, nParam, nData, beta,alpha);
255  for(i=0;i<nParam;i++){
256  for(j=0;j<nParam;j++){
257  alpsav[i][j]=alpha[i*nParam+j];
258  }
259  }
260  det=invrt_matrix(alpha, nParam);
261 
262  if(det==0)return 1;
263  for(i=0;i<nParam;i++){
264  if(alpha[i*nParam+i]>=0.){
265  perr[i]=sqrt(alpha[i*nParam+i]);
266  }else{
267  perr[i]=alpha[i*nParam+i];
268  }
269  }
270 
271  return 0;
272  }
pionana::beta
double beta(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:59
alpha
double alpha
Definition: doAna.cpp:15
gshf::MarqFitAlg::fgauss
void fgauss(const float yd[], const float p[], const int npar, const int ndat, std::vector< float > &res)
Definition: MarqFitAlg.cxx:8
test.p
p
Definition: test.py:223
gshf::MarqFitAlg::dgauss
void dgauss(const float p[], const int npar, const int ndat, std::vector< float > &dydp)
Definition: MarqFitAlg.cxx:22
gshf::MarqFitAlg::invrt_matrix
float invrt_matrix(std::vector< float > &alphaf, const int npar)
Definition: MarqFitAlg.cxx:130
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:51
float gshf::MarqFitAlg::cal_xi2 ( const std::vector< float > &  res,
const int  ndat 
)
private

Definition at line 40 of file MarqFitAlg.cxx.

40  {
41  int i;
42  float xi2;
43  xi2=0.;
44  for(i=0;i<ndat;i++){
45  xi2+=res[i]*res[i];
46  }
47  return xi2;
48  }
void gshf::MarqFitAlg::dgauss ( const float  p[],
const int  npar,
const int  ndat,
std::vector< float > &  dydp 
)
private

Definition at line 22 of file MarqFitAlg.cxx.

22  {
23  #if defined WITH_OPENMP
24  //#pragma GCC ivdep
25  #pragma omp simd
26  #endif
27  for(int i=0;i<ndat;i++){
28  for(int j=0;j<npar;j+=3){
29  const float xmu=float(i)-p[j+1];
30  const float xmu_sg=xmu/p[j+2];
31  const float xmu_sg2=xmu_sg*xmu_sg;
32  dydp[i*npar+j] = std::exp(-0.5*xmu_sg2);
33  dydp[i*npar+j+1]=p[j]*dydp[i*npar+j]*xmu_sg/p[j+2];
34  dydp[i*npar+j+2]=dydp[i*npar+j+1]*xmu_sg;
35  }
36  }
37  }
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
test.p
p
Definition: test.py:223
void gshf::MarqFitAlg::fgauss ( const float  yd[],
const float  p[],
const int  npar,
const int  ndat,
std::vector< float > &  res 
)
private

Definition at line 8 of file MarqFitAlg.cxx.

8  {
9  #if defined WITH_OPENMP
10  #pragma omp simd
11  #endif
12  for(int i=0;i<ndat;i++){
13  float yf=0.;
14  for(int j=0;j<npar;j+=3){
15  yf = yf + p[j]*std::exp(-0.5*std::pow((float(i)-p[j+1])/p[j+2],2));
16  }
17  res[i]=yd[i]-yf;
18  }
19  }
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
test.p
p
Definition: test.py:223
float gshf::MarqFitAlg::invrt_matrix ( std::vector< float > &  alphaf,
const int  npar 
)
private

Definition at line 130 of file MarqFitAlg.cxx.

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

Definition at line 274 of file MarqFitAlg.cxx.

275  {
276  std::vector<float> plimmin(nParam,std::numeric_limits<float>::lowest());
277  std::vector<float> plimmax(nParam,std::numeric_limits<float>::max());
278  return mrqdtfit(lambda, p, &plimmin[0], &plimmax[0], y, nParam, nData, chiSqr, dchiSqr);
279  }
gshf::MarqFitAlg::mrqdtfit
int mrqdtfit(float &lambda, float p[], float y[], const int nParam, const int nData, float &chiSqr, float &dchiSqr)
Definition: MarqFitAlg.cxx:274
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
test.p
p
Definition: test.py:223
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 281 of file MarqFitAlg.cxx.

282  {
283  int j;
284  float nu,rho,lzmlh,amax,chiSq0;
285 
286  std::vector<float> res(nData);
287  std::vector<float> beta(nParam);
288  std::vector<float> dp(nParam);
289  std::vector<float> alpsav(nParam);
290  std::vector<float> psav(nParam);
291  std::vector<float> dydp(nData*nParam);
292  std::vector<float> alpha(nParam*nParam);
293 
294  bool haslimits = false;
295  for(j=0;j<nParam;j++){
296  if (plimmin[j]>std::numeric_limits<float>::lowest() || plimmax[j]<std::numeric_limits<float>::max()) {
297  haslimits = true;
298  break;
299  }
300  }
301 
302  fgauss(y, p, nParam, nData, res);
303  chiSq0=cal_xi2(res, nData);
304  dgauss(p, nParam, nData, dydp);
305  setup_matrix(res, dydp, nParam, nData, beta, alpha);
306  if(lambda<0.){
307  amax=-999.;
308  for(j = 0; j < nParam; j++){
309  if(alpha[j*nParam+j]>amax)amax=alpha[j*nParam+j];
310  }
311  lambda=0.001*amax;
312  }
313  for(j = 0; j < nParam; j++){
314  alpsav[j]=alpha[j*nParam+j];
315  alpha[j*nParam+j]=alpsav[j]+lambda;
316  }
317  solve_matrix(beta, alpha, nParam, dp);
318 
319  nu=2.;
320  rho=-1.;
321 
322  do{
323  for(j=0;j<nParam;j++){
324  psav[j] = p[j];
325  p[j] = p[j] + dp[j];
326  }
327  fgauss(y, p, nParam, nData, res);
328  chiSqr = cal_xi2(res, nData);
329  if (haslimits) {
330  for(j=0;j<nParam;j++){
331  if (p[j]<=plimmin[j] || p[j]>=plimmax[j]) chiSqr*=10000;//penalty for going out of limits!
332  }
333  }
334 
335  lzmlh=0.;
336  for(j=0;j<nParam;j++){
337  lzmlh+=dp[j]*(lambda*dp[j]+beta[j]);
338  }
339  rho=2.*(chiSq0-chiSqr)/lzmlh;
340  if (rho<0.){
341  for (j=0;j<nParam;j++)p[j]=psav[j];
342  chiSqr=chiSq0;
343  lambda = nu*lambda;
344  nu=2.*nu;
345  for(j = 0; j < nParam; j++){
346  alpha[j*nParam+j]=alpsav[j]+lambda;
347  }
348  solve_matrix(beta, alpha, nParam, dp);
349  }
350  } while(rho<0.);
351  lambda=lambda*fmax(0.333333,1.-pow(2.*rho-1.,3));
352  dchiSqr=chiSqr-chiSq0;
353  return 0;
354 
355  }
pionana::beta
double beta(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:59
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:84
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
alpha
double alpha
Definition: doAna.cpp:15
gshf::MarqFitAlg::fgauss
void fgauss(const float yd[], const float p[], const int npar, const int ndat, std::vector< float > &res)
Definition: MarqFitAlg.cxx:8
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
test.p
p
Definition: test.py:223
gshf::MarqFitAlg::cal_xi2
float cal_xi2(const std::vector< float > &res, const int ndat)
Definition: MarqFitAlg.cxx:40
gshf::MarqFitAlg::dgauss
void dgauss(const float p[], const int npar, const int ndat, std::vector< float > &dydp)
Definition: MarqFitAlg.cxx:22
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:51
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 51 of file MarqFitAlg.cxx.

52  {
53  int i,j,k;
54 
55  /* ... Calculate beta */
56  #if defined WITH_OPENMP
57  #pragma omp simd
58  //#pragma GCC ivdep
59  #endif
60  for(j=0;j<npar;j++){
61  beta[j]=0.0;
62  for(i=0;i<ndat;i++){
63  beta[j]+=res[i]*dydp[i*npar+j];
64  }
65  }
66 
67  /* ... Calculate alpha */
68  #if defined WITH_OPENMP
69  #pragma omp simd
70  //#pragma GCC ivdep
71  #endif
72  for (j = 0; j < npar; j++){
73  for (k = j; k < npar; k++){
74  alpha[j*npar+k]=0.0;
75  for(i=0;i<ndat;i++){
76  alpha[j*npar+k]+=dydp[i*npar+j]*dydp[i*npar+k];
77  }
78  if(k!=j)alpha[k*npar+j]=alpha[j*npar+k];
79  }
80  }
81  }
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 84 of file MarqFitAlg.cxx.

85  {
86  int i,j,k,imax;
87  float hmax,hsav;
88 
89  std::vector<std::vector<float>> h(npar, std::vector<float>(npar+1,0));
90 
91  /* ... set up augmented N x N+1 matrix */
92  for(i=0;i<npar;i++){
93  h[i][npar]=beta[i];
94  for(j=0;j<npar;j++){
95  h[i][j]=alpha[i*npar+j];
96  }
97  }
98 
99  /* ... diagonalize N x N matrix but do only terms required for solution */
100  for(i=0;i<npar;i++){
101  hmax=h[i][i];
102  imax=i;
103  for(j=i+1;j<npar;j++){
104  if(h[j][i]>hmax){
105  hmax=h[j][i];
106  imax=j;
107  }
108  }
109  if(imax!=i){
110  for(k=0;k<=npar;k++){
111  hsav=h[i][k];
112  h[i][k]=h[imax][k];
113  h[imax][k]=hsav;
114  }
115  }
116  for(j=0;j<npar;j++){
117  if(j==i)continue;
118  for(k=i;k<npar;k++){
119  h[j][k+1]-=h[i][k+1]*h[j][i]/h[i][i];
120  }
121  }
122  }
123  /* ... scale (N+1)'th column with factor which normalizes the diagonal */
124  for(i=0;i<npar;i++){
125  dp[i]=h[i][npar]/h[i][i];
126  }
127 
128  }
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11