|
| 1 | +#include "RcppArmadillo.h" |
| 2 | +#include "roptim.h" |
| 3 | +using namespace roptim; |
| 4 | + |
| 5 | +// [[Rcpp::depends(RcppArmadillo)]] |
| 6 | +// [[Rcpp::depends(roptim)]] |
| 7 | + |
| 8 | +std::vector<size_t> CantorExpansion(size_t n, std::vector<size_t> s) { |
| 9 | + std::vector<size_t> out(s.size()); |
| 10 | + std::vector<size_t>::iterator it; |
| 11 | + it = s.begin(); |
| 12 | + it = s.insert(it, 1); |
| 13 | + size_t G[s.size()]; |
| 14 | + std::partial_sum(s.begin(), s.end(), G, std::multiplies<size_t>()); |
| 15 | + size_t k; |
| 16 | + while(n > 0) { |
| 17 | + k = 1; |
| 18 | + while(G[k] <= n) { |
| 19 | + k++; |
| 20 | + } |
| 21 | + out[k - 1] = n / G[k - 1]; |
| 22 | + n = n % G[k - 1]; |
| 23 | + } |
| 24 | + return out; |
| 25 | +} |
| 26 | + |
| 27 | +arma::mat grid(const size_t d) { |
| 28 | + std::array<double, 3> x = {0.01, 0.5, 0.99}; |
| 29 | + size_t p = pow((size_t)3, d); |
| 30 | + arma::mat out(d, p); |
| 31 | + std::vector<size_t> threes(d, 3); |
| 32 | + for(size_t n = 0; n < p; n++) { |
| 33 | + std::vector<size_t> indices = CantorExpansion(n, threes); |
| 34 | + for(size_t i = 0; i < d; i++) { |
| 35 | + out.at(i, n) = x[indices[i]]; |
| 36 | + } |
| 37 | + } |
| 38 | + return out; |
| 39 | +} |
| 40 | + |
| 41 | +arma::vec logit(const arma::vec& u) { |
| 42 | + return arma::log(u / (1.0 - u)); |
| 43 | +} |
| 44 | + |
| 45 | +double dlogit(double u) { |
| 46 | + return 1.0 / (u * (1.0 - u)); |
| 47 | +} |
| 48 | + |
| 49 | +arma::vec ldlogit(const arma::vec& u) { |
| 50 | + return -arma::log(u % (1.0-u)); |
| 51 | +} |
| 52 | + |
| 53 | +arma::vec ldlogis(const arma::vec& x) { |
| 54 | + return x - 2.0 * arma::log1p(arma::exp(x)); |
| 55 | +} |
| 56 | + |
| 57 | +arma::vec dldlogis(const arma::vec& x) { |
| 58 | + return 1.0 - 2.0 / (1.0 + arma::exp(-x)); |
| 59 | +} |
| 60 | + |
| 61 | +double log_f(const arma::vec& u, const arma::mat& P, const arma::vec& b) { |
| 62 | + const arma::vec x = P * logit(u) + b; |
| 63 | + return arma::sum(ldlogis(x)) + arma::sum(ldlogit(u)); |
| 64 | +} |
| 65 | + |
| 66 | +double dlog_f(const double ui, const arma::vec& Pi, const arma::vec& y) { |
| 67 | + return dlogit(ui) * arma::sum(Pi % y) + (2.0 * ui - 1.0) / (ui * (1.0 - ui)); |
| 68 | +} |
| 69 | + |
| 70 | +class Logf : public Functor { |
| 71 | + public: |
| 72 | + arma::mat P; |
| 73 | + arma::vec b; |
| 74 | + double operator()(const arma::vec& u) override { return log_f(u, P, b); } |
| 75 | + void Gradient(const arma::vec& u, arma::vec& gr) override { |
| 76 | + const size_t d = P.n_cols; |
| 77 | + gr = arma::zeros<arma::vec>(d); |
| 78 | + const arma::vec y = dldlogis(P * logit(u) + b); |
| 79 | + for(size_t i = 0; i < d; i++) { |
| 80 | + gr(i) = dlog_f(u[i], P.col(i), y); |
| 81 | + } |
| 82 | + } |
| 83 | +}; |
| 84 | + |
| 85 | +class uLogf1 : public Functor { |
| 86 | +public: |
| 87 | + arma::mat P; |
| 88 | + arma::vec b; |
| 89 | + arma::vec mu; |
| 90 | + size_t j; |
| 91 | + double operator()(const arma::vec& u) override { |
| 92 | + const size_t d = P.n_cols; |
| 93 | + return -log_f(u, P, b) - (d+2) * log(mu[j] - u.at(j)); |
| 94 | + } |
| 95 | + void Gradient(const arma::vec& u, arma::vec& gr) override { |
| 96 | + const size_t d = P.n_cols; |
| 97 | + gr = arma::zeros<arma::vec>(d); |
| 98 | + const arma::vec y = dldlogis(P * logit(u) + b); |
| 99 | + for(size_t i = 0; i < d; i++) { |
| 100 | + if(i == j){ |
| 101 | + gr(i) = -dlog_f(u[i], P.col(i), y) + (d+2) / (mu[i] - u.at(i)); |
| 102 | + }else{ |
| 103 | + gr(i) = -dlog_f(u[i], P.col(i), y); |
| 104 | + } |
| 105 | + } |
| 106 | + } |
| 107 | +}; |
| 108 | + |
| 109 | +class uLogf2 : public Functor { |
| 110 | +public: |
| 111 | + arma::mat P; |
| 112 | + arma::vec b; |
| 113 | + arma::vec mu; |
| 114 | + size_t j; |
| 115 | + double operator()(const arma::vec& u) override { |
| 116 | + const size_t d = P.n_cols; |
| 117 | + return log_f(u, P, b) + (d+2) * log(u.at(j) - mu[j]); |
| 118 | + } |
| 119 | + void Gradient(const arma::vec& u, arma::vec& gr) override { |
| 120 | + const size_t d = P.n_cols; |
| 121 | + gr = arma::zeros<arma::vec>(d); |
| 122 | + const arma::vec y = dldlogis(P * logit(u) + b); |
| 123 | + for(size_t i = 0; i < d; i++) { |
| 124 | + if(i == j){ |
| 125 | + gr(i) = dlog_f(u[i], P.col(i), y) - (d+2) / (mu[i] - u.at(i)); |
| 126 | + }else{ |
| 127 | + gr(i) = dlog_f(u[i], P.col(i), y); |
| 128 | + } |
| 129 | + } |
| 130 | + } |
| 131 | +}; |
| 132 | + |
| 133 | +Rcpp::List get_umax0(const arma::mat& P, const arma::vec& b, arma::vec init) { |
| 134 | + double eps = sqrt(std::numeric_limits<double>::epsilon()); |
| 135 | + Logf logf; |
| 136 | + logf.P = P; |
| 137 | + logf.b = b; |
| 138 | + Roptim<Logf> opt("L-BFGS-B"); |
| 139 | + opt.control.trace = 0; |
| 140 | + opt.control.maxit = 10000; |
| 141 | + opt.control.fnscale = -1.0; // maximize |
| 142 | + //opt.control.factr = 1.0; |
| 143 | + opt.set_hessian(false); |
| 144 | + arma::vec lwr = arma::zeros(init.size()) + eps; |
| 145 | + arma::vec upr = arma::ones(init.size()) - eps; |
| 146 | + opt.set_lower(lwr); |
| 147 | + opt.set_upper(upr); |
| 148 | + opt.minimize(logf, init); |
| 149 | + if(opt.convergence() != 0){ |
| 150 | + Rcpp::Rcout << "-- umax -----------------------" << std::endl; |
| 151 | + opt.print(); |
| 152 | + } |
| 153 | + //Rcpp::Rcout << "-------------------------" << std::endl; |
| 154 | + // opt.print(); |
| 155 | + return Rcpp::List::create(Rcpp::Named("par") = opt.par(), |
| 156 | + Rcpp::Named("value") = opt.value()); |
| 157 | +} |
| 158 | + |
| 159 | +// [[Rcpp::export]] |
| 160 | +Rcpp::List get_umax(const arma::mat& P, const arma::vec& b) { |
| 161 | + const size_t d = P.n_cols; |
| 162 | + const arma::mat inits = grid(d); |
| 163 | + const size_t n = inits.n_cols; |
| 164 | + std::vector<arma::vec> pars(n); |
| 165 | + arma::vec values(n); |
| 166 | + for(size_t i = 0; i < n; i++) { |
| 167 | + const Rcpp::List L = get_umax0(P, b, inits.col(i)); |
| 168 | + const arma::vec par = L["par"]; |
| 169 | + pars[i] = par; |
| 170 | + // double value = L["value"]; |
| 171 | + values(i) = L["value"]; |
| 172 | + } |
| 173 | + const size_t imax = values.index_max(); |
| 174 | + return Rcpp::List::create( |
| 175 | + Rcpp::Named("mu") = pars[imax], |
| 176 | + Rcpp::Named("umax") = pow(exp(values(imax)), 2.0 / (2.0 + d))); |
| 177 | +} |
| 178 | + |
| 179 | +// [[Rcpp::export]] |
| 180 | +double get_vmin_i( |
| 181 | + const arma::mat& P, const arma::vec& b, const size_t i, const arma::vec& mu |
| 182 | +) { |
| 183 | + double eps = sqrt(std::numeric_limits<double>::epsilon()) / 3.0; |
| 184 | + uLogf1 ulogf1; |
| 185 | + ulogf1.P = P; |
| 186 | + ulogf1.b = b; |
| 187 | + ulogf1.j = i; |
| 188 | + ulogf1.mu = mu; |
| 189 | + Roptim<uLogf1> opt("L-BFGS-B"); |
| 190 | + opt.control.trace = 0; |
| 191 | + opt.control.maxit = 10000; |
| 192 | + //opt.control.fnscale = 1.0; // minimize |
| 193 | + //opt.control.factr = 1.0; |
| 194 | + opt.set_hessian(false); |
| 195 | + const size_t d = P.n_cols; |
| 196 | + arma::vec init = 0.5 * arma::ones(d); |
| 197 | + init.at(i) = mu.at(i) / 2.0; |
| 198 | + arma::vec lwr = arma::zeros(d) + eps; |
| 199 | + arma::vec upr = arma::ones(d); |
| 200 | + upr.at(i) = mu.at(i); |
| 201 | + opt.set_lower(lwr); |
| 202 | + opt.set_upper(upr - eps); |
| 203 | + opt.minimize(ulogf1, init); |
| 204 | + if(opt.convergence() != 0){ |
| 205 | + Rcpp::Rcout << "-- vmin -----------------------" << std::endl; |
| 206 | + opt.print(); |
| 207 | + } |
| 208 | + //Rcpp::Rcout << "-------------------------" << std::endl; |
| 209 | + return -exp(-opt.value() / (d+2)); |
| 210 | +} |
| 211 | + |
| 212 | +// [[Rcpp::export]] |
| 213 | +arma::vec get_vmin( |
| 214 | + const arma::mat& P, const arma::vec& b, const arma::vec& mu |
| 215 | +) { |
| 216 | + const size_t d = P.n_cols; |
| 217 | + arma::vec vmin(d); |
| 218 | + for(size_t i = 0; i < d; i++){ |
| 219 | + vmin.at(i) = get_vmin_i(P, b, i, mu); |
| 220 | + } |
| 221 | + return vmin; |
| 222 | +} |
| 223 | + |
| 224 | +double get_vmax_i( |
| 225 | + const arma::mat& P, const arma::vec& b, const size_t i, const arma::vec& mu |
| 226 | +) { |
| 227 | + double eps = sqrt(std::numeric_limits<double>::epsilon()) / 3.0; |
| 228 | + uLogf2 ulogf2; |
| 229 | + ulogf2.P = P; |
| 230 | + ulogf2.b = b; |
| 231 | + ulogf2.j = i; |
| 232 | + ulogf2.mu = mu; |
| 233 | + Roptim<uLogf2> opt("L-BFGS-B"); |
| 234 | + opt.control.trace = 0; |
| 235 | + opt.control.maxit = 10000; |
| 236 | + opt.control.fnscale = -1.0; // maximize |
| 237 | + //opt.control.factr = 1.0; |
| 238 | + opt.set_hessian(false); |
| 239 | + const size_t d = P.n_cols; |
| 240 | + arma::vec init = 0.5 * arma::ones(d); |
| 241 | + init.at(i) = (mu.at(i) + 1.0) / 2.0; |
| 242 | + arma::vec lwr = arma::zeros(d); |
| 243 | + lwr.at(i) = mu.at(i); |
| 244 | + arma::vec upr = arma::ones(d) - eps; |
| 245 | + opt.set_lower(lwr + eps); |
| 246 | + opt.set_upper(upr); |
| 247 | + opt.minimize(ulogf2, init); |
| 248 | + if(opt.convergence() != 0){ |
| 249 | + Rcpp::Rcout << "-- vmax -----------------------" << std::endl; |
| 250 | + opt.print(); |
| 251 | + } |
| 252 | + return exp(opt.value() / (d+2)); |
| 253 | +} |
| 254 | + |
| 255 | +arma::vec get_vmax( |
| 256 | + const arma::mat& P, const arma::vec& b, const arma::vec& mu |
| 257 | +) { |
| 258 | + const size_t d = P.n_cols; |
| 259 | + arma::vec vmax(d); |
| 260 | + for(size_t i = 0; i < d; i++){ |
| 261 | + vmax.at(i) = get_vmax_i(P, b, i, mu); |
| 262 | + } |
| 263 | + return vmax; |
| 264 | +} |
| 265 | + |
| 266 | +// [[Rcpp::export]] |
| 267 | +Rcpp::List get_bounds(const arma::mat& P, const arma::vec& b){ |
| 268 | + Rcpp::List L = get_umax(P, b); |
| 269 | + arma::vec mu = L["mu"]; |
| 270 | + double umax = L["umax"]; |
| 271 | + arma::vec vmin = get_vmin(P, b, mu); |
| 272 | + arma::vec vmax = get_vmax(P, b, mu); |
| 273 | + return Rcpp::List::create(Rcpp::Named("umax") = umax, |
| 274 | + Rcpp::Named("mu") = mu, |
| 275 | + Rcpp::Named("vmin") = vmin, |
| 276 | + Rcpp::Named("vmax") = vmax); |
| 277 | +} |
| 278 | + |
| 279 | + |
| 280 | +// std::uniform_real_distribution<double> runif(0.0, 1.0); |
| 281 | +// std::default_random_engine generator(seed); |
| 282 | +// runif(generator) |
| 283 | +std::default_random_engine generator; |
| 284 | +std::uniform_real_distribution<double> runif(0.0, 1.0); |
| 285 | + |
| 286 | +// [[Rcpp::export]] |
| 287 | +arma::mat rcd(const size_t n, const arma::mat& P, const arma::vec& b){ |
| 288 | + //, const size_t seed){ |
| 289 | +// std::default_random_engine generator(seed); |
| 290 | +// std::uniform_real_distribution<double> runif(0.0, 1.0); |
| 291 | + const size_t d = P.n_cols; |
| 292 | + arma::mat tout(d, n); |
| 293 | + const Rcpp::List bounds = get_bounds(P, b); |
| 294 | + const double umax = bounds["umax"]; |
| 295 | + const arma::vec mu = bounds["mu"]; |
| 296 | + const arma::vec vmin = bounds["vmin"]; |
| 297 | + const arma::vec vmax = bounds["vmax"]; |
| 298 | + size_t k = 0; |
| 299 | + while(k < n){ |
| 300 | + const double u = umax * runif(generator); |
| 301 | + arma::vec v(d); |
| 302 | + for(size_t i = 0; i < d; i++){ |
| 303 | + v.at(i) = vmin.at(i) + (vmax.at(i) - vmin.at(i)) * runif(generator); |
| 304 | + } |
| 305 | + const arma::vec x = v / sqrt(u) + mu; |
| 306 | + bool test = arma::all(x > 0.0) && arma::all(x < 1.0) && |
| 307 | + (d+2) * log(u) < 2.0 * log_f(x, P, b); |
| 308 | + if(test){ |
| 309 | + tout.col(k) = logit(x); |
| 310 | + k++; |
| 311 | + } |
| 312 | + } |
| 313 | + return tout.t(); |
| 314 | +} |
| 315 | + |
| 316 | + |
| 317 | +//////////////////////////////////////////////////////////////////////////////// |
| 318 | +double plogis(double x){ |
| 319 | + return 1.0/(1.0 + exp(-x)); |
| 320 | +} |
| 321 | + |
| 322 | +double qlogis(double u){ |
| 323 | + return log(u/(1.0-u)); |
| 324 | +} |
| 325 | + |
| 326 | +double MachineEps = std::numeric_limits<double>::epsilon(); |
| 327 | + |
| 328 | +double rlogis1(double x){ |
| 329 | + double b = plogis(x); |
| 330 | + if(b <= MachineEps){ |
| 331 | + return x; |
| 332 | + } |
| 333 | + std::uniform_real_distribution<double> ru(MachineEps, b); |
| 334 | + return qlogis(ru(generator)); |
| 335 | +} |
| 336 | + |
| 337 | +double rlogis2(double x){ |
| 338 | + double a = plogis(x); |
| 339 | + if(a == 1){ |
| 340 | + return x; |
| 341 | + } |
| 342 | + std::uniform_real_distribution<double> ru(a, 1); |
| 343 | + return qlogis(ru(generator)); |
| 344 | +} |
0 commit comments