remove relaxed fit

DESCRIPTION

@@ -13,6 +13,6 @@ Encoding: UTF-8
 LazyData: true
 Imports: 
     Rcpp
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.0
 LinkingTo: 
     Rcpp, RcppEigen

NAMESPACE

@@ -4,7 +4,6 @@ export(basil)
 export(basil_base)
 export(compute_dual_norm)
 export(compute_residual)
-export(coxph_MKL)
 export(get_dual_norm)
 export(get_parameter_matrix)
 export(get_residual)

R/RcppExports.R

@@ -1,17 +1,22 @@
-# Generated by using Rcpp::compileAttributes() -> do not edit by hand
-# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand Generator
+# token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-fit_aligned <- function(X, status, rankmin, rankmax, order_list, B0, lambda_1_all, lambda_2_all, pfac, step_size = 1.0, niter = 2000L, linesearch_beta = 1.1, eps = 5e-7, debug = FALSE) {
-    .Call(`_mrcox_fit_aligned`, X, status, rankmin, rankmax, order_list, B0, lambda_1_all, lambda_2_all, pfac, step_size, niter, linesearch_beta, eps, debug)
+fit_aligned <- function(X, status, rankmin, rankmax, order_list, B0, lambda_1_all, 
+    lambda_2_all, pfac, step_size = 1, niter = 2000L, linesearch_beta = 1.1, eps = 1e-05)
+    {
+    .Call(`_mrcox_fit_aligned`, X, status, rankmin, rankmax, order_list, B0, lambda_1_all, 
+        lambda_2_all, pfac, step_size, niter, linesearch_beta, eps)
 }
 
 #' @export
-compute_dual_norm <- function(grad, alpha, tol) {
+compute_dual_norm <- function(grad, alpha, tol)
+{
     .Call(`_mrcox_compute_dual_norm`, grad, alpha, tol)
 }
 
 #' @export
-compute_residual <- function(X, status, rankmin, rankmax, order_list, v) {
+compute_residual <- function(X, status, rankmin, rankmax, order_list, v)
+{
     .Call(`_mrcox_compute_residual`, X, status, rankmin, rankmax, order_list, v)
 }
 

R/basil.R

@@ -4,10 +4,8 @@
 #' @importFrom magrittr %>%
 #' @export
 basil_base <- function(genotype.pfile, phe.file, responsid, covs, nlambda, lambda.min.ratio, 
-    alpha, p.factor, configs, num_lambda_per_iter, num_to_add, max_num_to_add, fit.fun, relaxed.fit.fun)
+    alpha, p.factor, configs, num_lambda_per_iter, num_to_add, max_num_to_add, fit.fun)
     {
-    fit.relaxed <- if(is.null(relaxed.fit.fun)) FALSE else TRUE
-    relaxed.save <- if(fit.relaxed) list() else NULL
     time.start <- Sys.time()
     responsid <- as.character(responsid)
     ### Get ids specified by psam --------------------------------------
@@ -173,9 +171,6 @@ basil_base <- function(genotype.pfile, phe.file, responsid, covs, nlambda, lambd
     }
     snpnetLoggerTimeDiff(sprintf("End metric evaluations for basil iteration %d.", 
         iter), time.basilmetric.start, indent = 3)
-
-    Ctrain_relaxed = Ctrain
-    Cval_relaxed = Cval
 
     ### Compute residuals and gradient-------------------------------
     residuals <- matrix(residuals, nrow = length(phe_train$ID), ncol = K, dimnames = list(paste(phe_train$ID, 
@@ -238,9 +233,6 @@ basil_base <- function(genotype.pfile, phe.file, responsid, covs, nlambda, lambd
         printf("Current maximum valid index is: %d\n", max_valid_index)
         printf("Current validation C-Indices are:\n")
         print(Cval[, 1:max_valid_index])
-
-        printf("Current validation (relaxed) C-Indices are:\n")
-        print(Cval_relaxed[, 1:max_valid_index])
 
         if (length(features.to.discard) > 0)
         {
@@ -356,34 +348,10 @@ basil_base <- function(genotype.pfile, phe.file, responsid, covs, nlambda, lambd
             for (j in 1:local_valid)
             {
                 out[[max_valid_index + j]] <- result[[j]]
-                if(fit.relaxed)
-                {
-                    relaxed.save[[max_valid_index + j]] <- list()
-                }
-
                 for (i in 1:K)
                 {
                   beta <- result[[j]][, i]
                   ind <- match(current_response[i], responsid)
-
-                  # relaxed Lasso fit
-                  if(fit.relaxed)
-                  {
-                      nnz_entries = which(beta != 0)
-                      beta_local = beta[nnz_entries]
-                      X_local = X[, nnz_entries]
-                      relaxed_fit = relaxed.fit.fun(X_local, y_list[[i]], status_list[[i]], beta0=beta_local, standardize=F)
-                      beta_local = beta
-                      beta_local[nnz_entries] = relaxed_fit
-                      relaxed_Ctrain = cindex::CIndex(X %*% beta_local, 
-                                                y_list[[i]], status_list[[i]])
-                      relaxed_Cval <- cindex::CIndex(X_val %*% beta_local, phe_val[[responses[i]]], 
-                                        phe_val[[status[i]]])
-                      Ctrain_relaxed[ind, max_valid_index + j] <- relaxed_Ctrain
-                      Cval_relaxed[ind, max_valid_index + j] <- relaxed_Cval
-                      relaxed.save[[max_valid_index + j]][[current_response[i]]] <- relaxed_fit
-
-                  }
                   Ctrain[ind, max_valid_index + j] <- cindex::CIndex(X %*% beta, 
                     y_list[[i]], status_list[[i]])
                   cval_tmp <- cindex::CIndex(X_val %*% beta, phe_val[[responses[i]]], 
@@ -402,13 +370,6 @@ basil_base <- function(genotype.pfile, phe.file, responsid, covs, nlambda, lambd
                 iter), time.basilmetric.start, indent = 3)
             # Save temp result to files
             save_list <- list(Ctrain = Ctrain, Cval = Cval, beta = out)
-                        if(fit.relaxed)
-            {
-                relaxed.save[[max_valid_index + j]][['Ctrain']] <- Ctrain_relaxed
-                relaxed.save[[max_valid_index + j]][['Cval']] <- Cval_relaxed
-                save_list[['relaxed']] <- relaxed.save
-            }
-
             save(save_list, file = file.path(configs[["save.dir"]], paste0("saveresult", 
                 iter, ".RData")))
 
@@ -516,5 +477,5 @@ basil <- function(genotype.pfile, phe.file, responsid, covs = NULL, nlambda = 10
     {
     basil_base(genotype.pfile, phe.file, responsid, covs, nlambda, lambda.min.ratio, 
         alpha, p.factor, configs, num_lambda_per_iter, num_to_add, max_num_to_add, 
-        solve_aligned, coxph_MKL)
+        solve_aligned)
 }

R/mrcox.R

@@ -80,42 +80,3 @@ get_dual_norm <- function(grad, alpha, tol = 1e-10)
     names(dnorm) <- rownames(grad)
     return(dnorm)
 }
-
-#' @export
-coxph_MKL <- function(X, y, status, beta0=NULL, standardize=T)
-{
-    if(any(is.na(X)))
-    {
-        stop("We do not allow NAs in the predictor matrix")
-    }
-
-    if(any(is.na(y)))
-    {
-        stop("We do not allow NAs in the time vector")
-    }
-
-    if(any(is.na(status)))
-    {
-        stop("We do not allow NAs in the status vector")
-    }
-
-    if(standardize)
-    {
-        mu = apply(X, 2, mean)
-        sigma  = apply(X, 2, sd)
-        X = sweep(X, 2, mu, FUN='-')
-        X = sweep(X, 2, sigma, FUN='/')
-    }
-    if(!is.null(beta0))
-    {
-        beta0 = as.matrix(beta0, ncol(X), 1)
-    }
-    result = solve_aligned(X, list(y), list(status), c(0.0), c(0.0), p.fac = NULL, 
-                B0 = beta0)
-    result = as.vector(result[["result"]][[1]])
-    if(standardize)
-    {
-        result  = result / sigma
-    }
-    result
-}

src/Makevars

@@ -1,2 +1,11 @@
+UNAME_S := $(shell uname -s)
+ifeq ($(UNAME_S),Linux)
 PKG_CXXFLAGS = -DMKL_ILP64 -m64 -I${MKLROOT}/include -Wno-ignored-attributes $(SHLIB_OPENMP_CXXFLAGS) -O3 -march=native
-PKG_LIBS=-Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a ${MKLROOT}/lib/intel64/libmkl_gnu_thread.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lgomp -lpthread -lm -ldl
+PKG_LIBS=-Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a ${MKLROOT}/lib/intel64/libmkl_gnu_thread.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lgomp -lpthread -lm -ldl
+#PKG_CXXFLAGS= -Wno-ignored-attributes $(SHLIB_OPENMP_CXXFLAGS) -O3 -march=native 
+endif
+
+#ifeq ($(UNAME_S),Darwin)
+#PKG_CXXFLAGS = -Wno-ignored-attributes -Xpreprocessor $(SHLIB_OPENMP_CXXFLAGS)
+#endif
+

src/RcppExports.cpp

@@ -8,8 +8,8 @@
 using namespace Rcpp;
 
 // fit_aligned
-Rcpp::List fit_aligned(Rcpp::NumericMatrix X, Rcpp::NumericMatrix status, Rcpp::IntegerMatrix rankmin, Rcpp::IntegerMatrix rankmax, Rcpp::List order_list, Rcpp::NumericMatrix B0, Rcpp::NumericVector lambda_1_all, Rcpp::NumericVector lambda_2_all, VectorXd pfac, double step_size, int niter, double linesearch_beta, double eps, bool debug);
-RcppExport SEXP _mrcox_fit_aligned(SEXP XSEXP, SEXP statusSEXP, SEXP rankminSEXP, SEXP rankmaxSEXP, SEXP order_listSEXP, SEXP B0SEXP, SEXP lambda_1_allSEXP, SEXP lambda_2_allSEXP, SEXP pfacSEXP, SEXP step_sizeSEXP, SEXP niterSEXP, SEXP linesearch_betaSEXP, SEXP epsSEXP, SEXP debugSEXP) {
+Rcpp::List fit_aligned(Rcpp::NumericMatrix X, Rcpp::NumericMatrix status, Rcpp::IntegerMatrix rankmin, Rcpp::IntegerMatrix rankmax, Rcpp::List order_list, Rcpp::NumericMatrix B0, Rcpp::NumericVector lambda_1_all, Rcpp::NumericVector lambda_2_all, VectorXd pfac, double step_size, int niter, double linesearch_beta, double eps);
+RcppExport SEXP _mrcox_fit_aligned(SEXP XSEXP, SEXP statusSEXP, SEXP rankminSEXP, SEXP rankmaxSEXP, SEXP order_listSEXP, SEXP B0SEXP, SEXP lambda_1_allSEXP, SEXP lambda_2_allSEXP, SEXP pfacSEXP, SEXP step_sizeSEXP, SEXP niterSEXP, SEXP linesearch_betaSEXP, SEXP epsSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -26,8 +26,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< int >::type niter(niterSEXP);
     Rcpp::traits::input_parameter< double >::type linesearch_beta(linesearch_betaSEXP);
     Rcpp::traits::input_parameter< double >::type eps(epsSEXP);
-    Rcpp::traits::input_parameter< bool >::type debug(debugSEXP);
-    rcpp_result_gen = Rcpp::wrap(fit_aligned(X, status, rankmin, rankmax, order_list, B0, lambda_1_all, lambda_2_all, pfac, step_size, niter, linesearch_beta, eps, debug));
+    rcpp_result_gen = Rcpp::wrap(fit_aligned(X, status, rankmin, rankmax, order_list, B0, lambda_1_all, lambda_2_all, pfac, step_size, niter, linesearch_beta, eps));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -62,7 +61,7 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_mrcox_fit_aligned", (DL_FUNC) &_mrcox_fit_aligned, 14},
+    {"_mrcox_fit_aligned", (DL_FUNC) &_mrcox_fit_aligned, 13},
     {"_mrcox_compute_dual_norm", (DL_FUNC) &_mrcox_compute_dual_norm, 3},
     {"_mrcox_compute_residual", (DL_FUNC) &_mrcox_compute_residual, 6},
     {NULL, NULL, 0}

src/aligned.cpp

@@ -195,6 +195,21 @@ void prox_l2(Eigen::Map<Eigen::MatrixXd> & B, const double step_size, double lam
     B =  ((B_row_norm.array() - lambda_2*step_size*penalty_factor.array())/(B_row_norm.array())).matrix().asDiagonal() * B;
 }
 
+void update_parameters(MatrixXd & B, const MatrixXd & grad, const MatrixXd &v, const double step_size,
+                       double lambda_1, double lambda_2, const VectorXd & penalty_factor,
+                       VectorXd & B_row_norm)
+{
+    int K = grad.cols();
+    B.noalias() = v - step_size*grad;
+    // Apply proximal operator here:
+    //Soft-thresholding
+    B = ((B.cwiseAbs().colwise() - lambda_1*step_size*penalty_factor).array().max(0) * B.array().sign()).matrix();
+    // Group soft-thresholding
+    // should be called the pmax of B_row_norm  and lambda_2*step_size
+    B_row_norm.noalias() = B.rowwise().norm().cwiseMax(lambda_2*step_size*penalty_factor);
+    B =  ((B_row_norm.array() - lambda_2*step_size*penalty_factor.array())/(B_row_norm.array())).matrix().asDiagonal() * B;
+}
+
 
 // [[Rcpp::export]]
 Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
@@ -209,8 +224,7 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
                        double step_size = 1.0,
                        int niter=2000,
                        double linesearch_beta = 1.1,
-                       double eps=5e-7, // convergence criteria
-                       bool debug=false
+                       double eps=1e-5 // convergence criteria
                        )
 {
     int N = X.rows();
@@ -271,6 +285,8 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
                 prox_l2(B, step_size, lambda_2, pfac, B_row_norm);
 
 
+                //update_parameters(B, grad, v, step_size, lambda_1, lambda_2, pfac, B_row_norm);
+
                 double cox_val_next = prob.get_value_only(B.data());
 
                 if(!std::isfinite(cox_val_next)){
@@ -289,15 +305,12 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
             }
 
 
-            if(value_change < eps){
-                if(debug){
-                    std::cout << "convergence based on value change reached in " << i <<" iterations\n";
-                    std::cout << "current step size is " << step_size << std::endl;
-                    gettimeofday(&end, NULL);
-                    double delta  = ((end.tv_sec  - start.tv_sec) * 1000000u + end.tv_usec - start.tv_usec) / 1.e6;
-                    std::cout <<  "elapsed time is " << delta << " seconds" << std::endl;
-                }
-
+            if(value_change < 5e-7){
+                std::cout << "convergence based on value change reached in " << i <<" iterations\n";
+                std::cout << "current step size is " << step_size << std::endl;
+                gettimeofday(&end, NULL);
+                double delta  = ((end.tv_sec  - start.tv_sec) * 1000000u + end.tv_usec - start.tv_usec) / 1.e6;
+                std::cout <<  "elapsed time is " << delta << " seconds" << std::endl;
                 Rcpp::checkUserInterrupt();
                 break;
             }
@@ -308,19 +321,17 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
             weight_old = weight_new;
 
             if (i != 0 && i % 100 == 0){
-                if(debug){
-                    std::cout << "reached " << i << " iterations\n";
-                    gettimeofday(&end, NULL);
-                    double delta  = ((end.tv_sec  - start.tv_sec) * 1000000u + end.tv_usec - start.tv_usec) / 1.e6;
-                    std::cout <<  "elapsed time is " << delta  << " seconds" << std::endl;
-                    std::cout << "current step size is " << step_size << std::endl;
-                }
+                std::cout << "reached " << i << " iterations\n";
+                gettimeofday(&end, NULL);
+                double delta  = ((end.tv_sec  - start.tv_sec) * 1000000u + end.tv_usec - start.tv_usec) / 1.e6;
+                std::cout <<  "elapsed time is " << delta  << " seconds" << std::endl;
+                std::cout << "current step size is " << step_size << std::endl;
                 Rcpp::checkUserInterrupt();
             }
         }
         result[lam_ind] = Bfull;
         residual_result[lam_ind] = prob.Rget_residual(B.data());
-        // std::cout << "Solution for the " <<  lam_ind+1 << "th lambda pair is obtained\n";
+        std::cout << "Solution for the " <<  lam_ind+1 << "th lambda pair is obtained\n";
     }
     return Rcpp::List::create(Rcpp::Named("result") = result,
                               Rcpp::Named("residual") = residual_result);