Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/objects/models/pmModel_QGAUSS.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/objects/models/pmModel_QGAUSS.c (modified) (15 diffs)

branches/simtest_nebulous_branches/psModules/src/objects/models/pmModel_QGAUSS.c

@@ -1 +1 @@
 /******************************************************************************
  * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
- * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
+ * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
  * statements of their own.  The models use a psVector to represent the set of parameters, with
  * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
@@ -20 +20 @@
    *****************************************************************************/
 
+#include <stdio.h>
+#include <pslib.h>
+
+#include "pmMoments.h"
+#include "pmPeaks.h"
+#include "pmSource.h"
+#include "pmModel.h"
+#include "pmModel_QGAUSS.h"
+
 # define PM_MODEL_FUNC            pmModelFunc_QGAUSS
 # define PM_MODEL_FLUX            pmModelFlux_QGAUSS
@@ -28 +37 @@
 # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS
 # define PM_MODEL_FIT_STATUS      pmModelFitStatus_QGAUSS
+# define PM_MODEL_SET_LIMITS      pmModelSetLimits_QGAUSS
+
+# define ALPHA   2.250
+# define ALPHA_M 1.250
+
+// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
+// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
+// values need to be pixel coords
+
+// Lax parameter limits
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
+static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+
+// Moderate parameter limits
+// Tolerate a small divot (k < 0)
+static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
+static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+
+// Strict parameter limits
+// k = PAR_7 < 0 is very undesirable (big divot in the middle)
+static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
+static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
+
+// Parameter limits to use
+static float *paramsMinUse = paramsMinLax;
+static float *paramsMaxUse = paramsMaxLax;
+static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
+
+static bool limitsApply = true;         // Apply limits?
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -48 +86 @@
     assert (z >= 0);
 
-    psF32 zp = pow(z,1.25);
+    psF32 zp = pow(z,ALPHA_M);
     psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
 
@@ -59 +97 @@
         // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
         psF32 t = r1*r;
-        psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
+        psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp);
 
         dPAR[PM_PAR_SKY]  = +1.0;
@@ -79 +117 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
+
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
-    float beta_lim = 0, params_min = 0, params_max = 0;
-    float f1 = 0, f2 = 0, q1 = 0, q2 = 0;
+    if (!limitsApply) {
+        return true;
+    }
+    psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds");
 
     // we need to calculate the limits for SXY specially
+    float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
+        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2  = 0.5*sqrt (q1);
+        q2 = 0.5*sqrtf(q1);
     }
 
     switch (mode) {
-    case PS_MINIMIZE_BETA_LIMIT:
-        switch (nParam) {
-        case PM_PAR_SKY:
-            beta_lim = 1000;
-            break;
-        case PM_PAR_I0:
-            beta_lim = 3e6;
-            break;
-        case PM_PAR_XPOS:
-            beta_lim = 5;
-            break;
-        case PM_PAR_YPOS:
-            beta_lim = 5;
-            break;
-        case PM_PAR_SXX:
-            beta_lim = 1.0;
-            break;
-        case PM_PAR_SYY:
-            beta_lim = 1.0;
-            break;
-        case PM_PAR_SXY:
-            beta_lim =  0.5*q2;
-            break;
-        case PM_PAR_7:
-            beta_lim = 2.0;
-            break;
-        default:
-            psAbort("invalid parameter %d for beta test", nParam);
-        }
-        if (fabs(beta[nParam]) > fabs(beta_lim)) {
-            beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim);
-            psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
-                     nParam, beta[nParam], beta_lim);
-            return false;
-        }
-        return true;
-    case PS_MINIMIZE_PARAM_MIN:
-        switch (nParam) {
-        case PM_PAR_SKY:
-            params_min = -1000;
-            break;
-        case PM_PAR_I0:
-            params_min =   0.01;
-            break;
-        case PM_PAR_XPOS:
-            params_min =  -100;
-            break;
-        case PM_PAR_YPOS:
-            params_min =  -100;
-            break;
-        case PM_PAR_SXX:
-            params_min =   0.5;
-            break;
-        case PM_PAR_SYY:
-            params_min =   0.5;
-            break;
-        case PM_PAR_SXY:
-            params_min =  -q2;
-            break;
-        case PM_PAR_7:
-            params_min =   0.1;
-            break;
-        default:
-            psAbort("invalid parameter %d for param min test", nParam);
-        }
-        if (params[nParam] < params_min) {
-            params[nParam] = params_min;
-            psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
-                     nParam, params[nParam], params_min);
-            return false;
-        }
-        return true;
-    case PS_MINIMIZE_PARAM_MAX:
-        switch (nParam) {
-        case PM_PAR_SKY:
-            params_max =   1e5;
-            break;
-        case PM_PAR_I0:
-            params_max =   1e8;
-            break;
-        case PM_PAR_XPOS:
-            params_max =   1e4;
-            break;
-        case PM_PAR_YPOS:
-            params_max =   1e4;
-            break;
-        case PM_PAR_SXX:
-            params_max =   100;
-            break;
-        case PM_PAR_SYY:
-            params_max =   100;
-            break;
-        case PM_PAR_SXY:
-            params_max =  +q2;
-            break;
-        case PM_PAR_7:
-            params_max =  20.0;
-            break;
-        default:
-            psAbort("invalid parameter %d for param max test", nParam);
-        }
-        if (params[nParam] > params_max) {
-            params[nParam] = params_max;
-            psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
-                     nParam, params[nParam], params_max);
-            return false;
-        }
-        return true;
+      case PS_MINIMIZE_BETA_LIMIT: {
+          psAssert(beta, "Require beta to limit beta");
+          float limit = betaUse[nParam];
+          if (nParam == PM_PAR_SXY) {
+              limit *= q2;
+          }
+          if (fabs(beta[nParam]) > fabs(limit)) {
+              beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit);
+              psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g",
+                      nParam, beta[nParam], limit);
+              return false;
+          }
+          return true;
+      }
+      case PS_MINIMIZE_PARAM_MIN: {
+          psAssert(params, "Require parameters to limit parameters");
+          psAssert(paramsMinUse, "Require parameter limits to limit parameters");
+          float limit = paramsMinUse[nParam];
+          if (nParam == PM_PAR_SXY) {
+              limit *= q2;
+          }
+          if (params[nParam] < limit) {
+              params[nParam] = limit;
+              psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g",
+                      nParam, params[nParam], limit);
+              return false;
+          }
+          return true;
+      }
+      case PS_MINIMIZE_PARAM_MAX: {
+          psAssert(params, "Require parameters to limit parameters");
+          psAssert(paramsMaxUse, "Require parameter limits to limit parameters");
+          float limit = paramsMaxUse[nParam];
+          if (nParam == PM_PAR_SXY) {
+              limit *= q2;
+          }
+          if (params[nParam] > limit) {
+              params[nParam] = limit;
+              psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g",
+                      nParam, params[nParam], limit);
+              return false;
+          }
+          return true;
+      }
     default:
         psAbort("invalid choice for limits");
@@ -213 +191 @@
 
 // make an initial guess for parameters
+// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters
 bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
 {
@@ -237 +216 @@
     if (!isfinite(shape.sxy)) return false;
 
-    // XXX turn this off here for now PAR[PM_PAR_SKY]  = moments->Sky;
     PAR[PM_PAR_SKY]  = 0.0;
     PAR[PM_PAR_I0]   = peak->flux;
@@ -273 +251 @@
     float f1, f2;
     for (z = DZ; z < 50; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
@@ -290 +268 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -296 +274 @@
     psF32 *PAR = params->data.F32;
 
-    if (flux <= 0)
-        return (1.0);
-    if (PAR[PM_PAR_I0] <= 0)
-        return (1.0);
-    if (flux >= PAR[PM_PAR_I0])
-        return (1.0);
+    if (flux <= 0) return 1.0;
+    if (PAR[PM_PAR_I0] <= 0) return 1.0;
+    if (flux >= PAR[PM_PAR_I0]) return 1.0;
+    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
 
     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
@@ -317 +293 @@
 
     // choose a z value guaranteed to be beyond our limit
-    float z0 = pow((1.0 / limit), (1.0 / 2.25));
-    float z1 = (1.0 / limit) / PAR[PM_PAR_7];
-    z1 = PS_MAX (z0, z1);
-    z0 = 0.0;
-
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
-    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
+    float z0 = 0.0;
+    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
+    psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
+    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
+
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    // use Newton-Raphson to minimize f(z) - limit = 0
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+        float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+        float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
+
+        float f = 1.0 / q;
+        float dfdz = -dqdz * f / q;
+
+        dz = (f - limit) / dfdz;
+
+        // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+        z -= dz;
+        z = PS_MAX(z, 0.0);
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -444 +423 @@
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
-
-    psF32 dP;
     bool  status;
 
@@ -451 +428 @@
     psF32 *dPAR = model->dparams->data.F32;
 
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP = sqrt (dP);
-
     status = true;
-//    status &= (dP < 0.5);
     status &= (PAR[PM_PAR_I0] > 0);
     status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
 
     return status;
+}
+
+
+void PM_MODEL_SET_LIMITS(pmModelLimitsType type)
+{
+    switch (type) {
+      case PM_MODEL_LIMITS_NONE:
+        paramsMinUse = NULL;
+        paramsMaxUse = NULL;
+        limitsApply = true;
+        break;
+      case PM_MODEL_LIMITS_IGNORE:
+        paramsMinUse = NULL;
+        paramsMaxUse = NULL;
+        limitsApply = false;
+        break;
+      case PM_MODEL_LIMITS_LAX:
+        paramsMinUse = paramsMinLax;
+        paramsMaxUse = paramsMaxLax;
+        limitsApply = true;
+        break;
+      case PM_MODEL_LIMITS_MODERATE:
+        paramsMinUse = paramsMinModerate;
+        paramsMaxUse = paramsMaxModerate;
+        limitsApply = true;
+        break;
+      case PM_MODEL_LIMITS_STRICT:
+        paramsMinUse = paramsMinStrict;
+        paramsMaxUse = paramsMaxStrict;
+        limitsApply = true;
+        break;
+      default:
+        psAbort("Unrecognised model limits type: %x", type);
+    }
+    return;
 }
 
@@ -472 +478 @@
 # undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS
+# undef PM_MODEL_SET_LIMITS
+# undef ALPHA
+# undef ALPHA_M