merging with upstream

Author: jonathan-taylor

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "C-software"]
+	path = C-software
+	url = https://github.com/selective-inference/C-software.git

.travis.yml

@@ -12,5 +12,10 @@ warnings_are_errors: true
 before_install:
   - tlmgr install index # for texlive and vignette?
   - R -e 'install.packages(c("Rcpp", "intervals"), repos="http://cloud.r-project.org")'
+  - cd C-software
+  - git submodule init
+  - git submodule update
+  - cd ..
+  - make src
   - make Rcpp
   - cd selectiveInference

C-software

@@ -3,12 +3,15 @@ Rcpp:
 	- rm -f selectiveInference/R/RcppExports.R
 	Rscript -e "library(Rcpp); Rcpp::compileAttributes('selectiveInference')"
 
-install: Rcpp
+install: Rcpp src
 	R CMD INSTALL selectiveInference
 
-build: 
+build: src 
 	R CMD build selectiveInference
 
-check: Rcpp build
+src:
+	cp C-software/src/* selectiveInference/src
+
+check: Rcpp build 
 	R CMD build selectiveInference
 	R CMD check selectiveInference_1.2.2.tar.gz # fix this to be a script variable

Makefile

@@ -9,7 +9,7 @@ Maintainer: Rob Tibshirani <[email protected]>
 Depends:
     glmnet,
     intervals,
-    survival
+    survival,
 Suggests:
     Rmpfr
 Description: New tools for post-selection inference, for use with forward

selectiveInference/DESCRIPTION

@@ -43,5 +43,6 @@ importFrom("stats", dnorm, lsfit, pexp, pnorm, predict,
            qnorm, rnorm, sd, uniroot, dchisq, model.matrix, pchisq)
 importFrom("stats", "coef", "df", "lm", "pf")
 importFrom("stats", "glm", "residuals", "vcov")
+importFrom("stats", "rbinom", "rexp")
 importFrom("Rcpp", "sourceCpp")
 importFrom("distr", "Norm", "DExp")

selectiveInference/NAMESPACE

@@ -319,16 +319,19 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
                            nsample, 
                            rows, 
  		           verbose=FALSE, 
-		           mu=NULL,             # starting value of mu
+		           bound=NULL,             # starting value of bound
    			   linesearch=TRUE,     # do a linesearch?
    		           scaling_factor=1.5,  # multiplicative factor for linesearch
 			   max_active=NULL,     # how big can active set get?
 			   max_try=10,          # how many steps in linesearch?
 			   warn_kkt=FALSE,      # warn if KKT does not seem to be satisfied?
-			   max_iter=100,        # how many iterations for each optimization problem
+			   max_iter=50,         # how many iterations for each optimization problem
+                           kkt_stop=TRUE,       # stop based on KKT conditions?
+                           parameter_stop=TRUE, # stop based on relative convergence of parameter?
+			   objective_stop=TRUE, # stop based on relative decrease in objective?
                            kkt_tol=1.e-4,       # tolerance for the KKT conditions
                            parameter_tol=1.e-4, # tolerance for relative convergence of parameter
-			   objective_tol=1.e-8  # tolerance for relative decrease in objective
+			   objective_tol=1.e-4  # tolerance for relative decrease in objective
                            ) {
 
 
@@ -339,8 +342,8 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
   p = ncol(Xinfo);
   M = matrix(0, length(rows), p);
 
-  if (is.null(mu)) {
-      mu = (1/sqrt(nsample)) * qnorm(1-(0.1/(p^2)))
+  if (is.null(bound)) {
+      bound = (1/sqrt(nsample)) * qnorm(1-(0.1/(p^2)))
   }
 
   xperc = 0;
@@ -356,13 +359,16 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
     output = debiasingRow(Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
                           is_wide,
                           row,
-                          mu,
+                          bound,
                           linesearch=linesearch,
                           scaling_factor=scaling_factor,
                           max_active=max_active,
 			  max_try=max_try,
 			  warn_kkt=FALSE,
 			  max_iter=max_iter,
+			  kkt_stop=kkt_stop,
+			  parameter_stop=parameter_stop,
+			  objective_stop=objective_stop,
 			  kkt_tol=kkt_tol,
 			  parameter_tol=parameter_tol,
 			  objective_tol=objective_tol)
@@ -387,16 +393,23 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
 debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
                          is_wide, 
                          row, 
-                         mu, 
+                         bound, 
 	                 linesearch=TRUE,     # do a linesearch?
 		         scaling_factor=1.5,  # multiplicative factor for linesearch
 		         max_active=NULL,     # how big can active set get?
 			 max_try=10,          # how many steps in linesearch?
 			 warn_kkt=FALSE,      # warn if KKT does not seem to be satisfied?
-			 max_iter=100,        # how many iterations for each optimization problem
+			 max_iter=50,         # how many iterations for each optimization problem
+                         kkt_stop=TRUE,       # stop based on KKT conditions?
+                         parameter_stop=TRUE, # stop based on relative convergence of parameter?
+                         objective_stop=TRUE, # stop based on relative decrease in objective?
                          kkt_tol=1.e-4,       # tolerance for the KKT conditions
 			 parameter_tol=1.e-4, # tolerance for relative convergence of parameter
+<<<<<<< HEAD
 			 objective_tol=1.e-8  # tolerance for relative decrease in objective
+=======
+			 objective_tol=1.e-4  # tolerance for relative decrease in objective
+>>>>>>> 232760d6aef5182e040b82e30555f4af5ad6803c
                          ) {
 
   p = ncol(Xinfo)
@@ -405,9 +418,11 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
       max_active = min(nrow(Xinfo), ncol(Xinfo))
   }
 
+   
   # Initialize variables 
 
   soln = rep(0, p)
+  soln = as.numeric(soln)
   ever_active = rep(0, p)
   ever_active[1] = row      # 1-based
   ever_active = as.integer(ever_active)
@@ -423,11 +438,15 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
 
   last_output = NULL
 
+  if (is_wide) {
+     Xsoln = as.numeric(rep(0, nrow(Xinfo)))
+  }
+
   while (counter_idx < max_try) {
 
       if (!is_wide) {
           result = solve_QP(Xinfo, # this is non-neg-def matrix
-                            mu, 
+                            bound, 
                             max_iter, 
                             soln, 
                             linear_func, 
@@ -438,6 +457,7 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
                             objective_tol, 
 			    parameter_tol,
                             max_active,
+<<<<<<< HEAD
 			    FALSE,        # objective_stop
 			    FALSE,        # kkt_stop
 			    TRUE)         # param_stop
@@ -446,6 +466,15 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
           result = solve_QP_wide(Xinfo, # this is a design matrix
                                  rep(mu, p),  # vector of Lagrange multipliers
 				 0,           # ridge_term 
+=======
+			    kkt_stop,
+			    objective_stop,
+			    parameter_stop)
+      } else {
+          result = solve_QP_wide(Xinfo,                      # this is a design matrix
+                                 as.numeric(rep(bound, p)),  # vector of Lagrange multipliers
+				 0,                          # ridge_term 
+>>>>>>> 232760d6aef5182e040b82e30555f4af5ad6803c
                                  max_iter, 
                                  soln, 
                                  linear_func, 
@@ -457,9 +486,15 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
                                  objective_tol, 
 				 parameter_tol,
                                  max_active,
+<<<<<<< HEAD
 				 FALSE,       # objective_stop
 				 FALSE,       # kkt_stop
 				 TRUE)        # param_stop
+=======
+				 kkt_stop,
+				 objective_stop,	
+				 parameter_stop)
+>>>>>>> 232760d6aef5182e040b82e30555f4af5ad6803c
 
       }
 
@@ -483,13 +518,13 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
          if ((iter < (max_iter+1)) && (counter_idx > 1)) { 
            break;      # we've found a feasible point and solved the problem            
          }
-         mu = mu * scaling_factor;
+         bound = bound * scaling_factor;
       } else {         # trying to drop the bound parameter further
          if ((iter == (max_iter + 1)) && (counter_idx > 1)) {
             result = last_output; # problem seems infeasible because we didn't solve it
    	    break;                # so we revert to previously found solution
          }
-         mu = mu / scaling_factor;
+         bound = bound / scaling_factor;
       }
 
       # If the active set has grown to a certain size
@@ -515,7 +550,8 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
   } 
 
   return(list(soln=result$soln,
-              kkt_check=result$kkt_check))
+              kkt_check=result$kkt_check,
+	      gradient=result$gradient))
 
 }
 

selectiveInference/R/funs.fixed.R

@@ -73,7 +73,7 @@ b1= -(mydiag(sign_bhat)%*%MM)%*%sign_bhat*lambda
        vup[jj]=junk$vup
        sd[jj]=junk$sd
 
-      junk2=TG.interval(bbar, A1, b1, vj, MM, alpha)
+      junk2=TG.interval(bbar, A1, b1, vj, MM, alpha, flip=(sign_bhat[jj]==-1))
        ci[jj,]=junk2$int
        tailarea[jj,] = junk2$tailarea
 

selectiveInference/R/funs.fixedCox.R

@@ -96,7 +96,7 @@ fixedLogitLassoInf=function(x,y,beta,lambda,alpha=.1, type=c("partial"), tol.bet
    vup[jj]=junk$vup
        sd[jj]=junk$sd
 
-     junk2=TG.interval(bbar, A1, b1, vj, MM,alpha=alpha)
+     junk2=TG.interval(bbar, A1, b1, vj, MM,alpha=alpha, flip=(sign_bhat[jj+1]==-1))
 
      ci[jj,]=junk2$int
      tailarea[jj,] = junk2$tailarea

selectiveInference/R/funs.fixedLogit.R

@@ -3,19 +3,19 @@
 #
 # min 1/2 || y - \beta_0 - X \beta ||_2^2 + \lambda || \beta ||_1 - \omega^T\beta + \frac{\epsilon}{2} \|\beta\|^2_2
 
-fit_randomized_lasso = function(X, 
-                                y, 
-                                lam, 
-                                noise_scale, 
-                                ridge_term, 
-                                noise_type=c('gaussian', 'laplace'),
-                                max_iter=100,        # how many iterations for each optimization problem
-                                kkt_tol=1.e-4,       # tolerance for the KKT conditions
-                                parameter_tol=1.e-8, # tolerance for relative convergence of parameter
-                                objective_tol=1.e-8, # tolerance for relative decrease in objective
-                                objective_stop=FALSE,
-                                kkt_stop=TRUE,
-                                param_stop=TRUE)
+randomizedLASSO = function(X, 
+                           y, 
+                           lam, 
+                           noise_scale, 
+                           ridge_term, 
+                           noise_type=c('gaussian', 'laplace'),
+                           max_iter=100,        # how many iterations for each optimization problem
+                           kkt_tol=1.e-4,       # tolerance for the KKT conditions
+                           parameter_tol=1.e-8, # tolerance for relative convergence of parameter
+                           objective_tol=1.e-8, # tolerance for relative decrease in objective
+                           objective_stop=FALSE,
+                           kkt_stop=TRUE,
+                           param_stop=TRUE)
 {
 
     n = nrow(X); p = ncol(X)
@@ -24,12 +24,11 @@ fit_randomized_lasso = function(X,
 
     if (noise_scale > 0) {
         if (noise_type == 'gaussian') {
-            D = Norm(mean=0, sd=noise_scale)
+            perturb_ = rnorm(p) * noise_scale
         }
         else if (noise_type == 'laplace') {
-            D = DExp(rate = 1 / noise_scale) # D is a Laplace distribution with rate = 1.
+            perturb_ = rexp(p) * (2 * rbinom(p, 1, 0.5) - 1) * noise_scale
         }
-        perturb_ = distr::r(D)(p)
     } else {
         perturb_ = rep(0, p)
     }
@@ -38,19 +37,21 @@ fit_randomized_lasso = function(X,
     if (length(lam) == 1) {
        lam = rep(lam, p)
     }
+
     if (length(lam) != p) {
        stop("Lagrange parameter should be single float or of length ncol(X)")
     }    
 
     soln = rep(0, p)
     Xsoln = rep(0, n)
-    linear_func = (- t(X) %*% y - perturb_)
+    linear_func = (- t(X) %*% y - perturb_) / n
+
     gradient = 1. * linear_func
     ever_active = rep(0, p)
     nactive = as.integer(0)
 
     result = solve_QP_wide(X,                  # design matrix
-    	                   lam,                # vector of Lagrange multipliers
+    	                   lam / n,            # vector of Lagrange multipliers
 		           ridge_term / n,     # ridge_term 
                            max_iter, 
                            soln, 
@@ -66,5 +67,65 @@ fit_randomized_lasso = function(X,
 		           objective_stop,     # objective_stop
 			   kkt_stop,           # kkt_stop
 			   param_stop)         # param_stop
-    return(result)
+    
+    sign_soln = sign(result$soln)
+
+    unpenalized = lam == 0
+    active = (!unpenalized) & (sign_soln != 0)
+    inactive = (!unpenalized) & (sign_soln == 0)
+
+    unpenalized_set = which(unpenalized)
+    active_set = which(active)
+    inactive_set = which(inactive)
+
+    # affine transform for optimization variables
+
+    E = c(unpenalized_set, active_set)
+    I = inactive_set
+    X_E = X[,E]
+    X_I = X[,I]
+    L_E = t(X) %*% X[,E]
+
+    coef_term = L_E
+    coef_term = coef_term %*% diag(c(rep(1, sum(unpenalized)), sign_soln[active]))  # coefficients are non-negative
+    coef_term[active,] = coef_term[active,] + ridge_term * diag(rep(1, sum(active)))  # ridge term
+
+    subgrad_term = matrix(0, p, sum(inactive)) # for subgrad
+    for (i in 1:sum(inactive)) {
+        subgrad_term[inactive_set[i], i] = 1
+    }
+
+    linear_term = cbind(coef_term,
+                        subgrad_term)
+
+    offset_term = rep(0, p)
+    offset_term[active] = lam[active] * sign_soln[active]
+
+    opt_transform = list(linear_term=linear_term,
+                         offset_term=offset_term)
+
+    # affine transform for internal (data) variables
+    # for now just use parametric in terms of
+    # (\bar{\beta}_E, X_{-E}^T(y-X_E\bar{\beta}_E)
+    # 
+    # we have to reconstruct -X^TY from this pair
+    #
+
+    active_term = -L_E                           # for \bar{\beta}_E
+
+    inactive_term = -subgrad_term
+    linear_term = cbind(active_term,
+                        inactive_term)
+    offset_term = rep(0, p)
+    internal_transform = list(linear_term = linear_term,
+                              offset_term = offset_term)
+
+    return(list(active_set = active_set,
+                inactive_set = inactive_set,
+                unpenalized_set = unpenalized_set,
+                sign_soln = sign_soln,
+                optimization_transform = opt_transform,
+                internal_transform = internal_transform
+                ))
+
 }