TMB Documentation  v1.9.1
convenience.hpp
Go to the documentation of this file.
1 // Copyright (C) 2013-2015 Kasper Kristensen
2 // License: GPL-2
3 
12 template <class Type>
14  if (x.size() != fac.size()) Rf_error("x and fac must have equal length.");
15  int nlevels = 0;
16  for (int i = 0; i < fac.size(); i++)
17  if (fac[i] >= nlevels) nlevels = fac[i] + 1;
18  vector<vector<Type> > ans(nlevels);
19  vector<int> lngt(nlevels);
20  lngt.setZero();
21  for (int i = 0; i < fac.size(); i++) lngt[fac[i]]++;
22  for (int i = 0; i < nlevels; i++) ans[i].resize(lngt[i]);
23  lngt.setZero();
24  for (int i = 0; i < fac.size(); i++) {
25  ans[fac[i]][lngt[fac[i]]] = x[i];
26  lngt[fac[i]]++;
27  }
28  return ans;
29 }
30 
32 template <template <class> class Vector, class Type>
33 Type sum(Vector<Type> x) {
34  return x.sum();
35 }
36 
41 template <class Type>
43  return A * x.matrix();
44 }
45 
47 template <class Type>
48 vector<Type> operator*(Eigen::SparseMatrix<Type> A, vector<Type> x) {
49  return (A * x.matrix()).array();
50 }
51 
57 template <class Type>
58 Type pnorm_approx(Type x) {
59  Type a = 993. / 880.;
60  Type b = 89. / 880.;
61  x = x / sqrt(Type(2));
62  return Type(.5) * tanh((a + b * x * x) * x) + Type(.5);
63 }
65 
66 
71 template <class Type>
72 Type qnorm_approx(Type x) {
73  Type a = 993. / 880.;
74  Type b = 89. / 880.;
75  Type y = .5 * (log(x) - log(1 - x));
76  Type p = a / b;
77  Type q = -y / b;
78  Type Delta0 = -3 * p;
79  Type Delta1 = 27 * q;
80  Type C = pow(.5 * Delta1 + .5 * sqrt(pow(Delta1, 2) - 4 * pow(Delta0, 3)),
81  Type(1) / Type(3));
82  return -(C + Delta0 / C) * sqrt(Type(2)) / Type(3);
83 }
85 
86 
91 template <class Type>
93  int n = x.size();
94  vector<Type> ans(n - 1);
95  for (int i = 0; i < n - 1; i++) ans[i] = x[i + 1] - x[i];
96  return ans;
97 }
98 
103 template<class Type>
104 Type logit(Type x){
105  return log(x/(Type(1.0)-x));
106 }
108 
109 
113 template<class Type>
114 Type invlogit(Type x){
115  return Type(1.0)/(Type(1.0)+exp(-x));
116 }
118 
119 
128 template<class Type>
129 Type matern(Type u, Type phi, Type kappa){
130  Type x = CppAD::CondExpEq(u, Type(0), Type(1), u / phi); /* Avoid NaN when u=0 */
131  Type ans = 1.0 / ( exp(lgamma(kappa)) * pow(2, kappa - 1.0) ) * pow(x, kappa) * besselK(x, kappa);
132  return CppAD::CondExpEq(u, Type(0), Type(1), ans);
133 }
134 
136 template<class Type>
137 Type squeeze(Type u){
138  Type eps = std::numeric_limits<double>::epsilon();
139  u = (1.0 - eps) * (u - .5) + .5;
140  return u;
141 }
142 
144 template <class Type>
145 Type max(const vector<Type> &x)
146 {
147  Type res = x[0];
148  for(int i=0; i < x.size(); i++){
149  res = CppAD::CondExpGt(res, x[i], res, x[i]);
150  }
151  return res;
152 }
153 
155 template <class Type>
156 Type min(const vector<Type> &x)
157 {
158  Type res = x[0];
159  for(int i = 0; i < x.size(); i++){
160  res = CppAD::CondExpLt(res, x[i], res, x[i]);
161  }
162  return res;
163 }
164 
180 template<class Type>
181 Type logspace_add(Type logx, Type logy) {
182  if ( !CppAD::Variable(logx) && logx == Type(-INFINITY) )
183  return logy;
184  if ( !CppAD::Variable(logy) && logy == Type(-INFINITY) )
185  return logx;
186  CppAD::vector<Type> tx(3);
187  tx[0] = logx;
188  tx[1] = logy;
189  tx[2] = 0; // order
190  return atomic::logspace_add(tx)[0];
191 }
192 
204 template<class Type>
205 Type logspace_sub(Type logx, Type logy) {
206  CppAD::vector<Type> tx(3);
207  tx[0] = logx;
208  tx[1] = logy;
209  tx[2] = 0; // order
210  return atomic::logspace_sub(tx)[0];
211 }
Vector class used by TMB.
Definition: vector.hpp:17
Type logit(Type x)
#define VECTORIZE1_t(FUN)
Vectorize 1-argument functions.
Definition: Vectorize.hpp:63
Type pnorm_approx(Type x)
Approximate normal cumulative distribution function, similar to R&#39;s pnorm (one-argument case only)...
Definition: convenience.hpp:58
vector< Type > diff(vector< Type > x)
Definition: convenience.hpp:92
Type invlogit(Type x)
vector< Type > operator*(matrix< Type > A, vector< Type > x)
Definition: convenience.hpp:42
Type lgamma(Type x)
Logarithm of gamma function (following R argument convention).
Definition: lgamma.hpp:12
Type min(const vector< Type > &x)
Matrix class used by TMB.
Definition: vector.hpp:101
Type qnorm_approx(Type x)
Approximate inverse normal cumulative distribution function, similar to R&#39;s qnorm (one-argument case ...
Definition: convenience.hpp:72
Type squeeze(Type u)
Type besselK(Type x, Type nu)
besselK function (same as besselK from R).
Array class used by TMB.
Definition: array.hpp:22
vector< vector< Type > > split(vector< Type > x, vector< int > fac)
Similar to R&#39;s split function: split(x,fac) devides x into groups defined by fac .
Definition: convenience.hpp:13
Type sum(Vector< Type > x)
Definition: convenience.hpp:33
Type matern(Type u, Type phi, Type kappa)
Type logspace_add(Type logx, Type logy)
Addition in log-space.
Type logspace_sub(Type logx, Type logy)
Subtraction in log-space.
Type max(const vector< Type > &x)
License: GPL v2