Changeset 41175

Timestamp:

Nov 27, 2019, 12:07:48 PM (6 years ago)

Author:

eugene

Message:

tweaking pmSourceFitModel tests

Location:

trunk/psModules/test/objects

Files:

• Makefile.am (modified) (1 diff)
• tap_pmSourceFitModel.c (modified) (11 diffs)
• tap_pmSourceFitModelBasic.c (copied) (copied from branches/eam_branches/ipp-20191011/psModules/test/objects/tap_pmSourceFitModelBasic.c )

trunk/psModules/test/objects/Makefile.am

@@ -12 +12 @@
 
 TEST_PROGS = \
-        tap_pmPeaks \
-        tap_pmGrowthCurve \
-        tap_pmMoments \
-        tap_pmSource \
-        tap_pmModel \
-        tap_pmModelUtils \
-        tap_pmModelClass \
-        tap_pmModel_CentralPixel \
-        tap_pmModel_CentralPixel_v2 \
-        tap_pmModel_SET_FWHM \
-        tap_pmPSF \
-        tap_pmTrend2D \
-        tap_pmResiduals \
-        tap_pmSourceExtendedPars \
-        tap_pmSourceSky \
-        tap_pmSourceMoments \
-        tap_pmSourceIO_PS1_DEV_0 \
-        tap_pmSourceIO_PS1_DEV_1 \
-        tap_pmSourceIO_SMPDATA \
-        tap_pmSourceContour \
-        tap_pmSourceUtils \
-        tap_pmSourceFitSet \
-        tap_pmPSF_IO \
-        tap_pmSourceIO_SMPDATA
+        tap_pmSourceFitModel
+
+#       tap_pmPeaks \
+#       tap_pmGrowthCurve \
+#       tap_pmMoments \
+#       tap_pmSource \
+#       tap_pmModel \
+#       tap_pmModelUtils \
+#       tap_pmModelClass \
+#       tap_pmModel_CentralPixel \
+#       tap_pmModel_CentralPixel_v2 \
+#       tap_pmModel_SET_FWHM \
+#       tap_pmPSF \
+#       tap_pmTrend2D \
+#       tap_pmResiduals \
+#       tap_pmSourceExtendedPars \
+#       tap_pmSourceSky \
+#       tap_pmSourceMoments \
+#       tap_pmSourceIO_PS1_DEV_0 \
+#       tap_pmSourceIO_PS1_DEV_1 \
+#       tap_pmSourceIO_SMPDATA \
+#       tap_pmSourceContour \
+#       tap_pmSourceUtils \
+#       tap_pmSourceFitSet \
+#       tap_pmPSF_IO \
+#       tap_pmSourceIO_SMPDATA
 
 if BUILD_TESTS

trunk/psModules/test/objects/tap_pmSourceFitModel.c

@@ -7 +7 @@
 #include "pstap.h"
 
+# define SKY_NOISE 10.0
+# define READ_NOISE 0.0
+
 bool fitModels (psRandom *seed, float flux, float radius, float sigma);
 bool fitModelFlux (psRandom *seed, float flux, float radius, float sigma);
+bool fitModelFluxOne (psRandom *seed, float flux, float radius, float sigma);
 
 // tests to check accuracy of fitted models for a range of fit radii, sigma, and flux.
@@ -18 +22 @@
 int main (void)
 {
-    pmModelGroupInit ();
-    pmSourceFitModelInit (15, 0.01, 1.0, true);
+    pmModelClassInit ();
+    // pmSourceFitModelInit (15, 0.01, 1.0, true);
 
     // psTraceSetLevel ("psModules.objects.pmSourceFitModel", 10);
-    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 10);
+    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
 
     plan_tests(240);
@@ -28 +32 @@
     // build a gauss-deviate vector (mean = 0.0, sigma = 1.0)
     psRandom *seed = psRandomAllocSpecific (PS_RANDOM_TAUS, 0);
+
+    // try a single model
+    // fitModelFluxOne (seed, 10000.0, 10.0, 2.0);
+    // return exit_status();
 
     static float radius[] = {3.0, 5.0, 7.0, 10.0, 15.0, 25.0};
@@ -36 +44 @@
     bool status = fitModels (seed, 10000.0, 10.0, 2.0);
     skip_start (!status, 240, "*** BASIC MODEL FITTING FAILS! *** : skipping related tests");
+    // exit (1);
 
     for (int i = 0; i < sizeof(sigma)/sizeof(float); i++) {
@@ -74 +83 @@
 
     float signal = 2*M_PI*sigma*sigma*flux;
-    float noise = sqrt(signal + 4*M_PI*sigma*sigma*(100 + PS_SQR(5)));
+    float noise = sqrt(signal + 4*M_PI*sigma*sigma*(SKY_NOISE*SKY_NOISE + PS_SQR(READ_NOISE)));
     float dMag = noise / signal;
     float dPos = sigma * dMag;
@@ -84 +93 @@
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     psVectorStats (stats, par1, NULL, NULL, 0);
-    ok ((stats->sampleStdev/dMag < 2.0), "Io ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dMag);
+    ok ((stats->sampleStdev/dMag < 1.5), "Io ref/fit stdev: %e : %e sigma",  stats->sampleStdev, stats->sampleStdev/dMag);
     psVectorStats (stats, par2, NULL, NULL, 0);
-    ok ((stats->sampleStdev/dPos < 2.0), "Xo ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dPos);
+    ok ((stats->sampleStdev/dPos < 1.5), "Xo ref-fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dPos);
     psVectorStats (stats, par3, NULL, NULL, 0);
-    ok ((stats->sampleStdev/dPos < 2.0), "Yo ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dPos);
+    ok ((stats->sampleStdev/dPos < 1.5), "Yo ref-fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dPos);
     psVectorStats (stats, par4, NULL, NULL, 0);
-    ok ((stats->sampleStdev/dMag < 2.0), "Sx ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dMag);
+    ok ((stats->sampleStdev/dMag < 1.5), "Sx ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dMag);
     psVectorStats (stats, par5, NULL, NULL, 0);
-    ok ((stats->sampleStdev/dMag < 2.0), "Sy ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dMag);
+    ok ((stats->sampleStdev/dMag < 1.5), "Sy ref/fit stdev: %e : %e sigma", stats->sampleStdev, stats->sampleStdev/dMag);
 
     psFree (par1);
@@ -107 +116 @@
 bool fitModelFlux (psRandom *seed, float flux, float radius, float sigma)
 {
+
+  psImageMaskType maskVal = 0x01;
 
     psVector *rnd = psVectorAlloc (1000, PS_TYPE_F32);
@@ -128 +139 @@
     source->modelEXT->params->data.F32[6] = 0;
 
-    source->pixels = psImageAlloc (100, 100, PS_TYPE_F32);
-    source->weight = psImageAlloc (100, 100, PS_TYPE_F32);
-    source->mask   = psImageAlloc (100, 100, PS_TYPE_U8);
+    source->pixels   = psImageAlloc (100, 100, PS_TYPE_F32);
+    source->variance = psImageAlloc (100, 100, PS_TYPE_F32);
+    source->maskObj  = psImageAlloc (100, 100, PS_TYPE_IMAGE_MASK);
     psImageInit (source->pixels, 0.0);
-    psImageInit (source->weight, 0.0);
-    psImageInit (source->mask, 0);
+    psImageInit (source->variance, 0.0);
+    psImageInit (source->maskObj, 0);
 
     // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
-    pmModelAdd (source->pixels, source->mask, source->modelEXT, PM_MODEL_OP_FULL);
+    pmModelAdd (source->pixels, source->maskObj, source->modelEXT, PM_MODEL_OP_FULL, maskVal);
     int npix = 0;
     for (int j = 0; j < source->pixels->numRows; j++) {
         for (int i = 0; i < source->pixels->numCols; i++) {
             float flux = source->pixels->data.F32[j][i];
-            float var = flux + 100 + PS_SQR(5);
+            float var = flux + PS_SQR(SKY_NOISE) + PS_SQR(READ_NOISE);
             source->pixels->data.F32[j][i] += rnd->data.F32[npix]*sqrt(var);
-            source->weight->data.F32[j][i] = var;
+            source->variance->data.F32[j][i] = var;
             npix ++;
             if (npix == rnd->n)
@@ -154 +165 @@
     // psFitsClose (fits);
 
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode          = PM_SOURCE_FIT_PSF;
+    fitOptions->covarFactor   = 1.0;
+
     // save the original model, modify params
     pmModel *guess = pmModelCopy (source->modelEXT);
-    guess->params->data.F32[1] *= 0.9;
-    guess->params->data.F32[2] += 1.0;
-    guess->params->data.F32[3] -= 1.0;
-    guess->params->data.F32[4] *= 0.9;
-    guess->params->data.F32[5] *= 0.9;
-
-    psImageKeepCircle (source->mask, 50, 50, radius, "OR", PM_MASK_MARK);
-    bool status = pmSourceFitModel (source, guess, PM_SOURCE_FIT_EXT);
+    if (fitOptions->mode == PM_SOURCE_FIT_PSF) {
+      guess->params->data.F32[1] *= 0.9;
+      guess->params->data.F32[2] += 1.0;
+      guess->params->data.F32[3] -= 1.0;
+    } 
+    if (fitOptions->mode == PM_SOURCE_FIT_EXT) {
+      guess->params->data.F32[1] *= 0.9;
+      guess->params->data.F32[4] *= 0.9;
+      guess->params->data.F32[5] *= 0.9;
+    }
+
+    // use maskVal to exclude pixels outside the circle
+    psImageKeepCircle (source->maskObj, 50, 50, radius, "OR", maskVal);
+
+    bool status = pmSourceFitModel (source, guess, fitOptions, maskVal);
+    // fprintf (stderr, "Io: %8.1f %8.1f\n", source->modelEXT->params->data.F32[1], guess->params->data.F32[1]);
     if (!status) {
         psFree (rnd);
@@ -170 +193 @@
         return false;
     }
-    psImageKeepCircle (source->mask, 50, 50, radius, "AND", PS_NOT_U8(PM_MASK_MARK));
-
-    par1->data.F32[par1->n] = (source->modelEXT->params->data.F32[1] / guess->params->data.F32[1]);
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(maskVal));
+
+    // par1->data.F32[par1->n] = (guess->params->data.F32[1]);
+    par1->data.F32[par1->n] = (guess->params->data.F32[1] / source->modelEXT->params->data.F32[1]);
     par2->data.F32[par2->n] = (source->modelEXT->params->data.F32[2] - guess->params->data.F32[2]);
     par3->data.F32[par3->n] = (source->modelEXT->params->data.F32[3] - guess->params->data.F32[3]);
@@ -187 +211 @@
     psFree (source);
     psFree (guess);
+    psFree (fitOptions);
 
     return true;
 }
+
+bool fitModelFluxOne (psRandom *seed, float flux, float radius, float sigma)
+{
+
+  psImageMaskType maskVal = 0x01;
+
+  psVector *rnd = psVectorAlloc (1000, PS_TYPE_F32);
+  for (int i = 0; i < rnd->n; i++) {
+    rnd->data.F32[i] = psRandomGaussian (seed);
+  }
+
+  // construct a model
+  pmSource *source = pmSourceAlloc ();
+  source->moments = pmMomentsAlloc ();
+
+  pmModelType type = pmModelClassGetType ("PS_MODEL_GAUSS");
+  source->modelEXT = pmModelAlloc (type);
+
+  source->modelEXT->params->data.F32[0] = 0;
+  source->modelEXT->params->data.F32[1] = flux;
+  source->modelEXT->params->data.F32[2] = 50;
+  source->modelEXT->params->data.F32[3] = 50;
+  source->modelEXT->params->data.F32[4] = 2.0*sqrt(sigma);
+  source->modelEXT->params->data.F32[5] = 2.0*sqrt(sigma);
+  source->modelEXT->params->data.F32[6] = 0;
+
+  source->pixels   = psImageAlloc (100, 100, PS_TYPE_F32);
+  source->variance = psImageAlloc (100, 100, PS_TYPE_F32);
+  source->maskObj  = psImageAlloc (100, 100, PS_TYPE_IMAGE_MASK);
+  psImageInit (source->pixels, 0.0);
+  psImageInit (source->variance, 0.0);
+  psImageInit (source->maskObj, 0);
+
+  // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
+  pmModelAdd (source->pixels, source->maskObj, source->modelEXT, PM_MODEL_OP_FULL, maskVal);
+  int npix = 0;
+  for (int j = 0; j < source->pixels->numRows; j++) {
+    for (int i = 0; i < source->pixels->numCols; i++) {
+      float flux = source->pixels->data.F32[j][i];
+      float var = flux + 100 + PS_SQR(5);
+      source->pixels->data.F32[j][i] += rnd->data.F32[npix]*sqrt(var);
+      source->variance->data.F32[j][i] = var;
+      npix ++;
+      if (npix == rnd->n)
+        npix = 0;
+    }
+  }
+
+  // psFits *fits = psFitsOpen ("test.fits", "w");
+  // psFitsWriteImage (fits, NULL, source->pixels, 0, NULL);
+  // psFitsClose (fits);
+
+  // EXT vs PSF fitting:
+  // EXT fits Io, Sxx, Sxy, Syy
+  // PSF fits Io, Xo, Yo
+  pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+  fitOptions->mode = PM_SOURCE_FIT_EXT;
+  // fitOptions->mode = PM_SOURCE_FIT_PSF;
+  fitOptions->covarFactor   = 1.0;
+
+  // save the original model, modify params
+  pmModel *guess = pmModelCopy (source->modelEXT);
+  if (fitOptions->mode == PM_SOURCE_FIT_PSF) {
+    guess->params->data.F32[1] *= 0.9;
+    guess->params->data.F32[2] += 1.0;
+    guess->params->data.F32[3] -= 1.0;
+  } 
+  if (fitOptions->mode == PM_SOURCE_FIT_EXT) {
+    guess->params->data.F32[1] *= 0.9;
+    guess->params->data.F32[4] *= 0.9;
+    guess->params->data.F32[5] *= 0.9;
+  }
+
+  // use maskVal to exclude pixels outside the circle
+  psImageKeepCircle (source->maskObj, 50, 50, radius, "OR", maskVal);
+
+  // I need to use psTrace to get verbosity from the fit...
+  bool status = pmSourceFitModel (source, guess, fitOptions, maskVal);
+  if (!status) {
+    psFree (rnd);
+    psFree (source);
+    psFree (guess);
+    psFree (fitOptions);
+    fprintf (stderr, "failed\n");
+    return false;
+  }
+  psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(maskVal));
+
+  psFree (rnd);
+  psFree (source);
+  psFree (guess);
+  psFree (fitOptions);
+
+  return true;
+}