Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitprofile.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitprofile.c	(revision 40753)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitprofile.c	(revision 40754)
@@ -4,12 +4,16 @@
 typedef enum {MODE_UPPER, MODE_LOWER} SigmoidMode;
 
-void sigmoid_fit (Vector *xprofile, Vector *fprofile, Vector *oprofile, double *par, int Npar, SigmoidMode mode, double limit);
+void sigmoid_fit (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar, FitFunc *FUNC);
+void sigmoid_fit_irls (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int NptsAll, opihi_flt *par, int Npar, SigmoidMode mode, double limit);
 
 opihi_flt fsigmoidDeimos (opihi_flt, opihi_flt *, int, opihi_flt *);
 opihi_flt rsigmoidDeimos (opihi_flt, opihi_flt *, int, opihi_flt *);
+
+double deimos_weight_cauchy (double x);
 
 // use these globals to carry into sigmoid_fit
 int VERBOSE = FALSE;
 int QUIET   = FALSE;
+int MaxIterations = 10;
 
 int deimos_fitprofile (int argc, char **argv) {
@@ -29,4 +33,5 @@
   Vector *xprofile = NULL;
   Vector *fprofile = NULL;
+  Vector *wprofile = NULL;
 
   QUIET = FALSE;
@@ -44,13 +49,21 @@
     remove_argument (N, &argc, argv);
   }
-
-  if (argc != 5) goto usage;
+  MaxIterations = 10;
+  if ((N = get_argument (argc, argv, "-max-iterations"))) {
+    remove_argument (N, &argc, argv);
+    MaxIterations = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 6) goto usage;
   if ((xprofile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((fprofile = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((wprofile = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if (xprofile->Nelements != fprofile->Nelements) goto usage;
+  if (xprofile->Nelements != wprofile->Nelements) goto usage;
 
   // supplied guess for boundaries
-  double Xs = atof(argv[3]);
-  double Xe = atof(argv[4]);
+  double Xs = atof(argv[4]);
+  double Xe = atof(argv[5]);
 
   // output vectors of fits
@@ -94,6 +107,9 @@
   par[2] = Sky;
   par[3] = BckL;
-
-  sigmoid_fit (xprofile, fprofile, Lprofile, par, Npar, MODE_LOWER, Xmid);
+  sigmoid_fit_irls (xprofile->elements.Flt, fprofile->elements.Flt, wprofile->elements.Flt, Lprofile->elements.Flt, xprofile->Nelements, par, Npar, MODE_LOWER, Xmid);
+  set_variable ("sigmoid_left_center", par[0]);
+  set_variable ("sigmoid_left_sigma",  par[1]);
+  set_variable ("sigmoid_left_top",    par[2]);
+  set_variable ("sigmoid_left_bottom", par[3]);
 
   par[0] = Xe;
@@ -101,11 +117,14 @@
   par[2] = Sky;
   par[3] = BckR;
-
-  sigmoid_fit (xprofile, fprofile, Rprofile, par, Npar, MODE_UPPER, Xmid);
+  sigmoid_fit_irls (xprofile->elements.Flt, fprofile->elements.Flt, wprofile->elements.Flt, Rprofile->elements.Flt, xprofile->Nelements, par, Npar, MODE_UPPER, Xmid);
+  set_variable ("sigmoid_right_center", par[0]);
+  set_variable ("sigmoid_right_sigma",  par[1]);
+  set_variable ("sigmoid_right_top",    par[2]);
+  set_variable ("sigmoid_right_bottom", par[3]);
  
   return TRUE;
   
 usage:
-  gprint (GP_ERR, "USAGE: deimos fitprofile xprof fprof Xs Xe\n");
+  gprint (GP_ERR, "USAGE: deimos fitprofile xprof fprof wprof Xs Xe\n");
   return FALSE;
 }
@@ -173,34 +192,85 @@
 }
 
-void sigmoid_fit (Vector *xprofile, Vector *fprofile, Vector *oprofile, double *par, int Npar, SigmoidMode mode, double limit) {
-
-  int Nelements = xprofile->Nelements;
-
-  ALLOCATE_PTR (xtemp, opihi_flt, Nelements);
-  ALLOCATE_PTR (ytemp, opihi_flt, Nelements);
-  ALLOCATE_PTR (wtemp, opihi_flt, Nelements);
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+
+// wprofile is the true poisson error
+void sigmoid_fit_irls (opihi_flt *xprofile, opihi_flt *fprofile, opihi_flt *wprofile, opihi_flt *oprofile, int NptsAll, opihi_flt *par, int Npar, SigmoidMode mode, double limit) {
+
+  FitFunc *FUNC = (mode == MODE_LOWER) ? fsigmoidDeimos : rsigmoidDeimos;
+
+  ALLOCATE_PTR (par_old, opihi_flt, Npar);
+
+  // allocate these as global static in this file?
+  ALLOCATE_PTR (xtemp, opihi_flt, NptsAll);
+  ALLOCATE_PTR (ytemp, opihi_flt, NptsAll);
+  ALLOCATE_PTR (wtemp, opihi_flt, NptsAll);
+  ALLOCATE_PTR (otemp, opihi_flt, NptsAll);
 
   // extract the points on the right side and fit
-  // XXX error model is pure poisson (allow external errors)
-
   int Npts = 0;
-  for (int i = 0; i < Nelements; i++) {
-    int keep = (mode == MODE_UPPER) ^ (xprofile->elements.Flt[i] < limit);
+  for (int i = 0; i < NptsAll; i++) {
+    int keep = (mode == MODE_UPPER) ^ (xprofile[i] < limit);
     if (keep) {
-      xtemp[Npts] = xprofile->elements.Flt[i];
-      ytemp[Npts] = fprofile->elements.Flt[i];
-      wtemp[Npts] = 1.0 / (fabs(fprofile->elements.Flt[i] + 0.01));
+      xtemp[Npts] = xprofile[i];
+      ytemp[Npts] = fprofile[i];
+      wtemp[Npts] = (wprofile[i] == 0.0) ? 1.0 : wprofile[i]; // crude
       Npts ++;
     }
   }
 
-  FitFunc *FUNC = (mode == MODE_LOWER) ? fsigmoidDeimos : rsigmoidDeimos;
-
-  double ochisq = mrqinit (xtemp, ytemp, wtemp, Npts, par, Npar, FUNC, VERBOSE);
+  // define initial weights as 1.0
+  ALLOCATE_PTR (wt, opihi_flt, Npts);
+  for (int i = 0; i < Npts; i++) { 
+    wt[i] = 1.0 / SQ(wtemp[i]);
+  }
+
+  // initial fit with ordinary least-squares
+  // otemp is the fitted value of the function at xtemp
+  sigmoid_fit (xtemp, ytemp, wt, otemp, Npts, par, Npar, FUNC);
+
+  int converged = FALSE;
+  for (int iterations = 0; !converged && (iterations < MaxIterations); iterations++) {
+
+    for (int i = 0; i < Npts; i++) {
+      wt[i] = deimos_weight_cauchy ((ytemp[i] - otemp[i]) / wtemp[i]) / SQ(wtemp[i]);
+    }
+
+    sigmoid_fit (xtemp, ytemp, wt, otemp, Npts, par, Npar, FUNC);
+
+    // save this solution
+    for (int i = 0; i < Npar; i++) {
+      par_old[i] = par[i];
+    }
+
+    converged = TRUE;
+    for (int i = 0; i < Npts; i++) {
+      if ((fabs(par[i] - par_old[i]) > FIT_TOLERANCE * fabs(par[i])) && 
+	  (fabs(par[i] - par_old[i]) > FLT_TOLERANCE))
+	converged = FALSE;
+    }
+  }
+
+  for (int i = 0; i < NptsAll; i++) {
+    oprofile[i] = FUNC (xprofile[i], par, Npar, NULL);
+  }
+
+  FREE (xtemp);
+  FREE (ytemp);
+  FREE (wtemp);
+  FREE (otemp);
+
+  FREE (par_old);
+  FREE (wt);
+}
+
+void sigmoid_fit (opihi_flt *xprofile, opihi_flt *yprofile, opihi_flt *wprofile, opihi_flt *oprofile, int Npts, opihi_flt *par, int Npar, FitFunc *FUNC) {
+
+  double ochisq = mrqinit (xprofile, yprofile, wprofile, Npts, par, Npar, FUNC, VERBOSE);
   double dchisq = ochisq + 2*Npts;
 
   int Niter;
   for (Niter = 0; (Niter < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); Niter++) {
-    double chisq = mrqmin (xtemp, ytemp, wtemp, Npts, par, Npar, FUNC, VERBOSE);
+    double chisq = mrqmin (xprofile, yprofile, wprofile, Npts, par, Npar, FUNC, VERBOSE);
     dchisq = ochisq - chisq;
     ochisq = chisq;
@@ -210,24 +280,15 @@
 
   if (!QUIET) gprint (GP_LOG, "sigmoid @ %f (%f) : top %f bottom %f\n", par[0], par[1], par[2], par[3]);
-  for (int i = 0; i < Nelements; i++) {
-    oprofile[0].elements.Flt[i] = FUNC (xprofile[0].elements.Flt[i], par, Npar, NULL);
-  }
-
-  if (mode == MODE_LOWER) {
-    set_variable ("sigmoid_left_center", par[0]);
-    set_variable ("sigmoid_left_sigma",  par[1]);
-    set_variable ("sigmoid_left_top",    par[2]);
-    set_variable ("sigmoid_left_bottom", par[3]);
-  } else {
-    set_variable ("sigmoid_right_center", par[0]);
-    set_variable ("sigmoid_right_sigma",  par[1]);
-    set_variable ("sigmoid_right_top",    par[2]);
-    set_variable ("sigmoid_right_bottom", par[3]);
-  }
+
+  for (int i = 0; i < Npts; i++) {
+    oprofile[i] = FUNC (xprofile[i], par, Npar, NULL);
+  }
+
   mrqfree (Npar);
-
-  FREE (xtemp);
-  FREE (ytemp);
-  FREE (wtemp);
-
 }  
+
+double deimos_weight_cauchy (double x) {
+  double r = x / 2.385;
+  return (1.0 / (1.0 + SQ(r)));
+}
+