FEDRA emulsion software from the OPERA Collaboration
minfc.hh
Go to the documentation of this file.
1#ifndef __RMINFC
2#define __RMINFC
3// **********************************************************************
4//
5// source:
6//
7// type: source code
8//
9// created: 26.07.2001
10//
11// author: Thorsten Glebe
12// HERA-B Collaboration
13// Max-Planck-Institut fuer Kernphysik
14// Saupfercheckweg 1
15// 69117 Heidelberg
16// Germany
17// E-mail: T.Glebe@mpi-hd.mpg.de
18//
19// Description: Function which finds the minimum of a function with one
20// variable. The "golden section search" is applied. The method
21// uses a fixed number of function evaluations.
22//
23// changes:
24// 26 Jul 2001 (TG) creation from CERNLIB minfc.F file (f2c + handoptimization)
25//
26// **********************************************************************
27//#include <cmath>
28#include "Functions.hh"
29
43template <class T, class P>
44bool Rminfc( T (*f)(P), P a, P b, double eps, double delta, P& x, P& y) {
45
46 // Local variables
47 static P c, d, h;
48 static int n;
49 static P v, w, fv, fw;
50 bool lge = true, llt = true;
51
52
53 n = -1;
54 if (a != b) {
55 n = round(log(fabs(a - b) / eps) * 2.08);
56 }
57 c = a;
58 d = b;
59 if (a > b) {
60 c = b;
61 d = a;
62 }
63
64 do {
65 h = d - c;
66 if (llt == true) {
67 v = c + h * 0.3819660112501051;
68 fv = (*f)(v);
69 }
70 if (lge == true) {
71 w = c + h * 0.6180339887498949;
72 fw = (*f)(w);
73 }
74 if (fv < fw) {
75 llt = true;
76 lge = false;
77 d = w;
78 w = v;
79 fw = fv;
80 } else {
81 llt = false;
82 lge = true;
83 c = v;
84 v = w;
85 fv = fw;
86 }
87 --n;
88 } while (n >= 0);
89
90 x = (c + d) * .5;
91 y = (*f)(x);
92 return (fabs(x - a) > delta) && (fabs(x - b) > delta);
93} // Rminfc
94
95
111template <class C, class T, class P>
112bool RminfcM(const C& part, T (C::*f)(P) const, P a, P b, double eps, double delta, P& x, P& y) {
113
114 // Local variables
115 static P c, d, h;
116 static int n;
117 static P v, w, fv, fw;
118 bool lge = true, llt = true;
119
120
121 n = -1;
122 if (a != b) {
123 n = round(log(fabs(a - b) / eps) * 2.08);
124 }
125 c = a;
126 d = b;
127 if (a > b) {
128 c = b;
129 d = a;
130 }
131
132 do {
133 h = d - c;
134 if (llt == true) {
135 v = c + h * 0.3819660112501051;
136 fv = (part.*f)(v);
137 }
138 if (lge == true) {
139 w = c + h * 0.6180339887498949;
140 fw = (part.*f)(w);
141 }
142 if (fv < fw) {
143 llt = true;
144 lge = false;
145 d = w;
146 w = v;
147 fw = fv;
148 } else {
149 llt = false;
150 lge = true;
151 c = v;
152 v = w;
153 fv = fw;
154 }
155 --n;
156 } while (n >= 0);
157
158 x = (c + d) * .5;
159 y = (part.*f)(x);
160 return (fabs(x - a) > delta) && (fabs(x - b) > delta);
161} // Rminfc
162
163#endif
round
int round(const T &x)
Definition: Functions.hh:83
d
void d()
Definition: RecDispEX.C:381
f
FILE * f
Definition: RecDispMC.C:150
log
FILE * log
fabs
Expr< UnaryOp< Fabs< T >, Expr< A, T, D >, T >, T, D > fabs(const Expr< A, T, D > &rhs)
Definition: UnaryOperators.hh:96
a
void a()
Definition: check_aligned.C:59
Rminfc
bool Rminfc(T(*f)(P), P a, P b, double eps, double delta, P &x, P &y)
Definition: minfc.hh:44
RminfcM
bool RminfcM(const C &part, T(C::*f)(P) const, P a, P b, double eps, double delta, P &x, P &y)
Definition: minfc.hh:112
w
void w(int rid=2, int nviews=2)
Definition: test.C:27