add code for bootstrap

NAMESPACE

@@ -42,6 +42,7 @@ export(listOpenCLDevices)
 export(meanLinearPredictor)
 export(mse)
 export(readCyclopsData)
+export(runBootstrap)
 export(setOpenCLDevice)
 export(simulateCyclopsData)
 export(splitTime)

R/ModelFit.R

@@ -784,6 +784,21 @@ getSEs <- function(object, covariates) {
     ses
 }
 
+
+#' @title Run Bootstrap for Cyclops model parameter
+#'
+#' @param object    A fitted Cyclops model object
+#' @param outFileName     Character: Output file name
+#' @param treatmentId     Character: variable to output
+#' @param replicates      Numeric: number of bootstrap samples
+#'
+#' @export
+runBootstrap <- function(object, outFileName, treatmentId, replicates) {
+    .checkInterface(object$cyclopsData, testOnly = TRUE)
+    bs <- .cyclopsRunBootstrap(object$cyclopsData$cyclopsInterfacePtr, outFileName, treatmentId, replicates)
+    bs
+}
+
 #' @title Confidence intervals for Cyclops model parameters
 #'
 #' @description

R/RcppExports.R

@@ -93,6 +93,10 @@
     .Call(`_Cyclops_cyclopsFitModel`, inRcppCcdInterface)
 }
 
+.cyclopsRunBootstrap <- function(inRcppCcdInterface, outFileName, treatmentId, replicates) {
+    .Call(`_Cyclops_cyclopsRunBootstrap`, inRcppCcdInterface, outFileName, treatmentId, replicates)
+}
+
 .cyclopsLogModel <- function(inRcppCcdInterface) {
     .Call(`_Cyclops_cyclopsLogModel`, inRcppCcdInterface)
 }

src/RcppCyclopsInterface.cpp

@@ -518,6 +518,32 @@ List cyclopsFitModel(SEXP inRcppCcdInterface) {
 	return list;
 }
 
+// [[Rcpp::export(".cyclopsRunBootstrap")]]
+List cyclopsRunBootstrap(SEXP inRcppCcdInterface, const std::string& outFileName, std::string& treatmentId, int replicates) {
+    using namespace bsccs;
+
+    XPtr<RcppCcdInterface> interface(inRcppCcdInterface);
+    interface->getArguments().doBootstrap = true;
+    interface->getArguments().outFileName = outFileName;
+    interface->getArguments().replicates = replicates;
+
+    // Save parameter point-estimates
+    std::vector<double> savedBeta;
+    for (int j = 0; j < interface->getCcd().getBetaSize(); ++j) {
+        savedBeta.push_back(interface->getCcd().getBeta(j));
+    } // TODO Handle above work in interface.runBootstrap
+    double timeUpdate = interface->runBoostrap(savedBeta, treatmentId);
+
+    interface->diagnoseModel(0.0, 0.0);
+
+    List list = List::create(
+        Rcpp::Named("interface")=interface,
+        Rcpp::Named("timeFit")=timeUpdate
+    );
+    RcppCcdInterface::appendRList(list, interface->getResult());
+    return list;
+}
+
 // [[Rcpp::export(".cyclopsLogModel")]]
 List cyclopsLogModel(SEXP inRcppCcdInterface) {
 	using namespace bsccs;

src/RcppCyclopsInterface.h

@@ -63,8 +63,8 @@ class RcppCcdInterface : public CcdInterface {
     	return CcdInterface::runCrossValidation(ccd, modelData);
     }
 
-    double runBoostrap(std::vector<double>& savedBeta) {
-    	return CcdInterface::runBoostrap(ccd, modelData, savedBeta);
+    double runBoostrap(std::vector<double>& savedBeta, std::string& treatmentId) {
+    	return CcdInterface::runBoostrap(ccd, modelData, savedBeta, treatmentId);
     }
 
     void logResultsToFile(const std::string& fileName, bool withASE);

src/RcppExports.cpp

@@ -305,6 +305,20 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cyclopsRunBootstrap
+List cyclopsRunBootstrap(SEXP inRcppCcdInterface, const std::string& outFileName, std::string& treatmentId, int replicates);
+RcppExport SEXP _Cyclops_cyclopsRunBootstrap(SEXP inRcppCcdInterfaceSEXP, SEXP outFileNameSEXP, SEXP treatmentIdSEXP, SEXP replicatesSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< SEXP >::type inRcppCcdInterface(inRcppCcdInterfaceSEXP);
+    Rcpp::traits::input_parameter< const std::string& >::type outFileName(outFileNameSEXP);
+    Rcpp::traits::input_parameter< std::string& >::type treatmentId(treatmentIdSEXP);
+    Rcpp::traits::input_parameter< int >::type replicates(replicatesSEXP);
+    rcpp_result_gen = Rcpp::wrap(cyclopsRunBootstrap(inRcppCcdInterface, outFileName, treatmentId, replicates));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cyclopsLogModel
 List cyclopsLogModel(SEXP inRcppCcdInterface);
 RcppExport SEXP _Cyclops_cyclopsLogModel(SEXP inRcppCcdInterfaceSEXP) {
@@ -820,6 +834,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_Cyclops_cyclopsSetControl", (DL_FUNC) &_Cyclops_cyclopsSetControl, 23},
     {"_Cyclops_cyclopsRunCrossValidationl", (DL_FUNC) &_Cyclops_cyclopsRunCrossValidationl, 1},
     {"_Cyclops_cyclopsFitModel", (DL_FUNC) &_Cyclops_cyclopsFitModel, 1},
+    {"_Cyclops_cyclopsRunBootstrap", (DL_FUNC) &_Cyclops_cyclopsRunBootstrap, 4},
     {"_Cyclops_cyclopsLogModel", (DL_FUNC) &_Cyclops_cyclopsLogModel, 1},
     {"_Cyclops_cyclopsInitializeModel", (DL_FUNC) &_Cyclops_cyclopsInitializeModel, 4},
     {"_Cyclops_listOpenCLDevices", (DL_FUNC) &_Cyclops_listOpenCLDevices, 0},

src/cyclops/CcdInterface.cpp

@@ -697,21 +697,29 @@ SelectorType CcdInterface::getDefaultSelectorTypeOrOverride(SelectorType selecto
 double CcdInterface::runBoostrap(
 		CyclicCoordinateDescent *ccd,
 		AbstractModelData *modelData,
-		std::vector<double>& savedBeta) {
+		std::vector<double>& savedBeta,
+		std::string& treatmentId) {
 	struct timeval time1, time2;
 	gettimeofday(&time1, NULL);
 
 	auto selectorType = getDefaultSelectorTypeOrOverride(
 		arguments.crossValidation.selectorType, modelData->getModelType());
 
-	BootstrapSelector selector(arguments.replicates, modelData->getPidVectorSTL(),
+	vector<int> ids;
+	if (selectorType == SelectorType::BY_ROW) {
+		std::cout << "runBoostrap SelectorType::BY_ROW \n";
+		ids.resize(modelData->getNumberOfRows());
+		std::iota(ids.begin(), ids.end(), 0);
+	}
+	BootstrapSelector selector(arguments.replicates, selectorType == SelectorType::BY_ROW ? ids : modelData->getPidVectorSTL(),
 			selectorType, arguments.seed, logger, error);
 	BootstrapDriver driver(arguments.replicates, modelData, logger, error);
 
 	driver.drive(*ccd, selector, arguments);
 	gettimeofday(&time2, NULL);
 
-	driver.logResults(arguments, savedBeta, ccd->getConditionId());
+//	driver.logResults(arguments, savedBeta, ccd->getConditionId());
+	driver.logHR(arguments, savedBeta, treatmentId);
 	return calculateSeconds(time1, time2);
 }
 

src/cyclops/CcdInterface.h

@@ -197,7 +197,8 @@ class CcdInterface {
     double runBoostrap(
             CyclicCoordinateDescent *ccd,
             AbstractModelData *modelData,
-            std::vector<double>& savedBeta);
+            std::vector<double>& savedBeta,
+	    std::string& treatmentId);
 
     void setDefaultArguments();
 

src/cyclops/drivers/BootstrapDriver.cpp

@@ -131,4 +131,60 @@ void BootstrapDriver::logResults(const CCDArguments& arguments, std::vector<doub
 	outLog.close();
 }
 
+void BootstrapDriver::logHR(const CCDArguments& arguments, std::vector<double>& savedBeta, std::string treatmentId) {
+
+//	int j = J-1;
+	int tId = 0;
+	while (modelData->getColumnLabel(tId) != treatmentId) tId++;
+
+	ofstream outLog(arguments.outFileName.c_str());
+	if (!outLog) {
+        std::ostringstream stream;
+		stream << "Unable to open log file: " << arguments.bsFileName;
+		error->throwError(stream);
+	}
+
+	string sep(","); // TODO Make option
+
+	// Raw estimates
+//	for(rvector::iterator it = estimates[j]->begin(); it != estimates[j]->end(); ++it) outLog << *it << endl;
+
+	// Stats
+	outLog << "Drug_concept_id" << sep << "Condition_concept_id" << sep <<
+		"score" << sep << "standard_error" << sep << "bs_mean" << sep << "bs_lower" << sep <<
+		"bs_upper" << sep << "bs_prob0" << endl;
+
+	for (int j = tId; j < J; ++j) {
+	outLog << modelData->getColumnLabel(j) <<
+		sep << treatmentId << sep;
+
+	double mean = 0.0;
+	double var = 0.0;
+	double prob0 = 0.0;
+	for (rvector::iterator it = estimates[j]->begin(); it != estimates[j]->end(); ++it) {
+		mean += *it;
+		var += *it * *it;
+		if (*it == 0.0) {
+			prob0 += 1.0;
+		}
+	}
+
+	double size = static_cast<double>(estimates[j]->size());
+	mean /= size;
+	var = (var / size) - (mean * mean);
+	prob0 /= size;
+
+	sort(estimates[j]->begin(), estimates[j]->end());
+	int offsetLower = static_cast<int>(size * 0.025);
+	int offsetUpper = static_cast<int>(size * 0.975);
+
+	double lower = *(estimates[j]->begin() + offsetLower);
+	double upper = *(estimates[j]->begin() + offsetUpper);
+
+	outLog << savedBeta[j] << sep;
+	outLog << std::sqrt(var) << sep << mean << sep << lower << sep << upper << sep << prob0 << endl;
+	}
+	outLog.close();
+}
+
 } // namespace

src/cyclops/drivers/BootstrapDriver.h

@@ -39,6 +39,8 @@ class BootstrapDriver : public AbstractDriver {
 
 	void logResults(const CCDArguments& arguments, std::vector<double>& savedBeta, std::string conditionId);
 
+void logHR(const CCDArguments& arguments, std::vector<double>& savedBeta, std::string treatmentId);
+
 private:
 	const int replicates;
 	AbstractModelData* modelData;

src/cyclops/drivers/BootstrapSelector.cpp

@@ -56,18 +56,20 @@ void BootstrapSelector::permute() {
 
 	// Get non-excluded indices
 	int N_new = indicesIncluded.size();
-	if (type == SelectorType::BY_PID) {
+//	if (type == SelectorType::BY_PID) {
 	    std::uniform_int_distribution<int> uniform(0, N_new - 1);
 		for (int i = 0; i < N_new; i++) {
             int ind =  uniform(prng);
 			int draw = indicesIncluded[ind];
 			selectedSet.insert(draw);
 		}
+/*
 	} else {
         std::ostringstream stream;
         stream << "BootstrapSelector::permute is not yet implemented.";
         error->throwError(stream);
 	}
+*/
 }
 
 void BootstrapSelector::getWeights(int batch, std::vector<double>& weights) {
@@ -80,16 +82,18 @@ void BootstrapSelector::getWeights(int batch, std::vector<double>& weights) {
 		return;
 	}
 
-	if (type == SelectorType::BY_PID) {
+//	if (type == SelectorType::BY_PID) {
 		for (size_t k = 0; k < K; k++) {
 			int count = selectedSet.count(ids.at(k));
 			weights[k] = static_cast<real>(count);
 		}
+/*
 	} else {
         std::ostringstream stream;
         stream << "BootstrapSelector::getWeights is not yet implemented.";
         error->throwError(stream);
 	}
+*/
 }
 
 void BootstrapSelector::getComplement(std::vector<double>& weights) {