Changeset 6848

Timestamp:

Apr 12, 2006, 8:18:46 PM (20 years ago)

Author:

magnier

Message:

additions to objects for further flexibility with options

Location:

branches/rel10_ifa/psModules/src

Files:

• objects/Makefile.am (modified) (2 diffs)
• objects/pmPSF.c (modified) (6 diffs)
• objects/pmPSF.h (modified) (4 diffs)
• objects/pmPSFtry.c (modified) (7 diffs)
• objects/pmPSFtry.h (modified) (3 diffs)
• objects/pmSourceFitModel.c (modified) (7 diffs)
• objects/pmSourceFitModel.h (modified) (3 diffs)
• objects/pmSourceFitSet.c (modified) (3 diffs)
• objects/pmSourceFitSet.h (modified) (2 diffs)
• objects/pmSourceIO_RAW.c (modified) (3 diffs)
• pslib/psAdditionals.h (modified) (1 diff)
• psmodules.h (modified) (1 diff)

branches/rel10_ifa/psModules/src/objects/Makefile.am

@@ -12 +12 @@
      pmSourceContour.c \
      pmSourceFitModel.c \
-     pmSourceFitSet.c \
      pmSourcePhotometry.c \
      pmSourceIO.c \
@@ -42 +41 @@
      pmSourceContour.h \
      pmSourceFitModel.h \
-     pmSourceFitSet.h \
      pmSourcePhotometry.h \
      pmSourceIO.h \

branches/rel10_ifa/psModules/src/objects/pmPSF.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-10 20:21:36 $
+ *  @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -80 +80 @@
  Object Normalization
  *****************************************************************************/
-pmPSF *pmPSFAlloc (pmModelType type)
+pmPSF *pmPSFAlloc (pmModelType type, bool poissonErrors)
 {
     int Nparams;
@@ -92 +92 @@
     psf->skyBias  = 0.0;
     psf->skySat   = 0.0;
+    psf->poissonErrors = poissonErrors;
 
     // the ApTrend components are (x, y, r2rflux, flux)
@@ -97 +98 @@
     pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
 
-    if (PM_PSF_POISSON_WEIGHTS) {
+    if (psf->poissonErrors) {
         psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     } else {
@@ -319 +320 @@
     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
 
+    psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors);
+
     int nPar = pmModelParameterCount (psf->type)    ;
     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
@@ -347 +350 @@
     pmModelType type = pmModelSetType (modelName);
 
-    pmPSF *psf = pmPSFAlloc (type);
+    bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS");
+    if (!status)
+        poissonErrors = true;
+
+    pmPSF *psf = pmPSFAlloc (type, poissonErrors);
 
     int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");

branches/rel10_ifa/psModules/src/objects/pmPSF.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.22.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-10 20:21:36 $
+ *  @version $Revision: 1.1.22.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -15 +15 @@
 # ifndef PM_PSF_H
 # define PM_PSF_H
-
-# define PM_PSF_POISSON_WEIGHTS 1
 
 /** pmPSF data structure
@@ -44 +42 @@
     int nPSFstars;   ///< number of stars used to measure PSF
     int nApResid;   ///< number of stars used to measure ApResid
+    bool poissonErrors;
 }
 pmPSF;
@@ -65 +64 @@
  */
 pmPSF *pmPSFAlloc(
-    pmModelType type                    ///< Add comment
+    pmModelType type,   // type of model for PSF
+    bool poissonErrors   ///< use poissonian errors or not?
 );
 

branches/rel10_ifa/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-10 20:21:36 $
+ *  @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -50 +50 @@
 
 // allocate a pmPSFtry based on the desired sources and the model (identified by name)
-pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName)
+pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors)
 {
 
@@ -59 +59 @@
     // XXX probably need to increment ref counter
     type           = pmModelSetType (modelName);
-    test->psf      = pmPSFAlloc (type);
+    test->psf      = pmPSFAlloc (type, poissonErrors);
     test->sources  = psMemIncrRefCounter(sources);
     test->modelEXT = psArrayAlloc (sources->n);
@@ -88 +88 @@
 // mask values indicate the reason the source was rejected:
 
-pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS)
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors)
 {
     bool status;
@@ -98 +98 @@
     int Npsf = 0;
 
-    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName);
+    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors);
 
     // stage 1:  fit an independent model (freeModel) to all sources
@@ -113 +113 @@
 
         psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED);
-        status = pmSourceFitModel (source, model, false);
+        status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT);
         psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED);
 
@@ -147 +147 @@
 
         psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED);
-        status = pmSourceFitModel (source, modelPSF, true);
+        status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF);
 
         // skip poor fits

branches/rel10_ifa/psModules/src/objects/pmPSFtry.h

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-17 17:13:42 $
+ *  @version $Revision: 1.3.4.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -80 +80 @@
 pmPSFtry *pmPSFtryAlloc(
     psArray *stars,                     ///< Add comment.
-    char *modelName                     ///< Add comment.
+    char *modelName,                     ///< Add comment.
+    bool poissonErrors   // use poissonian or constant errors?
 );
 
@@ -94 +95 @@
     psArray *sources,                   ///< Add comment.
     char *modelName,                    ///< Add comment.
-    float radius                        ///< Add comment.
+    float radius,                     ///< Add comment.
+    bool poissonErrors   // use poissonian or constant errors?
 );
 

branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-10 20:21:36 $
+ *  @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28 +28 @@
 #include "pmSourceFitModel.h"
 
-// save a static values so they may be set externally
+// 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;
-
-bool pmSourceFitModelInit (float nIter, float tol)
+static bool PM_PSF_POISSON_WEIGHTS = true;
+
+bool pmSourceFitModelInit (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;
 }
@@ -42 +45 @@
 bool pmSourceFitModel (pmSource *source,
                        pmModel *model,
-                       const bool PSF)
+                       pmSourceFitMode mode)
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -65 +68 @@
     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
 
-    int nParams = PSF ? 4 : params->n;
     psF32 dSky = source->moments->dSky;
 
@@ -110 +112 @@
     y->n = nPix;
     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();
+
+    // set parameter mask based on fitting mode
+    paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+    psVectorInit (paramMask, 1);
+
+    int nParams = 0;
+    switch (mode) {
+    case PM_SOURCE_FIT_NORM:
+        // NORM-only model fits only source normalization (Io)
+        nParams = 1;
+        paramMask->data.U8[1] = 0;
+        break;
+    case PM_SOURCE_FIT_PSF:
+        // PSF model only fits x,y,Io
+        nParams = 3;
+        paramMask->data.U8[1] = 0;
+        paramMask->data.U8[2] = 0;
+        paramMask->data.U8[3] = 0;
+        break;
+    case PM_SOURCE_FIT_EXT:
+        // EXT model fits all params (except sky)
+        nParams = params->n - 1;
+        psVectorInit (paramMask, 0);
+        paramMask->data.U8[0] = 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");
+        break;
+    default:
+        psAbort ("pmSourceFitModel", "invalid fitting mode");
+    }
+    constrain->paramMask = paramMask;
+
     if (nPix <  nParams + 1) {
         psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
@@ -117 +165 @@
         psFree (y);
         psFree (yErr);
+        psFree (myMin);
+        psFree (constrain);
+        psFree (paramMask);
         return(false);
     }
-
-    // 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();
-
-    // PSF model only fits x,y,Io, EXT model fits all (both skip sky)
-    paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-    paramMask->data.U8[0] = 1;
-    if (PSF) {
-        for (int i = 1; i < 4; i++) {
-            paramMask->data.U8[i] = 0;
-        }
-        for (int i = 4; i < paramMask->n; i++) {
-            paramMask->data.U8[i] = 1;
-        }
-    } else {
-        for (int i = 1; i < paramMask->n; i++) {
-            paramMask->data.U8[i] = 0;
-        }
-    }
-    constrain->paramMask = paramMask;
 
     // Set the parameter range checks
@@ -210 +236 @@
     return(rc);
 }
+
+/********************* 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.
+// 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;
+
+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
+*/
+
+bool pmSourceFitSet (pmSource *source,
+                     psArray *modelSet,
+                     pmSourceFitMode mode)
+{
+    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 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
+    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
+    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];
+            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(__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]);
+        }
+    }
+
+    // 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 (constrain->paramMask, 1);
+        for (int i = 0; i < nSrc; i++) {
+            constrain->paramMask->data.U8[1 + i*nPar] = 0;
+        }
+        break;
+    case PM_SOURCE_FIT_PSF:
+        // PSF model only fits x,y,Io
+        nParams = 3*nSrc;
+        psVectorInit (constrain->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;
+        }
+        break;
+    case PM_SOURCE_FIT_EXT:
+        // 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");
+        break;
+    default:
+        psAbort ("pmSourceFitModel", "invalid fitting mode");
+    }
+
+    // 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);
+}
+

branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-09 03:14:23 $
+ *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -13 +13 @@
 # define PM_SOURCE_FIT_MODEL_H
 
+typedef enum {
+    PM_SOURCE_FIT_NORM,
+    PM_SOURCE_FIT_PSF,
+    PM_SOURCE_FIT_EXT,
+    PM_SOURCE_FIT_PSF_AND_SKY,
+    PM_SOURCE_FIT_EXT_AND_SKY
+} pmSourceFitMode;
+
 bool pmSourceFitModelInit(
     float nIter,   ///< max number of allowed iterations
-    float tol      ///< convergence criterion
+    float tol,      ///< convergence criterion
+    bool poissonErrors   // use poisson errors for fits?
 );
 
 /** pmSourceFitModel()
  *
- * 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.
+ * 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.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image
+ * data. The function returns TRUE on success or FALSE on failure.
  *
  */
@@ -30 +38 @@
     pmSource *source,   ///< The input pmSource
     pmModel *model,   ///< model to be fitted
-    const bool PSF   ///< Treat model as PSF or EXT?
+    pmSourceFitMode mode  ///< define parameters to be fitted
+);
+
+
+// initialize data for a group of object models
+bool pmModelFitSetInit (pmModelType type);
+
+// clear data for a group of object models
+void pmModelFitSetClear (void);
+
+// function used to fit a group of object models
+psF32 pmModelFitSet(psVector *deriv,
+                    const psVector *params,
+                    const psVector *x);
+
+/** 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.pixels and source.mask entries. 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
+    pmSourceFitMode mode  ///< define parameters to be fitted
 );
 

branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-10 20:21:36 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -33 +33 @@
 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)
 {
@@ -97 +113 @@
 */
 
+/*
 # define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
 # define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
+*/
 
 bool pmSourceFitSet (pmSource *source,

branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h

@@ -3 +3 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-07 06:33:35 $
+ *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 
 void pmModelFitSetClear (void);
+
+bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors);
 
 /** pmSourceFitSet()

branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-01 02:47:20 $
+ *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-04-13 06:18:46 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -91 +91 @@
             fprintf (f, "%9.6f ", dPAR[j]);
         }
-        fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
+        fprintf (f, ": %8.4f %2d %#5x %7.3f %7.3f  %7.1f %7.2f %4d %2d\n",
                  source[0].apMag, source[0].type, source[0].mode,
                  log10(model[0].chisq/model[0].nDOF),
+                 log10(model[0].chisqNorm/model[0].nDOF),
                  source[0].moments->SN,
                  model[0].radiusTMP,
@@ -146 +147 @@
             fprintf (f, "%9.6f ", dPAR[j]);
         }
-        fprintf (f, ": %7.4f  %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n",
+        fprintf (f, ": %7.4f  %2d %#5x %7.3f %7.3f  %7.1f %7.2f %4d %2d\n",
                  source->apMag,
                  source[0].type, source[0].mode,
                  log10(model[0].chisq/model[0].nDOF),
+                 log10(model[0].chisqNorm/model[0].nDOF),
                  source[0].moments->SN,
                  model[0].radiusTMP,

branches/rel10_ifa/psModules/src/pslib/psAdditionals.h

@@ -31 +31 @@
 
 // strip whitespace from head and tail of string
-int          psStringStrip (char *string);
+int psStringStrip (char *string);
 
 // write out header with NAXIS=0

branches/rel10_ifa/psModules/src/psmodules.h

@@ -63 +63 @@
 # include <pmSourceSky.h>
 # include <pmSourceFitModel.h>
-# include <pmSourceFitSet.h>
 # include <pmSourceContour.h>
 # include <pmGrowthCurve.h>