Isis 3.0 Application Source Code Reference |
Home |
00001 /** 00002 * @file 00003 * $Revision: 1.2 $ 00004 * $Date: 2008/02/08 21:31:01 $ 00005 * $Id: NonLinearLSQ.cpp,v 1.2 2008/02/08 21:31:01 kbecker Exp $ 00006 * 00007 * Unless noted otherwise, the portions of Isis written by the USGS are 00008 * public domain. See individual third-party library and package descriptions 00009 * for intellectual property information, user agreements, and related 00010 * information. 00011 * 00012 * Although Isis has been used by the USGS, no warranty, expressed or 00013 * implied, is made by the USGS as to the accuracy and functioning of such 00014 * software and related material nor shall the fact of distribution 00015 * constitute any such warranty, and no responsibility is assumed by the 00016 * USGS in connection therewith. 00017 * 00018 * For additional information, launch 00019 * $ISISROOT/doc//documents/Disclaimers/Disclaimers.html 00020 * in a browser or see the Privacy & Disclaimers page on the Isis website, 00021 * http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on 00022 * http://www.usgs.gov/privacy.html. 00023 */ 00024 #include <string> 00025 #include <vector> 00026 #include <numeric> 00027 #include <iostream> 00028 #include <sstream> 00029 00030 #include <gsl/gsl_errno.h> 00031 #include <gsl/gsl_blas.h> 00032 #include <gsl/gsl_multifit_nlin.h> 00033 00034 #include "NonLinearLSQ.h" 00035 #include "iException.h" 00036 00037 using namespace std; 00038 00039 namespace Isis { 00040 00041 00042 int NonLinearLSQ::curvefit() { 00043 00044 size_t n(nSize()); 00045 size_t p(nParms()); 00046 00047 // Initialize the solver function information 00048 _nlsqPointer d = { this }; 00049 gsl_multifit_function_fdf mf; 00050 mf.f = &f; 00051 mf.df = &df; 00052 mf.fdf = &fdf; 00053 mf.n = n; 00054 mf.p = p; 00055 mf.params = &d; 00056 00057 const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder; 00058 gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc(T, n, p); 00059 00060 _fitParms = guess(); 00061 gsl_vector *x = NlsqTogsl(_fitParms); 00062 gsl_matrix *covar = gsl_matrix_alloc(p, p); 00063 gsl_multifit_fdfsolver_set(s, &mf, x); 00064 00065 _nIters = 0; 00066 checkIteration(_nIters, gslToNlsq(s->x), NLVector(p,999.0), 00067 gsl_blas_dnrm2(s->f), GSL_CONTINUE); 00068 00069 00070 do { 00071 _nIters++; 00072 00073 _status = gsl_multifit_fdfsolver_iterate(s); 00074 _fitParms = gslToNlsq(s->x); 00075 00076 gsl_multifit_covar(s->J, 0.0, covar); 00077 _uncert = getUncertainty(covar); 00078 00079 _status = checkIteration(_nIters, _fitParms, _uncert, gsl_blas_dnrm2(s->f), 00080 _status); 00081 if ( _status ) { break; } 00082 if(!doContinue()) { break; } 00083 00084 _status = gsl_multifit_test_delta(s->dx, s->x, absErr(), relErr()); 00085 } while ((_status == GSL_CONTINUE) && (_nIters < _maxIters)); 00086 00087 // Clean up 00088 gsl_multifit_fdfsolver_free(s); 00089 gsl_matrix_free(covar); 00090 00091 return (_status); 00092 } 00093 00094 00095 void NonLinearLSQ::Terminate(const std::string &message) { 00096 _userMessage = message; 00097 _userTerminated = true; 00098 _status = GSL_SUCCESS; 00099 return; 00100 } 00101 00102 void NonLinearLSQ::Abort(const std::string &reason) { 00103 _userMessage = reason; 00104 _userTerminated = true; 00105 _status = GSL_FAILURE; 00106 return; 00107 } 00108 00109 00110 int NonLinearLSQ::f(const gsl_vector *x, void *params, gsl_vector *fx) { 00111 NonLinearLSQ *nlsq = (static_cast<_nlsqPointer *> (params))->nlsq; 00112 int n = nlsq->nSize(); 00113 NLVector fxv = nlsq->f_x(nlsq->gslToNlsq(x)); 00114 for (int i = 0 ; i < n ; i++ ) { 00115 gsl_vector_set(fx, i, fxv[i]); 00116 } 00117 return (GSL_SUCCESS); 00118 } 00119 00120 int NonLinearLSQ::df(const gsl_vector *x, void *params, gsl_matrix *J) { 00121 NonLinearLSQ *nlsq = (static_cast<_nlsqPointer *> (params))->nlsq; 00122 int n = nlsq->nSize(); 00123 int p = nlsq->nParms(); 00124 00125 NLMatrix m = nlsq->df_x(nlsq->gslToNlsq(x)); 00126 00127 for (int i = 0 ; i < n ; i++ ) { 00128 for (int j = 0 ; j < p ; j++ ) { 00129 gsl_matrix_set(J, i, j, m[i][j]); 00130 } 00131 } 00132 return (GSL_SUCCESS); 00133 } 00134 00135 int NonLinearLSQ::fdf(const gsl_vector *x, void *params, gsl_vector *fx, 00136 gsl_matrix *J) { 00137 f(x,params,fx); 00138 df(x,params,J); 00139 return (GSL_SUCCESS); 00140 } 00141 00142 00143 NonLinearLSQ::NLVector NonLinearLSQ::getUncertainty(const gsl_matrix *covar) 00144 const { 00145 NLVector unc(covar->size1); 00146 for (size_t i = 0 ; i < covar->size1 ; i++ ) { 00147 unc[i] = sqrt(gsl_matrix_get(covar, i, i)); 00148 } 00149 return (unc); 00150 } 00151 00152 NonLinearLSQ::NLVector NonLinearLSQ::gslToNlsq(const gsl_vector *v) const { 00153 size_t n = v->size; 00154 NLVector Nv(n); 00155 for (size_t i = 0 ; i < n ; i++) { 00156 Nv[i] = gsl_vector_get(v, i); 00157 } 00158 return (Nv); 00159 } 00160 00161 NonLinearLSQ::NLMatrix NonLinearLSQ::gslToNlsq(const gsl_matrix *m) const { 00162 size_t nrows = m->size1; 00163 size_t ncols = m->size2; 00164 NLMatrix Nm(nrows, ncols); 00165 for (size_t i = 0 ; i < nrows ; i++) { 00166 for (size_t j = 0 ; j < ncols ; j++) { 00167 Nm[i][j] = gsl_matrix_get(m,i,j); 00168 } 00169 } 00170 return (Nm); 00171 } 00172 00173 gsl_vector *NonLinearLSQ::NlsqTogsl(const NonLinearLSQ::NLVector &v, 00174 gsl_vector *gv) const { 00175 if (gv == 0) { 00176 gv = gsl_vector_alloc(v.dim()); 00177 } 00178 else if (gv->size != (size_t) v.dim()) { 00179 ostringstream mess; 00180 mess << "Size of NL vector (" << v.dim() << ") not same as GSL vector (" 00181 << gv->size << ")"; 00182 throw iException::Message(iException::Programmer, mess.str(), _FILEINFO_); 00183 } 00184 00185 for (int i = 0 ; i < v.dim() ; i++) { 00186 gsl_vector_set(gv, i, v[i]); 00187 } 00188 return (gv); 00189 } 00190 00191 gsl_matrix *NonLinearLSQ::NlsqTogsl(const NonLinearLSQ::NLMatrix &m, 00192 gsl_matrix *gm) const { 00193 if (gm == 0) { 00194 gm = gsl_matrix_alloc(m.dim1(), m.dim2()); 00195 } 00196 else if ((gm->size1 != (size_t) m.dim1()) && 00197 (gm->size2 != (size_t) m.dim2()) ) { 00198 ostringstream mess; 00199 mess << "Size of NL matrix (" << m.dim1() << "," << m.dim2() 00200 << ") not same as GSL matrix (" << gm->size1 << "," << gm->size2 00201 << ")"; 00202 throw iException::Message(iException::Programmer, mess.str(), _FILEINFO_); 00203 } 00204 00205 for (int i = 0 ; i < m.dim1() ; i++) { 00206 for (int j = 0 ; j < m.dim2() ; j++) { 00207 gsl_matrix_set(gm, i, j, m[i][j]); 00208 } 00209 } 00210 return (gm); 00211 } 00212 00213 } // namespace ISIS