return derivatives optionally with getProfileLikelihood

R/ModelFit.R

@@ -942,13 +942,22 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
     prof
 }
 
-.initAdaptiveProfile <- function(object, parm, bounds, includePenalty) {
+.initAdaptiveProfile <- function(object, parm, bounds, includePenalty, returnDerivatives) {
     # If an MLE was found, let's not throw that bit of important information away:
     if (object$return_flag == "SUCCESS" &&
         coef(object)[as.character(parm)] > bounds[1] &&
         coef(object)[as.character(parm)] < bounds[2]) {
-        profile <- tibble(point = coef(object)[as.character(parm)],
-                          value = fixedGridProfileLogLikelihood(object, parm, coef(object)[as.character(parm)], includePenalty)$value)
+        eval <- fixedGridProfileLogLikelihood(object, parm, coef(object)[as.character(parm)], includePenalty,
+                                              returnDerivatives)
+        point <- coef(object)[as.character(parm)]
+        if (returnDerivatives) {
+            profile <- tibble(point = point,
+                              value = eval$value,
+                              derivative = eval$derivative)
+        } else {
+            profile <- tibble(point = point,
+                              value = eval$value)
+        }
     } else {
         profile <- tibble()
     }
@@ -967,6 +976,7 @@ confint.cyclopsFit <- function(object, parm, level = 0.95, #control,
 #' @param tolerance Absolute tolerance allowed for adaptive profiling
 #' @param initialGridSize Initial grid size for adaptive profiling
 #' @param includePenalty    Logical: Include regularized covariate penalty in profile
+#' @param returnDerivatives Return derivative of log-likelihood wrt the parameter `param` at each evaluation point
 #'
 #' @return
 #' A data frame containing the profile log likelihood. Returns NULL when the adaptive profiling fails
@@ -979,7 +989,8 @@ getCyclopsProfileLogLikelihood <- function(object,
                                            bounds = NULL,
                                            tolerance = 1E-3,
                                            initialGridSize = 10,
-                                           includePenalty = TRUE) {
+                                           includePenalty = TRUE,
+                                           returnDerivatives = FALSE) {
     maxResets <- 10
     if (!xor(is.null(x), is.null(bounds))) {
         stop("Must provide either `x` or `bounds`, but not both.")
@@ -990,7 +1001,7 @@ getCyclopsProfileLogLikelihood <- function(object,
         if (length(bounds) != 2 || bounds[1] >= bounds[2]) {
             stop("Must provide bounds[1] < bounds[2]")
         }
-        profile <- .initAdaptiveProfile(object, parm, bounds, includePenalty)
+        profile <- .initAdaptiveProfile(object, parm, bounds, includePenalty, returnDerivatives)
 
         # Start with sparse grid:
         grid <- seq(bounds[1], bounds[2], length.out = initialGridSize)
@@ -999,7 +1010,7 @@ getCyclopsProfileLogLikelihood <- function(object,
         priorMaxMaxError <- Inf
         resetsPerformed <- 0
         while (length(grid) != 0) {
-            ll <- fixedGridProfileLogLikelihood(object, parm, grid, includePenalty)
+            ll <- fixedGridProfileLogLikelihood(object, parm, grid, includePenalty, returnDerivatives)
             profile <- bind_rows(profile, ll) %>% arrange(.data$point)
             invalid <- is.nan(profile$value) | is.infinite(profile$value)
             if (any(invalid)) {
@@ -1063,21 +1074,21 @@ getCyclopsProfileLogLikelihood <- function(object,
             grid <- (profile$point[exceed] + profile$point[exceed + 1]) / 2
         }
     } else { # Use x
-        profile <- fixedGridProfileLogLikelihood(object, parm, x, includePenalty)
+        profile <- fixedGridProfileLogLikelihood(object, parm, x, includePenalty, returnDerivatives)
     }
 
     return(profile)
 }
 
-fixedGridProfileLogLikelihood <- function(object, parm, x, includePenalty) {
+fixedGridProfileLogLikelihood <- function(object, parm, x, includePenalty, returnDerivatives) {
 
     .checkInterface(object$cyclopsData, testOnly = TRUE)
     parm <- .checkCovariates(object$cyclopsData, parm)
     threads <- object$threads
 
     if (getNumberOfCovariates(object$cyclopsData) == 1 || length(x) == 1) {
         grid <- .cyclopsGetProfileLikelihood(object$cyclopsData$cyclopsInterfacePtr, parm, x,
-                                             threads, includePenalty)
+                                             threads, includePenalty, returnDerivatives)
     } else {
         # Partition sequence
         y <- sort(x)
@@ -1087,9 +1098,9 @@ fixedGridProfileLogLikelihood <- function(object, parm, x, includePenalty) {
 
         # Execute: TODO chunk and repeat until ill-conditioned
         gridLower <- .cyclopsGetProfileLikelihood(object$cyclopsData$cyclopsInterfacePtr, parm, lower,
-                                                  threads, includePenalty)
+                                                  threads, includePenalty, returnDerivatives)
         gridUpper <- .cyclopsGetProfileLikelihood(object$cyclopsData$cyclopsInterfacePtr, parm, upper,
-                                                  threads, includePenalty)
+                                                  threads, includePenalty, returnDerivatives)
         # Merge
         grid <- rbind(gridLower, gridUpper)
         grid <- grid[order(grid$point),]

R/RcppExports.R

@@ -81,8 +81,8 @@
     invisible(.Call(`_Cyclops_cyclopsSetParameterizedPrior`, inRcppCcdInterface, priorTypeName, priorFunction, startingParameters, excludeNumeric))
 }
 
-.cyclopsGetProfileLikelihood <- function(inRcppCcdInterface, inCovariate, points, threads, includePenalty) {
-    .Call(`_Cyclops_cyclopsGetProfileLikelihood`, inRcppCcdInterface, inCovariate, points, threads, includePenalty)
+.cyclopsGetProfileLikelihood <- function(inRcppCcdInterface, inCovariate, points, threads, includePenalty, returnDerivatives) {
+    .Call(`_Cyclops_cyclopsGetProfileLikelihood`, inRcppCcdInterface, inCovariate, points, threads, includePenalty, returnDerivatives)
 }
 
 .cyclopsProfileModel <- function(inRcppCcdInterface, sexpCovariates, threads, threshold, override, includePenalty) {
@@ -117,8 +117,8 @@
     .Call(`_Cyclops_cyclopsRunBootstrap`, inRcppCcdInterface, outFileName, treatmentId, replicates)
 }
 
-.cyclopsGetLogLikelihoodGradient <- function(inRcppCcdInterface) {
-    .Call(`_Cyclops_cyclopsGetLogLikelihoodGradient`, inRcppCcdInterface)
+.cyclopsGetLogLikelihoodGradient <- function(inRcppCcdInterface, index) {
+    .Call(`_Cyclops_cyclopsGetLogLikelihoodGradient`, inRcppCcdInterface, index)
 }
 
 .cyclopsLogModel <- function(inRcppCcdInterface) {

src/RcppCyclopsInterface.cpp

@@ -392,7 +392,8 @@ void cyclopsSetParameterizedPrior(SEXP inRcppCcdInterface,
 DataFrame cyclopsGetProfileLikelihood(SEXP inRcppCcdInterface,
                                  SEXP inCovariate,
                                  const std::vector<double> points,
-                                 int threads, bool includePenalty) {
+                                 int threads, bool includePenalty,
+                                 bool returnDerivatives) {
     using namespace bsccs;
     XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
 
@@ -402,12 +403,27 @@ DataFrame cyclopsGetProfileLikelihood(SEXP inRcppCcdInterface,
     //const IdType covariate = as<IdType>(inCovariate);
 
     std::vector<double> values(points.size());
-    interface->evaluateProfileModel(covariate, points, values, threads, includePenalty);
 
-    return DataFrame::create(
-        Rcpp::Named("point") = points,
-        Rcpp::Named("value") = values
-    );
+    if (!returnDerivatives) {
+        interface->evaluateProfileModel(covariate, points, values, nullptr, threads, includePenalty);
+
+        return DataFrame::create(
+            Rcpp::Named("point") = points,
+            Rcpp::Named("value") = values
+        );
+
+    } else {
+
+        std::vector<double> derivatives(points.size());
+
+        interface->evaluateProfileModel(covariate, points, values, &derivatives, threads, includePenalty);
+
+        return DataFrame::create(
+            Rcpp::Named("point") = points,
+            Rcpp::Named("value") = values,
+            Rcpp::Named("derivative") = derivatives
+        );
+    }
 }
 
 // [[Rcpp::export(".cyclopsProfileModel")]]
@@ -597,27 +613,42 @@ List cyclopsRunBootstrap(SEXP inRcppCcdInterface, const std::string& outFileName
 }
 
 // [[Rcpp::export(".cyclopsGetLogLikelihoodGradient")]]
-NumericVector cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterface) {
+NumericVector cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterface, int index) {
     using namespace bsccs;
 
     XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
 
     auto& ccd = interface->getCcd();
     auto& data = interface->getModelData();
 
+    // const auto offset = data.getHasOffsetCovariate();
+    const int offset = 0;
 
-    const auto offset = data.getHasOffsetCovariate();
-    const auto length = ccd.getBetaSize() - offset;
-
-    NumericVector gradient(length);
-
-    for (int i = 0; i < length; ++i) {
-        gradient[i] = ccd.getLogLikelihoodGradient(i + offset);
-    }
-
-    return gradient;
+    return ccd.getLogLikelihoodGradient(index + offset);
 }
 
+// // [[Rcpp::export(".cyclopsGetLogLikelihoodGradient")]]
+// NumericVector cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterface) {
+//     using namespace bsccs;
+//
+//     XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+//
+//     auto& ccd = interface->getCcd();
+//     auto& data = interface->getModelData();
+//
+//
+//     const auto offset = data.getHasOffsetCovariate();
+//     const auto length = ccd.getBetaSize() - offset;
+//
+//     NumericVector gradient(length);
+//
+//     for (int i = 0; i < length; ++i) {
+//         gradient[i] = ccd.getLogLikelihoodGradient(i + offset);
+//     }
+//
+//     return gradient;
+// }
+
 // [[Rcpp::export(".cyclopsLogModel")]]
 List cyclopsLogModel(SEXP inRcppCcdInterface) {
 	using namespace bsccs;

src/RcppCyclopsInterface.h

@@ -55,8 +55,10 @@ class RcppCcdInterface : public CcdInterface {
     double evaluateProfileModel(const IdType covariate,
                                 const std::vector<double>& points,
                                 std::vector<double>& values,
+                                std::vector<double>* derivatives,
                                 int threads, bool includePenalty) {
-        return CcdInterface::evaluateProfileModel(ccd, modelData, covariate, points, values, threads, includePenalty);
+        return CcdInterface::evaluateProfileModel(ccd, modelData, covariate, points, values,
+                                                  derivatives, threads, includePenalty);
     }
 
     double runCrossValidation() {

src/RcppExports.cpp

@@ -245,8 +245,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // cyclopsGetProfileLikelihood
-DataFrame cyclopsGetProfileLikelihood(SEXP inRcppCcdInterface, SEXP inCovariate, const std::vector<double> points, int threads, bool includePenalty);
-RcppExport SEXP _Cyclops_cyclopsGetProfileLikelihood(SEXP inRcppCcdInterfaceSEXP, SEXP inCovariateSEXP, SEXP pointsSEXP, SEXP threadsSEXP, SEXP includePenaltySEXP) {
+DataFrame cyclopsGetProfileLikelihood(SEXP inRcppCcdInterface, SEXP inCovariate, const std::vector<double> points, int threads, bool includePenalty, bool returnDerivatives);
+RcppExport SEXP _Cyclops_cyclopsGetProfileLikelihood(SEXP inRcppCcdInterfaceSEXP, SEXP inCovariateSEXP, SEXP pointsSEXP, SEXP threadsSEXP, SEXP includePenaltySEXP, SEXP returnDerivativesSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -255,7 +255,8 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const std::vector<double> >::type points(pointsSEXP);
     Rcpp::traits::input_parameter< int >::type threads(threadsSEXP);
     Rcpp::traits::input_parameter< bool >::type includePenalty(includePenaltySEXP);
-    rcpp_result_gen = Rcpp::wrap(cyclopsGetProfileLikelihood(inRcppCcdInterface, inCovariate, points, threads, includePenalty));
+    Rcpp::traits::input_parameter< bool >::type returnDerivatives(returnDerivativesSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsGetProfileLikelihood(inRcppCcdInterface, inCovariate, points, threads, includePenalty, returnDerivatives));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -375,13 +376,14 @@ BEGIN_RCPP
 END_RCPP
 }
 // cyclopsGetLogLikelihoodGradient
-NumericVector cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterface);
-RcppExport SEXP _Cyclops_cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterfaceSEXP) {
+NumericVector cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterface, int index);
+RcppExport SEXP _Cyclops_cyclopsGetLogLikelihoodGradient(SEXP inRcppCcdInterfaceSEXP, SEXP indexSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
-    rcpp_result_gen = Rcpp::wrap(cyclopsGetLogLikelihoodGradient(inRcppCcdInterface));
+    Rcpp::traits::input_parameter< int >::type index(indexSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsGetLogLikelihoodGradient(inRcppCcdInterface, index));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -897,7 +899,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsSetPrior", (DL_FUNC) &_Cyclops_cyclopsSetPrior, 6},
     {"_Cyclops_cyclopsTestParameterizedPrior", (DL_FUNC) &_Cyclops_cyclopsTestParameterizedPrior, 4},
     {"_Cyclops_cyclopsSetParameterizedPrior", (DL_FUNC) &_Cyclops_cyclopsSetParameterizedPrior, 5},
-    {"_Cyclops_cyclopsGetProfileLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetProfileLikelihood, 5},
+    {"_Cyclops_cyclopsGetProfileLikelihood", (DL_FUNC) &_Cyclops_cyclopsGetProfileLikelihood, 6},
     {"_Cyclops_cyclopsProfileModel", (DL_FUNC) &_Cyclops_cyclopsProfileModel, 6},
     {"_Cyclops_cyclopsPredictModel", (DL_FUNC) &_Cyclops_cyclopsPredictModel, 1},
     {"_Cyclops_cyclopsSetControl", (DL_FUNC) &_Cyclops_cyclopsSetControl, 23},
@@ -906,7 +908,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsClearCacheForJava", (DL_FUNC) &_Cyclops_cyclopsClearCacheForJava, 0},
     {"_Cyclops_cyclopsFitModel", (DL_FUNC) &_Cyclops_cyclopsFitModel, 1},
     {"_Cyclops_cyclopsRunBootstrap", (DL_FUNC) &_Cyclops_cyclopsRunBootstrap, 4},
-    {"_Cyclops_cyclopsGetLogLikelihoodGradient", (DL_FUNC) &_Cyclops_cyclopsGetLogLikelihoodGradient, 1},
+    {"_Cyclops_cyclopsGetLogLikelihoodGradient", (DL_FUNC) &_Cyclops_cyclopsGetLogLikelihoodGradient, 2},
     {"_Cyclops_cyclopsLogModel", (DL_FUNC) &_Cyclops_cyclopsLogModel, 1},
     {"_Cyclops_cyclopsInitializeModel", (DL_FUNC) &_Cyclops_cyclopsInitializeModel, 4},
     {"_Cyclops_listGPUDevices", (DL_FUNC) &_Cyclops_listGPUDevices, 0},

src/cyclops/CcdInterface.cpp

@@ -250,6 +250,14 @@ struct OptimizationProfile {
 		return y;
 	}
 
+	double derivative() {
+	    double d = ccd.getLogLikelihoodGradient(index);
+	    if (includePenalty) {
+	        d += ccd.getLogPriorGradient(index);
+	    }
+	    return d;
+	}
+
 	double getMaximum() {
 		return threshold;
 	}
@@ -533,6 +541,7 @@ double CcdInterface::evaluateProfileModel(CyclicCoordinateDescent *ccd, Abstract
                                           const IdType covariate,
                                           const std::vector<double>& points,
                                           std::vector<double>& values,
+                                          std::vector<double>* derivatives,
                                           int inThreads,
                                           bool includePenalty) {
 
@@ -563,13 +572,18 @@ double CcdInterface::evaluateProfileModel(CyclicCoordinateDescent *ccd, Abstract
         x0s[j] = ccd->getBeta(j);
     }
 
-    auto evaluate = [this, index, includePenalty](const double point,
+    bool doDerivative = derivatives != nullptr;
+
+    auto evaluate = [this, index, includePenalty, doDerivative](const double point,
                                                   CyclicCoordinateDescent* ccd) {
 
         OptimizationProfile eval(*ccd, arguments, index,
                                  0.0, 0.0, includePenalty);
 
-        return eval.objective(point);
+        double objective = eval.objective(point);
+        double derivative = doDerivative ? eval.derivative() : 0.0;
+
+        return std::pair<double,double>(objective, derivative);
     };
 
     ccd->makeDirty(); // Reset internal memory
@@ -584,16 +598,24 @@ double CcdInterface::evaluateProfileModel(CyclicCoordinateDescent *ccd, Abstract
 
     if (nThreads == 1) {
         for (int i = 0; i < points.size(); ++i) {
-            values[i] = evaluate(points[i], ccd);
+            auto pair = evaluate(points[i], ccd);
+            values[i] = pair.first;
+            if (doDerivative) {
+                (*derivatives)[i] = pair.second;
+            }
         }
     } else {
         auto scheduler = TaskScheduler<IncrementableIterator<size_t>>(
             IncrementableIterator<size_t>(0), //boost::make_counting_iterator(0),
             IncrementableIterator<size_t>(points.size()), //boost::make_counting_iterator(static_cast<int>(points.size())),
             nThreads);
 
-        auto oneTask = [&evaluate, &scheduler, &ccdPool, &points, &values](unsigned long task) {
-            values[task] = evaluate(points[task], ccdPool[scheduler.getThreadIndex(task)]);
+        auto oneTask = [&evaluate, &scheduler, &ccdPool, &points, &values, &derivatives, doDerivative](unsigned long task) {
+            auto pair = evaluate(points[task], ccdPool[scheduler.getThreadIndex(task)]);
+            values[task] = pair.first;
+            if (doDerivative) {
+                (*derivatives)[task] = pair.second;
+            }
         };
 
         // Run all tasks in parallel

src/cyclops/CcdInterface.h

@@ -188,6 +188,7 @@ class CcdInterface {
             const IdType covariate,
             const std::vector<double>& points,
             std::vector<double>& values,
+            std::vector<double>* derivatives,
             int threads,
             bool includePenalty);