Changeset 6481 for trunk/psphot/src/pmSourceFitSet.c

Timestamp:

Feb 23, 2006, 6:16:11 PM (20 years ago)

Author:

eugene

Message:

fixed some errors with double sources, added second-pass linear PSF fit, multiple-depths

File:

• trunk/psphot/src/pmSourceFitSet.c (modified) (7 diffs)

trunk/psphot/src/pmSourceFitSet.c

@@ -74 +74 @@
     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);
-    psVector *paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL;
+
+    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];
+    constrain->paramMask->data.U8[0]   = 0;
     for (int i = 0; i < nSrc; i++) {
       model = modelSet->data[i];
@@ -85 +105 @@
         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) {
-            paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1;
+            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]);
+        }
     }
 
@@ -126 +158 @@
     }
     if (nPix <  nParams + 1) {
-        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
-        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        psTrace (__func__, 4, "insufficient valid pixels\n");
+        psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__);
         model->status = PM_MODEL_BADARGS;
         psFree (x);
@@ -134 +166 @@
         psFree (params);
         psFree (dparams);
-        psFree (paramMask);
+        psFree(constrain->paramMask);
+        psFree(constrain->paramMin);
+        psFree(constrain->paramMax);
+        psFree(constrain->paramDelta);
+        psFree(constrain);
         return(false);
     }
@@ -143 +179 @@
     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
                             PM_SOURCE_FIT_MODEL_TOLERANCE);
-    psMinConstrain *constrain = psMinConstrainAlloc();
-    constrain->paramMask = paramMask;
-
-    // Set the parameter range checks
-    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-    modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
 
-    psTrace (".pmObjects.pmSourceFitSet", 5, "fitting function\n");
+    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 ((paramMask != NULL) && paramMask->data.U8[i])
+        if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i])
             continue;
         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
@@ -163 +193 @@
 
     // get the Gauss-Newton distance for fixed model parameters
-    if (paramMask != NULL) {
+    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 (!paramMask->data.U8[i])
+            if (!constrain->paramMask->data.U8[i])
                 continue;
             dparams->data.F32[i] = delta->data.F64[i];
@@ -220 +250 @@
 
     rc = (onPic && fitStatus);
-    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
+    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);
 }