Skip to content

Commit c940537

Browse files
jonathan-taylorjonathan-taylor
committed
added laplace densities
1 parent c175a8b commit c940537

File tree

2 files changed

+145
-13
lines changed

2 files changed

+145
-13
lines changed

‎selectiveInference/src/Rcpp-randomized.cpp

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,3 +63,66 @@ Rcpp::NumericVector log_density_gaussian_conditional_(double noise_scale,
6363

6464
return(result);
6565
}
66+
67+
// [[Rcpp::export]]
68+
Rcpp::NumericVector log_density_laplace_(double noise_scale, // Scale of randomization
69+
Rcpp::NumericMatrix internal_linear, // A_D -- linear part for data
70+
Rcpp::NumericMatrix internal_state, // D -- data state -- matrix of shape (nopt, npts)
71+
Rcpp::NumericMatrix optimization_linear, // A_O -- linear part for optimization variables
72+
Rcpp::NumericMatrix optimization_state, // O -- optimization state -- matrix of shape (ninternal, npts)
73+
Rcpp::NumericVector offset) { // h -- offset in affine transform -- "p" dimensional
74+
75+
int npt = internal_state.ncol(); // Function is vectorized
76+
if (optimization_state.ncol() != npt) { // Assuming each column is an internal or opt state because arrays are column major
77+
Rcpp::stop("Number of optimization samples should equal the number of (internally represented) data.");
78+
}
79+
80+
int ndim = optimization_linear.nrow();
81+
if (internal_linear.nrow() != ndim) {
82+
Rcpp::stop("Dimension of optimization range should be the same as the dimension of the data range.");
83+
}
84+
int ninternal = internal_linear.ncol();
85+
int noptimization = optimization_linear.ncol();
86+
87+
Rcpp::NumericVector result(npt);
88+
89+
int ipt;
90+
for (ipt=0; ipt<npt; ipt++) {
91+
result[ipt] = log_density_laplace(noise_scale,
92+
ndim,
93+
ninternal,
94+
noptimization,
95+
(double *) internal_linear.begin(),
96+
((double *) internal_state.begin() + ipt * ninternal),
97+
(double *) optimization_linear.begin(),
98+
((double *) optimization_state.begin() + ipt * noptimization),
99+
(double *) offset.begin());
100+
}
101+
102+
return(result);
103+
}
104+
105+
// [[Rcpp::export]]
106+
Rcpp::NumericVector log_density_laplace_conditional_(double noise_scale, // Scale of randomization
107+
Rcpp::NumericMatrix optimization_linear, // A_O -- linear part for optimization variables
108+
Rcpp::NumericMatrix optimization_state, // O -- optimization state -- matrix of shape (ninternal, npts)
109+
Rcpp::NumericVector offset) { // h -- offset in affine transform -- "p" dimensional
110+
111+
int npt = optimization_state.ncol(); // Function is vectorized
112+
int ndim = optimization_linear.nrow();
113+
int noptimization = optimization_linear.ncol();
114+
115+
Rcpp::NumericVector result(npt);
116+
117+
int ipt;
118+
for (ipt=0; ipt<npt; ipt++) {
119+
result[ipt] = log_density_laplace_conditional(noise_scale,
120+
ndim,
121+
noptimization,
122+
(double *) optimization_linear.begin(),
123+
((double *) optimization_state.begin() + ipt * noptimization),
124+
(double *) offset.begin());
125+
}
126+
127+
return(result);
128+
}

‎tests/randomized/test_randomized.R

Lines changed: 82 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -12,25 +12,94 @@ smoke_test = function() {
1212
}
1313
A = smoke_test()
1414

15-
density_test = function() {
15+
gaussian_density_test = function() {
1616

17+
noise_scale = 10.
1718
random_lasso = smoke_test()
1819
p = nrow(random_lasso$internal_transform$linear_term)
1920
internal_state = matrix(rnorm(p * 20), p, 20)
2021
optimization_state = matrix(rnorm(p * 20), p, 20)
2122
offset = rnorm(p)
2223

23-
selectiveInference:::log_density_gaussian_(10.,
24-
random_lasso$internal_transform$linear_term,
25-
internal_state,
26-
random_lasso$optimization_transform$linear_term,
27-
optimization_state,
28-
offset)
29-
30-
selectiveInference:::log_density_gaussian_conditional_(10.,
31-
random_lasso$optimization_transform$linear_term,
32-
optimization_state,
33-
offset)
24+
V1 = selectiveInference:::log_density_gaussian_(noise_scale,
25+
random_lasso$internal_transform$linear_term,
26+
internal_state,
27+
random_lasso$optimization_transform$linear_term,
28+
optimization_state,
29+
offset)
30+
A1 = random_lasso$internal_transform$linear_term
31+
A2 = random_lasso$optimization_transform$linear_term
32+
arg = A1 %*% internal_state + A2 %*% optimization_state + offset
33+
V2 = -apply(arg^2, 2, sum) / (2 * noise_scale^2)
34+
print(sqrt(sum((V1-V2)^2) / sum(V1^2)))
35+
36+
U1 = selectiveInference:::log_density_gaussian_conditional_(noise_scale,
37+
random_lasso$optimization_transform$linear_term,
38+
optimization_state,
39+
offset)
40+
arg = A2 %*% optimization_state + offset
41+
U2 = -apply(arg^2, 2, sum) / (2 * noise_scale^2)
42+
print(sqrt(sum((U1-U2)^2) / sum(U1^2)))
43+
44+
# test that a single column matrix works -- numeric should not
45+
46+
print(selectiveInference:::log_density_gaussian_conditional_(noise_scale,
47+
random_lasso$optimization_transform$linear_term,
48+
optimization_state[,1,drop=FALSE],
49+
offset))
50+
print(selectiveInference:::log_density_gaussian_(noise_scale,
51+
random_lasso$internal_transform$linear_term,
52+
internal_state[,1,drop=FALSE],
53+
random_lasso$optimization_transform$linear_term,
54+
optimization_state[,1,drop=FALSE],
55+
offset))
56+
57+
}
58+
59+
gaussian_density_test()
60+
61+
laplace_density_test = function() {
62+
63+
noise_scale = 10.
64+
random_lasso = smoke_test()
65+
p = nrow(random_lasso$internal_transform$linear_term)
66+
internal_state = matrix(rnorm(p * 20), p, 20)
67+
optimization_state = matrix(rnorm(p * 20), p, 20)
68+
offset = rnorm(p)
69+
70+
V1 = selectiveInference:::log_density_laplace_(noise_scale,
71+
random_lasso$internal_transform$linear_term,
72+
internal_state,
73+
random_lasso$optimization_transform$linear_term,
74+
optimization_state,
75+
offset)
76+
A1 = random_lasso$internal_transform$linear_term
77+
A2 = random_lasso$optimization_transform$linear_term
78+
arg = A1 %*% internal_state + A2 %*% optimization_state + offset
79+
V2 = -apply(abs(arg), 2, sum) / noise_scale
80+
print(sqrt(sum((V1-V2)^2) / sum(V1^2)))
81+
82+
U1 = selectiveInference:::log_density_laplace_conditional_(noise_scale,
83+
random_lasso$optimization_transform$linear_term,
84+
optimization_state,
85+
offset)
86+
arg = A2 %*% optimization_state + offset
87+
U2 = -apply(abs(arg), 2, sum) / noise_scale
88+
print(sqrt(sum((U1-U2)^2) / sum(U1^2)))
89+
90+
# test that a single column matrix works -- numeric should not
91+
92+
print(selectiveInference:::log_density_laplace_conditional_(noise_scale,
93+
random_lasso$optimization_transform$linear_term,
94+
optimization_state[,1,drop=FALSE],
95+
offset))
96+
print(selectiveInference:::log_density_laplace_(noise_scale,
97+
random_lasso$internal_transform$linear_term,
98+
internal_state[,1,drop=FALSE],
99+
random_lasso$optimization_transform$linear_term,
100+
optimization_state[,1,drop=FALSE],
101+
offset))
102+
34103
}
35104

36-
density_test()
105+
laplace_density_test()

0 commit comments

Comments
 (0)