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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173 | #ifndef C_OPTIM_H
#define C_OPTIM_H
#include <Rcpp.h>
namespace rstpm2 {
typedef double optimfn(int, double *, void *);
typedef void optimgr(int, double *, double *, void *);
/* type of pointer to the target and gradient functions for Nlm */
typedef void (*fcn_p)(int, double *, double *, void *);
/* type of pointer to the hessian functions for Nlm */
typedef void (*d2fcn_p)(int, int, double *, double *, void *);
double min(double a, double b);
double max(double a, double b);
double bound(double x, double lower, double upper);
/**
Adapt a function object (functor) for NelderMead and BFGS
**/
template<class T>
double adapt_functor(int n, double * beta, void * par) {
T * model = (T *) par;
Rcpp::NumericVector x(beta,beta+n);
return model->operator()(x);
}
/**
Adapt an objective function for NelderMead and BFGS
**/
template<class T>
double adapt_objective(int n, double * beta, void * par) {
T * model = (T *) par;
Rcpp::NumericVector x(beta,beta+n);
return model->objective(x);
}
/**
Adapt a gradient function for BFGS
**/
template<class T>
void adapt_gradient(int n, double * beta, double * grad, void * par) {
T * model = (T *) par;
Rcpp::NumericVector x(beta,beta+n);
Rcpp::NumericVector vgrad = model->gradient(x);
for (int i=0; i<n; ++i) grad[i] = vgrad[i];
}
class NelderMead {
public:
NelderMead(int trace = 0, int maxit = 500,
double abstol = - INFINITY,
double reltol = 1.0e-8,
double alpha = 1.0, double beta = 0.5, double gamma = 2.0,
double epshess = 6.055454e-06, bool hessianp = true);
virtual void optim(optimfn fn, Rcpp::NumericVector init, void * ex);
template<class T>
void optim(Rcpp::NumericVector init, T object) {
optim(&adapt_objective<T>,init,(void *) &object);
}
virtual Rcpp::NumericMatrix calc_hessian(optimfn fn, void * ex);
int n, trace, maxit, fail, fncount;
double abstol, reltol, alpha, beta, gamma, Fmin, epshess;
bool hessianp;
Rcpp::NumericVector coef;
Rcpp::NumericMatrix hessian;
};
class BFGS {
public:
BFGS(int trace = 0, int maxit = 100,
double abstol = - INFINITY,
double reltol = 1.0e-8, int report = 10, double epshess = 1.0e-8, bool hessianp = true);
virtual void optim(optimfn fn, optimgr gr, Rcpp::NumericVector init, void * ex);
virtual void optim(int n, optimfn fn, optimgr gr, double * init, void * ex);
virtual double calc_objective(optimfn fn, Rcpp::NumericVector coef, void * ex);
virtual double calc_objective(optimfn fn, void * ex);
virtual Rcpp::NumericMatrix calc_hessian(optimgr gr, void * ex);
template<class T>
void optim(Rcpp::NumericVector init, T object) {
optim(&adapt_objective<T>,&adapt_gradient<T>,init,(void *) &object);
}
int n, trace, maxit, report, fncount, grcount, fail;
double abstol, reltol, Fmin, epshess;
bool hessianp;
Rcpp::NumericVector coef;
Rcpp::NumericMatrix hessian;
};
class Nlm {
public:
Nlm(double fscale = 1.0, // nlm()
int method = 2, // cf. nlm: method=1
int iexp = 1, // nlm()
int msg = 9, // nlm()
int ndigit = 12, // nlm()
int itnlim = 50, // nlm()
int iagflg = 1, // nlm()
int iahflg = 0, // nlm()
double dlt = 1.0, // nlm
double gradtl = 1.0e-6, // nlm()
double stepmx = 0.0, // set to -1.0 to get nlm()'s behaviour
double steptl = 1.0e-6, // nlm()
double epshess = 6.055454e-06,
int itrmcd = 0,
int itncnt = 0,
bool hessianp = true);
void optim(fcn_p fcn, fcn_p d1fcn, Rcpp::NumericVector init, void * state);
void optim(fcn_p fcn, Rcpp::NumericVector init, void * state);
double calc_objective(fcn_p fn, Rcpp::NumericVector coef, void * ex);
double calc_objective(fcn_p fn, void * ex);
Rcpp::NumericMatrix calc_hessian(fcn_p gr, void * ex);
void set_print_level(int);
double fscale;
int method;
int iexp;
int msg;
int ndigit;
int itnlim;
int iagflg;
int iahflg;
double dlt;
double gradtl;
double stepmx;
double steptl;
double epshess;
int itrmcd;
int itncnt;
bool hessianp;
Rcpp::NumericVector coef;
Rcpp::NumericMatrix hessian;
};
typedef double (*Brent_fminfn)(double, void *);
double Brent_fmin(double ax, double bx, double (*f)(double, void *),
void *info, double tol);
double R_zeroin2( /* An estimate of the root */
double ax, /* Left border | of the range */
double bx, /* Right border| the root is seeked*/
double fa, double fb, /* f(a), f(b) */
double (*f)(double x, void *info), /* Function under investigation */
void *info, /* Add'l info passed on to f */
double *Tol, /* Acceptable tolerance */
int *Maxit); /* Max # of iterations */
/**
Adapt a function object (functor) to work with Brent_fmin()
**/
template<class T, class X>
double Brent_fmin_functor(X x, void * par) {
T * model = (T *) par;
return model->operator()(x);
}
/**
Use Brent_fmin with a function object (functor)
**/
template<class T>
double BrentFmin(double a, double b, T obj, double eps = 1.0e-8) {
return Brent_fmin(a,b,&Brent_fmin_functor<T,double>,(void *) &obj,eps);
}
Rcpp::NumericMatrix qr_q(const Rcpp::NumericMatrix& X, double tol = 1E-12);
} // anonymous rstpm2
#endif /* c_optim_h */
|