NonLinearLSQ.h

Go to the documentation of this file.

#ifndef NonLinearLSQ_h
#define NonLinearLSQ_h

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>

#include "tnt_array1d.h"
#include "tnt_array1d_utils.h"
#include "tnt_array2d.h"
#include "tnt_array2d_utils.h"

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>

#include "IException.h"

namespace Isis {

  class NonLinearLSQ {
    public:
      typedef TNT::Array1D<double> NLVector;
      typedef TNT::Array2D<double> NLMatrix;

      //  Constructors and Destructor
      NonLinearLSQ() : _fitParms(), _uncert(), _nIters(0), _maxIters(50),
                       _status(0), _userTerminated(false), _userMessage() {}

      virtual ~NonLinearLSQ() { }

      virtual int nSize() const = 0;
      virtual int nParms() const = 0;

      void setMaxIters(int m) { _maxIters = m; }

      int maxIters() const { return (_maxIters); }

      virtual NLVector guess()  = 0;
      virtual NLVector f_x(const NLVector &x) = 0;
      virtual NLMatrix df_x(const NLVector &x) = 0;

      virtual double absErr() const { return (1.0E-4); }
      virtual double relErr() const { return (1.0E-4); }

      int curvefit();
      inline int status() const { return (_status); }
      inline bool success() const { return (_status == GSL_SUCCESS); }
      inline bool success(int status) const { return (status == GSL_SUCCESS); }
      inline std::string statusstr() const { 
        return (std::string(gsl_strerror(_status)));
      }
      inline std::string statusstr(int status) const { 
        return (std::string(gsl_strerror(status)));
      }

      virtual int checkIteration(const int Iter, const NLVector &fitcoefs, 
                                 const NLVector &uncerts, double cplxconj,
                                 int Istatus) { 
        return (Istatus);
      }

      inline NLVector coefs() const { return (_fitParms); }
      inline NLVector uncert() const { return (_uncert); }
      inline int nIterations() const { return (_nIters); }

    protected:
      void Terminate(const std::string &message = "");
      void Abort(const std::string &reason = "");

      bool doContinue() const { return (!_userTerminated); }


    private:
      NLVector    _fitParms;
      NLVector    _uncert;
      int         _nIters;
      int         _maxIters;
      int         _status;
      bool        _userTerminated;
      std::string _userMessage;

      struct _nlsqPointer {
          NonLinearLSQ *nlsq;
      };
      static int f(const gsl_vector *x, void *params, gsl_vector *fx);
      static int df(const gsl_vector *x, void *params, gsl_matrix *J);
      static int fdf(const gsl_vector *x, void *params, gsl_vector *fx, 
                     gsl_matrix *J);

      NLVector gslToNlsq(const gsl_vector *v) const;
      NLMatrix gslToNlsq(const gsl_matrix *m) const;
      gsl_vector *NlsqTogsl(const NLVector &v, gsl_vector *gv = 0) const;
      gsl_matrix *NlsqTogsl(const NLMatrix &m, gsl_matrix *gm = 0) const;
      NLVector getUncertainty(const gsl_matrix *m) const;
  };

}     // namespace Isis
#endif