Changeset 41124 for branches/eam_branches/ipp-20191011

Timestamp:

Nov 10, 2019, 9:46:19 AM (7 years ago)

Author:

eugene

Message:

some updates to pmSourceFitModel tests

File:

• branches/eam_branches/ipp-20191011/psModules/test/objects/tap_pmSourceFitModel.c (modified) (11 diffs)

branches/eam_branches/ipp-20191011/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.
@@ -22 +26 @@
 
     // 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);
+    // exit (1);
 
     for (int i = 0; i < sizeof(sigma)/sizeof(float); i++) {
@@ -75 +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;
@@ -85 +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);
@@ -144 +152 @@
         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->variance->data.F32[j][i] = var;
@@ -157 +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;
-
-    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode          = PM_SOURCE_FIT_EXT;
-    fitOptions->covarFactor   = 1.0;
+    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
@@ -173 +186 @@
 
     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);
@@ -181 +195 @@
     psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(maskVal));
 
-    par1->data.F32[par1->n] = (source->modelEXT->params->data.F32[1] / guess->params->data.F32[1]);
+    // 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]);
@@ -196 +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;
+}