$$\phantom{******** If you see this on the webpage then the browser could not locate *********}$$
$$\phantom{******** jsMath/easy/load.js or the variable root is set wrong in this file *********}$$
$$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$$ $$\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$$ $$\newcommand{\A}{{\cal A}}$$ $$\newcommand{\W}{{\cal W}}$$

/*

*/

#ifdef _WIN32
  #include <windows.h>
  #include <time.h>
  #include <float.h>
  #define sleep Sleep
#else
  #include <stdlib.h>
  #define BOOL int
  #define TRUE 1
  #define FALSE 0
#endif

#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include <math.h>
#include <string.h>
#include <vector>
using std::vector;

#include "lmmin_fit.h"
#include "callbk_reportstream.h"
#include "lmmin.h"

#define UNSET_VALUE 86738.82583

class qub_lmmin_state
{
public:
  lmmin_CurveFunc curve;
  int Nparam;
  double *params;
  double *lo, *hi;
  int *can_fit;
  int Ndata;
  float *yy, *ww;
  lmmin_IterFunc on_iter;
  callbk_reportfun on_report;
  void *obj;
  int max_iter;
  double toler;
  float *ff;
  double *ssr;
  int iter;
  vector<int> useCpar;
  int Ncurvepar;
  bool stop;

  qub_lmmin_state(lmmin_CurveFunc _curve, int _Nparam, double *_params, double *_lo, double *_hi, int *_can_fit,
			      int _Ndata, float *_yy, float *_ww, lmmin_IterFunc _on_iter, callbk_reportfun _on_report, void *_obj,
			      int _max_iter, double _toler, float *_ff, double *_ssr)
    : curve(_curve), Nparam(_Nparam), params(_params), lo(_lo), hi(_hi), can_fit(_can_fit),
      Ndata(_Ndata), yy(_yy), ww(_ww), on_iter(_on_iter), on_report(_on_report), obj(_obj),
      max_iter(_max_iter), toler(_toler), ff(_ff), ssr(_ssr), iter(0), stop(false)
  {}
};

void qub_lmmin_eval( double* par, int QUBFAST_VAR_NOT_USED m_dat, double* fvec, 
		     qub_lmmin_state *self, int *info )
{
  if ( self->stop ) {
    *info = -1;
    return;
  }
  bool outOfBounds = false;
  for (int i=self->Ncurvepar-1; i>=0; --i) {
    double lim = self->lo[ self->useCpar[i] ];
    if ( (fabs(lim - UNSET_VALUE) > 1e-8) && (lim > par[i]) )
      outOfBounds = true;
    lim = self->hi[ self->useCpar[i] ];
    if ( (fabs(lim - UNSET_VALUE) > 1e-8) && (lim < par[i]) )
      outOfBounds = true;
    if ( outOfBounds )
      break;
    self->params[ self->useCpar[i] ] = par[i];
  }
  if ( outOfBounds ) {
    for (int i=self->Ndata-1; i>=0; --i)
      fvec[i] = 1e10;
    return;
  }
  self->curve(self->obj, self->params, self->ff);
  double ssr = 0.0;
  for (int i=self->Ndata-1; i>=0; --i) {
    double resid = self->ww[i] * (self->yy[i] - self->ff[i]);
    fvec[i] = resid;
    ssr += resid * resid;
  }
  * self->ssr = ssr;
}

void qub_lmmin_print( QUBFAST_VAR_NOT_USED int n_par, QUBFAST_VAR_NOT_USED double* par, QUBFAST_VAR_NOT_USED int m_dat, QUBFAST_VAR_NOT_USED double* fvec, 
		      qub_lmmin_state *self, QUBFAST_VAR_NOT_USED int iflag, QUBFAST_VAR_NOT_USED int iter, QUBFAST_VAR_NOT_USED int nfev )
{
  self->iter = iter;
  if ( self->on_iter && ! self->on_iter(self->obj, self->params, self->iter) )
    self->stop = true;
}


QUBFAST_API int qub_lmmin_fit(lmmin_CurveFunc curve, int Nparam, double *params, double *lo, double *hi, int *can_fit,
			      int Ndata, float *yy, float *ww, lmmin_IterFunc on_iter, callbk_reportfun on_report, void *obj,
			      int max_iter, double toler, float *ff, double *ssr)
{
  *ssr = 0.0;
  qub_lmmin_state self(curve, Nparam, params, lo, hi, can_fit, Ndata, yy, ww, on_iter, on_report, obj, max_iter, toler, ff, ssr);
  lm_control_type control;
  control.maxcall = max_iter;
  control.stepbound = 100.0;
  control.epsilon = toler;
  control.ftol = toler;
  control.xtol = toler;
  control.gtol = toler;
  control.curvature = NULL;
  
  for (int i=0; i<Nparam; ++i)
    if ( can_fit[i] )
      self.useCpar.push_back(i);
  self.Ncurvepar = (int) self.useCpar.size();
  
  vector<double> curvepars( self.Ncurvepar );
  for (int i=0; i<self.Ncurvepar; ++i)
    curvepars[i] = params[self.useCpar[i]];

  lm_minimize( Ndata, self.Ncurvepar, &(curvepars[0]),
	       (lm_evaluate_ftype*)qub_lmmin_eval, (lm_print_ftype*)qub_lmmin_print, &self, &control );
  
  for (int i=self.Ncurvepar-1; i>=0; --i)
    params[ self.useCpar[i] ] = curvepars[i];

  // make sure the ssr and fit curve are for final params:
  vector<double> fvec(Ndata);
  int info;
  qub_lmmin_eval( &(curvepars[0]), Ndata, &(fvec[0]), &self, &info);

#if defined(_WIN32)
	unsigned int fstat = _clearfp(); // just in case
	_fpreset();
#endif  
  
  return self.iter;
}