Skip to content

Commit 334e007

Browse files
jonathan-taylorjonathan-taylor
authored
Merge pull request #32 from jonathan-taylor/master
small cleanup of fixedLassoPoly, travis script got travis tests passing again
2 parents 782290e + eb2a5e4 commit 334e007

File tree

5 files changed

+73
-73
lines changed

5 files changed

+73
-73
lines changed

‎Makefile

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ Rcpp:
44
Rscript -e "library(Rcpp); Rcpp::compileAttributes('selectiveInference')"
55

66
install: Rcpp
7-
R CMD install selectiveInference
7+
R CMD INSTALL selectiveInference
88

99
build:
1010
R CMD build selectiveInference

‎selectiveInference/R/RcppExports.R

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,19 @@
1-
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
1+
# This file was generated by Rcpp::compileAttributes
22
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
33

44
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)
66
}
77

88
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)
1010
}
1111

1212
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)
1414
}
1515

1616
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)
1818
}
1919

‎selectiveInference/R/funs.fixed.R

Lines changed: 45 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -80,12 +80,9 @@ fixedLassoInf <- function(x, y, beta,
8080
warning(paste("Solution beta does not satisfy the KKT conditions",
8181
"(to within specified tolerances)"))
8282

83-
tol.coef = tol.beta * sqrt(n^2 / colSums(x^2))
84-
# print(tol.coef)
85-
vars = which(abs(beta) > tol.coef)
86-
# vars = abs(beta) > tol.coef
87-
# print(beta)
88-
# print(vars)
83+
tol.coef = tol.beta * sqrt(n / colSums(x^2))
84+
vars = which(abs(beta) > tol.coef)
85+
8986
if(sum(vars)==0){
9087
cat("Empty model",fill=T)
9188
return()
@@ -97,10 +94,17 @@ fixedLassoInf <- function(x, y, beta,
9794
"'thresh' parameter, for a more accurate convergence."))
9895

9996
# Get lasso polyhedral region, of form Gy >= u
100-
logical.vars=rep(FALSE,p)
101-
logical.vars[vars]=TRUE
102-
if (type == 'full') out = fixedLassoPoly(x,y,lambda,beta,logical.vars,inactive=TRUE)
103-
else out = fixedLassoPoly(x,y,lambda,beta,logical.vars)
97+
98+
logical.vars=rep(FALSE,p)
99+
logical.vars[vars]=TRUE
100+
101+
if (type == 'full') {
102+
out = fixedLassoPoly(x, y, lambda, beta, logical.vars, inactive=TRUE)
103+
}
104+
else {
105+
out = fixedLassoPoly(x, y, lambda, beta, logical.vars)
106+
}
107+
104108
A = out$A
105109
b = out$b
106110

@@ -233,34 +237,43 @@ logical.vars[vars]=TRUE
233237

234238
fixedLassoPoly =
235239
function(X, y, lambda, beta, active, inactive = FALSE) {
236-
Xa = X[,active,drop=F]
237-
Xac = X[,!active,drop=F]
238-
Xai = pinv(crossprod(Xa))
239-
Xap = Xai %*% t(Xa)
240-
241-
za = sign(beta[active])
242-
if (length(za)>1) dz = diag(za)
243-
if (length(za)==1) dz = matrix(za,1,1)
240+
241+
XA = X[, active, drop=FALSE]
242+
XI = X[, !active, drop=FALSE]
243+
XAi = pinv(crossprod(XA))
244+
XAp = XAi %*% t(XA)
245+
Ir = t(XI) %*% t(XAp) # matrix in the "irrepresentable" condition
246+
247+
if(length(lambda)>1) {
248+
lambdaA= lambda[active]
249+
lambdaI = lambda[!active]
250+
} else {
251+
lambdaA = rep(lambda, sum(active))
252+
lambdaI = rep(lambda, sum(!active))
253+
}
254+
255+
penalized = lambdaA != 0
256+
signA = sign(beta[active])
257+
active_subgrad = signA * lambdaA
258+
if (length(signA)>1) sign_diag = diag(signA)
259+
if (length(signA)==1) sign_diag = matrix(signA, 1, 1)
244260

245261
if (inactive) { # should we include the inactive constraints?
246-
R = diag(1,nrow(Xa)) - Xa %*% Xap # R is residual forming matrix of selected model
262+
RA = diag(rep(1, nrow(XA))) - XA %*% XAp # RA is residual forming matrix of selected model
247263

248264
A = rbind(
249-
1/lambda * t(Xac) %*% R,
250-
-1/lambda * t(Xac) %*% R,
251-
-dz %*% Xap
265+
t(XI) %*% RA,
266+
-t(XI) %*% RA,
267+
-(sign_diag %*% XAp)[penalized,] # no constraints for unpenalized
252268
)
253-
lambda2=lambda
254-
if(length(lambda)>1) lambda2=lambda[active]
269+
255270
b = c(
256-
1 - t(Xac) %*% t(Xap) %*% za,
257-
1 + t(Xac) %*% t(Xap) %*% za,
258-
-lambda2 * dz %*% Xai %*% za)
271+
lambdaI - Ir %*% active_subgrad,
272+
lambdaI + Ir %*% active_subgrad,
273+
-(sign_diag %*% XAi %*% active_subgrad)[penalized])
259274
} else {
260-
A = -dz %*% Xap
261-
lambda2=lambda
262-
if(length(lambda)>1) lambda2=lambda[active]
263-
b = -lambda2 * dz %*% Xai %*% za
275+
A = -(sign_diag %*% XAp)[penalized,] # no constraints for unpenalized
276+
b = -(sign_diag %*% XAi %*% active_subgrad)[penalized]
264277
}
265278

266279
return(list(A=A, b=b))
@@ -362,7 +375,6 @@ debiasingRow = function (Xinfo, # could be X or t(X) %*% X / n dep
362375
# Initialize variables
363376

364377
soln = rep(0, p)
365-
Xsoln = rep(0, n)
366378
ever_active = rep(0, p)
367379
ever_active[1] = row # 1-based
368380
ever_active = as.integer(ever_active)
@@ -393,6 +405,7 @@ debiasingRow = function (Xinfo, # could be X or t(X) %*% X / n dep
393405
objective_tol,
394406
max_active)
395407
} else {
408+
Xsoln = rep(0, nrow(Xinfo))
396409
result = solve_QP_wide(Xinfo, # this is a design matrix
397410
mu,
398411
max_iter,

‎selectiveInference/man/debiasingMatrix.Rd

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -108,7 +108,7 @@ set.seed(10)
108108
n = 50
109109
p = 100
110110
X = matrix(rnorm(n * p), n, p)
111-
S = t(X) %*% X / n
111+
S = t(X) \%*\% X / n
112112
M = debiasingMatrix(S, FALSE, n, c(1,3,5))
113113
M2 = debiasingMatrix(X, TRUE, n, c(1,3,5))
114114
max(M - M2)

‎selectiveInference/src/RcppExports.cpp

Lines changed: 21 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Generated by using Rcpp::compileAttributes() -> do not edit by hand
1+
// This file was generated by Rcpp::compileAttributes
22
// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
33

44
#include <Rcpp.h>
@@ -7,10 +7,10 @@ using namespace Rcpp;
77

88
// solve_QP
99
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);
10-
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) {
10+
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) {
1111
BEGIN_RCPP
12-
Rcpp::RObject rcpp_result_gen;
13-
Rcpp::RNGScope rcpp_rngScope_gen;
12+
Rcpp::RObject __result;
13+
Rcpp::RNGScope __rngScope;
1414
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Sigma(SigmaSEXP);
1515
Rcpp::traits::input_parameter< double >::type bound(boundSEXP);
1616
Rcpp::traits::input_parameter< int >::type maxiter(maxiterSEXP);
@@ -22,16 +22,16 @@ BEGIN_RCPP
2222
Rcpp::traits::input_parameter< double >::type kkt_tol(kkt_tolSEXP);
2323
Rcpp::traits::input_parameter< double >::type objective_tol(objective_tolSEXP);
2424
Rcpp::traits::input_parameter< int >::type max_active(max_activeSEXP);
25-
rcpp_result_gen = Rcpp::wrap(solve_QP(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active));
26-
return rcpp_result_gen;
25+
__result = Rcpp::wrap(solve_QP(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active));
26+
return __result;
2727
END_RCPP
2828
}
2929
// solve_QP_wide
3030
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);
31-
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) {
31+
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) {
3232
BEGIN_RCPP
33-
Rcpp::RObject rcpp_result_gen;
34-
Rcpp::RNGScope rcpp_rngScope_gen;
33+
Rcpp::RObject __result;
34+
Rcpp::RNGScope __rngScope;
3535
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type X(XSEXP);
3636
Rcpp::traits::input_parameter< double >::type bound(boundSEXP);
3737
Rcpp::traits::input_parameter< int >::type maxiter(maxiterSEXP);
@@ -44,49 +44,36 @@ BEGIN_RCPP
4444
Rcpp::traits::input_parameter< double >::type kkt_tol(kkt_tolSEXP);
4545
Rcpp::traits::input_parameter< double >::type objective_tol(objective_tolSEXP);
4646
Rcpp::traits::input_parameter< int >::type max_active(max_activeSEXP);
47-
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));
48-
return rcpp_result_gen;
47+
__result = Rcpp::wrap(solve_QP_wide(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active));
48+
return __result;
4949
END_RCPP
5050
}
5151
// update1_
5252
Rcpp::List update1_(Rcpp::NumericMatrix Q2, Rcpp::NumericVector w, int m, int k);
53-
RcppExport SEXP _selectiveInference_update1_(SEXP Q2SEXP, SEXP wSEXP, SEXP mSEXP, SEXP kSEXP) {
53+
RcppExport SEXP selectiveInference_update1_(SEXP Q2SEXP, SEXP wSEXP, SEXP mSEXP, SEXP kSEXP) {
5454
BEGIN_RCPP
55-
Rcpp::RObject rcpp_result_gen;
56-
Rcpp::RNGScope rcpp_rngScope_gen;
55+
Rcpp::RObject __result;
56+
Rcpp::RNGScope __rngScope;
5757
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Q2(Q2SEXP);
5858
Rcpp::traits::input_parameter< Rcpp::NumericVector >::type w(wSEXP);
5959
Rcpp::traits::input_parameter< int >::type m(mSEXP);
6060
Rcpp::traits::input_parameter< int >::type k(kSEXP);
61-
rcpp_result_gen = Rcpp::wrap(update1_(Q2, w, m, k));
62-
return rcpp_result_gen;
61+
__result = Rcpp::wrap(update1_(Q2, w, m, k));
62+
return __result;
6363
END_RCPP
6464
}
6565
// downdate1_
6666
Rcpp::List downdate1_(Rcpp::NumericMatrix Q1, Rcpp::NumericMatrix R, int j0, int m, int n);
67-
RcppExport SEXP _selectiveInference_downdate1_(SEXP Q1SEXP, SEXP RSEXP, SEXP j0SEXP, SEXP mSEXP, SEXP nSEXP) {
67+
RcppExport SEXP selectiveInference_downdate1_(SEXP Q1SEXP, SEXP RSEXP, SEXP j0SEXP, SEXP mSEXP, SEXP nSEXP) {
6868
BEGIN_RCPP
69-
Rcpp::RObject rcpp_result_gen;
70-
Rcpp::RNGScope rcpp_rngScope_gen;
69+
Rcpp::RObject __result;
70+
Rcpp::RNGScope __rngScope;
7171
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type Q1(Q1SEXP);
7272
Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type R(RSEXP);
7373
Rcpp::traits::input_parameter< int >::type j0(j0SEXP);
7474
Rcpp::traits::input_parameter< int >::type m(mSEXP);
7575
Rcpp::traits::input_parameter< int >::type n(nSEXP);
76-
rcpp_result_gen = Rcpp::wrap(downdate1_(Q1, R, j0, m, n));
77-
return rcpp_result_gen;
76+
__result = Rcpp::wrap(downdate1_(Q1, R, j0, m, n));
77+
return __result;
7878
END_RCPP
7979
}
80-
81-
static const R_CallMethodDef CallEntries[] = {
82-
{"_selectiveInference_solve_QP", (DL_FUNC) &_selectiveInference_solve_QP, 11},
83-
{"_selectiveInference_solve_QP_wide", (DL_FUNC) &_selectiveInference_solve_QP_wide, 12},
84-
{"_selectiveInference_update1_", (DL_FUNC) &_selectiveInference_update1_, 4},
85-
{"_selectiveInference_downdate1_", (DL_FUNC) &_selectiveInference_downdate1_, 5},
86-
{NULL, NULL, 0}
87-
};
88-
89-
RcppExport void R_init_selectiveInference(DllInfo *dll) {
90-
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
91-
R_useDynamicSymbols(dll, FALSE);
92-
}

0 commit comments

Comments
 (0)