1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
/**
\file linmin.hpp
\brief linear search
\author Junhua Gu
*/
#ifndef LINMIN_HPP
#define LINMIN_HPP
#define OPT_HEADER
#include "mnbrak.hpp"
#include "brent.hpp"
#include <cmath>
#include <core/opt_traits.hpp>
namespace opt_utilities
{
template <typename rT,typename pT>
class func_adaptor
:public func_obj<rT,rT>
{
private:
const pT p1,xi1;
const func_obj<rT,pT>* pfoo;
func_adaptor(){}
func_adaptor(const func_adaptor&)
:func_obj<rT,rT>(),p1(),xi1(),pfoo(NULL_PTR)
{}
public:
/*
void set_origin(pT& p2)
{
p1=p2;
}
void set_direction(pT& xi2)
{
xi1=xi2;
}
void set_func_obj(func_obj<rT,pT>& foo)
{
pfoo=&foo;
}*/
public:
func_adaptor(const pT& _p,const pT& _xi,const func_obj<rT,pT>& pf)
:p1(_p),xi1(_xi),pfoo(&pf)
{}
private:
func_obj<rT,rT>* do_clone()const
{
return new func_adaptor(*this);
}
rT do_eval(const rT& x)
{
//assert(p1.size()==xi1.size());
pT xt;
opt_eq(xt,p1);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(size_t i=0;i<get_size(xt);++i)
{
//get_element(xt,i)+=x*get_element((pT)xi1,i);
set_element(xt,i,
get_element(xt,i)+x*get_element((pT)xi1,i));
//get_element((pT)xi1,i);
}
return const_cast<func_obj<rT,pT>&>(*pfoo).eval(xt);
//return x;
}
};
template<typename rT,typename pT>
void linmin(pT& p,pT& xi,rT& fret,func_obj<rT,pT>& func)
{
// assert(p.size()==10);
//assert(xi.size()==10);
func_adaptor<rT,pT> fadpt(p,xi,func);
int j=0;
const rT TOL=std::sqrt(std::numeric_limits<rT>::epsilon());
rT xx=0,xmin=0,fx=0,fb=0,fa=0,bx=0,ax=0;
int n=(int)get_size(p);
ax=0.;
xx=1.;
mnbrak(ax,xx,bx,fa,fx,fb,fadpt);
//cout<<xx<<endl;
fret=brent(ax,xx,bx,fadpt,TOL,xmin);
//cout<<xmin<<endl;
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(j=0;j<n;++j)
{
//get_element(xi,j)*=xmin;
set_element(xi,j,
get_element(xi,j)*xmin);
//get_element(p,j)+=get_element(xi,j);
set_element(p,j,
get_element(p,j)+get_element(xi,j));
}
// delete xicom_p;
//delete pcom_p;
}
}
#endif
|
