Changeset 4630 for trunk/psphot/src/LocalSky.c

Timestamp:

Jul 27, 2005, 3:26:57 PM (21 years ago)

Author:

eugene

Message:

various tests and changes

File:

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

trunk/psphot/src/LocalSky.c

@@ -62 +62 @@
 bool pmSourceFitModel_EAM(psSource *source,
                           psModel *model,
-                          const bool PSF,
-                          float SOFT)
+                          const bool PSF)
 {
     PS_PTR_CHECK_NULL(source, false);
@@ -74 +73 @@
     psBool onPic     = true;
     psBool rc        = true;
-    psF32  Ro;
+    psF32  Ro, ymodel;
 
 
@@ -81 +80 @@
 
     psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
+    psModelLimits modelLimits = psModelLimits_GetFunction (model->type);
 
     psVector *params = model->params;
@@ -86 +86 @@
     psVector *paramMask = NULL;
 
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
     int nParams = PSF ? params->n - 4 : params->n;
-    psF32 Xo = params->data.F32[2];
-    psF32 Yo = params->data.F32[3];
+    psF32 So = params->data.F32[0];
 
     // find the number of valid pixels
@@ -122 +125 @@
                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
 
-                // XXX just for test purposes, use the raw radius.  a better
-                // choice is to ask what is the value of z and scale by that
-                Ro = hypot ((Xo - coord->data.F32[0]), (Yo - coord->data.F32[1]));
-
-                // XXX enhance the noise (doubled if R = 10 pix, etc. note the square)
-                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * (1 + PS_SQR(Ro/SOFT)));
+                ymodel = modelFunc (NULL, model->params, coord);
+                
+                // this test enhances the noise based on deviation from the model flux
+                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro);
                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then 
                 //           the minimization function calculates sq()
@@ -147 +149 @@
         paramMask->data.U8[i] = 1;
       }
-    }       
-
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    }  
+
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < params->n; i++) {
+        covar->data.F64[0][i] = beta_lim->data.F32[i];
+        covar->data.F64[1][i] = params_min->data.F32[i];
+        covar->data.F64[2][i] = params_max->data.F32[i];
+    }
 
     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
@@ -326 +334 @@
 }
 
+bool pmModelFitStatus (psModel *model) {
+
+    bool status;
+
+    psModelFitStatusFunc statusFunc = psModelFitStatusFunc_GetFunction (model->type);
+    status = statusFunc (model);
+
+    return (status);
+}
+
+// XXX EAM this implementation of MinLM includes limits on params & dparams
+psBool psMinimizeLMChi2_EAM(psMinimization *min,
+                            psImage *covar,
+                            psVector *params,
+                            const psVector *paramMask,
+                            const psArray *x,
+                            const psVector *y,
+                            const psVector *yErr,
+                            psMinimizeLMChi2Func func)
+{
+    PS_PTR_CHECK_NULL(min, NULL);
+    PS_VECTOR_CHECK_NULL(params, NULL);
+    PS_VECTOR_CHECK_EMPTY(params, NULL);
+    PS_PTR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_EMPTY(y, NULL);
+    PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
+    PS_PTR_CHECK_NULL(func, NULL);
+
+    // this function has test and current values for several things
+    // the current best value is in lower case
+    // the next guess value is in upper case
+
+    // allocate internal arrays (current vs Guess)
+    psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
+    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
+    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32);
+    psVector *dy     = NULL;
+    psF64 Chisq = 0.0;
+    psF64 lambda = 0.001;
+
+    psVector *beta_lim = NULL;
+    psVector *param_min = NULL;
+    psVector *param_max = NULL;
+
+    // if we are provided a covar image, we expect to find these three vectors in first three rows
+    if (covar != NULL) {
+        beta_lim  = psVectorAlloc (params->n, PS_TYPE_F32);
+        param_min = psVectorAlloc (params->n, PS_TYPE_F32);
+        param_max = psVectorAlloc (params->n, PS_TYPE_F32);
+        for (int i = 0; i < params->n; i++) {
+            beta_lim->data.F32[i] = covar->data.F64[0][i];
+            param_min->data.F32[i] = covar->data.F64[1][i];
+            param_max->data.F32[i] = covar->data.F64[2][i];
+        }
+        psImageRecycle (covar, params->n, params->n, PS_TYPE_F64);
+    }
+        
+
+    // why is this needed here??? the initial guess on params is provided by the user
+    Params = psVectorCopy (Params, params, PS_TYPE_F32);
+
+    // the user provides the error or NULL.  we need to convert
+    // to appropriate weights
+    dy = psVectorAlloc (y->n, PS_TYPE_F32);
+    if (yErr != NULL) {
+        for (int i = 0; i < dy->n; i++) {
+            if (yErr->data.F32[i] == 0.0) {
+                dy->data.F32[i] = 1.0;
+                // mask this?  bad pixel, obviously...
+            } else {
+                dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
+            }
+        }
+    } else {
+        for (int i = 0; i < dy->n; i++) {
+            dy->data.F32[i] = 1.0;
+        }
+    }
+
+    // calculate initial alpha and beta, set chisq (min->value)
+    min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(min->value)) {
+        min->iter = min->maxIter;
+        return (false);
+    }
+    # ifndef PS_NO_TRACE
+    // dump some useful info if trace is defined
+    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+        p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
+        p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
+        p_psVectorPrint (psTraceGetDestination(), params, "params guess");
+    }
+    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
+        p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+    }
+    # endif /* PS_NO_TRACE */
+
+    // iterate until the tolerance is reached, or give up
+    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+
+        // set a new guess for Alpha, Beta, Params
+        p_psMinLM_GuessABP_EAM (Alpha, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, lambda);
+
+        // measure linear model prediction
+        psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda);
+
+        # ifndef PS_NO_TRACE
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
+            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
+            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
+        }
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
+            p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+        }
+        # endif /* PS_NO_TRACE */
+
+        // calculate Chisq for new guess, update Alpha & Beta
+        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func);
+
+        // XXX EAM alternate convergence criterion:
+        // compare the delta (min->value - Chisq) with the
+        // expected delta from the linear model (dLinear)
+        // accept new guess (if improvement), or increase lambda
+        psF64 rho = (min->value - Chisq) / dLinear;
+
+        psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
+        # ifndef PS_NO_TRACE
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
+            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
+            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
+        }
+        # endif /* PS_NO_TRACE */
+
+        /* if (Chisq < min->value) {  */
+        if (rho > 0.0) {  
+            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
+            min->value = Chisq;
+            alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
+            beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
+            params = psVectorCopy (params, Params, PS_TYPE_F32);
+            lambda *= 0.1;
+        } else {
+            lambda *= 10.0;
+        }
+        min->iter ++;
+    }
+    psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+
+    // construct & return the covariance matrix (if requested)
+    if (covar != NULL) {
+      p_psMinLM_GuessABP_EAM (covar, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, 0.0);
+    }      
+
+    // free the internal temporary data
+    psFree (alpha);
+    psFree (Alpha);
+    psFree (beta);
+    psFree (Beta);
+    psFree (Params);
+    psFree (dy);
+
+    if (min->iter == min->maxIter) {
+      return (false);
+    } 
+    return (true);
+}
+
+// XXX EAM : can we use static copies of LUv, LUm, A?
+psBool p_psMinLM_GuessABP_EAM (psImage  *Alpha,
+                               psVector *Beta,
+                               psVector *Params,
+                               const psImage  *alpha,
+                               const psVector *beta,
+                               const psVector *params,
+                               const psVector *paramMask,
+                               const psVector *beta_lim,
+                               const psVector *params_min,
+                               const psVector *params_max,
+                               psF64 lambda)
+{
+
+    # define USE_LU_DECOMP 1
+    # if (USE_LU_DECOMP)
+    psVector *LUv = NULL;
+    psImage  *LUm = NULL;
+    psImage  *A   = NULL;
+    psF32    det;
+
+    // LU decomposition version
+    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n");
+
+    // set new guess values (creates matrix A)
+    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    // solve A*beta = Beta (Alpha = 1/A)
+    // these operations do not modify the input values (creates LUm, LUv)
+    LUm   = psMatrixLUD (NULL, &LUv, A);
+    Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
+    Alpha = psMatrixInvert (Alpha, A, &det);
+
+    # else
+    // gauss-jordan version
+    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version");
+
+    // set new guess values (creates matrix A)
+    Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
+    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    psGaussJordan (Alpha, Beta);
+    # endif
+
+    // apply Beta to get new Params values
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+        // continue;
+        // compare Beta to beta limits
+        if (beta_lim != NULL) {
+            if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) {
+                Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]);
+            }
+        }
+        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+        // compare new params to param limits
+        if (params_max != NULL) {
+            Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]);
+        }
+        if (params_min != NULL) {
+            Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]);
+        }
+    }
+
+    # if (USE_LU_DECOMP)
+    psFree (A);
+    psFree (LUm);
+    psFree (LUv);
+    # endif
+
+    return true;
+}
+