[Teem-commit] teem/src/ten grads.c,1.1,1.2 ten.h,1.77,1.78 tendGrads.c,1.3,1.4

[<kin...@us...>]

Update of /cvsroot/teem/teem/src/ten
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7475

Modified Files:
	grads.c ten.h tendGrads.c 
Log Message:
added second phase of tenGradientDistribute in which randomized signs of gradients are tried to minimize the length of the mean gradient vector

Index: grads.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/grads.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** grads.c	27 Feb 2004 10:28:07 -0000	1.1
--- grads.c	27 Feb 2004 12:50:36 -0000	1.2
***************
*** 31,36 ****
      ret->drag = 0.01;
      ret->dt = 0.05;
-     ret->minVelocity = 0.00001;
      ret->jitter = 0.05;
      ret->srand = AIR_TRUE;
      ret->snap = 0;
--- 31,38 ----
      ret->drag = 0.01;
      ret->dt = 0.05;
      ret->jitter = 0.05;
+     ret->minVelocity = 0.000001;
+     ret->minMean = 0.0001;
+     ret->minMeanImprovement = 0.00001;
      ret->srand = AIR_TRUE;
      ret->snap = 0;
***************
*** 241,250 ****
  }
  
  /*
  ******** tenGradientDistribute
  **
! ** takes the given list of gradients, normalizes their lengths, optionally jitters
! ** their positions, does point repulsion, and then selects a combination of
! ** directions with minimum vector sum.
  */
  int
--- 243,300 ----
  }
  
+ double
+ _tenGradientParty(double *grad, int num) {
+   double mean[3];
+   int ii;
+   
+   ELL_3V_SET(mean, 0, 0, 0);
+   for (ii=0; ii<num; ii++) {
+     if (airRandInt(2)) {
+       ELL_3V_SCALE(grad + 3*ii, -1, grad + 3*ii);
+     }
+     ELL_3V_INCR(mean, grad + 3*ii);
+   }
+   ELL_3V_SCALE(mean, 1.0/num, mean);
+   return ELL_3V_LEN(mean);
+ }
+ 
+ int
+ tenGradientMeanMinimize(Nrrd *nout, Nrrd *nin, tenGradientParm *tgparm) {
+   char me[]="tenGradientMeanMinimize", err[AIR_STRLEN_MED];
+   int num;
+   double *grad, len, lastLen, improv;
+ 
+   if (!nout || tenGradientCheck(nin, nrrdTypeUnknown, 2)) {
+     sprintf(err, "%s: got NULL pointer or invalid input", me);
+     biffAdd(TEN, err); return 1;
+   }
+   if (nrrdConvert(nout, nin, nrrdTypeDouble)) {
+     sprintf(err, "%s: can't initialize output with input", me);
+     biffMove(TEN, err, NRRD); return 1;
+   }
+   num = nout->axis[1].size;
+   grad = (double*)(nout->data);
+    
+   lastLen = _tenGradientParty(grad, num);
+   do {
+     do {
+       len = _tenGradientParty(grad, num);
+     } while (len > lastLen);
+     improv = lastLen - len;
+     lastLen = len;
+     fprintf(stderr, "%s: improvement: %g  (mean length = %g)\n",
+ 	    me, improv, len);
+   } while (improv > tgparm->minMeanImprovement
+ 	   && len > tgparm->minMean);
+   
+   return 0;
+ }
+ 
  /*
  ******** tenGradientDistribute
  **
! ** takes the given list of gradients, normalizes their lengths,
! ** optionally jitters their positions, does point repulsion, and then
! ** selects a combination of directions with minimum vector sum.
  */
  int
***************
*** 344,359 ****
        }
      } else {
!       if (!(iter % 25)) {
! 	fprintf(stderr, "%s: %d: meanVelocity = %g\n", me, iter, meanVelocity);
        }
      }
    }
  
!   /* select output gradients
!      - either do the trivial thing
!      - or run through 2^(#grads) to find one with maximum covariance
!      print covariance matrix eigenvalues
!   */
!   nrrdCopy(nout, npos);
  
    airMopOkay(mop); 
--- 394,409 ----
        }
      } else {
!       if (!(iter % 1000)) {
! 	fprintf(stderr, "%s: iteration = %d: meanVelocity = %g\n",
! 		me, iter, meanVelocity);
        }
      }
    }
  
!   fprintf(stderr, "%s: optimizing balance ... \n", me);
!   if (tenGradientMeanMinimize(nout, npos, tgparm)) {
!     sprintf(err, "%s: failed to minimize vector sum of gradients", me);
!     biffAdd(TEN, err); return 1;
!   }
  
    airMopOkay(mop); 
***************
*** 372,376 ****
    }
    if (!( num >= 3 )) {
!     sprintf(err, "%s: can generate at minimum 3 gradient directions (not %d)", me, num);
      biffAdd(TEN, err); return 1;
    }
--- 422,427 ----
    }
    if (!( num >= 3 )) {
!     sprintf(err, "%s: can generate minimum of 3 gradient directions "
! 	    "(not %d)", me, num);
      biffAdd(TEN, err); return 1;
    }
***************
*** 388,392 ****
    return 0;
  }
- 
- 
- 
--- 439,440 ----

Index: ten.h
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/ten.h,v
retrieving revision 1.77
retrieving revision 1.78
diff -C2 -d -r1.77 -r1.78
*** ten.h	27 Feb 2004 10:23:56 -0000	1.77
--- ten.h	27 Feb 2004 12:50:37 -0000	1.78
***************
*** 378,383 ****
      drag, 
      dt,
      minVelocity,
!     jitter; 
    int srand, snap, single, minIteration, maxIteration;
  } tenGradientParm;
--- 378,385 ----
      drag, 
      dt,
+     jitter,
      minVelocity,
!     minMean,
!     minMeanImprovement;
    int srand, snap, single, minIteration, maxIteration;
  } tenGradientParm;
***************
*** 400,405 ****
  TEEM_API int tenGradientRandom(Nrrd *ngrad, int num, int srand);
  TEEM_API int tenGradientJitter(Nrrd *nout, Nrrd *nin, double dist);
! TEEM_API int tenGradientDistribute(Nrrd *nout, Nrrd *nin, tenGradientParm *tgparm);
! TEEM_API int tenGradientGenerate(Nrrd *nout, int num, tenGradientParm *tgparm);
  
  /* enumsTen.c */
--- 402,411 ----
  TEEM_API int tenGradientRandom(Nrrd *ngrad, int num, int srand);
  TEEM_API int tenGradientJitter(Nrrd *nout, Nrrd *nin, double dist);
! TEEM_API int tenGradientMeanMinimize(Nrrd *nout, Nrrd *nin,
! 				     tenGradientParm *tgparm);
! TEEM_API int tenGradientDistribute(Nrrd *nout, Nrrd *nin,
! 				   tenGradientParm *tgparm);
! TEEM_API int tenGradientGenerate(Nrrd *nout, int num,
! 				 tenGradientParm *tgparm);
  
  /* enumsTen.c */

Index: tendGrads.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendGrads.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** tendGrads.c	27 Feb 2004 10:41:03 -0000	1.3
--- tendGrads.c	27 Feb 2004 12:50:37 -0000	1.4
***************
*** 21,33 ****
  #include "privateTen.h"
  
! #define INFO "Calculate gradient directions for DWI acquisition"
  char *_tend_gradsInfoL =
    (INFO
     ", based on a simulation of anti-podal point pairs repelling each other "
!    "on the unit sphere surface. This can either distribute more uniformly a given "
!    "set of gradients, or it can make a new distribution from scratch.  A more "
!    "clever implementation could decrease drag with time, as the solution converges, "
!    "to get closer to the minimum energy configuration faster.  In the mean time, you "
!    "can run a second pass on the output of the first pass, using lower drag. ");
  
  int
--- 21,36 ----
  #include "privateTen.h"
  
! #define INFO "Calculate balanced gradient directions for DWI acquisition"
  char *_tend_gradsInfoL =
    (INFO
     ", based on a simulation of anti-podal point pairs repelling each other "
!    "on the unit sphere surface. This can either distribute more uniformly "
!    "a given set of gradients, or it can make a new distribution from scratch. "
!    "A more clever implementation could decrease drag with time, as the "
!    "solution converges, to get closer to the minimum energy configuration "
!    "faster.  In the mean time, you can run a second pass on the output of "
!    "the first pass, using lower drag. A second phase of the algorithm "
!    "tries random signs in gradient directions in trying to find an optimally "
!    "balanced set of directions.");
  
  int
***************
*** 57,63 ****
    hestOptAdd(&hopt, "dt", "dt", airTypeDouble, 1, 1, &(tgparm->dt), "0.05",
  	     "time increment in solver");
!   hestOptAdd(&hopt, "drag", "drag", airTypeDouble, 1, 1, &(tgparm->drag), "0.01",
! 	     "viscous drag, to keep things stable");
!   hestOptAdd(&hopt, "charge", "charge", airTypeDouble, 1, 1, &(tgparm->charge), "0.2",
  	     "amount of charge on each particle");
    hestOptAdd(&hopt, "single", NULL, airTypeInt, 0, 0, &(tgparm->single), NULL,
--- 60,67 ----
    hestOptAdd(&hopt, "dt", "dt", airTypeDouble, 1, 1, &(tgparm->dt), "0.05",
  	     "time increment in solver");
!   hestOptAdd(&hopt, "drag", "drag", airTypeDouble, 1, 1, &(tgparm->drag),
! 	     "0.0005", "viscous drag, to keep things stable");
!   hestOptAdd(&hopt, "charge", "charge", airTypeDouble, 1, 1,
! 	     &(tgparm->charge), "0.2",
  	     "amount of charge on each particle");
    hestOptAdd(&hopt, "single", NULL, airTypeInt, 0, 0, &(tgparm->single), NULL,
***************
*** 69,77 ****
  	     "positions should be saved out.  By default (not using this "
  	     "option), there is no such snapshot behavior");
!   hestOptAdd(&hopt, "jitter", "jitter", airTypeDouble, 1, 1, &(tgparm->jitter), "0.05",
  	     "amount by which to perturb points when given an input nrrd");
    hestOptAdd(&hopt, "maxiter", "# iters", airTypeInt, 1, 1,
  	     &(tgparm->maxIteration), "1000000",
  	     "max number of iterations for which to run the simulation");
    hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, "-",
  	     "file to write output nrrd to");
--- 73,97 ----
  	     "positions should be saved out.  By default (not using this "
  	     "option), there is no such snapshot behavior");
!   hestOptAdd(&hopt, "jitter", "jitter", airTypeDouble, 1, 1,
! 	     &(tgparm->jitter), "0.05",
  	     "amount by which to perturb points when given an input nrrd");
+   hestOptAdd(&hopt, "minvelo", "vel", airTypeDouble, 1, 1, 
+ 	     &(tgparm->minVelocity), "0.00001",
+ 	     "low threshold on mean velocity of repelling points, "
+ 	     "at which point repulsion phase of algorithm terminates. ");
    hestOptAdd(&hopt, "maxiter", "# iters", airTypeInt, 1, 1,
  	     &(tgparm->maxIteration), "1000000",
  	     "max number of iterations for which to run the simulation");
+   hestOptAdd(&hopt, "minimprov", "delta", airTypeDouble, 1, 1, 
+ 	     &(tgparm->minMeanImprovement), "0.00001",
+ 	     "in the second phase of the algorithm, "
+ 	     "when stochastically optimizing the balance of gradients, "
+ 	     "the (small) improvement in length of mean gradient "
+ 	     "which triggers termination (as further improvements "
+ 	     "are unlikely. ");
+   hestOptAdd(&hopt, "minmean", "len", airTypeDouble, 1, 1,
+ 	     &(tgparm->minMean), "0.0001",
+ 	     "if length of mean gradient falls below this, finish "
+ 	     "the balancing phase");
    hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, "-",
  	     "file to write output nrrd to");