1 /*
2 * util.c
3 *
4 *
5 * Part of TREE-PUZZLE 5.2 (July 2004)
6 *
7 * (c) 2003-2004 by Heiko A. Schmidt, Korbinian Strimmer, and Arndt von Haeseler
8 * (c) 1999-2003 by Heiko A. Schmidt, Korbinian Strimmer,
9 * M. Vingron, and Arndt von Haeseler
10 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
11 *
12 * All parts of the source except where indicated are distributed under
13 * the GNU public licence. See http://www.opensource.org for details.
14 *
15 * ($Id$)
16 *
17 */
18
19 #ifdef HAVE_CONFIG_H
20 # include <config.h>
21 #endif
22
23 /*
24 * Use Scalable Parallel Random Number Generator
25 * See sprng.cs.fsu.edu for details
26 * from National Center for Supercomputing Applications (NCSA)
27 * See www.ncsa.uiuc.edu
28 */
29 #ifdef NOSPRNG
30 # ifdef SPRNG
31 # undef SPRNG
32 # endif
33 #else
34 # define SPRNG 1
35 #endif
36
37 #ifdef SPRNG
38 # ifdef PARALLEL
39 # define USE_MPI
40 # ifndef EXTERN
41 # define EXTERN extern
42 # endif
43 # include "ppuzzle.h"
44 # endif
45 # if 0
46 # define SIMPLE_SPRNG
47 # endif
48 # include "sprng/sprng.h"
49 #endif
50
51
52
53 #include "util.h"
54
55 #define STDOUT stdout
56 #ifndef PARALLEL /* because printf() runs significantly faster */
57 /* than fprintf(stdout) on an Apple McIntosh */
58 /* (HS) */
59 # define FPRINTF printf
60 # define STDOUTFILE
61 #else
62 # define FPRINTF fprintf
63 # define STDOUTFILE STDOUT,
64 extern int PP_NumProcs;
65 extern int PP_Myid;
66 long int PP_randn;
67 long int PP_rand;
68 void PP_Finalize();
69 #endif
70
71
72 /**********************************
73 * memory allocation error handler
74 */
75
maerror(char * message)76 void maerror(char *message)
77 {
78 fprintf(STDOUT, "\n\n\nUnable to proceed (lack of memory: %s)\n\n", message);
79 fprintf(STDOUT, "Hint for Macintosh users:\n");
80 fprintf(STDOUT, "Use the <Get Info> command of the Finder to increase the memory partition!\n\n");
81 # if PARALLEL
82 PP_Finalize();
83 # endif
84 exit(1);
85 } /* maerror */
86
87
88 /******************************************************
89 * memory allocate double vectors, matrices, and cubes
90 */
91
new_dvector(int n)92 dvector new_dvector(int n)
93 {
94 dvector v;
95
96 v = (dvector) calloc((size_t) n, sizeof(double));
97 if (v == NULL) maerror("step 1 in new_dvector");
98
99 return v;
100 } /* new_dvector */
101
102
103 /******************/
104
new_dmatrix(int nrow,int ncol)105 dmatrix new_dmatrix(int nrow, int ncol)
106 {
107 int i;
108 dmatrix m;
109
110 m = (dmatrix) calloc((size_t) nrow, sizeof(dvector));
111 if (m == NULL) maerror("step 1 in in new_dmatrix");
112
113 *m = (dvector) calloc((size_t) (nrow * ncol), sizeof(double));
114 if (*m == NULL) maerror("step 2 in in new_dmatrix");
115
116 for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
117
118 return m;
119 } /* new_dmatrix */
120
121
122 /******************/
123
new_dcube(int ntri,int nrow,int ncol)124 dcube new_dcube(int ntri, int nrow, int ncol)
125 {
126 int i, j;
127 dcube c;
128
129 c = (dcube) calloc((size_t) ntri, sizeof(dmatrix));
130 if (c == NULL) maerror("step 1 in in new_dcube");
131
132 *c = (dmatrix) calloc((size_t) (ntri * nrow), sizeof(dvector));
133 if (*c == NULL) maerror("step 2 in in new_dcube");
134
135 **c = (dvector) calloc((size_t) (ntri * nrow * ncol), sizeof(double));
136 if (**c == NULL) maerror("step 3 in in new_dcube");
137
138 for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
139
140 for (i = 1; i < ntri; i++) {
141 c[i] = c[i-1] + nrow;
142 c[i][0] = c[i-1][0] + nrow * ncol;
143 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
144 }
145
146 return c;
147 } /* new_dcube */
148
149
150 /******************/
151
free_dvector(dvector v)152 void free_dvector(dvector v)
153 {
154 free((double *) v);
155 } /* free_dvector */
156
157
158 /******************/
159
free_dmatrix(dmatrix m)160 void free_dmatrix(dmatrix m)
161 {
162 free((double *) *m);
163 free((double *) m);
164 } /* free_dmatrix */
165
166
167 /******************/
168
free_dcube(dcube c)169 void free_dcube(dcube c)
170 {
171 free((double *) **c);
172 free((double *) *c);
173 free((double *) c);
174 } /* free_dcube */
175
176
177 /****************************************************
178 * memory allocate char vectors, matrices, and cubes
179 */
180
new_cvector(int n)181 cvector new_cvector(int n)
182 {
183 cvector v;
184
185 v = (cvector) calloc((size_t) n, sizeof(char));
186 if (v == NULL) maerror("step1 in new_cvector");
187
188 return v;
189 } /* new_cvector */
190
191
192 /******************/
193
new_cmatrix(int nrow,int ncol)194 cmatrix new_cmatrix(int nrow, int ncol)
195 {
196 int i;
197 cmatrix m;
198
199 m = (cmatrix) calloc((size_t) nrow, sizeof(cvector));
200 if (m == NULL) maerror("step 1 in new_cmatrix");
201
202 *m = (cvector) calloc((size_t) (nrow * ncol), sizeof(char));
203 if (*m == NULL) maerror("step 2 in new_cmatrix");
204
205 for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
206
207 return m;
208 } /* new_cmatrix */
209
210
211 /******************/
212
new_ccube(int ntri,int nrow,int ncol)213 ccube new_ccube(int ntri, int nrow, int ncol)
214 {
215 int i, j;
216 ccube c;
217
218 c = (ccube) calloc((size_t) ntri, sizeof(cmatrix));
219 if (c == NULL) maerror("step 1 in new_ccube");
220
221 *c = (cmatrix) calloc((size_t) (ntri * nrow), sizeof(cvector));
222 if (*c == NULL) maerror("step 2 in new_ccube");
223
224 **c = (cvector) calloc((unsigned) (ntri * nrow * ncol), sizeof(char));
225 if (**c == NULL) maerror("step 3 in new_ccube");
226
227 for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
228
229 for (i = 1; i < ntri; i++) {
230 c[i] = c[i-1] + nrow;
231 c[i][0] = c[i-1][0] + nrow * ncol;
232 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
233 }
234
235 return c;
236 } /* new_ccube */
237
238
239 /******************/
240
free_cvector(cvector v)241 void free_cvector(cvector v)
242 {
243 free((char *) v);
244 } /* free_cvector */
245
246
247 /******************/
248
free_cmatrix(cmatrix m)249 void free_cmatrix(cmatrix m)
250 {
251 free((char *) *m);
252 free((char *) m);
253 } /* free_cmatrix */
254
255
256 /******************/
257
free_ccube(ccube c)258 void free_ccube(ccube c)
259 {
260 free((char *) **c);
261 free((char *) *c);
262 free((char *) c);
263 } /* free_ccube */
264
265
266 /***************************************************
267 * memory allocate int vectors, matrices, and cubes
268 */
269
new_ivector(int n)270 ivector new_ivector(int n)
271 {
272 ivector v;
273
274 v = (ivector) calloc((size_t) n, sizeof(int));
275 if (v == NULL) maerror("step 1 in new_ivector");
276
277 return v;
278 } /* new_ivector */
279
280
281 /******************/
282
new_imatrix(int nrow,int ncol)283 imatrix new_imatrix(int nrow, int ncol)
284 {
285 int i;
286 imatrix m;
287
288 m = (imatrix) calloc((size_t) nrow, sizeof(ivector));
289 if (m == NULL) maerror("step 1 in new_imatrix");
290
291 *m = (ivector) calloc((size_t) (nrow * ncol), sizeof(int));
292 if (*m == NULL) maerror("step 2 in new_imatrix");
293
294 for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
295
296 return m;
297 } /* new_imatrix */
298
299
300 /******************/
301
new_icube(int ntri,int nrow,int ncol)302 icube new_icube(int ntri, int nrow, int ncol)
303 {
304 int i, j;
305 icube c;
306
307 c = (icube) calloc((size_t) ntri, sizeof(imatrix));
308 if (c == NULL) maerror("step 1 in new_icube");
309
310 *c = (imatrix) calloc((size_t) (ntri * nrow), sizeof(ivector));
311 if (*c == NULL) maerror("step 2 in new_icube");
312
313 **c = (ivector) calloc((size_t) (ntri * nrow * ncol), sizeof(int));
314 if (**c == NULL) maerror("step 3 in new_icube");
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 } /* new_icube */
326
327
328 /******************/
329
free_ivector(ivector v)330 void free_ivector(ivector v)
331 {
332 free((int *) v);
333 } /* free_ivector */
334
335
336 /******************/
337
free_imatrix(imatrix m)338 void free_imatrix(imatrix m)
339 {
340 free((int *) *m);
341 free((int *) m);
342 } /* free_imatrix */
343
344
345 /******************/
346
free_icube(icube c)347 void free_icube(icube c)
348 {
349 free((int *) **c);
350 free((int *) *c);
351 free((int *) c);
352 } /* free_icube */
353
354
355 /***************************************************
356 * memory allocate uli vectors, matrices, and cubes
357 */
358
new_ulivector(int n)359 ulivector new_ulivector(int n)
360 {
361 ulivector v;
362
363 v = (ulivector) calloc((size_t) n, sizeof(uli));
364 if (v == NULL) maerror("step 1 in new_ulivector");
365
366 return v;
367 } /* new_ulivector */
368
369
370 /******************/
371
new_ulimatrix(int nrow,int ncol)372 ulimatrix new_ulimatrix(int nrow, int ncol)
373 {
374 int i;
375 ulimatrix m;
376
377 m = (ulimatrix) calloc((size_t) nrow, sizeof(ulivector));
378 if (m == NULL) maerror("step 1 in new_ulimatrix");
379
380 *m = (ulivector) calloc((size_t) (nrow * ncol), sizeof(uli));
381 if (*m == NULL) maerror("step 2 in new_ulimatrix");
382
383 for (i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;
384
385 return m;
386 } /* new_ulimatrix */
387
388
389 /******************/
390
new_ulicube(int ntri,int nrow,int ncol)391 ulicube new_ulicube(int ntri, int nrow, int ncol)
392 {
393 int i, j;
394 ulicube c;
395
396 c = (ulicube) calloc((size_t) ntri, sizeof(ulimatrix));
397 if (c == NULL) maerror("step 1 in new_ulicube");
398
399 *c = (ulimatrix) calloc((size_t) (ntri * nrow), sizeof(ulivector));
400 if (*c == NULL) maerror("step 2 in new_ulicube");
401
402 **c = (ulivector) calloc((size_t) (ntri * nrow * ncol), sizeof(uli));
403 if (**c == NULL) maerror("step 3 in new_ulicube");
404
405 for (j = 1; j < nrow; j++) c[0][j] = c[0][j-1] + ncol;
406
407 for (i = 1; i < ntri; i++) {
408 c[i] = c[i-1] + nrow;
409 c[i][0] = c[i-1][0] + nrow * ncol;
410 for (j = 1; j < nrow; j++) c[i][j] = c[i][j-1] + ncol;
411 }
412
413 return c;
414 } /* new_ulicube */
415
416
417 /******************/
418
free_ulivector(ulivector v)419 void free_ulivector(ulivector v)
420 {
421 free((uli *) v);
422 } /* free_ulivector */
423
424
425 /******************/
426
free_ulimatrix(ulimatrix m)427 void free_ulimatrix(ulimatrix m)
428 {
429 free((uli *) *m);
430 free((uli *) m);
431 } /* free_ulimatrix */
432
433
434 /******************/
435
free_ulicube(ulicube c)436 void free_ulicube(ulicube c)
437 {
438 free((uli *) **c);
439 free((uli *) *c);
440 free((uli *) c);
441 } /* free_ulicube */
442
443
444 /******************************************************************************/
445 /* check time */
446 /******************************************************************************/
447 #define TIMECHECK_INTERVALL 300
448 /* check timer and generate message roughly every TIMECHECK_INTERVALL seconds */
checktime(time_t * starttime,time_t * lasttime,time_t * nowtime,uli done,uli sumtodo,int * mflag)449 void checktime(time_t *starttime, /* starttime of overall part */
450 time_t *lasttime, /* lasttime of print out */
451 time_t *nowtime, /* current time (taken inside) */
452 uli done, /* number of tasks done */
453 uli sumtodo, /* sum of task todo */
454 int *mflag) /* LF needed to end line? */
455 {
456 double tc2, mintogo, minutes, hours;
457
458 /* check timer */
459 time(nowtime);
460 if ( (*nowtime - *lasttime) > TIMECHECK_INTERVALL) {
461 /* every TIMECHECK_INTERVAL seconds */
462 /* percentage of completed tasks */
463 if (*mflag == 0) {
464 fprintf(STDOUT, "\n");
465 *mflag = 1;
466 }
467 tc2 = 100.0 * done/sumtodo;
468 mintogo = (100.0-tc2) *
469 (double) (*nowtime-*starttime)/60.0/tc2;
470 hours = floor(mintogo/60.0);
471 minutes = mintogo - 60.0*hours;
472 fprintf(STDOUT, "%2.2f%%", tc2);
473 fprintf(STDOUT, " completed (remaining");
474 fprintf(STDOUT, " time: %.0f", hours);
475 fprintf(STDOUT, " hours %.0f", minutes);
476 fprintf(STDOUT, " minutes)\n");
477 fflush(STDOUT);
478 *lasttime = *nowtime;
479 }
480 /* else {fprintf(stdout,"t%ld,", *nowtime - *lasttime); fflush(stdout);} */
481
482 } /* checktime */
483
484
485 #ifndef SPRNG
486
487 /******************************************************************************/
488 /* random numbers generator (Numerical recipes) */
489 /******************************************************************************/
490
491 /* variable */
492 long idum;
493
494 /* definitions */
495 #define IM1 2147483563
496 #define IM2 2147483399
497 #define AM (1.0/IM1)
498 #define IMM1 (IM1-1)
499 #define IA1 40014
500 #define IA2 40692
501 #define IQ1 53668
502 #define IQ2 52774
503 #define IR1 12211
504 #define IR2 3791
505 #define NTAB 32
506 #define NDIV (1+IMM1/NTAB)
507 #define EPS 1.2e-7
508 #define RNMX (1.0-EPS)
509
randomunitintervall()510 double randomunitintervall()
511 /* Long period (> 2e18) random number generator. Returns a uniform random
512 deviate between 0.0 and 1.0 (exclusive of endpoint values).
513
514 Source:
515 Press et al., "Numerical recipes in C", Cambridge University Press, 1992
516 (chapter 7 "Random numbers", ran2 random number generator) */
517 {
518 int j;
519 long k;
520 static long idum2=123456789;
521 static long iy=0;
522 static long iv[NTAB];
523 double temp;
524
525 if (idum <= 0) {
526 if (-(idum) < 1)
527 idum=1;
528 else
529 idum=-(idum);
530 idum2=(idum);
531 for (j=NTAB+7;j>=0;j--) {
532 k=(idum)/IQ1;
533 idum=IA1*(idum-k*IQ1)-k*IR1;
534 if (idum < 0)
535 idum += IM1;
536 if (j < NTAB)
537 iv[j] = idum;
538 }
539 iy=iv[0];
540 }
541 k=(idum)/IQ1;
542 idum=IA1*(idum-k*IQ1)-k*IR1;
543 if (idum < 0)
544 idum += IM1;
545 k=idum2/IQ2;
546 idum2=IA2*(idum2-k*IQ2)-k*IR2;
547 if (idum2 < 0)
548 idum2 += IM2;
549 j=iy/NDIV;
550 iy=iv[j]-idum2;
551 iv[j] = idum;
552 if (iy < 1)
553 iy += IMM1;
554 if ((temp=AM*iy) > RNMX)
555 return RNMX;
556 else
557 return temp;
558 } /* randomunitintervall */
559
560 #undef IM1
561 #undef IM2
562 #undef AM
563 #undef IMM1
564 #undef IA1
565 #undef IA2
566 #undef IQ1
567 #undef IQ2
568 #undef IR1
569 #undef IR2
570 #undef NTAB
571 #undef NDIV
572 #undef EPS
573 #undef RNMX
574
575
initrandom(int seed)576 int initrandom(int seed) /* RAND4 */
577 {
578 srand((unsigned) time(NULL));
579 if (seed < 0)
580 seed = rand();
581 idum=-(long) seed;
582 fflush(stderr);
583 # ifndef PARALLEL
584 fprintf(stderr,"Using RAND4 Random Number Generator\n");
585 # ifdef RANDVERBOSE1
586 fprintf(stderr, "!!! random seed set to %d !!!\n", seed);
587 # endif /* RANDVERBOSE1 */
588 # else /* PARALLEL */
589 {
590 int n;
591 for (n=0; n<PP_Myid; n++)
592 (void) randomunitintervall();
593 # ifdef RANDVERBOSE1
594 fprintf(stderr, "(%2d) !!! random seed set to %d, %dx drawn !!!\n", PP_Myid, seed, n);
595 # endif /* RANDVERBOSE1 */
596 }
597 # endif
598 return (seed);
599 } /* initrandom */
600
601 /******************/
602
603 #else /* USE_SPRNG */
604
605 /******************/
606
607 int *randstream;
608
initrandom(int seed)609 int initrandom(int seed)
610 {
611 srand((unsigned) time(NULL));
612 if (seed < 0)
613 seed = make_sprng_seed();
614 # ifndef PARALLEL
615 fprintf(stderr,"Using SPRNG -- Scalable Parallel Random Number Generator\n");
616 fflush(stderr);
617 randstream = init_sprng(0,1,seed,SPRNG_DEFAULT); /*init stream*/
618
619 # ifdef RANDVERBOSE1
620 fprintf(stderr, "!!! random seed set to %d !!!\n", seed);
621 printf(" Printing information about new stream\n");
622 print_sprng(randstream);
623 # endif /* RANDVERBOSE 1 */
624 # else /* PARALLEL */
625 if (PP_IamMaster) {
626 fprintf(stderr,"Using SPRNG -- Scalable Parallel Random Number Generator\n");
627 fflush(stderr);
628 }
629 /* MPI_Bcast(&seed, 1, MPI_UNSIGNED, PP_MyMaster, MPI_COMM_WORLD); */
630 randstream = init_sprng(PP_Myid,PP_NumProcs,seed,SPRNG_DEFAULT); /*initialize stream*/
631 # ifdef RANDVERBOSE1
632 fprintf(stderr, "(%2d) !!! random seed set to %d !!!\n", PP_Myid, seed);
633 printf(" Printing information about new stream\n");
634 print_sprng(randstream);
635 # endif /* RANDVERBOSE1 */
636 # endif /* PARALLEL */
637 return (seed);
638 } /* initrandom */
639
640
641 #endif /* USE_SPRNG */
642
643 /******************/
644
645
646 /* returns a random integer in the range [0; n - 1] */
randominteger(int n)647 int randominteger(int n)
648 {
649 # ifndef FIXEDINTRAND
650 # ifndef PARALLEL
651 # ifdef SPRNG
652 return (int) floor(sprng(randstream)*n);
653 # else /* NO_SPRNG */
654 int t;
655 t = (int) floor(randomunitintervall()*n);
656 return t;
657 # endif /* NO_SPRNG */
658 # else /* NOT PARALLEL */
659 # ifdef SPRNG
660 return (int) floor(sprng(randstream)*n);
661 # else /* NO_SPRNG */
662 int m;
663 for (m=1; m<PP_NumProcs; m++)
664 (void) randomunitintervall();
665 PP_randn+=(m-1); PP_rand++;
666 return (int) floor(randomunitintervall()*n);
667 # endif /* NO_SPRNG */
668 # endif /* NOT PARALLEL */
669 # else /* FIXEDINTRAND */
670 fprintf(stderr, "!!! fixed \"random\" integers for testing purposes !!!\n");
671 return (int)0;
672 # endif /* FIXEDINTRAND */
673 } /* randominteger */
674
675
676 /******************/
677
678 /* Draw s numbers from the set 0,1,2,..,t-1 and put them
679 into slist (every number can be drawn only one time) */
chooser(int t,int s,ivector slist)680 void chooser(int t, int s, ivector slist)
681 {
682 int i;
683 int randidx;
684 ivector isfree;
685
686 isfree = new_ivector(t);
687 for (i = 0; i < t; i++) isfree[i] = i;
688
689 for (i = t; i > t-s; i--) {
690 /* random integer in the range [0, still to do = i-1] */
691 randidx = randominteger(i);
692
693 /* integer from array of free ints to slist */
694 slist[t - i] = isfree[randidx];
695
696 /* reduce list by moving last free to last chosen */
697 isfree[randidx] = isfree[i-1];
698 }
699 free_ivector(isfree);
700 } /* chooser */
701
702
703
704 /******************/
705
706 /* Draw s numbers from the set 0,1,2,..,t-1 and put them
707 into slist (every number can be drawn only one time) */
fixedchooser(int t,int s,ivector slist)708 void fixedchooser(int t, int s, ivector slist)
709 {
710 int i;
711
712 for (i = 0; i < s; i++) {
713 /* random integer in the range [0, s-1] */
714 slist[i] = i;
715 }
716 } /* fixedchooser */
717
718
719 /******************************************************************************/
720 /* Saver variants of ANSI C routines */
721 /******************************************************************************/
722
723 /* a realloc function that works also on non-ANSI compilers */
myrealloc(void * p,size_t s)724 void *myrealloc(void *p, size_t s)
725 {
726 if (p == NULL) return malloc(s);
727 else return realloc(p, s);
728 } /* myrealloc */
729
730
731 /******************/
732
733 /* safer variant of gets */
734 /* Reads characters from stdin until a newline character or EOF
735 is received. The newline is not made part of the string.
736 If an error occurs a null string \0 is returned */
mygets()737 cvector mygets()
738 {
739 int c, n;
740 cvector str;
741
742 str = new_cvector(100);
743
744 n = 0;
745 c = getchar();
746 while ((c != '\n') && (c != '\r')
747 && (n < 99) && (c != EOF)
748 && !ferror(stdin))
749 {
750 str[n] = (char) c;
751 c = getchar();
752 n++;
753 }
754 if (c == EOF || ferror(stdin))
755 {
756 str[0] = '\0';
757 }
758 else
759 {
760 str[n] = '\0';
761 }
762
763 return str;
764 } /* mygets */
765
766
767
768 /******************************************************************************/
769 /* minimization of a function by Brents method (Numerical Recipes) */
770 /******************************************************************************/
771
772 double brent(double, double, double, double (*f )(double ), double, double *, double *, double, double, double, int *);
773
774
775 #define ITMAX 100
776 #define CGOLD 0.3819660
777 #define GOLD 1.618034
778 #define GLIMIT 100.0
779 #define TINY 1.0e-20
780 #define ZEPS 1.0e-10
781 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
782 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
783
784 /* Brents method in one dimension */
brent(double ax,double bx,double cx,double (* f)(double),double tol,double * foptx,double * f2optx,double fax,double fbx,double fcx,int * converged)785 double brent(
786 double ax, /* three x-values ax < bx < cx */
787 double bx, /* or ax > bx > cx */
788 double cx,
789 double (*f)(double), /* function to minimize, f(bx)=min(f(ax),f(bx),f(cx)) */
790 double tol, /* describes the precision */
791 double *foptx,
792 double *f2optx,
793 double fax, /* f(ax), already precomputed outside of brent() */
794 double fbx, /* f(bx), already precomputed outside of brent() */
795 double fcx, /* f(cx), already precomputed outside of brent() */
796 int *converged) /* not converged error flag */
797 {
798 int iter;
799 double a,b,d=0,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
800 double xw,wv,vx;
801 double e=0.0;
802
803 *converged=TRUE;
804
805 a=(ax < cx ? ax : cx); /* is the order correct? If not, exchange ax and cx */
806 b=(ax > cx ? ax : cx);
807
808 x=bx; /* store 'middle' x, and y=fx */
809 fx=fbx;
810
811 if (fax < fcx) {
812 w=ax;
813 fw=fax;
814 v=cx;
815 fv=fcx;
816 } else {
817 w=cx;
818 fw=fcx;
819 v=ax;
820 fv=fax;
821 }
822
823 for (iter=1;iter<=ITMAX;iter++) { /* iteration loop */
824 xm=0.5*(a+b);
825 tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
826 /* is precision OK? I.e., are a and b close enough. */
827 if (fabs(x-xm) <= (tol2-0.5*(b-a))) { /* yes -> finish */
828 *foptx = fx; /* return current f(x) */
829 if ((x==w) && (x==v)) {
830 /* to avoid FPE errors by zero-division */
831 v=x-ZEPS;
832 w=x+ZEPS;
833 fv=(*f)(v);
834 fw=(*f)(w);
835 }
836 xw = x-w;
837 wv = w-v;
838 vx = v-x;
839 *f2optx = 2.0*(fv*xw + fx*wv + fw*vx)/
840 (v*v*xw + x*x*wv + w*w*vx); /* return S.E. of f(x) */
841 #if 0
842 fprintf(stderr, "XXXX: x=%f f(x)=%f=%f S.E.=%f w=%f v=%f\n", x, *foptx, fx, *f2optx, w, v);
843 #endif
844 return x; /* return minimum x-value */
845 }
846 if (fabs(e) > tol1) { /* compute a parabolic fit */
847 r=(x-w)*(fx-fv);
848 q=(x-v)*(fx-fw);
849 p=(x-v)*q-(x-w)*r;
850 q=2.0*(q-r);
851 if (q > 0.0) p = -p;
852 q=fabs(q);
853 etemp=e;
854 e=d;
855 /* determine acceptability of the parabolic fit */
856 if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
857 /* compute golden section of the larger segment max(a-x,b-x) */
858 d=CGOLD*(e=(x >= xm ? a-x : b-x));
859 else {
860 /* accept the parabolic step */
861 d=p/q;
862 u=x+d;
863 if (u-a < tol2 || b-u < tol2)
864 d=SIGN(tol1,xm-x);
865 }
866 } else {
867 /* compute golden section of the larger segment max(a-x,b-x) */
868 d=CGOLD*(e=(x >= xm ? a-x : b-x));
869 }
870 u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
871 fu=(*f)(u);
872
873 /* new u better than old x ? */
874 if (fu <= fx) { /* yes -> accept it */
875 if (u >= x) a=x; else b=x;
876 SHFT(v,w,x,u) /* ...and use the new values */
877 SHFT(fv,fw,fx,fu)
878 } else {
879 if (u < x) a=u; else b=u;
880 if (fu <= fw || w == x) {
881 v=w;
882 w=u;
883 fv=fw;
884 fw=fu;
885 } else if (fu <= fv || v == x || v == w) {
886 v=u;
887 fv=fu;
888 }
889 }
890 } /* end loop */
891
892 /* the function should never get here, i.e., the search did not converge within ITMAX iterations */
893 *converged=FALSE;
894
895 *foptx = fx; /* return f(x) */
896 if ((x==w) && (x==v)) { /* regularize to avoid FPE errors by zero-division */
897 v=x-ZEPS;
898 w=x+ZEPS;
899 fv=(*f)(v);
900 fw=(*f)(w);
901 }
902 xw = x-w;
903 wv = w-v;
904 vx = v-x;
905 *f2optx = 2.0*(fv*xw + fx*wv + fw*vx)/
906 (v*v*xw + x*x*wv + w*w*vx); /* return S.E. of f(x) */
907 #if 0
908 fprintf(stderr, "YYYY: x=%f f(x)=%f=%f S.E.=%f w=%f v=%f\n", x, *foptx, fx, *f2optx, w, v);
909 #endif
910 return x; /* return minimum x-value */
911 } /* brent */
912 #undef ITMAX
913 #undef CGOLD
914 #undef ZEPS
915 #undef SHFT
916 #undef SIGN
917 #undef GOLD
918 #undef GLIMIT
919 #undef TINY
920
921
922 /******************/
923
924 /* one-dimensional minimization - as input a lower and an upper limit and a trial
925 value for the minimum is needed: xmin < xguess < xmax
926 the function and a fractional tolerance has to be specified
927 onedimenmin returns the optimal x value and the value of the function
928 and its second derivative at this point
929 */
onedimenmin(double xmin,double xguess,double xmax,double (* f)(double),double tol,double * fx,double * f2x)930 double onedimenmin(
931 double xmin, double xguess, double xmax, double (*f)(double),
932 double tol, double *fx, double *f2x)
933 {
934 double eps, optx, ax, bx, cx, fa, fb, fc;
935 int converged; /* not converged error flag */
936
937 /* first attempt to bracketize minimum */
938 eps = xguess*tol*50.0;
939 ax = xguess - eps;
940 if (ax < xmin) ax = xmin;
941 bx = xguess;
942 cx = xguess + eps;
943 if (cx > xmax) cx = xmax;
944
945 /* check if this works */
946 fa = (*f)(ax);
947 fb = (*f)(bx);
948 fc = (*f)(cx);
949
950 /* if it works use these borders else be conservative */
951 if ((fa < fb) || (fc < fb)) {
952 if (ax != xmin) fa = (*f)(xmin);
953 if (cx != xmax) fc = (*f)(xmax);
954 optx = brent(xmin, xguess, xmax, f, tol, fx, f2x, fa, fb, fc, &converged);
955 } else
956 optx = brent(ax, bx, cx, f, tol, fx, f2x, fa, fb, fc, &converged);
957
958 return optx; /* return optimal x */
959 } /* onedimenmin */
960
961
962 /******************/
963
964 /* two-dimensional minimization with borders and calculations of standard errors */
965 /* we optimize along basis vectors - not very optimal but it seems to work well */
twodimenmin(double tol,int active1,double min1,double * x1,double max1,double (* func1)(double),double * err1,int active2,double min2,double * x2,double max2,double (* func2)(double),double * err2)966 void twodimenmin(double tol,
967 int active1, double min1, double *x1, double max1, double (*func1)(double), double *err1,
968 int active2, double min2, double *x2, double max2, double (*func2)(double), double *err2)
969 {
970 int it, nump, change;
971 double x1old, x2old;
972 double fx, f2x;
973
974 it = 0;
975 nump = 0;
976
977 /* count number of parameters */
978 if (active1) nump++;
979 if (active2) nump++;
980
981 do { /* repeat until nothing changes any more */
982 it++;
983 change = FALSE;
984
985 /* optimize first variable */
986 if (active1) {
987
988 if ((*x1) <= min1) (*x1) = min1 + 0.2*(max1-min1);
989 if ((*x1) >= max1) (*x1) = max1 - 0.2*(max1-min1);
990 x1old = (*x1);
991 (*x1) = onedimenmin(min1, (*x1), max1, func1, tol, &fx, &f2x);
992 if ((*x1) < min1) (*x1) = min1;
993 if ((*x1) > max1) (*x1) = max1;
994 /* same tolerance as 1D minimization */
995 if (fabs((*x1) - x1old) > 3.3*tol) change = TRUE;
996
997 /* standard error */
998 f2x = fabs(f2x);
999 if (1.0/(max1*max1) < f2x) (*err1) = sqrt(1.0/f2x);
1000 else (*err1) = max1;
1001
1002 }
1003
1004 /* optimize second variable */
1005 if (active2) {
1006
1007 if ((*x2) <= min2) (*x2) = min2 + 0.2*(max2-min2);
1008 if ((*x2) >= max2) (*x2) = max2 - 0.2*(max2-min2);
1009 x2old = (*x2);
1010 (*x2) = onedimenmin(min2, (*x2), max2, func2, tol, &fx, &f2x);
1011 if ((*x2) < min2) (*x2) = min2;
1012 if ((*x2) > max2) (*x2) = max2;
1013 /* same tolerance as 1D minimization */
1014 if (fabs((*x2) - x2old) > 3.3*tol) change = TRUE;
1015
1016 /* standard error */
1017 f2x = fabs(f2x);
1018 if (1.0/(max2*max2) < f2x) (*err2) = sqrt(1.0/f2x);
1019 else (*err2) = max2;
1020
1021 }
1022
1023 if (nump == 1) return;
1024
1025 } while (it != MAXITS && change);
1026
1027 return;
1028 } /* twodimenmin */
1029
1030