Index: trunk/psModules/src/objects/pmSourceFitSet.c
===================================================================
--- trunk/psModules/src/objects/pmSourceFitSet.c	(revision 6945)
+++ 	(revision )
@@ -1,334 +1,0 @@
-/** @file  pmSourceFitSet.c
- *
- *  @author EAM, IfA
- *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- *
- */
-
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
-#include "pslib.h"
-#include "pmHDU.h"
-#include "pmFPA.h"
-#include "pmPeaks.h"
-#include "pmMoments.h"
-#include "pmGrowthCurve.h"
-#include "pmModel.h"
-#include "pmPSF.h"
-#include "pmSource.h"
-#include "pmModelGroup.h"
-#include "pmSourceFitSet.h"
-
-// 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 psVector *onePar;
-static psVector *oneDeriv;
-
-// save as static values so they may be set externally
-static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
-static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
-static bool PM_PSF_POISSON_WEIGHTS = true;
-
-bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors)
-{
-
-    PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;
-    PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
-
-    PM_PSF_POISSON_WEIGHTS = poissonErrors;
-
-    return true;
-}
-
-bool pmModelFitSetInit (pmModelType type)
-{
-
-    mFunc = pmModelFunc_GetFunction (type);
-    nPar  = pmModelParameterCount (type);
-
-    onePar = psVectorAlloc (nPar, PS_DATA_F32);
-    oneDeriv = psVectorAlloc (nPar, PS_DATA_F32);
-
-    return true;
-}
-
-void pmModelFitSetClear (void)
-{
-
-    psFree (onePar);
-    psFree (oneDeriv);
-    return;
-}
-
-psF32 pmModelFitSet(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) / (nPar - 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++) {
-            PAR[n] = pars[nOff + n];
-        }
-        if (deriv == NULL) {
-            value = mFunc (NULL, onePar, x);
-        } else {
-            value = mFunc (oneDeriv, onePar, x);
-            for (int n = 1; n < nPar; 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
-*/
-
-/*
-# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
-# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
-*/
-
-bool pmSourceFitSet (pmSource *source,
-                     psArray *modelSet,
-                     const bool PSF)
-{
-    psTrace(__func__, 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->peak, 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);
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-
-    // base values on first model
-    pmModel *model = modelSet->data[0];
-
-    // set the static variables
-    pmModelFitSetInit (model->type);
-
-    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);
-
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-
-    // define limits for single-source model
-    psVector *oneDelta;
-    psVector *oneParMin;
-    psVector *oneParMax;
-    modelLimits (&oneDelta, &oneParMin, &oneParMax);
-
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    constrain->paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
-
-    // all but the sky are allowed to vary independently (subject to PSF)
-    // set the parameters from the multiple models
-    // also, set limits based on single-source limits
-    params->data.F32[0] = model->params->data.F32[0];
-    constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0];
-    constrain->paramMin->data.F32[0]   = oneParMin->data.F32[0];
-    constrain->paramMax->data.F32[0]   = oneParMax->data.F32[0];
-
-    // XXX all models skip the sky
-    constrain->paramMask->data.U8[0]   = 1;
-    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];
-            dparams->data.F32[i*nPar + n] = model->dparams->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];
-            if (PSF) {
-                constrain->paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
-            }
-        }
-    }
-    psFree (oneDelta);
-    psFree (oneParMin);
-    psFree (oneParMax);
-
-    if (psTraceGetLevel(__func__) >= 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]);
-        }
-    }
-
-    // PSF model only fits first 4 parameters, EXT model fits all
-    int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1;
-
-    // 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 = psArrayAlloc(nPix);
-    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(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++;
-        }
-    }
-    if (nPix <  nParams + 1) {
-        psTrace (__func__, 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__);
-        model->status = PM_MODEL_BADARGS;
-        psFree (x);
-        psFree (y);
-        psFree (yErr);
-        psFree (params);
-        psFree (dparams);
-        psFree(constrain->paramMask);
-        psFree(constrain->paramMin);
-        psFree(constrain->paramMax);
-        psFree(constrain->paramDelta);
-        psFree(constrain);
-        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);
-
-    // parameter errors from the covariance matrix
-    for (int i = 0; i < dparams->n; i++) {
-        if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
-    }
-
-    // get the Gauss-Newton distance for fixed model parameters
-    if (constrain->paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
-        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!constrain->paramMask->data.U8[i])
-                continue;
-            dparams->data.F32[i] = delta->data.F64[i];
-        }
-        psFree (delta);
-    }
-
-    // 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++) {
-            model->params->data.F32[n] = params->data.F32[i*nPar + n];
-            model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
-        }
-        // 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->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE;
-
-        // 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);
-        if (!onPic) {
-            model->status = PM_MODEL_OFFIMAGE;
-        }
-    }
-
-    source->mode |= PM_SOURCE_MODE_FITTED;
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(constrain->paramMask);
-    psFree(constrain->paramMin);
-    psFree(constrain->paramMax);
-    psFree(constrain->paramDelta);
-    psFree(constrain);
-    psFree(params);
-    psFree(dparams);
-
-    // free static memory used by pmModelFitSet
-    pmModelFitSetClear ();
-
-    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(__func__, 3, "---- %s end : status %d ----\n", __func__, rc);
-    return(rc);
-}
-
Index: trunk/psModules/src/objects/pmSourceFitSet.h
===================================================================
--- trunk/psModules/src/objects/pmSourceFitSet.h	(revision 6945)
+++ 	(revision )
@@ -1,40 +1,0 @@
-/** @file  pmSourceFitset.h
- *
- *  @author EAM, IfA; GLG, MHPCC
- *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-17 18:01:05 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- *
- */
-
-# ifndef PM_SOURCE_FIT_SET_H
-# define PM_SOURCE_FIT_SET_H
-
-psF32 pmModelFitSet(psVector *deriv,
-                    const psVector *params,
-                    const psVector *x);
-
-bool pmModelFitSetInit (pmModelType type);
-
-void pmModelFitSetClear (void);
-
-bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors);
-
-/** pmSourceFitSet()
- *
- * Fit the requested model to the specified source. The starting guess for the
- * model is given by the input source.model parameter values. The pixels of
- * interest are specified by the source.pixelsand source.maskentries. This
- * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
- * on success or FALSE on failure.
- *
- */
-bool pmSourceFitSet(
-    pmSource *source,   ///< The input pmSource
-    psArray *modelSet,   ///< model to be fitted
-    const bool PSF   ///< Treat model as PSF or EXT?
-);
-
-# endif /* PM_SOURCE_FIT_MODEL_H */