LMFitter.cxx
Go to the documentation of this file.
1 
12 #ifdef _MSC_VER
13 #include "msdevstudio/MSconfig.h"
14 #endif
15 
16 #include "LMFitter.h"
17 
18 #include "NumLinAlg.h"
19 #include "StatedFCN.h"
20 
21 #include <algorithm>
22 
23 #include <climits>
24 #include <cmath>
25 #include <cassert>
26 
27 #ifdef ITERATOR_MEMBER_DEFECT
28 using namespace std;
29 #else
30 using std::abs;
31 using std::distance;
32 using std::swap;
33 using std::vector;
34 #endif
35 
36 using namespace hippodraw::Numeric;
37 
38 using namespace hippodraw;
39 
40 LMFitter::
41 LMFitter ( const char * name )
42  : Fitter ( name ),
43  m_chi_cutoff ( 0.000001 ),
44  m_start_lambda ( 0.001 ),
45  m_lambda_shrink_factor( 9.8 ),
46  m_lambda_expand_factor( 10.2 )
47 {
48 }
49 
50 
52 LMFitter ( const LMFitter & fitter )
53  : Fitter ( fitter ),
54  m_chi_cutoff ( fitter.m_chi_cutoff ),
55  m_start_lambda ( fitter.m_start_lambda ),
56  m_lambda_shrink_factor ( fitter.m_lambda_shrink_factor ),
57  m_lambda_expand_factor ( fitter.m_lambda_expand_factor )
58 {
59 }
60 
61 
62 Fitter *
64 clone ( ) const
65 {
66  return new LMFitter ( *this );
67 }
68 
72 {
73  m_fcn -> calcAlphaBeta ( m_alpha, m_beta );
74  unsigned int num_parms = m_beta.size();
75 
76  unsigned int j = 0; // for MS VC++
77  for ( ; j < num_parms; j++ ) {
78  for ( unsigned int k = 0; k < j; k++ ) {
79  m_alpha[k][j] = m_alpha[j][k];
80  }
81  }
82 
83  j = 0; // for MS VC++
84  for ( ; j < num_parms; j++ ) {
85  m_alpha[j][j] *= ( 1.0 + m_lambda );
86  }
87 }
88 
94 int LMFitter::calcCovariance ( std::vector < std::vector < double > >& cov )
95 {
96  m_lambda = 0;
97  calcAlpha ();
98 
99  // Invert the semi hessian to get the error covarience matrix
100  // Since we use cholesky factorisation we can conclude if the
101  // given matrix is positive definite or not. So the return value
102  // is EXIT_SUCCESS if m_alpha is Positive Definite, EXIT_FAILURE otherwise.
103  return invertMatrix ( m_alpha, cov );
104 }
105 
107 {
108  unsigned int num_parms = m_beta.size ();
109 
110  vector< int > ipiv ( num_parms, 0 );
111 
112  vector< int > indxr ( num_parms, -1 );
113  vector< int > indxc ( num_parms, -1 );
114 
115  unsigned int irow = UINT_MAX;
116  unsigned int icol = UINT_MAX;
117 
118  for ( unsigned int i = 0; i < num_parms; i++ ) {
119  double big = 0.0;
120 
121  for ( unsigned int j = 0; j < num_parms; j++ ) {
122  if ( ipiv[j] != 1 ) {
123 
124  for ( unsigned int k = 0; k < num_parms; k++ ) {
125  if ( ipiv[k] == 0 ) {
126  if ( abs ( m_alpha[j][k] ) >= big ) {
127  big = abs ( m_alpha[j][k] );
128  irow = j;
129  icol = k;
130  }
131  }
132  else if ( ipiv[k] > 1 ) {
133  return false;
134  }
135  }
136  }
137  }
138 
139  if ( irow == UINT_MAX ) { // something is wrong, can't do fit.
140  return false;
141  }
142 
143  ++ipiv[icol];
144  if ( irow != icol ) {
145  for ( unsigned int l = 0; l < num_parms; l++ ) {
146  swap ( m_alpha[irow][l], m_alpha[icol][l] );
147  }
148  swap ( m_beta[irow], m_beta[icol] );
149  }
150  indxr[i] = irow;
151  indxc[i] = icol;
152  if ( m_alpha[icol][icol] == 0.0 ) {
153  return false;
154  }
155  double pivinv = 1.0 / m_alpha[icol][icol];
156  m_alpha[icol][icol] = 1.0;
157 
158  for ( unsigned int l = 0; l < num_parms; l++ ) {
159  m_alpha[icol][l] *= pivinv;
160  }
161  m_beta[icol] *= pivinv;
162 
163  for ( unsigned int ll = 0; ll < num_parms; ll++ ) {
164  if ( ll != icol ) {
165  double dum = m_alpha[ll][icol];
166  m_alpha[ll][icol] = 0.0;
167 
168  for ( unsigned int l = 0; l < num_parms; l++ ) {
169  m_alpha[ll][l] -= m_alpha[icol][l] * dum;
170  }
171  m_beta[ll] -= m_beta[icol] * dum;
172  }
173  }
174  }
175 
176  for ( int l = num_parms - 1; l >= 0; l-- ) {
177  if ( indxr[l] != indxc[l] ) {
178 
179  for ( unsigned int k = 0; k < num_parms; k++ ) {
180  swap ( m_alpha[k][indxr[l]], m_alpha[k][indxc[l]] );
181  }
182  }
183  }
184  return true;
185 }
186 
188 {
189  calcAlpha ();
190  bool ok = solveSystem ();
191 
192  return ok;
193 }
194 
196 {
198 
199  int i = 0;
200  for ( ; i < m_max_iterations; i++ ) {
201 
202  double old_chisq = objectiveValue ();
203 
204  vector< double > old_parms;
205  m_fcn -> fillFreeParameters ( old_parms );
206 
207  bool ok = calcStep ();
208  assert ( old_parms.size() == m_beta.size() );
209 
210  vector< double > new_parms ( old_parms );
211  vector< double >::iterator pit = new_parms.begin ( );
212  vector< double >::iterator dit = m_beta.begin ( );
213 
214  while ( pit != new_parms.end () ) {
215  *pit++ += *dit++;
216  }
217  m_fcn -> setFreeParameters ( new_parms );
218 
219  double new_chisq = objectiveValue ();
220 
221  if ( abs ( old_chisq - new_chisq ) < m_chi_cutoff ) break;
222 
223  if ( new_chisq < old_chisq ) {
225  }
226  else {
228  m_fcn -> setFreeParameters ( old_parms );
229  }
230 
231  if ( ! ok ) return ok;
232  }
233 
234  return i < m_max_iterations;
235 }
236 
237 
238 void
240 setFCN ( StatedFCN * fcn )
241 {
242  Fitter::setFCN ( fcn );
243  fcn -> setNeedsDerivatives ( true );
244 }

Generated for HippoDraw Class Library by doxygen