Converting Optim (and it's Objective Faction) to Rcpp

Question

I have a section of code (see below) drawn from a larger model I am working on. This seems to be the bottle neck and I can see optim has been an issue for many others in the past. I know nothing about C++ and I can see this has been a recommended fix using the Rcpp package.

All my attempts to do this have failed so its time to bite the bullet and ask for some help and hopefully someone can help me. My optimisation code is as follows (with dummy data):

if(!require(TruncatedDistributions)){#_                                         
  install.packages("TruncatedDistributions", 
                   repos="http://R-Forge.R-project.org")
}#_____________________________________
library(TruncatedDistributions)

if(!require(Rcpp)){#___________________                                        
  install.packages("Rcpp")
}#_____________________________________
library(Rcpp)



# Inputs

Starting_Values    <- 0.15
labour_wage_elast  <- 0.45
wage               <- 0.15
employment         <- 152000
alpha              <- 1.5
sd                 <- 1.24
rowNum             <- 1
labour_value       <- 2531
max_workers        <- 352000
salary_lower_bound <- 0.01
salary_upper_bound <- 2


CLASSICAL_OPTIMISE_WAGE_EMPLOYEES <- function(Starting_Values,
                                              labour_wage_elast,
                                              wage,
                                              employment,
                                              alpha,
                                              sd,
                                              rowNum,
                                              labour_value,
                                              max_workers,
                                              salary_lower_bound,
                                              salary_upper_bound) {#~~~~~~~~~~~~    
  set.seed(rowNum)                                              
  
  # QL (Quantity of Labour)
  
  QL = ((Starting_Values[1] - wage) / wage * labour_wage_elast ) * employment + employment       
  
  QL = min(QL, max_workers)
  
  
  # Labour Value Function
  
  value = sum(rtgamma(seq(1:QL), 
                      shape = alpha, 
                      scale = (Starting_Values[1]  / sd ^ 2),
                      a = salary_lower_bound, 
                      b = salary_upper_bound))                                  
  
  difference  = abs(labour_value - value)                                       
  
  
  return(difference)
  
}#~~~~~~~~~ END of CLASSICAL_OPTIMISE_WAGE_EMPLOYEES  ~~~~~~~~~~~~~~~~~~~~~~~~~~




Estimates   <- optim(Starting_Values, CLASSICAL_OPTIMISE_WAGE_EMPLOYEES,
                     labour_wage_elast  = labour_wage_elast,
                     wage               = wage,
                     employment         = employment,
                     alpha              = alpha,
                     sd                 = sd, 
                     rowNum             = rowNum, 
                     labour_value       = labour_value, 
                     max_workers        = max_workers, 
                     salary_lower_bound = salary_lower_bound,
                     salary_upper_bound = salary_upper_bound,
                     method = "L-BFGS-B", lower = 0.001, 
                     upper = 2)

last_wage   <- wage 

wage        <- Estimates$par[1]

employment  =  ((wage - last_wage) / last_wage) * labour_wage_elast * employment + employment

Thank you very much in advance.

EDIT:

I have written the following:

Objective Function (which seems to work ok)

cppFunction('

//[[Rcpp::depends(Rcpp)]]
//[[Rcpp::depends(TruncatedDistributions)]]

  #include<cmath>
  #include<iostream>
  #include<Rcpp.h>
  #include<algorithm>
  #include<cstdlib>
  
  //[[Rcpp::export]]
  
 double classical_optimise_wage_employees(Rcpp::NumericVector Starting_Values,
                                         double labour_wage_elast,
                                         double wage,
                                         double employment,
                                         double alpha,
                                         double sd,
                                         int rowNum,
                                         double labour_value,
                                         double max_workers,
                                         double salary_lower_bound,
                                         double salary_upper_bound) {
                                         
  
  std::srand(rowNum);
  

  double QL = ((Starting_Values[0] - wage) / wage * labour_wage_elast ) * employment + employment;
  QL = std::min(QL, max_workers);
  
    
  Rcpp::Environment pkg = Rcpp::Environment::namespace_env("TruncatedDistributions");


  Rcpp::Function rtgamma_1 = Rcpp::as<Rcpp::Function>(pkg["rtgamma"]);
  

  NumericVector values = rtgamma_1(QL, alpha, Starting_Values[0] / pow(sd, 2) , salary_lower_bound, salary_upper_bound);
  
  double value = Rcpp::sum(values);
  

  return std::abs(labour_value - value);}')


# Calling it in R
classical_optimise_wage_employees(Starting_Values,
                                  labour_wage_elast,
                                  wage,
                                  employment,
                                  alpha,
                                  sd,
                                  rowNum,
                                  labour_value,
                                  max_workers,
                                  salary_lower_bound,
                                  salary_upper_bound)

And a function to use the optim function (which for some reason is not working)

cppFunction('

#include<Rcpp.h>
  
  //[[Rcpp::export]]
  
double solving_function(Rcpp::NumericVector Starting_Values,
                                         double labour_wage_elast,
                                         double wage,
                                         double employment,
                                         double alpha,
                                         double sd,
                                         int rowNum,
                                         double labour_value,
                                         double max_workers,
                                         double salary_lower_bound,
                                         double salary_upper_bound) {
                                         
  

            Rcpp::Environment stats("package:stats"); 
            Rcpp::Function optim = stats["optim"];
            
            
              Rcpp::List opt_results = optim(Rcpp::_["par"]    = Starting_Values,
                                             Rcpp::_["fn"]     = Rcpp::InternalFunction(&classical_optimise_wage_employees),
                                             Rcpp::_["method"] = "L-BFGS-B",
                                             Rcpp::_["labour_wage_elast"]  = labour_wage_elast,
                                             Rcpp::_["wage"]               = wage,
                                             Rcpp::_["employment"]         = employment,
                                             Rcpp::_["alpha"]              = alpha,
                                             Rcpp::_["sd"]                 = sd, 
                                             Rcpp::_["rowNum"]             = rowNum, 
                                             Rcpp::_["labour_value"]       = labour_value, 
                                             Rcpp::_["max_workers"]        = max_workers, 
                                             Rcpp::_["salary_lower_bound"] = salary_lower_bound,
                                             Rcpp::_["salary_upper_bound"] = salary_upper_bound,
                                             Rcpp::_["lower"] = 0.001,
                                             Rcpp::_["upper"] = 2);
            
            
            return opt_results[0];}')

I hope these edits help others to solve this problem. Thank you.