The vignette looks like a manuscript.
Was it published somewhere that I could cite?

Hi Yi,

I am getting my package ready for CRAN and get the following note when running the check on R-devel:

File ‘psychonetrics/libs/psychonetrics.so’:
  Found ‘_ZSt4cout’, possibly from ‘std::cout’ (C++)
    Object: ‘04_generalfit_psychonetrics_BFGS.o’

Compiled code should not call entry points which might terminate R nor
write to stdout/stderr instead of to the console, nor use Fortran I/O
nor system RNGs.

I don't see any reference to cout or so in the corresponding code though. I was wondering if this NOTE is perhaps due to the roptim import here?

Hi, I want to pass other arguments to the roptim in Rcpp but I don't know how to do that.
For example, I want to do the same thing in rcpp as the following in R's optim():

optim(par=r0, fn=fn, gr=gr, X=X, time=time, q=q, k=k, method = 'L-BFGS-B')

Do you know how could I do that? Thank you very much!

Hi,

I have a custom function for which there is no analytical gradient that I can find. Is it still possible to use this package to do optimization in Rcpp (without defining a gradient function)?

Thanks,
Teng

Hi, thank you for making this package which might potentially facilitate one of the projects that I am working on. There is one thing that I couldn't follow, which is to reproduce the example described in the README.md. I downloaded the demo.cpp and try to run the following R code. However, it gave me lots of error messages, which do not fit here.

library(Rcpp)
sourceCpp("~/demo.cpp") # you may need to change the directory
rosen_bfgs()

Please let me know if this is the right way to execute the rcpp code. Thank you!

Hi Yi Pan,

Thank you very much for your earlier help! I was encountering one more problem I hope you could help me with. In optimization, I require parameters to form positive definite matrices. However, sometimes the optimizer would go to a parameter value that would lead to non positive definite matrices. I was wondering if there is some way to specify when the parameters are in a non tolerable space? I tried building in a simple if statement in the fit function that returns NA_REAL or a very large number, but the first solution gives the message:

 L-BFGS-B needs finite values of 'fn'

And the second solution seems to work but leads to poor parameter estimates.

Best, Sacha

I think your package could be very handy for my problem; thanks for developing it.

I prototyped some R code, which works fine:

library(bigsnpr)
bigsnp <- snp_attachExtdata()
G <- bigsnp$genotypes

simu <- snp_simuPheno(G, 0.2, 500, alpha = 0.5)

log_var <- log(big_colstats(G, ind.col = simu$set)$var)
beta2 <- simu$effects^2

FUN <- function(x, log_var, beta2) {
  S <- 1 + x[[1]]; sigma2 <- x[[2]]
  S * sum(log_var) + length(log_var) * log(sigma2) + sum(beta2 / exp(S * log_var)) / sigma2
}

DER <- function(x, log_var, beta2) {
  S <- 1 + x[[1]]; sigma2 <- x[[2]]
  res1 <- sum(log_var) - sum(log_var * beta2 / exp(S * log_var)) / sigma2
  res2 <- length(log_var) / sigma2 - sum(beta2 / exp(S * log_var)) / sigma2^2
  c(res1, res2)
}

optim(par = c(-0.5, 0.2 / 500), fn = FUN, method = "L-BFGS-B", gr = DER,
      lower = c(-2, 0.2 / 5000), upper = c(1, 0.2 / 50),
      log_var = log_var, beta2 = beta2)

Then I tried implementing it in C++ using your package:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

#include <roptim.h>
// [[Rcpp::depends(roptim)]]

using namespace roptim;

class MLE : public Functor {
public:
  MLE(const arma::vec& a_, const arma::vec& b_) : a(a_), b(b_) {}

  double operator()(const arma::vec& par) override {

    double S = par[0] + 1;
    double sigma2 = par[1];

    int m = a.size();
    double sum_a = 0, sum_c = 0;
    for (int j = 0; j < m; j++) {
      sum_a += a[j];
      sum_c += b[j] / ::exp(S * a[j]);
    }

    return S * sum_a + m * ::log(sigma2) + sum_c / sigma2;
  }

  void Gradient(const arma::vec& par, arma::vec& gr) override {

    double S = par[0] + 1;
    double sigma2 = par[1];

    int m = a.size();
    double sum_a = 0, sum_c = 0, sum_ac = 0;
    for (int j = 0; j < m; j++) {
      double c_j = b[j] / ::exp(S * a[j]);
      sum_a += a[j];
      sum_c += c_j;
      sum_ac += a[j] * c_j;
    }

    gr[0] = sum_a - sum_ac / sigma2;
    gr[1] = (m - sum_c / sigma2) / sigma2;
  }

private:
  arma::vec a;
  arma::vec b;
};


// [[Rcpp::export]]
arma::vec test_MLE(const arma::vec& a_, const arma::vec& b_,
                   arma::vec par, const arma::vec& lower, const arma::vec& upper) {

  MLE mle(a_, b_);
  Roptim<MLE> opt("L-BFGS-B");
  opt.set_lower(lower);
  opt.set_upper(upper);
  opt.set_hessian(false);
  opt.control.maxit = 100;

  opt.minimize(mle, par);

  Rcpp::Rcout << "-------------------------" << std::endl;
  opt.print();

  return par;
}

Then in R:

Rcpp::sourceCpp('tmp-tests/proto-MLE-optim-C++.cpp')
x <- c(0.5, 0.2 / 500)
test_MLE(log_var, beta2, x,
         lower = c(-2, 0.2 / 5000), upper = c(1, 0.2 / 50))

but it crashes..

I have been trying to debug it, but I am not sure if it comes from a silly mistake on my part, or something I didn't understand about how to use your package.
Any idea?

Dear Yi Pan,

Thank you for the roptim package! It works very well and seems to be exactly what I was looking for. I am working on a developmental version of my R package psychonetrics which uses roptim in this c++ code. The package works well on Linux, but when trying to compile the package on Windows I get the following errors:

I am entirely at a loss what this error could be, which seems not to be related to the source code of my package at all. I tried commenting out lines 82 - 85 of the code that uses roptim, but this did not solve the problem.

Thank you very much for your time!