Changeset 6525 for trunk/psLib/src/math/psMinimizeLMM.c

Timestamp:

Mar 6, 2006, 10:32:36 AM (20 years ago)

Author:

gusciora

Message:

....

File:

• trunk/psLib/src/math/psMinimizeLMM.c (modified) (9 diffs)

trunk/psLib/src/math/psMinimizeLMM.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-03-06 20:32:19 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-03-06 20:32:36 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -239 +239 @@
     psF64 ymodel;
     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
+    deriv->n = deriv->nalloc;
 
     // zero alpha and beta for summing below
@@ -301 +302 @@
 NOTES: EAM: this is my re-implementation of MinLM
  
-XXX: EAM : constraints added 2006.02.15
+XXX: Must implement code to handle the constrain->min, ->max, ->delta members.
  
 XXX: Put the ASSERTS in.
@@ -328 +329 @@
     PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
     psVector *paramMask = NULL;
-    psVector *paramDelta = NULL;
-    psVector *paramMin = NULL;
-    psVector *paramMax = NULL;
     if (constrain != NULL) {
-        // XXX EAM : fill out the asserts
-        paramDelta = constrain->paramDelta;
-        paramMin = constrain->paramMin;
-        paramMax = constrain->paramMax;
         paramMask = constrain->paramMask;
         if (paramMask != NULL) {
             PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_U8, false);
             PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false);
+        }
+        if (constrain->paramMin != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: the constrain->paramMin vector is currently ignored.\n");
+        }
+        if (constrain->paramMax != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: the constrain->paramMax vector is currently ignored.\n");
+        }
+        if (constrain->paramDelta != NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: the constrain->paramDelta vector is currently ignored.\n");
         }
     }
@@ -368 +371 @@
     psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F64);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
+    beta->n = beta->nalloc;
+    Beta->n = Beta->nalloc;
+    Params->n = Params->nalloc;
     psVector *dy     = NULL;
     psF64 Chisq = 0.0;
     psF64 lambda = 0.001;
+    // XXX: Code this properly.  Don't use mustFree00.
+    psBool mustFree00 = false;
+    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) {
+        mustFree00 = true;
+        beta_lim  = psVectorAlloc(params->n, PS_TYPE_F32);
+        param_min = psVectorAlloc(params->n, PS_TYPE_F32);
+        param_max = psVectorAlloc(params->n, PS_TYPE_F32);
+        beta_lim->n = beta_lim->nalloc;
+        param_min->n = param_min->nalloc;
+        param_max->n = param_max->nalloc;
+        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
@@ -381 +409 @@
     } else {
         dy = psVectorAlloc(y->n, PS_TYPE_F32);
+        dy->n = dy->nalloc;
         psVectorInit(dy, 1.0);
     }
@@ -410 +439 @@
         // set a new guess for Alpha, Beta, Params
         p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask,
-                           paramDelta, paramMin, paramMax, lambda);
+                           beta_lim, param_min, param_max, lambda);
 
         // measure linear model prediction
@@ -467 +496 @@
     if (covar != NULL) {
         p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask,
-                           paramDelta, paramMin, paramMax, 0.0);
+                           beta_lim, param_min, param_max, 0.0);
     }
 
@@ -478 +507 @@
     if (yWt == NULL) {
         psFree(dy);
+    }
+    if (mustFree00 == true) {
+        psFree(beta_lim);
+        psFree(param_min);
+        psFree(param_max);
     }
     if (min->iter == min->maxIter) {