amir sampler working with subgrad condition false as well

Author: Jelena Markovic

selectiveInference/R/funs.randomized.R

@@ -374,14 +374,23 @@ randomizedLassoInf = function(X,
   }
   inactive_set = lasso_soln$inactive_set
 
-
   noise_scale = lasso_soln$noise_scale # set to default value in randomizedLasso
-
- if (condition_subgrad==TRUE){
+  
+  constraints = matrix(0,nactive,2)
+  constraints[,2] = Inf
+  if (condition_subgrad==TRUE){
    condl_lasso=conditional_density(noise_scale, lasso_soln)
    lasso_soln = condl_lasso$lasso_soln
-   reduced_opt_transform = condl_lasso$reduced_opt_transform
- } 
+   cur_opt_transform = condl_lasso$reduced_opt_transform
+   } else{
+   if (nactive<p){
+     subgrad_constraints = matrix(-lam, p-nactive, 2)
+     subgrad_constraints[,2]=lam
+     constraints = rbind(constraints, subgrad_constraints)
+   }
+    cur_opt_transform = list(linear_term = lasso_soln$optimization_transform$linear_term,
+                             offset_term = lasso_soln$optimization_transform$offset_term+lasso_soln$observed_raw)
+  }
 
   ndim = length(lasso_soln$observed_opt_state)
 
@@ -391,8 +400,9 @@ randomizedLassoInf = function(X,
   } else if (sampler == "A"){
     opt_samples = gaussian_sampler(noise_scale, 
                                  lasso_soln$observed_opt_state, 
-                                 reduced_opt_transform$linear_term,
-                                 reduced_opt_transform$offset_term,
+                                 cur_opt_transform$linear_term,
+                                 cur_opt_transform$offset_term,
+                                 constraints,
                                  nsamples=nsample)
     opt_sample = opt_samples[(burnin+1):nsample,]
   }
@@ -439,7 +449,6 @@ randomizedLassoInf = function(X,
         cur_linear = reduced_target_opt_linear[,2:ncol(reduced_target_opt_linear)]
         cur_offset = temp %*% opt_transform$offset_term
         cur_transform = list(linear_term = as.matrix(cur_linear), offset_term = cur_offset)
-
         raw = target_transform$linear_term * observed_target[i] + target_transform$offset_term
     } else {
         cur_transform = opt_transform

selectiveInference/R/sampler.R

@@ -74,7 +74,12 @@ log_concave_sampler = function(negative_log_density,
   return (samples)
 }
 
-gaussian_sampler = function(noise_scale, observed, linear_term, offset_term, nsamples){
+gaussian_sampler = function(noise_scale, 
+                            observed, 
+                            linear_term, 
+                            offset_term, 
+                            constraints,
+                            nsamples){
 
   negative_log_density = function(x) {
     recon = linear_term %*% x+offset_term
@@ -84,9 +89,6 @@ gaussian_sampler = function(noise_scale, observed, linear_term, offset_term, nsa
     recon = linear_term %*% x+offset_term
     return(t(linear_term)%*% recon/(noise_scale^2))
   }
-  dim = length(observed)
-  constraints = matrix(0,dim,2)
-  constraints[,2] = Inf
 
   return(log_concave_sampler(negative_log_density, 
                              grad_negative_log_density,

tests/randomized/test_instances.R

@@ -61,7 +61,7 @@ test_KKT=function(){
 
 
 
-collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=TRUE, lam=1.2){
+collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=FALSE, lam=1.2){
 
   rho=0.3
   sigma=1
@@ -76,7 +76,7 @@ collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=TRUE, l
                                                      lam=lam, 
                                                      sigma=sigma,
                                                      level=level, 
-                                                     sampler = "R",
+                                                     sampler = "A",
                                                      burnin=1000, 
                                                      nsample=5000, 
                                                      condition_subgrad=condition_subgrad)
@@ -104,7 +104,7 @@ collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=TRUE, l
 }
 
 set.seed(1)
-collect_results(n=100, p=2000, s=0, lam=2.5)
+collect_results(n=100, p=20, s=0, lam=1.2)
 #test_randomized_lasso()
 #test_KKT()