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