Index: trunk/psModules/src/objects/pmSourceFitModel.c
===================================================================
--- trunk/psModules/src/objects/pmSourceFitModel.c	(revision 10199)
+++ trunk/psModules/src/objects/pmSourceFitModel.c	(revision 10259)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-11-26 22:26:51 $
+ *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-11-29 02:41:31 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -36,5 +36,5 @@
 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
 static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
-static bool PM_PSF_POISSON_WEIGHTS = true;
+static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
 
 bool pmSourceFitModelInit (float nIter, float tol, float weight, bool poissonErrors)
@@ -44,5 +44,5 @@
     PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
     PM_SOURCE_FIT_MODEL_WEIGHT = weight;
-    PM_PSF_POISSON_WEIGHTS = poissonErrors;
+    PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = poissonErrors;
 
     return true;
@@ -53,34 +53,23 @@
                        pmSourceFitMode mode)
 {
-    psTrace("psModules.objects", 5, "---- %s() begin ----\n", __func__);
+    psTrace("psModules.objects", 5, "---- %s begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     PS_ASSERT_PTR_NON_NULL(source->mask, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
 
-    // XXX EAM : is it necessary for the mask & weight to exist?  the
-    //           tests below could be conditions (!NULL)
-
     psBool fitStatus = true;
     psBool onPic     = true;
     psBool rc        = true;
 
-    psVector *params = model->params;
-    psVector *dparams = model->dparams;
-    psVector *paramMask = NULL;
-
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-
-    psF32 dSky = source->moments->dSky;
-
     // maximum number of valid pixels
     psS32 nPix = source->pixels->numRows * source->pixels->numCols;
 
-    // construct the coordinate and value entries
+    // arrays to hold the data to be fitted
     psArray *x = psArrayAllocEmpty(nPix);
     psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
     psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
 
+    // fill in the coordinate and value entries
     nPix = 0;
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
@@ -102,13 +91,13 @@
             x->data[nPix] = (psPtr *) coord;
             y->data.F32[nPix] = source->pixels->data.F32[i][j];
-            // psMinimizeLMChi2 takes wt = 1/dY^2
-            if (PM_PSF_POISSON_WEIGHTS) {
+
+            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
+            // as weight to avoid the bias from systematic errors here we would just use the
+            // source sky variance
+            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
                 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
             } else {
-                yErr->data.F32[nPix] = 1.0 / dSky;
-            }
-            // XXX EAM : suggestion from RHL is to use the local sky as weight
-            // to avoid the bias from systematic errors
-            // here we would just use the source sky variance
+                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
+            }
             nPix++;
         }
@@ -118,16 +107,20 @@
     yErr->n = nPix;
 
-    // XXX EAM : the new minimization API supplies the constraints as a struct
-    // XXX the chisq if a fcn of source flux. the minimization criterion should
-    // probably  be scaled somehow by flux.  measure the trend?  it depends on
-    // the about of systematic error in the models themselves...
-    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
-    psMinConstrain *constrain = psMinConstrainAlloc();
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+    if (!modelFunc)
+        psAbort ("pmSourceFitModel", "invalid model function");
+    pmModelLimits checkLimits = pmModelLimits_GetFunction (model->type);
+    if (!checkLimits)
+        psAbort ("pmSourceFitModel", "invalid model limits function");
+
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+    constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+    constraint->checkLimits = checkLimits;
 
     // set parameter mask based on fitting mode
-    paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-    psVectorInit (paramMask, 1);
-
     int nParams = 0;
     switch (mode) {
@@ -135,67 +128,52 @@
         // NORM-only model fits only source normalization (Io)
         nParams = 1;
-        paramMask->data.U8[PM_PAR_I0] = 0;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.U8[PM_PAR_I0] = 0;
         break;
     case PM_SOURCE_FIT_PSF:
         // PSF model only fits x,y,Io
         nParams = 3;
-        paramMask->data.U8[PM_PAR_I0] = 0;
-        paramMask->data.U8[PM_PAR_XPOS] = 0;
-        paramMask->data.U8[PM_PAR_YPOS] = 0;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.U8[PM_PAR_I0] = 0;
+        constraint->paramMask->data.U8[PM_PAR_XPOS] = 0;
+        constraint->paramMask->data.U8[PM_PAR_YPOS] = 0;
         break;
     case PM_SOURCE_FIT_EXT:
         // EXT model fits all params (except sky)
         nParams = params->n - 1;
-        psVectorInit (paramMask, 0);
-        paramMask->data.U8[PM_PAR_SKY] = 1;
-        break;
-    case PM_SOURCE_FIT_PSF_AND_SKY:
-        nParams = 4;
-        psAbort ("pmSourceFitModel", "PSF + SKY not currently available");
-        break;
-    case PM_SOURCE_FIT_EXT_AND_SKY:
-        nParams = params->n;
-        psAbort ("pmSourceFitModel", "EXT + SKY not currently available");
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.U8[PM_PAR_SKY] = 1;
         break;
     default:
         psAbort ("pmSourceFitModel", "invalid fitting mode");
     }
-    constrain->paramMask = paramMask;
+    // force the floating parameters to fall within the contraint ranges
+    for (int i = 0; i < params->n; i++) {
+        checkLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        checkLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+    }
 
     if (nPix <  nParams + 1) {
         psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
-        psTrace ("psModules.objects", 5, "---- %s(false) end ----\n", __func__);
         model->status = PM_MODEL_BADARGS;
         psFree (x);
         psFree (y);
         psFree (yErr);
-        psFree (myMin);
-        psFree (constrain);
-        psFree (paramMask);
+        psFree (constraint);
         return(false);
     }
 
-    // Set the parameter range checks
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-    modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
-
-    // force the floating parameters to fall within the contraint ranges
-    for (int i = 0; i < params->n; i++) {
-        if (constrain->paramMask->data.U8[i])
+    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if (psTraceGetLevel("psModules.objects") >= 4) {
+            fprintf (stderr, "%f ", params->data.F32[i]);
+        }
+        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
             continue;
-        params->data.F32[i] = PS_MIN(constrain->paramMax->data.F32[i], params->data.F32[i]);
-        params->data.F32[i] = PS_MAX(constrain->paramMin->data.F32[i], params->data.F32[i]);
-    }
-
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-
-    psTrace (__func__, 5, "fitting function\n");
-
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc);
-    for (int i = 0; i < dparams->n; i++) {
-        psTrace ("psModules.objects", 4, "%f ", params->data.F32[i]);
-        if ((paramMask != NULL) && paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     }
     psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
@@ -208,19 +186,19 @@
     // get the Gauss-Newton distance for fixed model parameters
     // hold the fitted parameters fixed; mask sky which is not fitted at all
-    if (paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+    if (constraint->paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
         psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
         altmask->data.U8[0] = 1;
         for (int i = 1; i < dparams->n; i++) {
-            altmask->data.U8[i] = (paramMask->data.U8[i]) ? 0 : 1;
+            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
         }
         psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, modelFunc);
         for (int i = 0; i < dparams->n; i++) {
-            if (!paramMask->data.U8[i])
+            if (!constraint->paramMask->data.U8[i])
                 continue;
             // note that delta is the value *subtracted* from the parameter
             // to get the new guess.  for dparams to represent the direction
             // of motion, we need to take -delta
-            dparams->data.F32[i] = -delta->data.F64[i];
+            dparams->data.F32[i] = -delta->data.F32[i];
         }
         psFree (delta);
@@ -251,9 +229,5 @@
     psFree(myMin);
     psFree(covar);
-    psFree(constrain->paramMask);
-    psFree(constrain->paramMin);
-    psFree(constrain->paramMax);
-    psFree(constrain->paramDelta);
-    psFree(constrain);
+    psFree(constraint);
 
     rc = (onPic && fitStatus);
@@ -266,25 +240,27 @@
 // these static variables are used by one pass of pmSourceFitSet to store
 // data for a model set.  If we are going to make psphot thread-safe, these
-// will have to go in a structure of their own.
+// will have to go in a structure of their own or be allocated once per thread
 // sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2,
 // nPar = nSrc*(nOnePar - 1) + 1
-static pmModelFunc mFunc;
-static int nPar;
+static pmModelFunc oneModelFunc;
+static pmModelLimits oneCheckLimits;
 static psVector *onePar;
 static psVector *oneDeriv;
-
-bool pmModelFitSetInit (pmModelType type)
-{
-
-    mFunc = pmModelFunc_GetFunction (type);
-    nPar  = pmModelParameterCount (type);
-
-    onePar = psVectorAlloc (nPar, PS_DATA_F32);
-    oneDeriv = psVectorAlloc (nPar, PS_DATA_F32);
+static int nOnePar;
+
+bool pmSourceFitSetInit (pmModelType type)
+{
+
+    oneModelFunc = pmModelFunc_GetFunction (type);
+    oneCheckLimits = pmModelLimits_GetFunction (type);
+    nOnePar = pmModelParameterCount (type);
+
+    onePar = psVectorAlloc (nOnePar, PS_DATA_F32);
+    oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32);
 
     return true;
 }
 
-void pmModelFitSetClear (void)
+void pmSourceFitSetClear (void)
 {
 
@@ -294,7 +270,19 @@
 }
 
-psF32 pmModelFitSet(psVector *deriv,
-                    const psVector *params,
-                    const psVector *x)
+bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)
+{
+    // convert the value of nParam into corresponding single model parameter entry
+    // convert params into corresponding single model parameter pointer
+
+    int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1;
+    float *paramSingle = params + nParam - nParamSingle;
+    float *betaSingle = betas + nParam - nParamSingle;
+    bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle);
+    return status;
+}
+
+psF32 pmSourceFitSet_Function(psVector *deriv,
+                              const psVector *params,
+                              const psVector *x)
 {
 
@@ -308,17 +296,17 @@
     psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;
 
-    int nSrc = (params->n - 1) / (nPar - 1);
+    int nSrc = (params->n - 1) / (nOnePar - 1);
 
     PAR[0] = model = pars[0];
     for (int i = 0; i < nSrc; i++) {
-        int nOff = i*nPar - i;
-        for (int n = 1; n < nPar; n++) {
+        int nOff = i*nOnePar - i;
+        for (int n = 1; n < nOnePar; n++) {
             PAR[n] = pars[nOff + n];
         }
         if (deriv == NULL) {
-            value = mFunc (NULL, onePar, x);
+            value = oneModelFunc (NULL, onePar, x);
         } else {
-            value = mFunc (oneDeriv, onePar, x);
-            for (int n = 1; n < nPar; n++) {
+            value = oneModelFunc (oneDeriv, onePar, x);
+            for (int n = 1; n < nOnePar; n++) {
                 dpars[nOff + n] = dPAR[n];
             }
@@ -333,5 +321,5 @@
 
 /*
-i:         0           1               2
+i:       0                   1                 2
 n:         1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6
 i*6 + n: 0 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
@@ -342,7 +330,6 @@
                      pmSourceFitMode mode)
 {
-    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     PS_ASSERT_PTR_NON_NULL(source->mask, false);
@@ -353,10 +340,51 @@
     psBool rc        = true;
 
-    // base values on first model
+    // maximum number of valid pixels
+    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAllocEmpty(nPix);
+    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
+    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
+
+    // fill in the coordinate and value entries
+    nPix = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            // skip masked points
+            if (source->mask->data.U8[i][j]) {
+                continue;
+            }
+            // skip zero-weight points
+            if (source->weight->data.F32[i][j] == 0) {
+                continue;
+            }
+            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+            // Convert i/j to image space:
+            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+            x->data[nPix] = (psPtr *) coord;
+            y->data.F32[nPix] = source->pixels->data.F32[i][j];
+
+            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
+            // as weight to avoid the bias from systematic errors here we would just use the
+            // source sky variance
+            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
+                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+            } else {
+                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
+            }
+            nPix++;
+        }
+    }
+    x->n = nPix;
+    y->n = nPix;
+    yErr->n = nPix;
+
+    // base values on first model (all models must be identical)
     pmModel *model = modelSet->data[0];
 
-    // set the static variables
-    pmModelFitSetInit (model->type);
-
+    // determine number of model parameters
     int nSrc = modelSet->n;
     int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
@@ -366,24 +394,15 @@
     psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
 
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-
-    // define limits for a single-source model
-    psVector *oneDelta;
-    psVector *oneParMin;
-    psVector *oneParMax;
-    modelLimits (&oneDelta, &oneParMin, &oneParMax);
-
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    constrain->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
-    constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-
-    // set the parameter guesses and constraints for the multiple models
-    // first the values for the single sky parameter
+    // set the static variables
+    pmSourceFitSetInit (model->type);
+
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+    constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
+    constraint->checkLimits = pmSourceFitSet_CheckLimits;
+
+    // set the parameter guesses for the multiple models
+    // first the value for the single sky parameter
     params->data.F32[0] = model->params->data.F32[0];
-    constrain->paramMin->data.F32[0]   = oneParMin->data.F32[0];
-    constrain->paramMax->data.F32[0]   = oneParMax->data.F32[0];
-    constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0];
 
     // next, the values for the source parameters
@@ -392,16 +411,10 @@
         for (int n = 1; n < nPar + 1; n++) {
             params->data.F32[i*nPar + n] = model->params->data.F32[n];
-            constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n];
-            constrain->paramMin->data.F32[i*nPar + n]   = oneParMin->data.F32[n];
-            constrain->paramMax->data.F32[i*nPar + n]   = oneParMax->data.F32[n];
-        }
-    }
-    psFree (oneDelta);
-    psFree (oneParMin);
-    psFree (oneParMax);
+        }
+    }
 
     if (psTraceGetLevel("psModules.objects") >= 5) {
         for (int i = 0; i < params->n; i++) {
-            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]);
+            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
         }
     }
@@ -413,7 +426,7 @@
         // NORM-only model fits only source normalization (Io)
         nParams = nSrc;
-        psVectorInit (constrain->paramMask, 1);
+        psVectorInit (constraint->paramMask, 1);
         for (int i = 0; i < nSrc; i++) {
-            constrain->paramMask->data.U8[1 + i*nPar] = 0;
+            constraint->paramMask->data.U8[1 + i*nPar] = 0;
         }
         break;
@@ -421,9 +434,9 @@
         // PSF model only fits x,y,Io
         nParams = 3*nSrc;
-        psVectorInit (constrain->paramMask, 1);
+        psVectorInit (constraint->paramMask, 1);
         for (int i = 0; i < nSrc; i++) {
-            constrain->paramMask->data.U8[1 + i*nPar] = 0;
-            constrain->paramMask->data.U8[2 + i*nPar] = 0;
-            constrain->paramMask->data.U8[3 + i*nPar] = 0;
+            constraint->paramMask->data.U8[1 + i*nPar] = 0;
+            constraint->paramMask->data.U8[2 + i*nPar] = 0;
+            constraint->paramMask->data.U8[3 + i*nPar] = 0;
         }
         break;
@@ -431,14 +444,6 @@
         // EXT model fits all params (except sky)
         nParams = nPar*nSrc;
-        psVectorInit (constrain->paramMask, 0);
-        constrain->paramMask->data.U8[0] = 1;
-        break;
-    case PM_SOURCE_FIT_PSF_AND_SKY:
-        nParams = 1 + 3*nSrc;
-        psAbort ("pmSourceFitModel", "PSF + SKY not currently available");
-        break;
-    case PM_SOURCE_FIT_EXT_AND_SKY:
-        nParams = 1 + nPar*nSrc;
-        psAbort ("pmSourceFitModel", "EXT + SKY not currently available");
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.U8[0] = 1;
         break;
     default:
@@ -448,48 +453,8 @@
     // force the floating parameters to fall within the contraint ranges
     for (int i = 0; i < params->n; i++) {
-        if (constrain->paramMask->data.U8[i])
-            continue;
-        params->data.F32[i] = PS_MIN(constrain->paramMax->data.F32[i], params->data.F32[i]);
-        params->data.F32[i] = PS_MAX(constrain->paramMin->data.F32[i], params->data.F32[i]);
-    }
-
-    // maximum number of valid pixels
-    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
-
-    // local sky variance
-    psF32 dSky = source->moments->dSky;
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAllocEmpty(nPix);
-    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
-    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
-
-    nPix = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            // skip masked points
-            if (source->mask->data.U8[i][j]) {
-                continue;
-            }
-            // skip zero-weight points
-            if (source->weight->data.F32[i][j] == 0) {
-                continue;
-            }
-            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-
-            // Convert i/j to image space:
-            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-            x->data[nPix] = (psPtr *) coord;
-            y->data.F32[nPix] = source->pixels->data.F32[i][j];
-            // psMinimizeLMChi2 takes wt = 1/dY^2
-            if (PM_PSF_POISSON_WEIGHTS) {
-                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
-            } else {
-                yErr->data.F32[nPix] = 1.0 / dSky;
-            }
-            nPix++;
-        }
-    }
+        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+    }
+
     if (nPix <  nParams + 1) {
         psTrace (__func__, 4, "insufficient valid pixels\n");
@@ -501,25 +466,15 @@
         psFree (params);
         psFree (dparams);
-        psFree(constrain->paramMask);
-        psFree(constrain->paramMin);
-        psFree(constrain->paramMax);
-        psFree(constrain->paramDelta);
-        psFree(constrain);
-        pmModelFitSetClear ();
+        psFree(constraint);
+        pmSourceFitSetClear ();
         return(false);
     }
-    x->n = nPix;
-    y->n = nPix;
-    yErr->n = nPix;
-
-    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
-                            PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-
-    psTrace (__func__, 5, "fitting function\n");
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet);
+
+    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function);
     if (!fitStatus) {
-        // psError(PS_ERR_UNKNOWN, false, "Failed to fit model (%d)\n", nSrc);
         psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);
     }
@@ -527,25 +482,25 @@
     // parameter errors from the covariance matrix
     for (int i = 0; i < dparams->n; i++) {
-        if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i])
+        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
             continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
     }
 
     // get the Gauss-Newton distance for fixed model parameters
-    if (constrain->paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+    if (constraint->paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
         psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
         altmask->data.U8[0] = 1;
         for (int i = 1; i < dparams->n; i++) {
-            altmask->data.U8[i] = (constrain->paramMask->data.U8[i]) ? 0 : 1;
-        }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmModelFitSet);
+            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
+        }
+        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function);
         for (int i = 0; i < dparams->n; i++) {
-            if (!constrain->paramMask->data.U8[i])
+            if (!constraint->paramMask->data.U8[i])
                 continue;
             // note that delta is the value *subtracted* from the parameter
             // to get the new guess.  for dparams to represent the direction
             // of motion, we need to take -delta
-            dparams->data.F32[i] = -delta->data.F64[i];
+            dparams->data.F32[i] = -delta->data.F32[i];
         }
         psFree (delta);
@@ -558,7 +513,12 @@
         model->params->data.F32[0] = params->data.F32[0];
         for (int n = 1; n < nPar + 1; n++) {
+            if (psTraceGetLevel("psModules.objects") >= 4) {
+                fprintf (stderr, "%f ", params->data.F32[i*nPar + n]);
+            }
             model->params->data.F32[n] = params->data.F32[i*nPar + n];
             model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
         }
+        psTrace ("psModules.objects", 4, " src %d", i);
+
         // save the resulting chisq, nDOF, nIter
         // these are not unique for any one source
@@ -571,12 +531,13 @@
 
         // models can go insane: reject these
-        onPic &= (model->params->data.F32[2] >= source->pixels->col0);
-        onPic &= (model->params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
-        onPic &= (model->params->data.F32[3] >= source->pixels->row0);
-        onPic &= (model->params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->row0 + source->pixels->numRows);
         if (!onPic) {
             model->status = PM_MODEL_OFFIMAGE;
         }
     }
+    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
 
     source->mode |= PM_SOURCE_MODE_FITTED;
@@ -587,18 +548,14 @@
     psFree(myMin);
     psFree(covar);
-    psFree(constrain->paramMask);
-    psFree(constrain->paramMin);
-    psFree(constrain->paramMax);
-    psFree(constrain->paramDelta);
-    psFree(constrain);
+    psFree(constraint);
     psFree(params);
     psFree(dparams);
 
-    // free static memory used by pmModelFitSet
-    pmModelFitSetClear ();
+    // free static memory used by pmSourceFitSet
+    pmSourceFitSetClear ();
 
     rc = (onPic && fitStatus);
     psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
-    psTrace("psModules.objects", 3, "---- %s end : status %d ----\n", __func__, rc);
+    psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);
     return(rc);
 }