Index: trunk/psModules/src/objects/pmSourceFitModel.c
===================================================================
--- trunk/psModules/src/objects/pmSourceFitModel.c	(revision 13932)
+++ trunk/psModules/src/objects/pmSourceFitModel.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-21 22:58:11 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,8 +27,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceFitModel.h"
 
@@ -116,15 +116,8 @@
     psVector *dparams = model->dparams;
 
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    if (!modelFunc)
-        psAbort("invalid model function");
-    pmModelLimits checkLimits = pmModelLimits_GetFunction (model->type);
-    if (!checkLimits)
-        psAbort("invalid model limits function");
-
     // create the minimization constraints
     psMinConstraint *constraint = psMinConstraintAlloc();
     constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-    constraint->checkLimits = checkLimits;
+    constraint->checkLimits = model->modelLimits;
 
     // set parameter mask based on fitting mode
@@ -156,6 +149,6 @@
     // 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);
+        model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     }
 
@@ -174,5 +167,5 @@
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
 
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, modelFunc);
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if (psTraceGetLevel("psModules.objects") >= 4) {
@@ -201,5 +194,5 @@
             altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
         }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, modelFunc);
+        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, model->modelFunc);
         for (int i = 0; i < dparams->n; i++) {
             if (!constraint->paramMask->data.U8[i])
@@ -237,332 +230,2 @@
 }
 
-/********************* Source Model Set Functions ***************************/
-
-// 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 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 oneModelFunc;
-static pmModelLimits oneCheckLimits;
-static psVector *onePar;
-static psVector *oneDeriv;
-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 pmSourceFitSetClear (void)
-{
-
-    psFree (onePar);
-    psFree (oneDeriv);
-    return;
-}
-
-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)
-{
-
-    psF32 value;
-    psF32 model;
-
-    psF32 *PAR = onePar->data.F32;
-    psF32 *dPAR = oneDeriv->data.F32;
-
-    psF32 *pars = params->data.F32;
-    psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;
-
-    int nSrc = (params->n - 1) / (nOnePar - 1);
-
-    PAR[0] = model = pars[0];
-    for (int i = 0; i < nSrc; i++) {
-        int nOff = i*nOnePar - i;
-        for (int n = 1; n < nOnePar; n++) {
-            PAR[n] = pars[nOff + n];
-        }
-        if (deriv == NULL) {
-            value = oneModelFunc (NULL, onePar, x);
-        } else {
-            value = oneModelFunc (oneDeriv, onePar, x);
-            for (int n = 1; n < nOnePar; n++) {
-                dpars[nOff + n] = dPAR[n];
-            }
-        }
-        model += value;
-    }
-    if (deriv != NULL) {
-        dpars[0] = dPAR[0]*2.0;
-    }
-    return (model);
-}
-
-/*
-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
-*/
-
-bool pmSourceFitSet (pmSource *source,
-                     psArray *modelSet,
-                     pmSourceFitMode mode,
-                     psMaskType maskVal)
-{
-    psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-
-    // 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->maskObj->data.U8[i][j] & maskVal) {
-                continue;
-            }
-            // skip zero-weight points
-            if (source->weight->data.F32[i][j] == 0) {
-                continue;
-            }
-            // skip nan values in image
-            if (!isfinite(source->pixels->data.F32[i][j])) {
-                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];
-
-    // determine number of model parameters
-    int nSrc = modelSet->n;
-    int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
-
-    // define parameter vectors for source set
-    psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-
-    // 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];
-
-    // next, the values for the source parameters
-    for (int i = 0; i < nSrc; i++) {
-        model = modelSet->data[i];
-        for (int n = 1; n < nPar + 1; n++) {
-            params->data.F32[i*nPar + n] = model->params->data.F32[n];
-        }
-    }
-
-    if (psTraceGetLevel("psModules.objects") >= 5) {
-        for (int i = 0; i < params->n; i++) {
-            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
-        }
-    }
-
-    // set the parameter masks based on the fitting mode
-    int nParams = 0;
-    switch (mode) {
-    case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = nSrc;
-        psVectorInit (constraint->paramMask, 1);
-        for (int i = 0; i < nSrc; i++) {
-            constraint->paramMask->data.U8[1 + i*nPar] = 0;
-        }
-        break;
-    case PM_SOURCE_FIT_PSF:
-        // PSF model only fits x,y,Io
-        nParams = 3*nSrc;
-        psVectorInit (constraint->paramMask, 1);
-        for (int i = 0; i < nSrc; i++) {
-            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;
-    case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = nPar*nSrc;
-        psVectorInit (constraint->paramMask, 0);
-        constraint->paramMask->data.U8[0] = 1;
-        break;
-    default:
-        psAbort("invalid fitting mode");
-    }
-
-    // force the floating parameters to fall within the contraint ranges
-    for (int i = 0; i < params->n; i++) {
-        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");
-        psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        psFree (x);
-        psFree (y);
-        psFree (yErr);
-        psFree (params);
-        psFree (dparams);
-        psFree(constraint);
-        pmSourceFitSetClear ();
-        return(false);
-    }
-
-    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) {
-        psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);
-    }
-
-    // parameter errors from the covariance matrix
-    for (int i = 0; i < dparams->n; i++) {
-        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
-    }
-
-    // get the Gauss-Newton distance for fixed model parameters
-    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] = (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 (!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.F32[i];
-        }
-        psFree (delta);
-        psFree (altmask);
-    }
-
-    // assign back the parameters to the models
-    for (int i = 0; i < nSrc; i++) {
-        model = modelSet->data[i];
-        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
-        model->chisq = myMin->value;
-        model->nIter = myMin->iter;
-        model->nDOF  = y->n - nParams;
-
-        // set the model success or failure status
-        model->flags |= PM_MODEL_STATUS_FITTED;
-        if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
-
-        // models can go insane: reject these
-        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->flags |= PM_MODEL_STATUS_OFFIMAGE;
-    }
-    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
-
-    source->mode |= PM_SOURCE_MODE_FITTED;
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(constraint);
-    psFree(params);
-    psFree(dparams);
-
-    // 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", 5, "---- %s end (%d) ----\n", __func__, rc);
-    return(rc);
-}
-