Fixed 'y > yub' error by reducing default tolerances min_d{x,y} in

tension_balance() and x_of_xo().  Now that this is known to be a
problem, I need to come up with a more permanent fix...

git-svn-id: svn://abax.physics.drexel.edu/sawsim/trunk@11 865a22a6-13cc-4084-8aa6-876098d8aa20

diff --git a/TODO b/TODO
new file mode 100644 (file)
index 0000000..9eaf623
--- /dev/null
+++ b/TODO
@@ -0,0 +1 @@
+fix arbitrary length scales and solution tolerances somhow

diff --git a/sawsim.nw b/sawsim.nw
index 1c01ab79f700624e27ca1e742b383013fef87cb7..6ec8e2ba80fe71e246cadd2cc4d56d0430497e3c 100644 (file)
--- a/sawsim.nw
+++ b/sawsim.nw
@@ -1867,8 +1867,8 @@ double tension_balance(int num_tension_groups,
                       double x)
 {
   static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
-  double F, x_guess, xo, lb, ub;
-  double min_dx=1e-10, min_dy=1e-15;
+  double F, xo_guess, xo, lb, ub;
+  double min_dx=1e-15, min_dy=1e-16;
   int max_steps=100, i;
 
   assert(num_tension_groups > 0);
@@ -1879,9 +1879,9 @@ double tension_balance(int num_tension_groups,
   if (num_tension_groups == 1) { /* only one group, no balancing required */
     xo = x;
   } else {
-    //printf("balancing for x = %g with ", x);
-    //for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, x_xo_data.xi[i]);
-    //printf("\n");
+    //fprintf(stderr, "balancing for x = %g with ", x);
+    //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
+    //fprintf(stderr, "\n");
     if (last_x == -1) { /* new group setup, reset x_xo_data */
       /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
       if (x_xo_data.xi != NULL)
@@ -1895,15 +1895,15 @@ double tension_balance(int num_tension_groups,
         x_xo_data.xi[i] = -1.0;
     }
     if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
-      if (x == 0) x_guess = 1.0;
-      else        x_guess = x/num_tension_groups;
+      if (x == 0) xo_guess = length_scale;
+      else        xo_guess = x/num_tension_groups;
 #ifdef DEBUG
-      fprintf(stderr, "search bracket for %g with guess of %g\n", x, x_guess);
+      fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
 #endif
-      oneD_bracket(x_of_xo, &x_xo_data, x, x_guess, &lb, &ub);
+      oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
     } else { /* work off the last balanced point */
 #ifdef DEBUG
-      fprintf(stderr, "step off the old bracketing x %g + %g = %g\n", last_x, x-last_x, x);
+      fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
 #endif
       if (x > last_x) {
         lb = x_xo_data.xi[0];
@@ -1912,18 +1912,23 @@ double tension_balance(int num_tension_groups,
         lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
         ub = x_xo_data.xi[0];
       } else { /* x == last_x */
-        printf("not moving\n");
+        fprintf(stderr, "not moving\n");
         lb= x_xo_data.xi[0];
         ub= x_xo_data.xi[0];
       }
     }
-    //printf("lb %g,\tub %g\n", lb, ub);
+#ifdef DEBUG
+    fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
+           __FUNCTION__, x, lb, ub);
+#endif
     xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_dx, min_dy, max_steps, NULL);
-    if (num_tension_groups > 1 && 0) {
-      printf("balanced  for x = %g with ", x);
-      for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, x_xo_data.xi[i]);
-      printf("\n");
+#ifdef DEBUG
+    if (num_tension_groups > 1) {
+      fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
+      for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
+      fprintf(stderr, "\n");
     }
+#endif
   }
   F = (*tension_handlers[0])(xo, params[0]);
   return F;
@@ -1949,7 +1954,7 @@ double x_of_xo(double xo, void *pdata)
   x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
   double F, x=0, xi, xi_guess, lb, ub;
   int i;
-  double min_dx=1e-10, min_dy=1e-14;
+  double min_dx=1e-15, min_dy=1e-16;
   int max_steps=100;
   assert(data->n_groups > 0);
   data->xi[0] = xo;
@@ -1957,7 +1962,7 @@ double x_of_xo(double xo, void *pdata)
   x += xo;
   if (data->xi)  data->xi[0] = xo;
   for (i=1; i < data->n_groups; i++) {
-    if (data->xi[i] == -1) xi_guess = 1.0;
+    if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
     else                   xi_guess = data->xi[i];
     oneD_bracket(data->tension_handler[i], data->handler_data[i], F,  xi_guess, &lb, &ub);
     xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_dx, min_dy, max_steps, NULL);
@@ -2001,17 +2006,25 @@ double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
   if (ylb == y) { x=lb; yx=ylb; goto end; }
   if (yub == y) { x=ub; yx=yub; goto end; }
 
-  //printf("lb %g, x %g, ub %g\tylb %g, y %g, yub %g\n", lb, x, ub, ylb, y, yub);
+#ifdef DEBUG
+  fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
+#endif
   assert(ylb < y); assert(yub > y);
   for (i=0; i<max_steps; i++) {
     x = (lb+ub)/2.0;
     yx = (*f)(x, params);
+#ifdef DEBUG
+  fprintf(stderr, "%s:\tbracket: lb %g, x %g, ub %g\tylb %g, yx %g, yub %g\t(y %g)\n", __FUNCTION__, lb, x, ub, ylb, yx, yub, y);
+#endif
     if (ub-lb < min_dx || yub-ylb < min_dy || yx == y) break;
     else if (yx > y)  { ub=x; yub=yx; }
     else /* yx < y */ { lb=x; ylb=yx; }
   }
  end:
   if (pYx) *pYx = yx;
+#ifdef DEBUG
+  fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
+#endif
   return x;
 }
 @ 
@@ -2028,20 +2041,30 @@ void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double
   double yx, step, x=xguess;
   int i=0;
   yx = (*f)(x, params);
-  //fprintf(stdout, "bracketing %g, start at f(%g) = %g\n", y, x, yx);
+#ifdef DEBUG
+  fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
+#endif
   if (yx > y) {
     assert(x > 0.0);
     *ub = x;
     *lb = 0;
+#ifdef DEBUG
+    fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
+#endif
   } else {
     *lb = x;
-    if (x == 0) x = 0.5; /* guess a scale of 1.0 */
+    if (x == 0) x = length_scale/2.0; /* HACK */
     while (yx < y) {
       x *= 2.0;
       yx = (*f)(x, params);
-      //fprintf(stdout, "increasing to f(%g) = %g\n", x, yx);
+#ifdef DEBUG
+      fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
+#endif
     }
     *ub = x;
+#ifdef DEBUG
+    fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
+#endif
   }
   //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
 }
@@ -2050,9 +2073,13 @@ void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double
 \subsection{Balancing implementation}
 
 <<tension-balance.c>>=
+//#define DEBUG
 <<license comment>>
 <<tension balance includes>>
 #include "tension_balance.h"
+
+double length_scale = 1e-10; /* HACK */
+
 <<tension balance internal definitions>>
 <<tension balance functions>>
 @