Changeset 7764 for trunk/psLib/src/math/psMinimizePowell.c

Timestamp:

Jun 29, 2006, 3:56:22 PM (20 years ago)

Author:

Paul Price

Message:

Removing static variables. See bug 768.

File:

• trunk/psLib/src/math/psMinimizePowell.c (modified) (19 diffs)

trunk/psLib/src/math/psMinimizePowell.c

@@ -11 +11 @@
  *  NOTE: XXX: The SDR is silent about data types.  F32 is implemented here.
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-24 23:43:15 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-06-30 01:56:22 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -64 +64 @@
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
-/*****************************************************************************/
-
-/*****************************************************************************/
-/* GLOBAL VARIABLES                                                          */
-/*****************************************************************************/
-static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
-static psVector *PowellValue;
-static psVector *PowellError;
-/*****************************************************************************/
-/* FILE STATIC VARIABLES                                                     */
 /*****************************************************************************/
 
@@ -239 +229 @@
 psTrace(__func__, 4, "Final bracket (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); \
 psTrace(__func__, 4, "---- p_psDetermineBracket() end ----\n"); \
+psFree(tmp); \
 return(bracket); \
 
@@ -256 +247 @@
     psF32 fc = 0.0;
     psS32 iter = 0;
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
+    psVector *tmp = psVectorAlloc(params->n, PS_TYPE_F32);
     psBool boolLineIsNull = true;
     psF32 prevMin = 0.0;
@@ -268 +259 @@
         psTrace(__func__, 2, "p_psDetermineBracket() called with zero line vector.\n");
         psTrace(__func__, 4, "---- p_psDetermineBracket() end (NULL) ----\n");
+        psFree(tmp);
         return(NULL);
     }
@@ -374 +366 @@
     psF32 fn = 0.0;
     psF32 mul = 0.0;
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpc, params->n, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(tmpn, params->n, PS_TYPE_F32);
     psS32 i = 0;
     psS32 boolLineIsNull = true;
@@ -405 +393 @@
     }
     numIterations = 0;
+
+    psVector *tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
+
     while (numIterations < PS_LINEMIN_MAX_ITERATIONS) {
         numIterations++;
@@ -462 +456 @@
             psFree(bracket);
             psTrace(__func__, 4, "---- LineMin() end.a (%f) (%f) ----\n", mul, min->value);
+            psFree(tmpa);
+            psFree(tmpb);
+            psFree(tmpc);
+            psFree(tmpn);
             return(mul);
         }
@@ -472 +470 @@
 
     psFree(bracket);
+    psFree(tmpa);
+    psFree(tmpb);
+    psFree(tmpc);
+    psFree(tmpn);
     return(mul);
 }
@@ -503 +505 @@
     PS_ASSERT_PTR_NON_NULL(func, NULL);
     psS32 numDims = params->n;
-    PS_VECTOR_GEN_STATIC_RECYCLED(pQP, numDims, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(u, numDims, PS_TYPE_F32);
-    PS_VECTOR_GEN_STATIC_RECYCLED(Q, numDims, PS_TYPE_F32);
     psS32 i = 0;
     psS32 j = 0;
@@ -533 +532 @@
     }
     PS_ASSERT_VECTORS_SIZE_EQUAL(params, myParamMask, NULL);
+
+
+    psVector *pQP = psVectorAlloc(numDims, PS_TYPE_F32);
+    psVector *u   = psVectorAlloc(numDims, PS_TYPE_F32);
+    psVector *Q   = psVectorAlloc(numDims, PS_TYPE_F32);
 
     // 1: Set v[i] to be the unit vectors for each dimension in params
@@ -578 +582 @@
                             "Could not perform line minimization.  Returning FALSE.\n");
                     psFree(v);
+                    psFree(pQP);
+                    psFree(u);
+                    psFree(Q);
                     return(false);
                 }
@@ -620 +627 @@
                     "Could not perform line minimization.  Returning FALSE.\n");
             psFree(v);
+            psFree(pQP);
+            psFree(u);
+            psFree(Q);
             return(false);
         }
@@ -626 +636 @@
         if (dummyMin.value > currFuncVal) {
             psFree(v);
+            psFree(pQP);
+            psFree(u);
+            psFree(Q);
             min->iter = iterationNumber;
             min->value = currFuncVal;
@@ -664 +677 @@
         if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
             psFree(v);
+            psFree(pQP);
+            psFree(u);
+            psFree(Q);
             // XXX: Ensure that currFuncVal is the correct value to use here.
             min->value = currFuncVal;
@@ -674 +690 @@
 
     psFree(v);
+    psFree(pQP);
+    psFree(u);
+    psFree(Q);
     min->iter = iterationNumber;
     psTrace(__func__, 4, "---- psMinimizePowell() end (0) (false) ----\n");
@@ -695 +714 @@
     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
     PS_ASSERT_VECTOR_NON_EMPTY(params, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(PowellValue, NAN);
-    PS_ASSERT_VECTOR_NON_EMPTY(PowellValue, NAN);
     PS_ASSERT_PTR_NON_NULL(coords, NAN);
 
@@ -704 +721 @@
     psVector *tmp;
 
-    tmp = Chi2PowellFunc(params, coords);
-    if (PowellError == NULL) {
+    psVector *values = coords->data[coords->n];
+    psVector *errors = coords->data[coords->n + 1];
+    psMinimizeChi2PowellFunc *func = coords->data[coords->n + 2];
+
+    PS_ASSERT_VECTOR_NON_NULL(values, NAN);
+    PS_ASSERT_VECTOR_NON_EMPTY(values, NAN);
+    PS_ASSERT_VECTOR_TYPE(values, PS_TYPE_F32, NAN);
+    if (errors) {
+        PS_ASSERT_VECTOR_NON_NULL(errors, NAN);
+        PS_ASSERT_VECTOR_NON_EMPTY(errors, NAN);
+        PS_ASSERT_VECTOR_TYPE(errors, PS_TYPE_F32, NAN);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NAN);
+    }
+
+    tmp = (*func)(params, coords);
+
+    if (errors == NULL) {
         for (i=0;i<coords->n;i++) {
-            d = (tmp->data.F32[i] - PowellValue->data.F32[i]);
+            d = (tmp->data.F32[i] - values->data.F32[i]);
             chi2+= d * d;
         }
     } else {
         for (i=0;i<coords->n;i++) {
-            d = (tmp->data.F32[i] - PowellValue->data.F32[i]) / PowellError->data.F32[i];
+            d = (tmp->data.F32[i] - values->data.F32[i]) / errors->data.F32[i];
             chi2+= d * d;
         }
@@ -741 +773 @@
     psMinimizeChi2PowellFunc model)
 {
-    PowellValue = (psVector *) value;
-    PowellError = (psVector *) error;
-
-    Chi2PowellFunc = model;
-
-    return(psMinimizePowell(min, params, constrain->paramMask, coords, myPowellChi2Func));
+    // Generate extended version of coords array, so we can pass in extra data to the chi^2 function
+    psArray *newCoords = psArrayAlloc(coords->n + 3);
+    for (long i = 0; i < coords->n; i++) {
+        newCoords->data[i] = psMemIncrRefCounter(coords->data[i]);
+    }
+    newCoords->n = coords->n;           // We deceive everyone else as to the length
+    // Casting away const: I'm not going to hurt you, just want to increment your reference counter is all
+    newCoords->data[coords->n] = psMemIncrRefCounter((psVector*)value);
+    newCoords->data[coords->n + 1] = psMemIncrRefCounter((psVector*)error);
+    newCoords->data[coords->n + 2] = &model;
+
+    bool success = psMinimizePowell(min, params, constrain->paramMask, newCoords, myPowellChi2Func);
+
+    newCoords->data[coords->n - 1] = NULL; // We can't free the array with a function pointer on it
+    psFree(newCoords);
+    return success;
 }