forester/archive/RIO/others/puzzle_mod/src/util.c

   1 /*
   2  * util.c
   3  *
   4  *
   5  * Part of TREE-PUZZLE 5.0 (June 2000)
   6  *
   7  * (c) 1999-2000 by Heiko A. Schmidt, Korbinian Strimmer,
   8  *                  M. Vingron, and Arndt von Haeseler
   9  * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
  10  *
  11  * All parts of the source except where indicated are distributed under
  12  * the GNU public licence.  See http://www.opensource.org for details.
  13  */
  14
  15
  16 #include "util.h"
  17
  18 #define STDOUT     stdout
  19 #ifndef PARALLEL             /* because printf() runs significantly faster */
  20                              /* than fprintf(stdout) on an Apple McIntosh  */
  21                              /* (HS) */
  22 #       define FPRINTF    printf
  23 #       define STDOUTFILE
  24 #else
  25 #       define FPRINTF    fprintf
  26 #       define STDOUTFILE STDOUT,
  27         extern int PP_NumProcs;
  28         extern int PP_Myid;
  29         long int PP_randn;
  30         long int PP_rand;
  31 #endif
  32
  33
  34 /*
  35  * memory allocation error handler
  36  */
  37
  38 void maerror(char *message)
  39 {
  40         FPRINTF(STDOUTFILE "\n\n\nUnable to proceed (lack of memory: %s)\n\n", message);
  41         FPRINTF(STDOUTFILE "Hint for Macintosh users:\n");
  42         FPRINTF(STDOUTFILE "Use the <Get Info> command of the Finder to increase the memory partition!\n\n");
  43         exit(1);
  44 }
  45
  46
  47 /*
  48  * memory allocate double vectors, matrices, and cubes
  49  */
  50
  51 dvector new_dvector(int n)
  52 {
  53         dvector v;
  54
  55         v = (dvector) malloc((unsigned) (n * sizeof(double)));
  56         if (v == NULL) maerror("step 1 in new_dvector");
  57
  58         return v;
  59 }
  60
  61 dmatrix new_dmatrix(int nrow, int ncol)
  62 {
  63         int i;
  64         dmatrix m;
  65
  66         m = (dmatrix) malloc((unsigned) (nrow * sizeof(dvector)));
  67         if (m == NULL) maerror("step 1 in in new_dmatrix");
  68
  69         *m = (dvector) malloc((unsigned) (nrow * ncol * sizeof(double)));
  70         if (*m == NULL) maerror("step 2 in in new_dmatrix");
  71
  72         for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
  73
  74         return m;
  75 }
  76
  77 dcube new_dcube(int ntri, int nrow, int ncol)
  78 {
  79         int i, j;
  80         dcube c;
  81
  82         c = (dcube) malloc((unsigned) (ntri * sizeof(dmatrix)));
  83         if (c == NULL) maerror("step 1 in in new_dcube");
  84
  85         *c = (dmatrix) malloc((unsigned) (ntri * nrow * sizeof(dvector)));
  86         if (*c == NULL) maerror("step 2 in in new_dcube");
  87
  88         **c = (dvector) malloc((unsigned) (ntri * nrow * ncol * sizeof(double)));
  89         if (**c == NULL) maerror("step 3 in in new_dcube");
  90
  91         for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
  92
  93         for (i = 1; i < ntri; i++) {
  94                 c[i] = c[i-1] + nrow;
  95                 c[i][0] = c[i-1][0] + nrow * ncol;
  96                 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
  97         }
  98
  99         return c;
 100 }
 101
 102 void free_dvector(dvector v)
 103 {
 104         free((double *) v);
 105 }
 106
 107 void free_dmatrix(dmatrix m)
 108 {
 109         free((double *) *m);
 110         free((double *) m);
 111 }
 112
 113 void free_dcube(dcube c)
 114 {
 115         free((double *) **c);
 116         free((double *) *c);
 117         free((double *) c);
 118 }
 119
 120
 121 /*
 122  * memory allocate char vectors, matrices, and cubes
 123  */
 124
 125 cvector new_cvector(int n)
 126 {
 127         cvector v;
 128
 129         v = (cvector) malloc((unsigned)n * sizeof(char));
 130         if (v == NULL) maerror("step1 in new_cvector");
 131
 132         return v;
 133 }
 134
 135 cmatrix new_cmatrix(int nrow, int ncol)
 136 {
 137         int i;
 138         cmatrix m;
 139
 140         m = (cmatrix) malloc((unsigned) (nrow * sizeof(cvector)));
 141         if (m == NULL) maerror("step 1 in new_cmatrix");
 142
 143         *m = (cvector) malloc((unsigned) (nrow * ncol * sizeof(char)));
 144         if (*m == NULL) maerror("step 2 in new_cmatrix");
 145
 146         for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
 147
 148         return m;
 149 }
 150
 151 ccube new_ccube(int ntri, int nrow, int ncol)
 152 {
 153         int i, j;
 154         ccube c;
 155
 156         c = (ccube) malloc((unsigned) (ntri * sizeof(cmatrix)));
 157         if (c == NULL) maerror("step 1 in new_ccube");
 158
 159         *c = (cmatrix) malloc((unsigned) (ntri * nrow * sizeof(cvector)));
 160         if (*c == NULL) maerror("step 2 in new_ccube");
 161
 162         **c = (cvector) malloc((unsigned) (ntri * nrow * ncol * sizeof(char)));
 163         if (**c == NULL) maerror("step 3 in new_ccube");
 164
 165         for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
 166
 167         for (i = 1; i < ntri; i++) {
 168                 c[i] = c[i-1] + nrow;
 169                 c[i][0] = c[i-1][0] + nrow * ncol;
 170                 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
 171         }
 172
 173         return c;
 174 }
 175
 176 void free_cvector(cvector v)
 177 {
 178         free((char *) v);
 179 }
 180
 181 void free_cmatrix(cmatrix m)
 182 {
 183         free((char *) *m);
 184         free((char *) m);
 185 }
 186
 187 void free_ccube(ccube c)
 188 {
 189         free((char *) **c);
 190         free((char *) *c);
 191         free((char *) c);
 192 }
 193
 194
 195 /*
 196  * memory allocate int vectors, matrices, and cubes
 197  */
 198
 199 ivector new_ivector(int n)
 200 {
 201         ivector v;
 202
 203         v = (ivector) malloc((unsigned) (n * sizeof(int)));
 204         if (v == NULL) maerror("step 1 in new_ivector");
 205
 206         return v;
 207 }
 208
 209 imatrix new_imatrix(int nrow, int ncol)
 210 {
 211         int i;
 212         imatrix m;
 213
 214         m = (imatrix) malloc((unsigned) (nrow * sizeof(ivector)));
 215         if (m == NULL) maerror("step 1 in new_imatrix");
 216
 217         *m = (ivector) malloc((unsigned) (nrow * ncol * sizeof(int)));
 218         if (*m == NULL) maerror("step 2 in new_imatrix");
 219
 220         for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
 221
 222         return m;
 223 }
 224
 225 icube new_icube(int ntri, int nrow, int ncol)
 226 {
 227         int i, j;
 228         icube c;
 229
 230         c = (icube) malloc((unsigned) (ntri * sizeof(imatrix)));
 231         if (c == NULL) maerror("step 1 in new_icube");
 232
 233         *c = (imatrix) malloc((unsigned) (ntri * nrow * sizeof(ivector)));
 234         if (*c == NULL) maerror("step 2 in new_icube");
 235
 236         **c = (ivector) malloc((unsigned) (ntri * nrow * ncol * sizeof(int)));
 237         if (**c == NULL) maerror("step 3 in new_icube");
 238
 239         for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
 240
 241         for (i = 1; i < ntri; i++) {
 242                 c[i] = c[i-1] + nrow;
 243                 c[i][0] = c[i-1][0] + nrow * ncol;
 244                 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
 245         }
 246
 247         return c;
 248 }
 249
 250 void free_ivector(ivector v)
 251 {
 252         free((int *) v);
 253 }
 254
 255 void free_imatrix(imatrix m)
 256 {
 257         free((int *) *m);
 258         free((int *) m);
 259 }
 260
 261 void free_icube(icube c)
 262 {
 263         free((int *) **c);
 264         free((int *) *c);
 265         free((int *) c);
 266 }
 267
 268
 269 /*
 270  * memory allocate uli vectors, matrices, and cubes
 271  */
 272
 273 ulivector new_ulivector(int n)
 274 {
 275         ulivector v;
 276
 277         v = (ulivector) malloc((unsigned) (n * sizeof(uli)));
 278         if (v == NULL) maerror("step 1 in new_ulivector");
 279
 280         return v;
 281 }
 282
 283 ulimatrix new_ulimatrix(int nrow, int ncol)
 284 {
 285         int i;
 286         ulimatrix m;
 287
 288         m = (ulimatrix) malloc((unsigned) (nrow * sizeof(ulivector)));
 289         if (m == NULL) maerror("step 1 in new_ulimatrix");
 290
 291         *m = (ulivector) malloc((unsigned) (nrow * ncol * sizeof(uli)));
 292         if (*m == NULL) maerror("step 2 in new_ulimatrix");
 293
 294         for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
 295
 296         return m;
 297 }
 298
 299 ulicube new_ulicube(int ntri, int nrow, int ncol)
 300 {
 301         int i, j;
 302         ulicube c;
 303
 304         c = (ulicube) malloc((unsigned) (ntri * sizeof(ulimatrix)));
 305         if (c == NULL) maerror("step 1 in new_ulicube");
 306
 307         *c = (ulimatrix) malloc((unsigned) (ntri * nrow * sizeof(ulivector)));
 308         if (*c == NULL) maerror("step 2 in new_ulicube");
 309
 310         **c = (ulivector) malloc((unsigned) (ntri * nrow * ncol * sizeof(uli)));
 311         if (**c == NULL) maerror("step 3 in new_ulicube");
 312
 313         for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
 314
 315         for (i = 1; i < ntri; i++) {
 316                 c[i] = c[i-1] + nrow;
 317                 c[i][0] = c[i-1][0] + nrow * ncol;
 318                 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
 319         }
 320
 321         return c;
 322 }
 323
 324 void free_ulivector(ulivector v)
 325 {
 326         free((uli *) v);
 327 }
 328
 329 void free_ulimatrix(ulimatrix m)
 330 {
 331         free((uli *) *m);
 332         free((uli *) m);
 333 }
 334
 335 void free_ulicube(ulicube c)
 336 {
 337         free((uli *) **c);
 338         free((uli *) *c);
 339         free((uli *) c);
 340 }
 341
 342
 343 /******************************************************************************/
 344 /* random numbers generator  (Numerical recipes)                              */
 345 /******************************************************************************/
 346
 347 /* definitions */
 348 #define IM1 2147483563
 349 #define IM2 2147483399
 350 #define AM (1.0/IM1)
 351 #define IMM1 (IM1-1)
 352 #define IA1 40014
 353 #define IA2 40692
 354 #define IQ1 53668
 355 #define IQ2 52774
 356 #define IR1 12211
 357 #define IR2 3791
 358 #define NTAB 32
 359 #define NDIV (1+IMM1/NTAB)
 360 #define EPS 1.2e-7
 361 #define RNMX (1.0-EPS)
 362
 363 /* variable */
 364 long idum;
 365
 366 double randomunitintervall()
 367 /* Long period (> 2e18) random number generator. Returns a uniform random
 368    deviate between 0.0 and 1.0 (exclusive of endpoint values).
 369
 370    Source:
 371    Press et al., "Numerical recipes in C", Cambridge University Press, 1992
 372    (chapter 7 "Random numbers", ran2 random number generator) */
 373 {
 374         int j;
 375         long k;
 376         static long idum2=123456789;
 377         static long iy=0;
 378         static long iv[NTAB];
 379         double temp;
 380
 381         if (idum <= 0) {
 382                 if (-(idum) < 1)
 383                         idum=1;
 384                 else
 385                         idum=-(idum);
 386                 idum2=(idum);
 387                 for (j=NTAB+7;j>=0;j--) {
 388                         k=(idum)/IQ1;
 389                         idum=IA1*(idum-k*IQ1)-k*IR1;
 390                         if (idum < 0)
 391                                 idum += IM1;
 392                         if (j < NTAB)
 393                                 iv[j] = idum;
 394                 }
 395                 iy=iv[0];
 396         }
 397         k=(idum)/IQ1;
 398         idum=IA1*(idum-k*IQ1)-k*IR1;
 399         if (idum < 0)
 400                 idum += IM1;
 401         k=idum2/IQ2;
 402         idum2=IA2*(idum2-k*IQ2)-k*IR2;
 403         if (idum2 < 0)
 404                 idum2 += IM2;
 405         j=iy/NDIV;
 406         iy=iv[j]-idum2;
 407         iv[j] = idum;
 408         if (iy < 1)
 409                 iy += IMM1;
 410         if ((temp=AM*iy) > RNMX)
 411                 return RNMX;
 412         else
 413                 return temp;
 414 }
 415
 416 #undef IM1
 417 #undef IM2
 418 #undef AM
 419 #undef IMM1
 420 #undef IA1
 421 #undef IA2
 422 #undef IQ1
 423 #undef IQ2
 424 #undef IR1
 425 #undef IR2
 426 #undef NTAB
 427 #undef NDIV
 428 #undef EPS
 429 #undef RNMX
 430
 431 int initrandom(int seed)
 432 {
 433    srand((unsigned) time(NULL));
 434    if (seed < 0)
 435         seed = rand();
 436    idum=-(long) seed;
 437 #  ifdef PARALLEL
 438         {
 439         int n;
 440         for (n=0; n<PP_Myid; n++)
 441                 (void) randomunitintervall();
 442 #       ifdef PVERBOSE1
 443            fprintf(stderr, "(%2d) !!! random seed set to %d, %dx drawn !!!\n", PP_Myid, seed, n);
 444 #       endif
 445         }
 446 #  else
 447 #       ifdef PVERBOSE1
 448            fprintf(stderr, "!!! random seed set to %d !!!\n", seed);
 449 #       endif
 450 #  endif
 451    return (seed);
 452 }  /* initrandom */
 453
 454
 455 /* returns a random integer in the range [0; n - 1] */
 456 int randominteger(int n)
 457 {
 458         int t;
 459 #  ifndef FIXEDINTRAND
 460 #       ifndef PARALLEL
 461                 t = (int) floor(randomunitintervall()*n);
 462                 return t;
 463 #       else
 464                 int m;
 465                 for (m=1; m<PP_NumProcs; m++)
 466                         (void) randomunitintervall();
 467                 PP_randn+=(m-1); PP_rand++;
 468                 return (int) floor(randomunitintervall()*n);
 469 #       endif
 470 #  else
 471         fprintf(stderr, "!!! fixed \"random\" integers for testing purposes !!!\n");
 472         return (int)0;
 473 #  endif /* FIXEDINTRAND */
 474 }
 475
 476 /* Draw s numbers from the set 0,1,2,..,t-1 and put them
 477    into slist (every number can be drawn only one time) */
 478 void chooser(int t, int s, ivector slist)
 479 {
 480         int i, j, k, l;
 481         ivector isfree;
 482
 483         isfree = new_ivector(t);
 484         for (i = 0; i < t; i++) isfree[i] = TRUE;
 485         for (i = 0; i < s; i++) {
 486                 /* random integer in the range [0, t-1-i] */
 487                 j = randominteger(t-i);
 488                 k = -1;
 489                 l = -1;
 490                 do {
 491                         k++;
 492                         if (isfree[k] == TRUE) l++;
 493                 } while ( l != j);
 494                 slist[i] = k;
 495                 isfree[k] = FALSE;
 496         }
 497         free_ivector(isfree);
 498 }
 499
 500 /* a realloc function that works also on non-ANSI compilers */
 501 void *myrealloc(void *p, size_t s)
 502 {
 503         if (p == NULL) return malloc(s);
 504         else return realloc(p, s);
 505 }
 506
 507 /* safer variant of gets */
 508 /* Reads characters from stdin until a newline character or EOF
 509    is received.  The newline is not made part of the string.
 510    If an error occurs a null string \0 is returned */
 511 cvector mygets()
 512 {
 513         int c, n;
 514         cvector str;
 515
 516         str = new_cvector(100);
 517
 518         n = 0;
 519         c = getchar();
 520         while (c != '\n' && c != '\r' && n < 99 && c != EOF && !ferror(stdin))
 521         {
 522                 str[n] = (char) c;
 523                 c = getchar();
 524                 n++;
 525         }
 526         if (c == EOF || ferror(stdin))
 527         {
 528                 str[0] = '\0';
 529         }
 530         else
 531         {
 532                 str[n] = '\0';
 533         }
 534
 535         return str;
 536 }
 537
 538
 539
 540 /******************************************************************************/
 541 /* minimization of a function by Brents method (Numerical Recipes)            */
 542 /******************************************************************************/
 543
 544 double brent(double, double, double, double (*f )(double ), double, double *, double *, double, double, double);
 545
 546
 547 #define ITMAX 100
 548 #define CGOLD 0.3819660
 549 #define GOLD 1.618034
 550 #define GLIMIT 100.0
 551 #define TINY 1.0e-20
 552 #define ZEPS 1.0e-10
 553 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
 554 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
 555
 556 /* Brents method in one dimension */
 557 double brent(double ax, double bx, double cx, double (*f)(double), double tol,
 558         double *foptx, double *f2optx, double fax, double fbx, double fcx)
 559 {
 560         int iter;
 561         double a,b,d=0,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
 562         double xw,wv,vx;
 563         double e=0.0;
 564
 565         a=(ax < cx ? ax : cx);
 566         b=(ax > cx ? ax : cx);
 567         x=bx;
 568         fx=fbx;
 569         if (fax < fcx) {
 570                 w=ax;
 571                 fw=fax;
 572                 v=cx;
 573                 fv=fcx;
 574         } else {
 575                 w=cx;
 576                 fw=fcx;
 577                 v=ax;
 578                 fv=fax;
 579         }
 580         for (iter=1;iter<=ITMAX;iter++) {
 581                 xm=0.5*(a+b);
 582                 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
 583                 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
 584                         *foptx = fx;
 585                         xw = x-w;
 586                         wv = w-v;
 587                         vx = v-x;
 588                         *f2optx = 2.0*(fv*xw + fx*wv + fw*vx)/
 589                                 (v*v*xw + x*x*wv + w*w*vx);
 590                         return x;
 591                 }
 592                 if (fabs(e) > tol1) {
 593                         r=(x-w)*(fx-fv);
 594                         q=(x-v)*(fx-fw);
 595                         p=(x-v)*q-(x-w)*r;
 596                         q=2.0*(q-r);
 597                         if (q > 0.0) p = -p;
 598                         q=fabs(q);
 599                         etemp=e;
 600                         e=d;
 601                         if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
 602                                 d=CGOLD*(e=(x >= xm ? a-x : b-x));
 603                         else {
 604                                 d=p/q;
 605                                 u=x+d;
 606                                 if (u-a < tol2 || b-u < tol2)
 607                                         d=SIGN(tol1,xm-x);
 608                         }
 609                 } else {
 610                         d=CGOLD*(e=(x >= xm ? a-x : b-x));
 611                 }
 612                 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
 613                 fu=(*f)(u);
 614                 if (fu <= fx) {
 615                         if (u >= x) a=x; else b=x;
 616                         SHFT(v,w,x,u)
 617                         SHFT(fv,fw,fx,fu)
 618                 } else {
 619                         if (u < x) a=u; else b=u;
 620                         if (fu <= fw || w == x) {
 621                                 v=w;
 622                                 w=u;
 623                                 fv=fw;
 624                                 fw=fu;
 625                         } else if (fu <= fv || v == x || v == w) {
 626                                 v=u;
 627                                 fv=fu;
 628                         }
 629                 }
 630         }
 631         *foptx = fx;
 632         xw = x-w;
 633         wv = w-v;
 634         vx = v-x;
 635         *f2optx = 2.0*(fv*xw + fx*wv + fw*vx)/
 636                 (v*v*xw + x*x*wv + w*w*vx);
 637         return x;
 638 }
 639 #undef ITMAX
 640 #undef CGOLD
 641 #undef ZEPS
 642 #undef SHFT
 643 #undef SIGN
 644 #undef GOLD
 645 #undef GLIMIT
 646 #undef TINY
 647
 648 /* one-dimensional minimization - as input a lower and an upper limit and a trial
 649   value for the minimum is needed: xmin < xguess < xmax
 650   the function and a fractional tolerance has to be specified
 651   onedimenmin returns the optimal x value and the value of the function
 652   and its second derivative at this point
 653   */
 654 double onedimenmin(double xmin, double xguess, double xmax, double (*f)(double),
 655         double tol, double *fx, double *f2x)
 656 {
 657         double eps, optx, ax, bx, cx, fa, fb, fc;
 658
 659         /* first attempt to bracketize minimum */
 660         eps = xguess*tol*50.0;
 661         ax = xguess - eps;
 662         if (ax < xmin) ax = xmin;
 663         bx = xguess;
 664         cx = xguess + eps;
 665         if (cx > xmax) cx = xmax;
 666
 667         /* check if this works */
 668         fa = (*f)(ax);
 669         fb = (*f)(bx);
 670         fc = (*f)(cx);
 671
 672         /* if it works use these borders else be conservative */
 673         if ((fa < fb) || (fc < fb)) {
 674                 if (ax != xmin) fa = (*f)(xmin);
 675                 if (cx != xmax) fc = (*f)(xmax);
 676                 optx = brent(xmin, xguess, xmax, f, tol, fx, f2x, fa, fb, fc);
 677         } else
 678                 optx = brent(ax, bx, cx, f, tol, fx, f2x, fa, fb, fc);
 679
 680         return optx; /* return optimal x */
 681 }
 682
 683 /* two-dimensional minimization with borders and calculations of standard errors */
 684 /* we optimize along basis vectors - not very optimal but it seems to work well */
 685 void twodimenmin(double tol,
 686         int active1, double min1, double *x1, double max1, double (*func1)(double), double *err1,
 687         int active2, double min2, double *x2, double max2, double (*func2)(double), double *err2)
 688 {
 689         int it, nump, change;
 690         double x1old, x2old;
 691         double fx, f2x;
 692
 693         it = 0;
 694         nump = 0;
 695
 696         /* count number of parameters */
 697         if (active1) nump++;
 698         if (active2) nump++;
 699
 700         do { /* repeat until nothing changes any more */
 701                 it++;
 702                 change = FALSE;
 703
 704                 /* optimize first variable */
 705                 if (active1) {
 706
 707                         if ((*x1) <= min1) (*x1) = min1 + 0.2*(max1-min1);
 708                         if ((*x1) >= max1) (*x1) = max1 - 0.2*(max1-min1);
 709                         x1old = (*x1);
 710                         (*x1) = onedimenmin(min1, (*x1), max1, func1, tol, &fx, &f2x);
 711                         if ((*x1) < min1) (*x1) = min1;
 712                         if ((*x1) > max1) (*x1) = max1;
 713                         /* same tolerance as 1D minimization */
 714                         if (fabs((*x1) - x1old) > 3.3*tol) change = TRUE;
 715
 716                         /* standard error */
 717                         f2x = fabs(f2x);
 718                         if (1.0/(max1*max1) < f2x) (*err1) = sqrt(1.0/f2x);
 719                         else (*err1) = max1;
 720
 721                 }
 722
 723                 /* optimize second variable */
 724                 if (active2) {
 725
 726                         if ((*x2) <= min2) (*x2) = min2 + 0.2*(max2-min2);
 727                         if ((*x2) >= max2) (*x2) = max2 - 0.2*(max2-min2);
 728                         x2old = (*x2);
 729                         (*x2) = onedimenmin(min2, (*x2), max2, func2, tol, &fx, &f2x);
 730                         if ((*x2) < min2) (*x2) = min2;
 731                         if ((*x2) > max2) (*x2) = max2;
 732                         /* same tolerance as 1D minimization */
 733                         if (fabs((*x2) - x2old) > 3.3*tol) change = TRUE;
 734
 735                         /* standard error */
 736                         f2x = fabs(f2x);
 737                         if (1.0/(max2*max2) < f2x) (*err2) = sqrt(1.0/f2x);
 738                         else (*err2) = max2;
 739
 740                 }
 741
 742                 if (nump == 1) return;
 743
 744         } while (it != MAXITS && change);
 745
 746         return;
 747 }
 748