WIP: running but not agreeing with other solver

jonathan-taylor · jonathan-taylor · commit c26a06b24308 · 2017-09-26T21:29:58.000-07:00
diff --git a/selectiveInference/R/RcppExports.R b/selectiveInference/R/RcppExports.R
@@ -1,19 +1,19 @@
-# This file was generated by Rcpp::compileAttributes
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 solve_QP <- function(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active) {
-    .Call('selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
+    .Call('_selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
 }
 
 solve_QP_wide <- function(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active) {
-    .Call('selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)
+    .Call('_selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)
 }
 
 update1_ <- function(Q2, w, m, k) {
-    .Call('selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
+    .Call('_selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
 }
 
 downdate1_ <- function(Q1, R, j0, m, n) {
-    .Call('selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
+    .Call('_selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
 }
 
diff --git a/selectiveInference/src/RcppExports.cpp b/selectiveInference/src/RcppExports.cpp
@@ -1,4 +1,4 @@
-// This file was generated by Rcpp::compileAttributes
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
 // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 #include <Rcpp.h>
@@ -7,10 +7,10 @@ using namespace Rcpp;
 
 // solve_QP
 Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma, double bound, int maxiter, Rcpp::NumericVector theta, Rcpp::NumericVector linear_func, Rcpp::NumericVector gradient, Rcpp::IntegerVector ever_active, Rcpp::IntegerVector nactive, double kkt_tol, double objective_tol, int max_active);
-RcppExport SEXP selectiveInference_solve_QP(SEXP SigmaSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
+RcppExport SEXP _selectiveInference_solve_QP(SEXP SigmaSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Sigma(SigmaSEXP);
     Rcpp::traits::input_parameter< double >::type bound(boundSEXP);
     Rcpp::traits::input_parameter< int >::type maxiter(maxiterSEXP);
@@ -22,16 +22,16 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< double >::type kkt_tol(kkt_tolSEXP);
     Rcpp::traits::input_parameter< double >::type objective_tol(objective_tolSEXP);
     Rcpp::traits::input_parameter< int >::type max_active(max_activeSEXP);
-    __result = Rcpp::wrap(solve_QP(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(solve_QP(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active));
+    return rcpp_result_gen;
 END_RCPP
 }
 // solve_QP_wide
 Rcpp::List solve_QP_wide(Rcpp::NumericMatrix X, double bound, int maxiter, Rcpp::NumericVector theta, Rcpp::NumericVector linear_func, Rcpp::NumericVector gradient, Rcpp::NumericVector X_theta, Rcpp::IntegerVector ever_active, Rcpp::IntegerVector nactive, double kkt_tol, double objective_tol, int max_active);
-RcppExport SEXP selectiveInference_solve_QP_wide(SEXP XSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP X_thetaSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
+RcppExport SEXP _selectiveInference_solve_QP_wide(SEXP XSEXP, SEXP boundSEXP, SEXP maxiterSEXP, SEXP thetaSEXP, SEXP linear_funcSEXP, SEXP gradientSEXP, SEXP X_thetaSEXP, SEXP ever_activeSEXP, SEXP nactiveSEXP, SEXP kkt_tolSEXP, SEXP objective_tolSEXP, SEXP max_activeSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type X(XSEXP);
     Rcpp::traits::input_parameter< double >::type bound(boundSEXP);
     Rcpp::traits::input_parameter< int >::type maxiter(maxiterSEXP);
@@ -44,36 +44,49 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< double >::type kkt_tol(kkt_tolSEXP);
     Rcpp::traits::input_parameter< double >::type objective_tol(objective_tolSEXP);
     Rcpp::traits::input_parameter< int >::type max_active(max_activeSEXP);
-    __result = Rcpp::wrap(solve_QP_wide(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(solve_QP_wide(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active));
+    return rcpp_result_gen;
 END_RCPP
 }
 // update1_
 Rcpp::List update1_(Rcpp::NumericMatrix Q2, Rcpp::NumericVector w, int m, int k);
-RcppExport SEXP selectiveInference_update1_(SEXP Q2SEXP, SEXP wSEXP, SEXP mSEXP, SEXP kSEXP) {
+RcppExport SEXP _selectiveInference_update1_(SEXP Q2SEXP, SEXP wSEXP, SEXP mSEXP, SEXP kSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Q2(Q2SEXP);
     Rcpp::traits::input_parameter< Rcpp::NumericVector >::type w(wSEXP);
     Rcpp::traits::input_parameter< int >::type m(mSEXP);
     Rcpp::traits::input_parameter< int >::type k(kSEXP);
-    __result = Rcpp::wrap(update1_(Q2, w, m, k));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(update1_(Q2, w, m, k));
+    return rcpp_result_gen;
 END_RCPP
 }
 // downdate1_
 Rcpp::List downdate1_(Rcpp::NumericMatrix Q1, Rcpp::NumericMatrix R, int j0, int m, int n);
-RcppExport SEXP selectiveInference_downdate1_(SEXP Q1SEXP, SEXP RSEXP, SEXP j0SEXP, SEXP mSEXP, SEXP nSEXP) {
+RcppExport SEXP _selectiveInference_downdate1_(SEXP Q1SEXP, SEXP RSEXP, SEXP j0SEXP, SEXP mSEXP, SEXP nSEXP) {
 BEGIN_RCPP
-    Rcpp::RObject __result;
-    Rcpp::RNGScope __rngScope;
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Q1(Q1SEXP);
     Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type R(RSEXP);
     Rcpp::traits::input_parameter< int >::type j0(j0SEXP);
     Rcpp::traits::input_parameter< int >::type m(mSEXP);
     Rcpp::traits::input_parameter< int >::type n(nSEXP);
-    __result = Rcpp::wrap(downdate1_(Q1, R, j0, m, n));
-    return __result;
+    rcpp_result_gen = Rcpp::wrap(downdate1_(Q1, R, j0, m, n));
+    return rcpp_result_gen;
 END_RCPP
 }
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_selectiveInference_solve_QP", (DL_FUNC) &_selectiveInference_solve_QP, 11},
+    {"_selectiveInference_solve_QP_wide", (DL_FUNC) &_selectiveInference_solve_QP_wide, 12},
+    {"_selectiveInference_update1_", (DL_FUNC) &_selectiveInference_update1_, 4},
+    {"_selectiveInference_downdate1_", (DL_FUNC) &_selectiveInference_downdate1_, 5},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_selectiveInference(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}
diff --git a/selectiveInference/src/quadratic_program_wide.c b/selectiveInference/src/quadratic_program_wide.c
@@ -388,6 +388,7 @@ int solve_wide(double *X_ptr,              /* A design matrix */
 		       ncase,
 		       bound,
 		       kkt_tol) == 1) {
+      fprintf(stderr, "break1\n");
       break;
     }
 
@@ -426,6 +427,7 @@ int solve_wide(double *X_ptr,              /* A design matrix */
 		       ncase,
 		       bound,
 		       kkt_tol) == 1) {
+      fprintf(stderr, "break2\n");
       break;
     }
 					  
@@ -434,6 +436,7 @@ int solve_wide(double *X_ptr,              /* A design matrix */
     // fprintf(stderr, "here18 %d\n", *nactive_ptr);
 
     if (*nactive_ptr >= max_active) {
+      fprintf(stderr, "break3\n");
       break;
     }
 
@@ -454,6 +457,7 @@ int solve_wide(double *X_ptr,              /* A design matrix */
       // fprintf(stderr, "here8\n");
 
       if ((fabs(old_value - new_value) < objective_tol * fabs(new_value)) && (iter > 0)) {
+	fprintf(stderr, "break5 %f %f %f %d\n", old_value, new_value, objective_tol, iter);
 	break;
       }
       old_value = new_value;
diff --git a/test.R b/test.R
@@ -2,28 +2,46 @@ set.seed(43)
 
     n = 100                                                                                                                                                                              
     p = 200                                                                                                                                                                              
-    lam = 0.1                                                                                                                                                                             
+    lam = 0.2                                                                                                                                                                             
     X = matrix(rnorm(n*p), n, p)                                                                                                                                                         
     Y = rnorm(n)                                                                                                                                                                         
     library(selectiveInference)                                                                                                                                                          
     p = ncol(X)                                                                                                                                                                          
     soln_R = rep(0, p)                                                                                                                                                                   
     grad = -t(X) %*% Y / n                                                                                                                                                                   
-    ever_active = c(1, rep(0, p-1))                                                                                                                                                      
+    ever_active = as.integer(c(1, rep(0, p-1)))
     nactive = as.integer(1)                                                                                                                                                              
     kkt_tol = 1.e-12                                                                                                                                                                     
-    objective_tol = 1.e-12                                                                                                                                                               
+    objective_tol = 1.e-16                                                                                                                                                               
     maxiter = 500                                                                                                                                                                        
-    soln_R = selectiveInference:::solve_QP(t(X) %*% X / n, lam, maxiter, soln_R, -t(X) %*% Y / n, grad, ever_active, nactive, kkt_tol, objective_tol, p)$soln                                    
+    soln_R = selectiveInference:::solve_QP(t(X) %*% X / n, lam, maxiter, soln_R, -t(X) %*% Y / n, grad, ever_active, nactive, kkt_tol, objective_tol, p)
+    print('active')
+    print(nactive)
+    print(ever_active)
+    print(soln_R$ever_active)
+    soln_R = soln_R$soln
+    soln_R_old = soln_R
     print(soln_R)                                                                                                                                                                        
     Xtheta = rep(0, n)                                                                                                                                                                   
     nactive = as.integer(1)
-    ever_active = c(1, rep(0, p-1))                                                                                                                                                      
+    ever_active = as.integer(c(1, rep(0, p-1)))
     soln_R = rep(0, p)
     grad = - t(X) %*% Y / n
     # test wide solver                                                                                                                                                                   
-    soln_R_wide = selectiveInference:::solve_QP_wide(X, lam, maxiter, soln_R, -t(X) %*% Y / n, grad, Xtheta, ever_active, nactive, kkt_tol, objective_tol, p)
-    print(soln_R_wide)
-    print(soln_R)
-    print(soln_R_wide$soln)
-    print(soln_R_wide$soln - soln_R)
+    soln_R_wide = selectiveInference:::solve_QP_wide(X, lam, maxiter, soln_R*1., -t(X) %*% Y / n, grad, Xtheta, ever_active, nactive, kkt_tol, objective_tol, p)
+    print(nactive)
+    print(soln_R_wide$ever_active)
+
+    print('diff')
+    print(soln_R_wide$soln - soln_R_old)
+    print(soln_R_wide$gradient[soln_R_wide$ever_active])
+    print(max(abs(soln_R_wide$gradient[-soln_R_wide$ever_active])))
+    print(soln_R_wide$kkt_check)
+    print(soln_R_wide$iter)
+#    print(Xtheta - X %*% soln_R_wide$soln)
+#     print(Xtheta)
+
+     soln_R_wide = selectiveInference:::solve_QP_wide(X, 0.7 * lam, maxiter, soln_R_wide$soln, -t(X) %*% Y / n, grad, Xtheta, ever_active, nactive, kkt_tol, objective_tol, p)
+#     print(Xtheta - X %*% soln_R_wide$soln)
+#     print(soln_R_wide$soln)
+     soln_R_wide = selectiveInference:::solve_QP_wide(X, 0.5 * lam, maxiter, soln_R_wide$soln, -t(X) %*% Y / n, grad, Xtheta, ever_active, nactive, kkt_tol, objective_tol, p)

Original file line number	Diff line number	Diff line change
`@@ -1,19 +1,19 @@`
`1`		`-# This file was generated by Rcpp::compileAttributes`
	`1`	`+# Generated by using Rcpp::compileAttributes() -> do not edit by hand`
`2`	`2`	`# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393`
`3`	`3`
`4`	`4`	`solve_QP <- function(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active) {`
`5`		`- .Call('selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)`
	`5`	`+ .Call('_selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)`
`6`	`6`	`}`
`7`	`7`
`8`	`8`	`solve_QP_wide <- function(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active) {`
`9`		`- .Call('selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)`
	`9`	`+ .Call('_selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)`
`10`	`10`	`}`
`11`	`11`
`12`	`12`	`update1_ <- function(Q2, w, m, k) {`
`13`		`- .Call('selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)`
	`13`	`+ .Call('_selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)`
`14`	`14`	`}`
`15`	`15`
`16`	`16`	`downdate1_ <- function(Q1, R, j0, m, n) {`
`17`		`- .Call('selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)`
	`17`	`+ .Call('_selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)`
`18`	`18`	`}`
`19`	`19`
Original file line number	Diff line number	Diff line change
`@@ -388,6 +388,7 @@ int solve_wide(double X_ptr, / A design matrix */`
`388`	`388`	`ncase,`
`389`	`389`	`bound,`
`390`	`390`	`kkt_tol) == 1) {`
	`391`	`+ fprintf(stderr, "break1\n");`
`391`	`392`	`break;`
`392`	`393`	`}`
`393`	`394`
`@@ -426,6 +427,7 @@ int solve_wide(double X_ptr, / A design matrix */`
`426`	`427`	`ncase,`
`427`	`428`	`bound,`
`428`	`429`	`kkt_tol) == 1) {`
	`430`	`+ fprintf(stderr, "break2\n");`
`429`	`431`	`break;`
`430`	`432`	`}`
`431`	`433`
`@@ -434,6 +436,7 @@ int solve_wide(double X_ptr, / A design matrix */`
`434`	`436`	`// fprintf(stderr, "here18 %d\n", *nactive_ptr);`
`435`	`437`
`436`	`438`	`if (*nactive_ptr >= max_active) {`
	`439`	`+ fprintf(stderr, "break3\n");`
`437`	`440`	`break;`
`438`	`441`	`}`
`439`	`442`
`@@ -454,6 +457,7 @@ int solve_wide(double X_ptr, / A design matrix */`
`454`	`457`	`// fprintf(stderr, "here8\n");`
`455`	`458`
`456`	`459`	`if ((fabs(old_value - new_value) < objective_tol * fabs(new_value)) && (iter > 0)) {`
	`460`	`+ fprintf(stderr, "break5 %f %f %f %d\n", old_value, new_value, objective_tol, iter);`
`457`	`461`	`break;`
`458`	`462`	`}`
`459`	`463`	`old_value = new_value;`