1 ///////////////////////////////////////////////////////////////////////////// 2 // Name: lm_lsqr.h 3 // Purpose: Levenberg-Marquart nonlinear least squares 2-D curve fitting 4 // Author: John Labenski, mostly others (see below) 5 // Modified by: 6 // Created: 6/5/2002 7 // Copyright: (c) John Labenski, mostly others (see below) 8 // Licence: Public domain 9 ///////////////////////////////////////////////////////////////////////////// 10 11 /* 12 * Solves or minimizes the sum of squares of m nonlinear 13 * functions of n variables. 14 * 15 * From public domain Fortran version 16 * of Argonne National Laboratories MINPACK 17 * 18 * argonne national laboratory. minpack project. march 1980. 19 * burton s. garbow, kenneth e. hillstrom, jorge j. more 20 * 21 * C translation by Steve Moshier http://www.moshier.net/ 22 * 23 * C++ "translation" for use w/ wxWindows by John Labenski 24 */ 25 26 #ifndef _LM_LEASTSQUARE_H_ 27 #define _LM_LEASTSQUARE_H_ 28 29 #include "wx/plotctrl/plotdefs.h" 30 class WXDLLIMPEXP_PLOTCTRL wxPlotData; 31 class WXDLLIMPEXP_PLOTCTRL wxPlotFunction; 32 33 // When SetLM_LeastSquareProgressHandler is called with a non NULL handler it will be 34 // called when fitting a curve every SetLM_LeastSquareProgressHandlerTicks 35 // text is the function name 36 // current is the current iteration, max is max allowed iterations 37 // note: current may exceed max by a few iterations in some cases 38 39 // Usage: create a function like this 40 // void LM_LeastSquareProgressHandler(const wxString &text, int current, int max) 41 // { [ do stuff... for example update a progress dialog ] 42 // wxString str = text + wxString::Format(wxT("\nIteration # %d of %d"), current, max); 43 // int percent = wxMin(int(100.0*current/max), 99); // iterations may overflow! 44 // return s_progressDialog->Update(percent, str); } 45 // 46 // then call SetLM_LeastSquareProgressHandler( LM_LeastSquareProgressHandler ); 47 48 extern "C" { 49 typedef bool (*LM_LeastSquareProgressHandler_)(const wxString &WXUNUSED(text), 50 int WXUNUSED(current), 51 int WXUNUSED(max)); 52 extern void SetLM_LeastSquareProgressHandler( LM_LeastSquareProgressHandler_ handler ); 53 extern void SetLM_LeastSquareProgressHandlerTicks( int iterations ); 54 } 55 56 //============================================================================= 57 // LM_LeastSquare - Levenberg-Marquart nonlinear least squares 2-D curve fitting 58 // 59 // Fit the plot function to the plot data 60 // 61 // Notes : 62 // the plot function must have fewer or equal vars than the plotdata has points 63 // the plotfunction MUST! have the LAST variable as 'x' 64 // you can set the starting values by filling 'initial_vals', size = (plotFunc.GetNumberVars - 1) 65 // if initial_vars = NULL then they are all 0.1 66 // 67 // Sample usage : 68 // wxString message; 69 // // Create some plotData, in this case from a known function 70 // wxPlotData data(wxPlotFunction("2.5*x*x-3*x+5+3.3*log(x)+13*exp(15*x/(x+4))", "x",dummy), 0, 1E-4, 10000); 71 // // Create the plotFunc we want to fit to the data, note: x is last var 72 // wxPlotFunction func("a*x*x+b*x+c+d*log(x)+e*exp(f*x/(x+g))", "a,b,c,d,e,f,g,x", message); 73 // LM_LeastSquare lmLeastSquare; 74 // if (lmLeastSquare.Create(data, func)) { 75 // lmLeastSquare.Fit(NULL); 76 // or Fit(init, init_count) where double init[init_count] = { a, b, c, ... } 77 // for (int k=0; k<lmLeastSquare.GetNumberVariables(); k++) // print a,b,c,... 78 // wxPrintf(wxT("%s=%g; "), func.GetVariableName(k).c_str(), lmLeastSquare.GetVariable(k)); 79 // 80 //============================================================================= 81 82 class WXDLLIMPEXP_PLOTCTRL LM_LeastSquare 83 { 84 public: 85 LM_LeastSquare(); ~LM_LeastSquare()86 virtual ~LM_LeastSquare() { Destroy(); } 87 88 // Initialize everything, returns sucess, on failure GetResultMessage() 89 // you may call Create and then Fit on a single instance on this 90 // as many times as you like. 91 bool Create(const wxPlotData &plotData, const wxPlotFunction &plotFunc); 92 // Has this been sucessfully created and is ready to be Fit() Ok()93 bool Ok() const { return m_ok; } 94 95 // After creation fit the plotFunc's vars to the plotData, returns # iterations 96 // initial_vals are initial guesses for the variables which may be NULL 97 // specify the number of initial variables with init_count 98 int Fit(const double *initial_vals = NULL, int init_count = 0); 99 // returns true if this is currently fitting IsFitting()100 bool IsFitting() const { return m_fitting; } 101 // abort a currently running fit, this may take a cycle or two so check 102 // IsFitting to determine when it is done. AbortFitting()103 void AbortFitting() { m_abort_fitting = true; } 104 // was the last fit aborted, reset when Create or Fit is called again GetAbortFitting()105 bool GetAbortFitting() const { return m_abort_fitting; } 106 107 // If you don't cal Fit(some_vars, count) then the variables are all 108 // initialized with this value, default = 0.1 GetInitValue()109 double GetInitValue() const { return m_init_value; } SetInitValue(double init_val)110 void SetInitValue(double init_val) { m_init_value = init_val; } 111 112 // Get the number of evaluations performed to find best fit GetNumberIterations()113 int GetNumberIterations() const { return m_nfev; } 114 // Get the euclidean norm of errors between data and function points GetEuclideanNorm()115 double GetEuclideanNorm() const { return m_fnorm; } 116 // Get the number of variables, i.e. (plotFunc.GetNumberVars() - 1, x is excluded) GetNumberVariables()117 int GetNumberVariables() const { return m_n; } 118 // Get the evaluated variables, size is (plotFunc.GetNumberVars() - 1, x is excluded) GetVariables()119 double *GetVariables() const { return m_x; } 120 // Get a single evaluated variable, 0 to GetNumberVariables()-1 121 double GetVariable(int n); 122 // Get an informational message about the results 123 wxString GetResultMessage() const; 124 125 protected: 126 void ReInit(); // only after a call to destroy - reset the vars 127 void Destroy(); 128 129 wxPlotData *m_plotData; 130 wxPlotFunction *m_plotFunc; 131 double m_init_value; 132 wxString m_resultMsg; 133 bool m_ok; 134 bool m_fitting; 135 bool m_abort_fitting; 136 137 // this is the function to calculate the difference 138 virtual void fcn(int m, int n, double x[], double fvec[], int *iflag); 139 140 void lmdif( int m, int n, double x[], double fvec[], double ftol, 141 double xtol, double gtol, int maxfev, double epsfcn, 142 double diag[], int mode, double factor, int nprint, int *info, 143 int *nfev, double fjac[], int ldfjac, int ipvt[], double qtf[], 144 double wa1[], double wa2[], double wa3[], double wa4[]); 145 146 // implementation - you probably don't want to mess with these! 147 148 void lmpar(int n, double r[], int ldr, int ipvt[], 149 double diag[], double qtb[], double delta, double *par, 150 double x[], double sdiag[], double wa1[], double wa2[]); 151 152 void qrfac(int m, int n, double a[], int lda, int pivot, int ipvt[], 153 int lipvt, double rdiag[], double acnorm[], double wa[]); 154 155 void qrsolv(int n, double r[], int ldr, int ipvt[], double diag[], 156 double qtb[], double x[], double sdiag[], double wa[]); 157 158 double enorm(int n, double x[]); 159 160 void fdjac2(int m,int n, double x[], double fvec[], double fjac[], 161 int ldfjac, int *iflag, double epsfcn, double wa[]); 162 163 int m_n; // # of variables of plotFunc 164 int m_m; // # of functions = points in plotData 165 int m_info; // index of info message strings 166 double m_fnorm; // euclidean norm of errors 167 double m_eps; // resolution of arithmetic 168 double m_dwarf; // smallest nonzero number 169 int m_nfev; // # iterations completed 170 unsigned long m_nan; // # if times function evaluation had a NaN 171 double m_ftol; // relative error in the sum of the squares, if less done 172 double m_xtol; // relative error between two iterations, if less done 173 double m_gtol; // cosine of the angle between fvec and any column of the jacobian, if less done 174 double m_epsfcn; // step length for the forward-difference approximation 175 double m_factor; // initial step bound 176 double *m_vars; // variables + 1, where last is var 'x' for wxPlotFunction 177 double *m_x; // variables (size m_n) 178 double *m_fvec; // output of evaluated functions (size m_m) 179 double *m_diag; // multiplicative scale factors for the variables, see m_mode 180 int m_mode; // =1 the vars scaled internally. if 2, scaling specified by m_diag. 181 double *m_fjac; // output m by n array 182 int m_ldfjac; // the leading dimension of the array fjac >= m_m 183 double *m_qtf; // output array the first n elements of the vector (q transpose)*fvec 184 int *m_ipvt; // integer output array of length n 185 int m_maxfev; // maximum number of iterations to try 186 187 private: 188 void Init(); 189 }; 190 191 #endif // _LM_LEASTSQUARE_H_ 192