Skip to content

Commit 2254e03

Browse files
tylermorganwalltylermorganwall
committed
-Fix global variables warning
-Split cpp files into multiple parts to avoid large object compilation error on windows
1 parent 25ccea3 commit 2254e03

File tree

7 files changed

+149
-140
lines changed

7 files changed

+149
-140
lines changed

‎R/plot_correlations.R‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -198,3 +198,4 @@ plot_correlations = function(
198198
rownames(results) = colnames(mm)[-1]
199199
invisible(results)
200200
}
201+
globalVariables(c("x", "y"))

‎R/plot_fds.R‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -441,3 +441,4 @@ plot_fds = function(
441441
}
442442
return(varsorderedscaled)
443443
}
444+
globalVariables(c("variance"))

‎src/blocked_optimality_functions.cpp‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
2+
#include <RcppEigen.h>
3+
4+
//**********************************************************
5+
//Everything below is for generating blocked optimal designs
6+
//**********************************************************
7+
8+
double calculateBlockedDOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& gls) {
9+
return((currentDesign.transpose()*gls*currentDesign).partialPivLu().determinant());
10+
}
11+
12+
double calculateBlockedDOptimalityLog(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& gls) {
13+
Eigen::MatrixXd XtX = (currentDesign.transpose()*gls*currentDesign);
14+
return(exp(XtX.llt().matrixL().toDenseMatrix().diagonal().array().log().sum()));
15+
}
16+
17+
double calculateBlockedIOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& momentsMatrix,const Eigen::MatrixXd& gls) {
18+
return(((currentDesign.transpose()*gls*currentDesign).llt().solve(momentsMatrix)).trace());
19+
}
20+
21+
double calculateBlockedAOptimality(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
22+
return((currentDesign.transpose()*gls*currentDesign).partialPivLu().inverse().trace());
23+
}
24+
25+
double calculateBlockedAliasTrace(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& aliasMatrix,const Eigen::MatrixXd& gls) {
26+
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
27+
Eigen::MatrixXd A = XtX.llt().solve(currentDesign.transpose()*aliasMatrix);
28+
return((A.transpose() * A).trace());
29+
}
30+
31+
double calculateBlockedGOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& gls) {
32+
Eigen::MatrixXd results = currentDesign*(currentDesign.transpose()*gls*currentDesign).partialPivLu().solve(currentDesign.transpose())*gls;
33+
return(results.diagonal().maxCoeff());
34+
}
35+
36+
double calculateBlockedTOptimality(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
37+
return((currentDesign.transpose()*gls*currentDesign).trace());
38+
}
39+
40+
double calculateBlockedEOptimality(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
41+
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
42+
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(XtX);
43+
return(eigensolver.eigenvalues().minCoeff());
44+
}
45+
46+
double calculateBlockedDEff(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
47+
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
48+
return(pow(XtX.partialPivLu().determinant(), 1/currentDesign.cols()) / currentDesign.rows());
49+
}
50+
51+
double calculateBlockedDEffNN(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
52+
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
53+
return(pow(XtX.partialPivLu().determinant(), 1.0/currentDesign.cols()));
54+
}
55+
56+
double calculateBlockedAliasTracePseudoInv(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& aliasMatrix,const Eigen::MatrixXd& gls) {
57+
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
58+
Eigen::MatrixXd A = XtX.partialPivLu().solve(currentDesign.transpose()*aliasMatrix);
59+
return((A.transpose() * A).trace());
60+
}
61+
62+
bool isSingularBlocked(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
63+
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
64+
return(!XtX.colPivHouseholderQr().isInvertible());
65+
}
66+
67+
double calculateBlockedCustomOptimality(const Eigen::MatrixXd& currentDesign, Rcpp::Function customBlockedOpt, const Eigen::MatrixXd& gls) {
68+
return Rcpp::as<double>(customBlockedOpt(Rcpp::Named("currentDesign", currentDesign),Rcpp::Named("vInv", gls)));
69+
}

‎src/exported_optimality.cpp‎

Lines changed: 0 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -40,42 +40,5 @@ double AOptimality(const Eigen::MatrixXd& currentDesign) {
4040
return(100 * currentDesign.cols() / ((currentDesign.rows() * (XtX.partialPivLu().inverse())).trace()));
4141
}
4242

43-
// [[Rcpp::export]]
44-
double calculateAOptimalityPseudo(const Eigen::MatrixXd& currentDesign) {
45-
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
46-
return(XtX.completeOrthogonalDecomposition().pseudoInverse().trace());
47-
}
48-
49-
// [[Rcpp::export]]
50-
double IOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& momentsMatrix, const Eigen::MatrixXd& blockedVar) {
51-
Eigen::MatrixXd XtX = (currentDesign.transpose() * (blockedVar.householderQr().solve(currentDesign))).householderQr().solve(momentsMatrix);
52-
return(XtX.trace());
53-
}
54-
55-
56-
// [[Rcpp::export]]
57-
double calcAliasTrace(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& aliasMatrix) {
58-
Eigen::MatrixXd XtX = currentDesign.transpose() * currentDesign;
59-
Eigen::MatrixXd A = XtX.llt().solve(currentDesign.transpose() * aliasMatrix);
60-
return((A.transpose() * A).trace());
61-
}
62-
63-
// [[Rcpp::export]]
64-
Eigen::MatrixXd covarianceMatrixPseudo(const Eigen::MatrixXd& design) {
65-
Eigen::MatrixXd XtX = design.transpose() * design;
66-
return(XtX.completeOrthogonalDecomposition().pseudoInverse());
67-
}
68-
69-
// [[Rcpp::export]]
70-
Eigen::MatrixXd getPseudoInverse(const Eigen::MatrixXd& currentDesign) {
71-
return(currentDesign.completeOrthogonalDecomposition().pseudoInverse());
72-
}
73-
74-
// [[Rcpp::export]]
75-
double GEfficiency(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& candset) {
76-
Eigen::MatrixXd V = (currentDesign.transpose()*currentDesign).partialPivLu().inverse();
77-
Eigen::MatrixXd results = candset*V*candset.transpose();
78-
return(100*((double)currentDesign.cols()/(double)currentDesign.rows())*1/results.diagonal().maxCoeff());
79-
}
8043

8144

‎src/exported_optimality_2.cpp‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
// [[Rcpp::depends(RcppEigen)]]
2+
3+
#include <RcppEigen.h>
4+
using namespace Rcpp;
5+
6+
// [[Rcpp::export]]
7+
double calculateAOptimalityPseudo(const Eigen::MatrixXd& currentDesign) {
8+
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
9+
return(XtX.completeOrthogonalDecomposition().pseudoInverse().trace());
10+
}
11+
12+
// [[Rcpp::export]]
13+
double IOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& momentsMatrix, const Eigen::MatrixXd& blockedVar) {
14+
Eigen::MatrixXd XtX = (currentDesign.transpose() * (blockedVar.householderQr().solve(currentDesign))).householderQr().solve(momentsMatrix);
15+
return(XtX.trace());
16+
}
17+
18+
19+
// [[Rcpp::export]]
20+
double calcAliasTrace(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& aliasMatrix) {
21+
Eigen::MatrixXd XtX = currentDesign.transpose() * currentDesign;
22+
Eigen::MatrixXd A = XtX.llt().solve(currentDesign.transpose() * aliasMatrix);
23+
return((A.transpose() * A).trace());
24+
}
25+
26+
// [[Rcpp::export]]
27+
Eigen::MatrixXd covarianceMatrixPseudo(const Eigen::MatrixXd& design) {
28+
Eigen::MatrixXd XtX = design.transpose() * design;
29+
return(XtX.completeOrthogonalDecomposition().pseudoInverse());
30+
}
31+
32+
// [[Rcpp::export]]
33+
Eigen::MatrixXd getPseudoInverse(const Eigen::MatrixXd& currentDesign) {
34+
return(currentDesign.completeOrthogonalDecomposition().pseudoInverse());
35+
}
36+
37+
// [[Rcpp::export]]
38+
double GEfficiency(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& candset) {
39+
Eigen::MatrixXd V = (currentDesign.transpose()*currentDesign).partialPivLu().inverse();
40+
Eigen::MatrixXd results = candset*V*candset.transpose();
41+
return(100*((double)currentDesign.cols()/(double)currentDesign.rows())*1/results.diagonal().maxCoeff());
42+
}

‎src/optimalityfunctions.cpp‎

Lines changed: 0 additions & 103 deletions
Original file line numberDiff line numberDiff line change
@@ -1,36 +1,5 @@
11
#include <RcppEigen.h>
22

3-
double calculateDOptimality(const Eigen::MatrixXd& currentDesign) {
4-
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
5-
return(XtX.partialPivLu().determinant()); //works without partialPivLu()
6-
}
7-
8-
double calculateDOptimalityLog(const Eigen::MatrixXd& currentDesign) {
9-
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
10-
return(XtX.llt().matrixL().toDenseMatrix().diagonal().array().log().sum());
11-
}
12-
13-
double calculateIOptimality(const Eigen::MatrixXd& currentV, const Eigen::MatrixXd& momentsMatrix) {
14-
return((currentV * momentsMatrix).trace());
15-
}
16-
17-
18-
double calculateGOptimality(const Eigen::MatrixXd& currentV, const Eigen::MatrixXd& currentDesign) {
19-
Eigen::MatrixXd results = currentDesign*currentV*currentDesign.transpose();
20-
return(results.diagonal().maxCoeff());
21-
}
22-
23-
double calculateTOptimality(const Eigen::MatrixXd& currentDesign) {
24-
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
25-
return(XtX.trace());
26-
}
27-
28-
double calculateEOptimality(const Eigen::MatrixXd& currentDesign) {
29-
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
30-
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(XtX);
31-
return(eigensolver.eigenvalues().minCoeff());
32-
}
33-
343
double calculateAOptimality(const Eigen::MatrixXd& currentV) {
354
return(currentV.trace());
365
}
@@ -53,10 +22,6 @@ double calculateDEff(const Eigen::MatrixXd& currentDesign, double numbercols, do
5322
return(pow(XtX.partialPivLu().determinant(), 1/numbercols) / numberrows);
5423
}
5524

56-
double calculateDEffLog(const Eigen::MatrixXd& currentDesign, double numbercols, double numberrows) {
57-
return(exp(calculateDOptimalityLog(currentDesign)/numbercols) / numberrows);
58-
}
59-
6025
double calculateDEffNN(const Eigen::MatrixXd& currentDesign, double numbercols) {
6126
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
6227
return(pow(XtX.partialPivLu().determinant(), 1.0/numbercols));
@@ -107,71 +72,3 @@ void search_candidate_set(const Eigen::MatrixXd& V, const Eigen::MatrixXd& candi
10772
}
10873
}
10974
}
110-
111-
//**********************************************************
112-
//Everything below is for generating blocked optimal designs
113-
//**********************************************************
114-
115-
double calculateBlockedDOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& gls) {
116-
return((currentDesign.transpose()*gls*currentDesign).partialPivLu().determinant());
117-
}
118-
119-
double calculateBlockedDOptimalityLog(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& gls) {
120-
Eigen::MatrixXd XtX = (currentDesign.transpose()*gls*currentDesign);
121-
return(exp(XtX.llt().matrixL().toDenseMatrix().diagonal().array().log().sum()));
122-
}
123-
124-
double calculateBlockedIOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& momentsMatrix,const Eigen::MatrixXd& gls) {
125-
return(((currentDesign.transpose()*gls*currentDesign).llt().solve(momentsMatrix)).trace());
126-
}
127-
128-
double calculateBlockedAOptimality(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
129-
return((currentDesign.transpose()*gls*currentDesign).partialPivLu().inverse().trace());
130-
}
131-
132-
double calculateBlockedAliasTrace(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& aliasMatrix,const Eigen::MatrixXd& gls) {
133-
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
134-
Eigen::MatrixXd A = XtX.llt().solve(currentDesign.transpose()*aliasMatrix);
135-
return((A.transpose() * A).trace());
136-
}
137-
138-
double calculateBlockedGOptimality(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& gls) {
139-
Eigen::MatrixXd results = currentDesign*(currentDesign.transpose()*gls*currentDesign).partialPivLu().solve(currentDesign.transpose())*gls;
140-
return(results.diagonal().maxCoeff());
141-
}
142-
143-
double calculateBlockedTOptimality(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
144-
return((currentDesign.transpose()*gls*currentDesign).trace());
145-
}
146-
147-
double calculateBlockedEOptimality(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
148-
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
149-
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(XtX);
150-
return(eigensolver.eigenvalues().minCoeff());
151-
}
152-
153-
double calculateBlockedDEff(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
154-
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
155-
return(pow(XtX.partialPivLu().determinant(), 1/currentDesign.cols()) / currentDesign.rows());
156-
}
157-
158-
double calculateBlockedDEffNN(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
159-
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
160-
return(pow(XtX.partialPivLu().determinant(), 1.0/currentDesign.cols()));
161-
}
162-
163-
double calculateBlockedAliasTracePseudoInv(const Eigen::MatrixXd& currentDesign, const Eigen::MatrixXd& aliasMatrix,const Eigen::MatrixXd& gls) {
164-
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
165-
Eigen::MatrixXd A = XtX.partialPivLu().solve(currentDesign.transpose()*aliasMatrix);
166-
return((A.transpose() * A).trace());
167-
}
168-
169-
bool isSingularBlocked(const Eigen::MatrixXd& currentDesign,const Eigen::MatrixXd& gls) {
170-
Eigen::MatrixXd XtX = currentDesign.transpose()*gls*currentDesign;
171-
return(!XtX.colPivHouseholderQr().isInvertible());
172-
}
173-
174-
double calculateBlockedCustomOptimality(const Eigen::MatrixXd& currentDesign, Rcpp::Function customBlockedOpt, const Eigen::MatrixXd& gls) {
175-
return Rcpp::as<double>(customBlockedOpt(Rcpp::Named("currentDesign", currentDesign),Rcpp::Named("vInv", gls)));
176-
}
177-

‎src/optimalityfunctions2.cpp‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
#include <RcppEigen.h>
2+
3+
double calculateDOptimality(const Eigen::MatrixXd& currentDesign) {
4+
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
5+
return(XtX.partialPivLu().determinant()); //works without partialPivLu()
6+
}
7+
8+
double calculateDOptimalityLog(const Eigen::MatrixXd& currentDesign) {
9+
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
10+
return(XtX.llt().matrixL().toDenseMatrix().diagonal().array().log().sum());
11+
}
12+
13+
double calculateIOptimality(const Eigen::MatrixXd& currentV, const Eigen::MatrixXd& momentsMatrix) {
14+
return((currentV * momentsMatrix).trace());
15+
}
16+
17+
18+
double calculateGOptimality(const Eigen::MatrixXd& currentV, const Eigen::MatrixXd& currentDesign) {
19+
Eigen::MatrixXd results = currentDesign*currentV*currentDesign.transpose();
20+
return(results.diagonal().maxCoeff());
21+
}
22+
23+
double calculateTOptimality(const Eigen::MatrixXd& currentDesign) {
24+
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
25+
return(XtX.trace());
26+
}
27+
28+
double calculateEOptimality(const Eigen::MatrixXd& currentDesign) {
29+
Eigen::MatrixXd XtX = currentDesign.transpose()*currentDesign;
30+
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(XtX);
31+
return(eigensolver.eigenvalues().minCoeff());
32+
}
33+
34+
double calculateDEffLog(const Eigen::MatrixXd& currentDesign, double numbercols, double numberrows) {
35+
return(exp(calculateDOptimalityLog(currentDesign)/numbercols) / numberrows);
36+
}

0 commit comments

Comments
 (0)