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 | #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 */
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 negll function for NelderMead and BFGS
**/
template<class T>
double adapt_negll(int n, double * beta, void * par) {
T * model = (T *) par;
Rcpp::NumericVector x(beta,beta+n);
return model->negll(x);
}
/**
Adapt a grad_negll function for BFGS
**/
template<class T>
void adapt_grad_negll(int n, double * beta, double * grad, void * par) {
T * model = (T *) par;
Rcpp::NumericVector x(beta,beta+n);
//grad = model->grad_negll(x);
}
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_functor<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 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);
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()
int itrmcd = 0,
int itncnt = 0,
bool hessianp = true
);
void optim(fcn_p fcn, fcn_p d1fcn, Rcpp::NumericVector init, void * state); // assumes iahflg=0
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;
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);
/**
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);
}
} // anonymous rstpm2
#endif /* c_optim_h */
|