remove unused parameter, exit when prox does not converge

rivas-lab · Dec 1, 2020 · 1f31cc7 · 1f31cc7
1 parent d00dcea
commit 1f31cc7
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 18 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,22 +1,17 @@
-# 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, 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)
+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) {
+    .Call(`_mrcox_fit_aligned`, X, status, rankmin, rankmax, order_list, B0, lambda_1_all, lambda_2_all, pfac, step_size, niter, linesearch_beta)
 }
 
 #' @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)
 }
 
diff --git a/src/RcppExports.cpp b/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);
-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) {
+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);
+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) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -25,8 +25,7 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< double >::type step_size(step_sizeSEXP);
     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_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));
+    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));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -61,7 +60,7 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_mrcox_fit_aligned", (DL_FUNC) &_mrcox_fit_aligned, 13},
+    {"_mrcox_fit_aligned", (DL_FUNC) &_mrcox_fit_aligned, 12},
     {"_mrcox_compute_dual_norm", (DL_FUNC) &_mrcox_compute_dual_norm, 3},
     {"_mrcox_compute_residual", (DL_FUNC) &_mrcox_compute_residual, 6},
     {NULL, NULL, 0}

diff --git a/src/aligned.cpp b/src/aligned.cpp
@@ -223,8 +223,7 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
                        VectorXd pfac,
                        double step_size = 1.0,
                        int niter=2000,
-                       double linesearch_beta = 1.1,
-                       double eps=1e-5 // convergence criteria
+                       double linesearch_beta = 1.1
                        )
 {
     int N = X.rows();
@@ -257,6 +256,8 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
     double weight_old, weight_new;
     double rhs_ls;
     struct timeval start, end;
+    int num_iter;
+    bool nan_stop = false;
 
     Rcpp::List result(nlambda);
     Rcpp::List residual_result(nlambda);
@@ -291,6 +292,10 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
 
                 if(!std::isfinite(cox_val_next)){
                     stop = (step_size < 1e-9);
+                    if(stop){
+                        nan_stop = true;
+                        goto terminate_nan;
+                    }
                 } else {
                     rhs_ls = cox_val + (grad.array() * (B - v).array()).sum() + (B-v).squaredNorm()/(2*step_size);
                     stop = (cox_val_next <= rhs_ls);
@@ -299,6 +304,7 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
 
                 if (stop){
                     value_change = abs(cox_val_next - cox_val)/fmax(1.0, abs(cox_val));
+                    num_iter = i;
                     break;
                 }
                 step_size /= linesearch_beta;
@@ -325,6 +331,10 @@ Rcpp::List fit_aligned(Rcpp::NumericMatrix X,
         }
         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, number of iterations: " << num_iter <<std::endl;
+    }
+    if(nan_stop){
+        terminate_nan:Rf_warning("Proximal gradient did not converge for some lambdas, return finite results only\n");
     }
     return Rcpp::List::create(Rcpp::Named("result") = result,
                               Rcpp::Named("residual") = residual_result);