Changeset 34071 for branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c

Timestamp:

Jun 25, 2012, 5:15:08 PM (14 years ago)

Author:

eugene

Message:

some tests to check working analysis; add some additional stats to psphotFitSourcesLinear

File:

• branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c (modified) (10 diffs)

branches/eam_branches/ipp-20120601/psphot/src/psphotFitSourcesLinear.c

@@ -273 +273 @@
             source->mode |= PM_SOURCE_MODE_PSFMODEL;
         }           
-        
+
         psArrayAdd (fitSources, 100, source);
     }
@@ -308 +308 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-        psSparseMatrixElement (sparse, i, i, f);
+        float MM = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
+        psSparseMatrixElement (sparse, i, i, MM);
 
         // if we have used CONSTANT errors, then we need to re-calculate the value of the parameter error
-        if (fitVarMode == PM_SOURCE_PHOTFIT_CONST) {
+        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
             float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
-            errors->data.F32[i] = 1.0 / sqrt(f);
+            errors->data.F32[i] = 1.0 / sqrt(MM);
         }
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
-        psSparseVectorElement (sparse, i, f);
+        float FM = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
+        psSparseVectorElement (sparse, i, FM);
 
         // add the per-source variances (border region)
@@ -435 +435 @@
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
+    // mean stats on fit windows
+    float sumRadius = 0.0;
+    float sumPixels = 0.0;
+    float sumSource = 0.0;
+
     // measure chisq for each source
     // for (int i = 0; final && (i < fitSources->n); i++) {
@@ -440 +445 @@
         pmSource *source = fitSources->data[i];
         pmModel *model = pmSourceGetModel (NULL, source);
+
+        // accumulate fit windows statistics
+        sumRadius += model->fitRadius;
+        sumPixels += M_PI*PS_SQR(model->fitRadius);
+        sumSource += 1.0;
+
         if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
             model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
@@ -446 +457 @@
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
+
+    float meanRadius = sumRadius / sumSource;
+    float meanPixels = sumPixels / sumSource;
 
     // psFree (index);
@@ -455 +469 @@
     psFree (border);
 
-    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
+    psLogMsg ("psphot.ensemble", PS_LOG_WARN, "measure ensemble of PSFs (mean radius = %f pixels, mean area = %f pixels: %f sec\n", meanRadius, meanPixels, psTimerMark ("psphot.linear"));
 
     psphotVisualPlotChisq (sources);
@@ -466 +480 @@
     return true;
 }
-
-// XXX do we need this?
-// XXX disallow the simultaneous sky fit and remove this code...
 
 // Calculate the weight terms for the sky fit component of the matrix.  This function operates
@@ -558 +569 @@
     );
 
-bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
+ bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, pmFPAfile *file, int index, psMetadata *recipe, pmReadout *readout, psArray *sources) {
 
     bool status = false;
@@ -599 +610 @@
         return false;
     }
-
     psFree (varModel);
     psFree (varModelStdev);
@@ -612 +622 @@
     }
     if (gain > 2.0) { /* warn? */ }
+    // XXX we are not actually using the gain, but need to test it to avoid gcc pedantic warnings
 
     // insert all of the source models