src/tools/goal-seek.c

   1 /*
   2  * goal-seek.c:  A generic root finder.
   3  *
   4  * Author:
   5  *   Morten Welinder (terra@gnome.org)
   6  *
   7  */
   8
   9 #undef DEBUG_GOAL_SEEK
  10 #ifdef STANDALONE
  11 #define DEBUG_GOAL_SEEK
  12 #endif
  13
  14 #include <gnumeric-config.h>
  15 #include <numbers.h>
  16 #include <gnumeric.h>
  17 #include <tools/goal-seek.h>
  18 #include <gnm-random.h>
  19 #include <value.h>
  20 #include <cell.h>
  21 #include <sheet.h>
  22
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <limits.h>
  26
  27
  28 static gboolean
  29 update_data (gnm_float x, gnm_float y, GnmGoalSeekData *data)
  30 {
  31         if (!gnm_finite (y))
  32                 return FALSE;
  33
  34         if (y > 0) {
  35                 if (data->havexpos) {
  36                         if (data->havexneg) {
  37                                 /*
  38                                  * When we have pos and neg, prefer the new point only
  39                                  * if it makes the pos-neg x-interval smaller.
  40                                  */
  41                                 if (gnm_abs (x - data->xneg) < gnm_abs (data->xpos - data->xneg)) {
  42                                         data->xpos = x;
  43                                         data->ypos = y;
  44                                 }
  45                         } else if (y < data->ypos) {
  46                                 /* We have pos only and our neg y is closer to zero.  */
  47                                 data->xpos = x;
  48                                 data->ypos = y;
  49                         }
  50                 } else {
  51                         data->xpos = x;
  52                         data->ypos = y;
  53                         data->havexpos = TRUE;
  54                 }
  55                 return FALSE;
  56         } else if (y < 0) {
  57                 if (data->havexneg) {
  58                         if (data->havexpos) {
  59                                 /*
  60                                  * When we have pos and neg, prefer the new point only
  61                                  * if it makes the pos-neg x-interval smaller.
  62                                  */
  63                                 if (gnm_abs (x - data->xpos) < gnm_abs (data->xpos - data->xneg)) {
  64                                         data->xneg = x;
  65                                         data->yneg = y;
  66                                 }
  67                         } else if (-y < -data->yneg) {
  68                                 /* We have neg only and our neg y is closer to zero.  */
  69                                 data->xneg = x;
  70                                 data->yneg = y;
  71                         }
  72
  73                 } else {
  74                         data->xneg = x;
  75                         data->yneg = y;
  76                         data->havexneg = TRUE;
  77                 }
  78                 return FALSE;
  79         } else {
  80                 /* Lucky guess...  */
  81                 data->have_root = TRUE;
  82                 data->root = x;
  83 #ifdef DEBUG_GOAL_SEEK
  84                 g_print ("update_data: got root %.20" GNM_FORMAT_g "\n", x);
  85 #endif
  86                 return TRUE;
  87         }
  88 }
  89
  90
  91 /*
  92  * Calculate a reasonable approximation to the derivative of a function
  93  * in a single point.
  94  */
  95 static GnmGoalSeekStatus
  96 fake_df (GnmGoalSeekFunction f, gnm_float x, gnm_float *dfx, gnm_float xstep,
  97          GnmGoalSeekData *data, void *user_data)
  98 {
  99         gnm_float xl, xr, yl, yr;
 100         GnmGoalSeekStatus status;
 101
 102 #ifdef DEBUG_GOAL_SEEK
 103         g_print ("fake_df (x=%.20" GNM_FORMAT_g ", xstep=%.20" GNM_FORMAT_g ")\n",
 104                 x, xstep);
 105 #endif
 106
 107         xl = x - xstep;
 108         if (xl < data->xmin)
 109                 xl = x;
 110
 111         xr = x + xstep;
 112         if (xr > data->xmax)
 113                 xr = x;
 114
 115         if (xl == xr) {
 116 #ifdef DEBUG_GOAL_SEEK
 117                 g_print ("==> xl == xr\n");
 118 #endif
 119                 return GOAL_SEEK_ERROR;
 120         }
 121
 122         status = f (xl, &yl, user_data);
 123         if (status != GOAL_SEEK_OK) {
 124 #ifdef DEBUG_GOAL_SEEK
 125                 g_print ("==> failure at xl\n");
 126 #endif
 127                 return status;
 128         }
 129 #ifdef DEBUG_GOAL_SEEK
 130         g_print ("==> xl=%.20" GNM_FORMAT_g "; yl=%.20" GNM_FORMAT_g "\n",
 131                  xl, yl);
 132 #endif
 133
 134         status = f (xr, &yr, user_data);
 135         if (status != GOAL_SEEK_OK) {
 136 #ifdef DEBUG_GOAL_SEEK
 137                 g_print ("==> failure at xr\n");
 138 #endif
 139                 return status;
 140         }
 141 #ifdef DEBUG_GOAL_SEEK
 142         g_print ("==> xr=%.20" GNM_FORMAT_g "; yr=%.20" GNM_FORMAT_g "\n",
 143                  xr, yr);
 144 #endif
 145
 146         *dfx = (yr - yl) / (xr - xl);
 147 #ifdef DEBUG_GOAL_SEEK
 148         g_print ("==> %.20" GNM_FORMAT_g "\n", *dfx);
 149 #endif
 150         return gnm_finite (*dfx) ? GOAL_SEEK_OK : GOAL_SEEK_ERROR;
 151 }
 152
 153 void
 154 goal_seek_initialize (GnmGoalSeekData *data)
 155 {
 156         data->havexpos = data->havexneg = data->have_root = FALSE;
 157         data->xpos = data->xneg = data->root = gnm_nan;
 158         data->ypos = data->yneg = gnm_nan;
 159         data->xmin = -1e10;
 160         data->xmax = +1e10;
 161         data->precision = 1e-10;
 162 }
 163
 164
 165 /**
 166  * goal_seek_point:
 167  * @f: (scope call): object function
 168  * @data: #GnmGoalSeekData state
 169  * @user_data: user data for @f
 170  * @x0: root guess
 171  *
 172  * Seek a goal using a single point.
 173  *
 174  * Returns:
 175  */
 176 GnmGoalSeekStatus
 177 goal_seek_point (GnmGoalSeekFunction f, GnmGoalSeekData *data,
 178                  void *user_data, gnm_float x0)
 179 {
 180         GnmGoalSeekStatus status;
 181         gnm_float y0;
 182
 183         if (data->have_root)
 184                 return GOAL_SEEK_OK;
 185
 186 #ifdef DEBUG_GOAL_SEEK
 187         g_print ("goal_seek_point\n");
 188 #endif
 189
 190         if (x0 < data->xmin || x0 > data->xmax)
 191                 return GOAL_SEEK_ERROR;
 192
 193         status = f (x0, &y0, user_data);
 194         if (status != GOAL_SEEK_OK)
 195                 return status;
 196
 197         if (update_data (x0, y0, data))
 198                 return GOAL_SEEK_OK;
 199
 200         return GOAL_SEEK_ERROR;
 201 }
 202
 203
 204 static GnmGoalSeekStatus
 205 goal_seek_newton_polish (GnmGoalSeekFunction f, GnmGoalSeekFunction df,
 206                          GnmGoalSeekData *data, void *user_data,
 207                          gnm_float x0, gnm_float y0)
 208 {
 209         int iterations;
 210         gnm_float last_df0 = 1;
 211         gboolean try_newton = TRUE;
 212         gboolean try_square = x0 != 0 && gnm_abs (x0) < 1e10;
 213
 214 #ifdef DEBUG_GOAL_SEEK
 215         g_print ("goal_seek_newton_polish\n");
 216 #endif
 217
 218         for (iterations = 0; iterations < 20; iterations++) {
 219                 if (try_square) {
 220                         gnm_float x1 = x0 * gnm_abs (x0);
 221                         gnm_float y1, r;
 222                         GnmGoalSeekStatus status = f (x1, &y1, user_data);
 223                         if (status != GOAL_SEEK_OK)
 224                                 goto nomore_square;
 225
 226                         if (update_data (x1, y1, data))
 227                                 return GOAL_SEEK_OK;
 228
 229                         r = gnm_abs (y1 / y0);
 230                         if (r >= 1)
 231                                 goto nomore_square;
 232
 233                         x0 = x1;
 234 #ifdef DEBUG_GOAL_SEEK
 235                         g_print ("polish square: x0=%.20" GNM_FORMAT_g "\n",
 236                                  x0);
 237 #endif
 238                         if (r > 0.5)
 239                                 goto nomore_square;
 240
 241                         continue;
 242
 243                 nomore_square:
 244                         try_square = FALSE;
 245                 }
 246
 247                 if (try_newton) {
 248                         gnm_float df0, r, x1, y1;
 249                         GnmGoalSeekStatus status = df
 250                                 ? df (x0, &df0, user_data)
 251                                 : fake_df (f, x0, &df0, gnm_abs (x0) / 1e6, data, user_data);
 252                         if (status != GOAL_SEEK_OK || df0 == 0)
 253                                 df0 = last_df0;  /* Bogus */
 254                         else
 255                                 last_df0 = df0;
 256
 257                         x1 = x0 - y0 / df0;
 258                         if (x1 < data->xmin || x1 > data->xmax)
 259                                 goto nomore_newton;
 260
 261                         status = f (x1, &y1, user_data);
 262                         if (status != GOAL_SEEK_OK)
 263                                 goto nomore_newton;
 264
 265                         if (update_data (x1, y1, data))
 266                                 return GOAL_SEEK_OK;
 267
 268                         r = gnm_abs (y1 / y0);
 269                         if (r >= 1)
 270                                 goto nomore_newton;
 271
 272                         x0 = x1;
 273 #ifdef DEBUG_GOAL_SEEK
 274                         g_print ("polish Newton: x0=%.20" GNM_FORMAT_g "\n",
 275                                  x0);
 276 #endif
 277                         if (r > 0.5)
 278                                 goto nomore_newton;
 279
 280                         continue;
 281
 282                 nomore_newton:
 283                         try_newton = FALSE;
 284                 }
 285
 286                 /* Nothing left to try.  */
 287                 break;
 288         }
 289
 290         if (goal_seek_bisection (f, data, user_data) == GOAL_SEEK_OK)
 291                 return GOAL_SEEK_OK;
 292
 293         data->root = x0;
 294         data->have_root = TRUE;
 295         return GOAL_SEEK_OK;
 296 }
 297
 298
 299 /**
 300  * goal_seek_newton:
 301  * @f: (scope call): object function
 302  * @df: (scope call) (nullable): object function derivative
 303  * @data: #GnmGoalSeekData state
 304  * @user_data: user data for @f and @df
 305  * @x0: root guess
 306  *
 307  * Seek a goal (root) using Newton's iterative method.
 308  *
 309  * The supplied function must (should) be continously differentiable in
 310  * the supplied interval.  If @df is %NULL, this function will
 311  * estimate the derivative.
 312  *
 313  * This method will find a root rapidly provided the initial guess, x0,
 314  * is sufficiently close to the root.  (The number of significant digits
 315  * (asympotically) goes like i^2 unless the root is a multiple root in
 316  * which case it is only like c*i.)
 317  */
 318 GnmGoalSeekStatus
 319 goal_seek_newton (GnmGoalSeekFunction f, GnmGoalSeekFunction df,
 320                   GnmGoalSeekData *data, void *user_data, gnm_float x0)
 321 {
 322         int iterations;
 323         gnm_float precision = data->precision / 2;
 324         gnm_float last_df0 = 1;
 325         gnm_float step_factor = 1e-6;
 326
 327         if (data->have_root)
 328                 return GOAL_SEEK_OK;
 329
 330 #ifdef DEBUG_GOAL_SEEK
 331         g_print ("goal_seek_newton\n");
 332 #endif
 333
 334         for (iterations = 0; iterations < 100; iterations++) {
 335                 gnm_float x1, y0, df0, stepsize;
 336                 GnmGoalSeekStatus status;
 337                 gboolean flat;
 338
 339 #ifdef DEBUG_GOAL_SEEK
 340                 g_print ("x0 = %.20" GNM_FORMAT_g "   (i=%d)\n", x0, iterations);
 341 #endif
 342
 343                 /* Check whether we have left the valid interval.  */
 344                 if (x0 < data->xmin || x0 > data->xmax)
 345                         return GOAL_SEEK_ERROR;
 346
 347                 status = f (x0, &y0, user_data);
 348                 if (status != GOAL_SEEK_OK)
 349                         return status;
 350
 351 #ifdef DEBUG_GOAL_SEEK
 352                 g_print ("                                        y0 = %.20" GNM_FORMAT_g "\n", y0);
 353 #endif
 354
 355                 if (update_data (x0, y0, data))
 356                         return GOAL_SEEK_OK;
 357
 358                 if (df)
 359                         status = df (x0, &df0, user_data);
 360                 else {
 361                         gnm_float xstep;
 362
 363                         if (gnm_abs (x0) < 1e-10)
 364                                 if (data->havexneg && data->havexpos)
 365                                         xstep = gnm_abs (data->xpos - data->xneg) / 1e6;
 366                                 else
 367                                         xstep = (data->xmax - data->xmin) / 1e6;
 368                         else
 369                                 xstep = step_factor * gnm_abs (x0);
 370
 371                         status = fake_df (f, x0, &df0, xstep, data, user_data);
 372                 }
 373                 if (status != GOAL_SEEK_OK)
 374                         return status;
 375
 376                 /* If we hit a flat spot, we are in trouble.  */
 377                 flat = (df0 == 0);
 378                 if (flat) {
 379                         last_df0 /= 2;
 380                         if (gnm_abs (last_df0) <= GNM_MIN)
 381                                 return GOAL_SEEK_ERROR;
 382                         df0 = last_df0;  /* Might be utterly bogus.  */
 383                 } else
 384                         last_df0 = df0;
 385
 386                 if (data->havexpos && data->havexneg)
 387                         x1 = x0 - y0 / df0;
 388                 else
 389                         /*
 390                          * Overshoot slightly to prevent us from staying on
 391                          * just one side of the root.
 392                          */
 393                         x1 = x0 - 1.000001 * y0 / df0;
 394
 395                 stepsize = gnm_abs (x1 - x0) / (gnm_abs (x0) + gnm_abs (x1));
 396
 397 #ifdef DEBUG_GOAL_SEEK
 398                 g_print ("                                        df0 = %.20" GNM_FORMAT_g "\n", df0);
 399                 g_print ("                                        ss = %.20" GNM_FORMAT_g "\n", stepsize);
 400 #endif
 401
 402                 if (stepsize < precision) {
 403                         goal_seek_newton_polish (f, df, data, user_data, x0, y0);
 404                         return GOAL_SEEK_OK;
 405                 }
 406
 407                 if (flat && iterations > 0) {
 408                         /*
 409                          * Verify that we made progress using our
 410                          * potentially bogus df0.
 411                          */
 412                         gnm_float y1;
 413
 414                         if (x1 < data->xmin || x1 > data->xmax)
 415                                 return GOAL_SEEK_ERROR;
 416
 417                         status = f (x1, &y1, user_data);
 418                         if (status != GOAL_SEEK_OK)
 419                                 return status;
 420
 421 #ifdef DEBUG_GOAL_SEEK
 422                         g_print ("                                        y1 = %.20" GNM_FORMAT_g "\n", y1);
 423 #endif
 424                         if (gnm_abs (y1) >= 0.9 * gnm_abs (y0))
 425                                 return GOAL_SEEK_ERROR;
 426                 }
 427
 428                 if (stepsize < step_factor)
 429                         step_factor = stepsize;
 430
 431                 x0 = x1;
 432         }
 433
 434         return GOAL_SEEK_ERROR;
 435 }
 436
 437 /**
 438  * goal_seek_bisection:
 439  * @f: (scope call): object function
 440  * @data: #GnmGoalSeekData state.
 441  * @user_data: user data for @f.
 442  *
 443  * Seek a goal (root) using bisection methods.
 444  *
 445  * The supplied function must (should) be continous over the interval.
 446  *
 447  * Caller must have located a positive and a negative point.
 448  *
 449  * This method will find a root steadily using bisection to narrow the
 450  * interval in which a root lies.
 451  *
 452  * It alternates between mid-point-bisection (semi-slow, but guaranteed
 453  * progress), secant-bisection (usually quite fast, but sometimes gets
 454  * nowhere), and Ridder's Method (usually fast, harder to fool than
 455  * the secant method).
 456  */
 457 GnmGoalSeekStatus
 458 goal_seek_bisection (GnmGoalSeekFunction f, GnmGoalSeekData *data,
 459                      void *user_data)
 460 {
 461         int iterations;
 462         gnm_float stepsize;
 463         int newton_submethod = 0;
 464
 465         if (data->have_root)
 466                 return GOAL_SEEK_OK;
 467
 468 #ifdef DEBUG_GOAL_SEEK
 469         g_print ("goal_seek_bisection\n");
 470 #endif
 471
 472         if (!data->havexpos || !data->havexneg)
 473                 return GOAL_SEEK_ERROR;
 474
 475         stepsize = gnm_abs (data->xpos - data->xneg)
 476                 / (gnm_abs (data->xpos) + gnm_abs (data->xneg));
 477
 478         /* log_2 (10) = 3.3219 < 4.  */
 479         for (iterations = 0; iterations < 100 + GNM_DIG * 4; iterations++) {
 480                 gnm_float xmid, ymid;
 481                 GnmGoalSeekStatus status;
 482                 enum { M_SECANT, M_RIDDER, M_NEWTON, M_MIDPOINT } method;
 483
 484                 method = (iterations % 4 == 0)
 485                         ? M_RIDDER
 486                         : ((iterations % 4 == 2)
 487                            ? M_NEWTON
 488                            : M_MIDPOINT);
 489
 490         again:
 491                 switch (method) {
 492                 default:
 493                         abort ();
 494
 495                 case M_SECANT:
 496                         xmid = data->xpos - data->ypos *
 497                                 ((data->xneg - data->xpos) /
 498                                  (data->yneg - data->ypos));
 499                         break;
 500
 501                 case M_RIDDER: {
 502                         gnm_float det;
 503
 504                         xmid = (data->xpos + data->xneg) / 2;
 505                         status = f (xmid, &ymid, user_data);
 506                         if (status != GOAL_SEEK_OK)
 507                                 continue;
 508                         if (ymid == 0) {
 509                                 update_data (xmid, ymid, data);
 510                                 return GOAL_SEEK_OK;
 511                         }
 512
 513                         det = gnm_sqrt (ymid * ymid - data->ypos * data->yneg);
 514                         if (det == 0)
 515                                  /* This might happen with underflow, I guess. */
 516                                 continue;
 517
 518                         xmid += (xmid - data->xpos) * ymid / det;
 519                         break;
 520                 }
 521
 522                 case M_MIDPOINT:
 523                         xmid = (data->xpos + data->xneg) / 2;
 524                         break;
 525
 526                 case M_NEWTON: {
 527                         gnm_float x0, y0, xstep, df0;
 528
 529                         /* This method is only effective close-in.  */
 530                         if (stepsize > 0.1) {
 531                                 method = M_MIDPOINT;
 532                                 goto again;
 533                         }
 534
 535                         switch (newton_submethod++ % 4) {
 536                         case 0: x0 = data->xpos; y0 = data->ypos; break;
 537                         case 2: x0 = data->xneg; y0 = data->yneg; break;
 538                         default:
 539                         case 3:
 540                         case 1:
 541                                 x0 = (data->xpos + data->xneg) / 2;
 542
 543                                 status = f (x0, &y0, user_data);
 544                                 if (status != GOAL_SEEK_OK)
 545                                         continue;
 546                         }
 547
 548                         xstep = gnm_abs (data->xpos - data->xneg) / 1e6;
 549                         status = fake_df (f, x0, &df0, xstep, data, user_data);
 550                         if (status != GOAL_SEEK_OK)
 551                                 continue;
 552
 553                         if (df0 == 0)
 554                                 continue;
 555
 556                         /*
 557                          * Overshoot by 1% to prevent us from staying on
 558                          * just one side of the root.
 559                          */
 560                         xmid = x0 - 1.01 * y0 / df0;
 561                 }
 562                 }
 563
 564                 if ((xmid < data->xpos && xmid < data->xneg) ||
 565                     (xmid > data->xpos && xmid > data->xneg)) {
 566                         /* We left the interval.  */
 567                         xmid = (data->xpos + data->xneg) / 2;
 568                         method = M_MIDPOINT;
 569                 }
 570
 571                 status = f (xmid, &ymid, user_data);
 572                 if (status != GOAL_SEEK_OK)
 573                         continue;
 574
 575 #ifdef DEBUG_GOAL_SEEK
 576                 {
 577                         const char *themethod;
 578                         switch (method) {
 579                         case M_MIDPOINT: themethod = "midpoint"; break;
 580                         case M_RIDDER: themethod = "Ridder"; break;
 581                         case M_SECANT: themethod = "secant"; break;
 582                         case M_NEWTON: themethod = "Newton"; break;
 583                         default: themethod = "?";
 584                         }
 585
 586                         g_print ("xmid = %.20" GNM_FORMAT_g " (%s)\n", xmid, themethod);
 587                         g_print ("                                        ymid = %.20" GNM_FORMAT_g "\n", ymid);
 588                 }
 589 #endif
 590
 591                 if (update_data (xmid, ymid, data)) {
 592                         return GOAL_SEEK_OK;
 593                 }
 594
 595                 stepsize = gnm_abs (data->xpos - data->xneg)
 596                         / (gnm_abs (data->xpos) + gnm_abs (data->xneg));
 597
 598 #ifdef DEBUG_GOAL_SEEK
 599                 g_print ("                                          ss = %.20" GNM_FORMAT_g "\n", stepsize);
 600 #endif
 601
 602                 if (stepsize < GNM_EPSILON) {
 603                         if (data->yneg < ymid)
 604                                 ymid = data->yneg, xmid = data->xneg;
 605
 606                         if (data->ypos < ymid)
 607                                 ymid = data->ypos, xmid = data->xpos;
 608
 609                         data->have_root = TRUE;
 610                         data->root = xmid;
 611                         return GOAL_SEEK_OK;
 612                 }
 613         }
 614         return GOAL_SEEK_ERROR;
 615 }
 616
 617 #undef SECANT_P
 618 #undef RIDDER_P
 619
 620 /**
 621  * goal_seek_trawl_uniformly:
 622  * @f: (scope call): object function
 623  * @data: #GnmGoalSeekData state
 624  * @user_data: user data for @f
 625  * @xmin: lower search bound
 626  * @xmax: upper search bound
 627  * @points: number of points to try.
 628  */
 629 GnmGoalSeekStatus
 630 goal_seek_trawl_uniformly (GnmGoalSeekFunction f,
 631                            GnmGoalSeekData *data, void *user_data,
 632                            gnm_float xmin, gnm_float xmax,
 633                            int points)
 634 {
 635         int i;
 636
 637         if (data->have_root)
 638                 return GOAL_SEEK_OK;
 639
 640 #ifdef DEBUG_GOAL_SEEK
 641         g_print ("goal_seek_trawl_uniformly\n");
 642 #endif
 643
 644         if (xmin > xmax || xmin < data->xmin || xmax > data->xmax)
 645                 return GOAL_SEEK_ERROR;
 646
 647         for (i = 0; i < points; i++) {
 648                 gnm_float x, y;
 649                 GnmGoalSeekStatus status;
 650
 651                 if (data->havexpos && data->havexneg)
 652                         break;
 653
 654                 x = xmin + (xmax - xmin) * random_01 ();
 655                 status = f (x, &y, user_data);
 656                 if (status != GOAL_SEEK_OK)
 657                         /* We are not depending on the result, so go on.  */
 658                         continue;
 659
 660 #ifdef DEBUG_GOAL_SEEK
 661                 g_print ("x = %.20" GNM_FORMAT_g "\n", x);
 662                 g_print ("                                        y = %.20" GNM_FORMAT_g "\n", y);
 663 #endif
 664
 665                 if (update_data (x, y, data))
 666                         return GOAL_SEEK_OK;
 667         }
 668
 669         /* We were not (extremely) lucky, so we did not actually hit the
 670            root.  We report this as an error.  */
 671         return GOAL_SEEK_ERROR;
 672 }
 673
 674 /**
 675  * goal_seek_trawl_normally:
 676  * @f: (scope call): object function
 677  * @data: #GnmGoalSeekData state
 678  * @user_data: user data for @f
 679  * @mu: search mean
 680  * @sigma: search standard deviation
 681  * @points: number of points to try.
 682  */
 683 GnmGoalSeekStatus
 684 goal_seek_trawl_normally (GnmGoalSeekFunction f,
 685                           GnmGoalSeekData *data, void *user_data,
 686                           gnm_float mu, gnm_float sigma,
 687                           int points)
 688 {
 689         int i;
 690
 691         if (data->have_root)
 692                 return GOAL_SEEK_OK;
 693
 694 #ifdef DEBUG_GOAL_SEEK
 695         g_print ("goal_seek_trawl_normally\n");
 696 #endif
 697
 698         if (sigma <= 0 || mu < data->xmin || mu > data->xmax)
 699                 return GOAL_SEEK_ERROR;
 700
 701         for (i = 0; i < points; i++) {
 702                 gnm_float x, y;
 703                 GnmGoalSeekStatus status;
 704
 705                 if (data->havexpos && data->havexneg)
 706                         break;
 707
 708                 x = mu + sigma * random_normal ();
 709                 if (x < data->xmin || x > data->xmax)
 710                         continue;
 711
 712                 status = f (x, &y, user_data);
 713                 if (status != GOAL_SEEK_OK)
 714                         /* We are not depending on the result, so go on.  */
 715                         continue;
 716
 717 #ifdef DEBUG_GOAL_SEEK
 718                 g_print ("x = %.20" GNM_FORMAT_g "\n", x);
 719                 g_print ("                                        y = %.20" GNM_FORMAT_g "\n", y);
 720 #endif
 721
 722                 if (update_data (x, y, data))
 723                         return GOAL_SEEK_OK;
 724         }
 725
 726         /* We were not (extremely) lucky, so we did not actually hit the
 727            root.  We report this as an error.  */
 728         return GOAL_SEEK_ERROR;
 729 }
 730
 731 /**
 732  * gnm_goal_seek_eval_cell:
 733  * @x: x-value for which to evaluate
 734  * @y: (out): location to store result
 735  * @data: user data
 736  *
 737  * Returns: An status indicating whether evaluation went ok.
 738  */
 739 GnmGoalSeekStatus
 740 gnm_goal_seek_eval_cell (gnm_float x, gnm_float *y, gpointer data_)
 741 {
 742         GnmGoalSeekCellData const *data = data_;
 743         GnmValue *v = value_new_float (x);
 744
 745         gnm_cell_set_value (data->xcell, v);
 746         cell_queue_recalc (data->xcell);
 747         gnm_cell_eval (data->ycell);
 748
 749         if (data->ycell->value &&
 750             VALUE_IS_NUMBER (data->ycell->value)) {
 751                 *y = value_get_as_float (data->ycell->value) - data->ytarget;
 752                 if (gnm_finite (*y))
 753                         return GOAL_SEEK_OK;
 754         }
 755
 756         return GOAL_SEEK_ERROR;
 757 }
 758
 759 GnmGoalSeekStatus
 760 gnm_goal_seek_cell (GnmGoalSeekData *data,
 761                     GnmGoalSeekCellData *celldata)
 762 {
 763         GnmGoalSeekStatus status;
 764         gboolean hadold;
 765         gnm_float oldx;
 766         GnmValue *v;
 767
 768         hadold = !VALUE_IS_EMPTY_OR_ERROR (celldata->xcell->value);
 769         oldx = hadold ? value_get_as_float (celldata->xcell->value) : 0;
 770
 771         /* PLAN A: Newton's iterative method from initial or midpoint.  */
 772         {
 773                 gnm_float x0;
 774
 775                 if (hadold && oldx >= data->xmin && oldx <= data->xmax)
 776                         x0 = oldx;
 777                 else
 778                         x0 = (data->xmin + data->xmax) / 2;
 779
 780                 status = goal_seek_newton (gnm_goal_seek_eval_cell, NULL,
 781                                            data, celldata,
 782                                            x0);
 783                 if (status == GOAL_SEEK_OK)
 784                         goto DONE;
 785         }
 786
 787         /* PLAN B: Trawl uniformly.  */
 788         if (!data->havexpos || !data->havexneg) {
 789                 status = goal_seek_trawl_uniformly (gnm_goal_seek_eval_cell,
 790                                                     data, celldata,
 791                                                     data->xmin, data->xmax,
 792                                                     100);
 793                 if (status == GOAL_SEEK_OK)
 794                         goto DONE;
 795         }
 796
 797         /* PLAN C: Trawl normally from middle.  */
 798         if (!data->havexpos || !data->havexneg) {
 799                 gnm_float sigma, mu;
 800                 int i;
 801
 802                 sigma = MIN (data->xmax - data->xmin, 1e6);
 803                 mu = (data->xmax + data->xmin) / 2;
 804
 805                 for (i = 0; i < 5; i++) {
 806                         sigma /= 10;
 807                         status = goal_seek_trawl_normally (gnm_goal_seek_eval_cell,
 808                                                            data, celldata,
 809                                                            mu, sigma, 30);
 810                         if (status == GOAL_SEEK_OK)
 811                                 goto DONE;
 812                 }
 813         }
 814
 815         /* PLAN D: Trawl normally from left.  */
 816         if (!data->havexpos || !data->havexneg) {
 817                 gnm_float sigma, mu;
 818                 int i;
 819
 820                 sigma = MIN (data->xmax - data->xmin, 1e6);
 821                 mu = data->xmin;
 822
 823                 for (i = 0; i < 5; i++) {
 824                         sigma /= 10;
 825                         status = goal_seek_trawl_normally (gnm_goal_seek_eval_cell,
 826                                                            data, celldata,
 827                                                            mu, sigma, 20);
 828                         if (status == GOAL_SEEK_OK)
 829                                 goto DONE;
 830                 }
 831         }
 832
 833         /* PLAN E: Trawl normally from right.  */
 834         if (!data->havexpos || !data->havexneg) {
 835                 gnm_float sigma, mu;
 836                 int i;
 837
 838                 sigma = MIN (data->xmax - data->xmin, 1e6);
 839                 mu = data->xmax;
 840
 841                 for (i = 0; i < 5; i++) {
 842                         sigma /= 10;
 843                         status = goal_seek_trawl_normally (gnm_goal_seek_eval_cell,
 844                                                            data, celldata,
 845                                                            mu, sigma, 20);
 846                         if (status == GOAL_SEEK_OK)
 847                                 goto DONE;
 848                 }
 849         }
 850
 851         /* PLAN F: Newton iteration with uniform net of starting points.  */
 852         if (!data->havexpos || !data->havexneg) {
 853                 int i;
 854                 const int N = 10;
 855
 856                 for (i = 1; i <= N; i++) {
 857                         gnm_float x0 =  data->xmin +
 858                                 (data->xmax - data->xmin) / (N + 1) * i;
 859
 860                         status = goal_seek_newton (gnm_goal_seek_eval_cell, NULL,
 861                                                    data, celldata,
 862                                                    x0);
 863                         if (status == GOAL_SEEK_OK)
 864                                 goto DONE;
 865                 }
 866         }
 867
 868         /* PLAN Z: Bisection.  */
 869         {
 870                 status = goal_seek_bisection (gnm_goal_seek_eval_cell,
 871                                               data, celldata);
 872                 if (status == GOAL_SEEK_OK)
 873                         goto DONE;
 874         }
 875
 876  DONE:
 877         if (status == GOAL_SEEK_OK)
 878                 v = value_new_float (data->root);
 879         else if (hadold)
 880                 v = value_new_float (oldx);
 881         else
 882                 v = value_new_empty ();
 883         sheet_cell_set_value (celldata->xcell, v);
 884
 885         return status;
 886 }