8000 speed up graflex again by suzannejin · Pull Request #50 · tpq/propr · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

speed up graflex again #50

New issue

Have a question about this project? Sign up for a free Gi 8000 tHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Sep 13, 2024
Merged

speed up graflex again #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c2e9bfa
add more tests to theta
suzannejin Sep 12, 2024
7a9cb0c
add more test graflex
suzannejin Sep 13, 2024
bb95ae6
update graflex
suzannejin Sep 13, 2024
51945ea
add more tests graflex
suzannejin Sep 13, 2024
1a06c77
update news.md
suzannejin Sep 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,8 @@
## propr 5.1.1
---------------------
* Speed up `graflex` related functions


## propr 5.1.0
---------------------
* Allowed `updateCutoffs` function to compute the FDR values for negative cutoffs
Expand Down
20 changes: 6 additions & 14 deletions R/RcppExports.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,16 +105,12 @@ ctzRcpp <- function(X) {
.Call(`_propr_ctzRcpp`, X)
}

get_lower_triangle <- function(mat) {
.Call(`_propr_get_lower_triangle`, mat)
get_triangle <- function(mat) {
.Call(`_propr_get_triangle`, mat)
}

shuffle_and_get_lower_triangle <- function(mat) {
.Call(`_propr_shuffle_and_get_lower_triangle`, mat)
}

binTab <- function(A, G) {
.Call(`_propr_binTab`, A, G)
get_triangle_from_index <- function(mat, index) {
.Call(`_propr_get_triangle_from_index`, mat, index)
}

getOR <- function(A, G) {
Expand All @@ -125,12 +121,8 @@ permuteOR <- function(A, Gstar, p = 100L) {
.Call(`_propr_permuteOR`, A, Gstar, p)
}

getFDR_over <- function(actual, permuted) {
.Call(`_propr_getFDR_over`, actual, permuted)
}

getFDR_under <- function(actual, permuted) {
.Call(`_propr_getFDR_under`, actual, permuted)
getFDR <- function(actual, permuted) {
.Call(`_propr_getFDR`, actual, permuted)
}

graflex <- function(A, G, p = 100L) {
Expand Down
81 changes: 28 additions & 53 deletions src/RcppExports.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,97 +317,74 @@ BEGIN_RCPP
return rcpp_result_gen;
END_RCPP
}
// get_lower_triangle
IntegerVector get_lower_triangle(IntegerMatrix& mat);
RcppExport SEXP _propr_get_lower_triangle(SEXP matSEXP) {
// get_triangle
IntegerVector get_triangle(const IntegerMatrix& mat);
RcppExport SEXP _propr_get_triangle(SEXP matSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< IntegerMatrix& >::type mat(matSEXP);
rcpp_result_gen = Rcpp::wrap(get_lower_triangle(mat));
Rcpp::traits::input_parameter< const IntegerMatrix& >::type mat(matSEXP);
rcpp_result_gen = Rcpp::wrap(get_triangle(mat));
return rcpp_result_gen;
END_RCPP
}
// shuffle_and_get_lower_triangle
IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat);
RcppExport SEXP _propr_shuffle_and_get_lower_triangle(SEXP matSEXP) {
// get_triangle_from_index
IntegerVector get_triangle_from_index(const IntegerMatrix& mat, const IntegerVector& index);
RcppExport SEXP _propr_get_triangle_from_index(SEXP matSEXP, SEXP indexSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< IntegerMatrix& >::type mat(matSEXP);
rcpp_result_gen = Rcpp::wrap(shuffle_and_get_lower_triangle(mat));
return rcpp_result_gen;
END_RCPP
}
// binTab
NumericMatrix binTab(IntegerVector& A, IntegerVector& G);
RcppExport SEXP _propr_binTab(SEXP ASEXP, SEXP GSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< IntegerVector& >::type A(ASEXP);
Rcpp::traits::input_parameter< IntegerVector& >::type G(GSEXP);
rcpp_result_gen = Rcpp::wrap(binTab(A, G));
Rcpp::traits::input_parameter< const IntegerMatrix& >::type mat(matSEXP);
Rcpp::traits::input_parameter< const IntegerVector& >::type index(indexSEXP);
rcpp_result_gen = Rcpp::wrap(get_triangle_from_index(mat, index));
return rcpp_result_gen;
END_RCPP
}
// getOR
NumericMatrix getOR(IntegerVector& A, IntegerVector& G);
NumericVector getOR(const IntegerVector& A, const IntegerVector& G);
RcppExport SEXP _propr_getOR(SEXP ASEXP, SEXP GSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< IntegerVector& >::type A(ASEXP);
Rcpp::traits::input_parameter< IntegerVector& >::type G(GSEXP);
Rcpp::traits::input_parameter< const IntegerVector& >::type A(ASEXP);
Rcpp::traits::input_parameter< const IntegerVector& >::type G(GSEXP);
rcpp_result_gen = Rcpp::wrap(getOR(A, G));
return rcpp_result_gen;
END_RCPP
}
// permuteOR
NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p);
NumericMatrix permuteOR(const IntegerMatrix& A, const IntegerVector& Gstar, int p);
RcppExport SEXP _propr_permuteOR(SEXP ASEXP, SEXP GstarSEXP, SEXP pSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< IntegerMatrix& >::type A(ASEXP);
Rcpp::traits::input_parameter< IntegerVector& >::type Gstar(GstarSEXP);
Rcpp::traits::input_parameter< const IntegerMatrix& >::type A(ASEXP);
Rcpp::traits::input_parameter< const IntegerVector& >::type Gstar(GstarSEXP);
Rcpp::traits::input_parameter< int >::type p(pSEXP);
rcpp_result_gen = Rcpp::wrap(permuteOR(A, Gstar, p));
return rcpp_result_gen;
END_RCPP
}
// getFDR_over
double getFDR_over(double actual, NumericVector permuted);
RcppExport SEXP _propr_getFDR_over(SEXP actualSEXP, SEXP permutedSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< double >::type actual(actualSEXP);
Rcpp::traits::input_parameter< NumericVector >::type permuted(permutedSEXP);
rcpp_result_gen = Rcpp::wrap(getFDR_over(actual, permuted));
return rcpp_result_gen;
END_RCPP
}
// getFDR_under
double getFDR_under(double actual, NumericVector permuted);
RcppExport SEXP _propr_getFDR_under(SEXP actualSEXP, SEXP permutedSEXP) {
// getFDR
List getFDR(double actual, const NumericVector& permuted);
RcppExport SEXP _propr_getFDR(SEXP actualSEXP, SEXP permutedSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< double >::type actual(actualSEXP);
Rcpp::traits::input_parameter< NumericVector >::type permuted(permutedSEXP);
rcpp_result_gen = Rcpp::wrap(getFDR_under(actual, permuted));
Rcpp::traits::input_parameter< const NumericVector& >::type permuted(permutedSEXP);
rcpp_result_gen = Rcpp::wrap(getFDR(actual, permuted));
return rcpp_result_gen;
END_RCPP
}
// graflex
NumericMatrix graflex(IntegerMatrix& A, IntegerMatrix& G, int p);
NumericVector graflex(const IntegerMatrix& A, const IntegerMatrix& G, int p);
RcppExport SEXP _propr_graflex(SEXP ASEXP, SEXP GSEXP, SEXP pSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< IntegerMatrix& >::type A(ASEXP);
Rcpp::traits::input_parameter< IntegerMatrix& >::type G(GSEXP);
Rcpp::traits::input_parameter< const IntegerMatrix& >::type A(ASEXP);
Rcpp::traits::input_parameter< const IntegerMatrix& >::type G(GSEXP);
Rcpp::traits::input_parameter< int >::type p(pSEXP);
rcpp_result_gen = Rcpp::wrap(graflex(A, G, p));
return rcpp_result_gen;
Expand Down Expand Up @@ -539,13 +516,11 @@ static const R_CallMethodDef CallEntries[] = {
{"_propr_count_less_equal_than", (DL_FUNC) &_propr_count_less_equal_than, 2},
{"_propr_count_greater_equal_than", (DL_FUNC) &_propr_count_greater_equal_than, 2},
{"_propr_ctzRcpp", (DL_FUNC) &_propr_ctzRcpp, 1},
{"_propr_get_lower_triangle", (DL_FUNC) &_propr_get_lower_triangle, 1},
{"_propr_shuffle_and_get_lower_triangle", (DL_FUNC) &_propr_shuffle_and_get_lower_triangle, 1},
{"_propr_binTab", (DL_FUNC) &_propr_binTab, 2},
{"_propr_get_triangle", (DL_FUNC) &_propr_get_triangle, 1},
{"_propr_get_triangle_from_index", (DL_FUNC) &_propr_get_triangle_from_index, 2},
{"_propr_getOR", (DL_FUNC) &_propr_getOR, 2},
{"_propr_permuteOR", (DL_FUNC) &_propr_permuteOR, 3},
{"_propr_getFDR_over", (DL_FUNC) &_propr_getFDR_over, 2},
{"_propr_getFDR_under", (DL_FUNC) &_propr_getFDR_under, 2},
{"_propr_getFDR", (DL_FUNC) &_propr_getFDR, 2},
{"_propr_graflex", (DL_FUNC) &_propr_graflex, 3},
{"_propr_lr2vlr", (DL_FUNC) &_propr_lr2vlr, 1},
{"_propr_lr2phi", (DL_FUNC) &_propr_lr2phi, 1},
Expand Down
124 changes: 56 additions & 68 deletions src/graflex.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,133 +5,121 @@ using namespace Rcpp;

// Function to extract the triangle of a square and symmetric IntegerMatrix
// [[Rcpp::export]]
IntegerVector get_lower_triangle(IntegerMatrix& mat) {
int nrow = mat.nrow();
IntegerVector get_triangle(const IntegerMatrix& mat) {
int ncol = mat.ncol();
int n = ncol * (ncol - 1) / 2;

IntegerVector triangle(n);

int k = 0;
for (int j = 0; j < ncol; ++j) {
for (int i = j+1; i < nrow; ++i) {
for (int i = j+1; i < ncol; ++i) {
triangle[k++] = mat(i, j);
}
}

return triangle;
}

// Function to shuffle a square and symmetric IntegerMatrix, and get the lower triangle
// Function to get the triangle of a matrix based on the given indeces
// [[Rcpp::export]]
IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat) {
int nrow = mat.nrow();
IntegerVector get_triangle_from_index(const IntegerMatrix& mat, const IntegerVector& index) {
int ncol = mat.ncol();
int n = ncol * (ncol - 1) / 2;

IntegerVector shuffled_triangle(n);
IntegerVector index = sample(ncol, ncol, false) - 1;
IntegerVector triangle(n);

int k = 0;
for (int i = 0; i < ncol; ++i) {
int index_i = index[i];
for (int j = 0; j < i; ++j) {
int index_j = index[j];
shuffled_triangle[k++] = mat(index_i, index_j);
for (int j = 0; j < ncol; ++j) {
for (int i = j+1; i < ncol; ++i) {
triangle[k++] = mat(index[i], index[j]);
}
}

return shuffled_triangle;
return triangle;
}

// Function to calculate the contingency table
// [[Rcpp::export]]
NumericMatrix binTab(IntegerVector& A, IntegerVector& G) {
NumericMatrix tab(2, 2);
NumericVector getOR(const IntegerVector& A, const IntegerVector& G) {
int n = A.size();

tab(0, 0) = sum((1 - A) * (1 - G)); // not in A and not in G
tab(0, 1) = sum((1 - A) * G); // not in A but in G
tab(1, 0) = sum(A * (1 - G)); // in A but not G
tab(1, 1) = sum(A * G); // in A and G

return tab;
}

// Function to calculate the odds ratio, and parse all information into a matrix
// [[Rcpp::export]]
NumericMatrix getOR(IntegerVector& A, IntegerVector& G) {

NumericMatrix tab = binTab(A, G);
double odds_ratio = (tab(0, 0) * tab(1, 1)) / (tab(0, 1) * tab(1, 0));

NumericMatrix or_table(1, 8);
or_table(0, 0) = tab(0, 0);
or_table(0, 1) = tab(0, 1);
or_table(0, 2) = tab(1, 0);
or_table(0, 3) = tab(1, 1);
or_table(0, 4) = odds_ratio;
or_table(0, 5) = log(odds_ratio);
// calculate the contingency table
int a = 0, b = 0, c = 0, d = 0;
for (int i = 0; i < n; ++i) {
if (A[i] == 0) {
if (G[i] == 0) ++a; // not in A and not in G
else ++b; // not in A but in G
} else {
if (G[i] == 0) ++c; // in A but not in G
else ++d; // in A and in G
}
}

// calculate the odds ratio
double odds_ratio = static_cast<double>(a * d) / (b * c);

return or_table;
return NumericVector::create(
a, b, c, d, odds_ratio, std::log(odds_ratio), R_NaN, R_NaN
);
}

// Function to calculate the odds ratio and other relevant info for each permutation
// [[Rcpp::export]]
NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p = 100) {

NumericMatrix permuteOR(const IntegerMatrix& A, const IntegerVector& Gstar, int p = 100) {
int n = A.ncol();
NumericMatrix or_table(p, 8);

// calculate odds ratio for each permutation
for (int i = 0; i < p; ++i) {
IntegerVector Astar = shuffle_and_get_lower_triangle(A);
IntegerVector idx = sample(n, n, false) - 1;
IntegerVector Astar = get_triangle_from_index(A, idx);
or_table(i, _) = getOR(Astar, Gstar);
}

return(or_table);
}

// [[Rcpp::export]]
double getFDR_over(double actual, NumericVector permuted) {
double n = permuted.size();
double fdr = 0.0;
for (int i = 0; i < n; ++i) {
if (permuted[i] >= actual) {
fdr += 1.0;
}
}
fdr /= n;
return fdr;
}
List getFDR(double actual, const NumericVector& permuted) {
int n = permuted.size();
int count_over = 0;
int count_under = 0;

// [[Rcpp::export]]
double getFDR_under(double actual, NumericVector permuted) {
double n = permuted.size();
double fdr = 0.0;
// Count values above and below the actual value
for (int i = 0; i < n; ++i) {
if (permuted[i] <= actual) {
fdr += 1.0;
}
double current = permuted[i];
if (current >= actual) ++count_over;
if (current <= actual) ++count_under;
}
fdr /= n;
return fdr;

// Calculate FDR for both "over" and "under"
double fdr_over = static_cast<double>(count_over) / n;
double fdr_under = static_cast<double>(count_under) / n;

// Return both FDR values as a named list
return List::create(
Named("over") = fdr_over,
Named("under") = fdr_under
);
}

// Function to calculate the odds ratio and FDR, given the adjacency matrix A and the knowledge graph G
// [[Rcpp::export]]
NumericMatrix graflex(IntegerMatrix& A, IntegerMatrix& G, int p = 100) {
NumericVector graflex(const IntegerMatrix& A, const IntegerMatrix& G, int p = 100) {

// get the actual odds ratio
IntegerVector Astar = get_lower_triangle(A);
IntegerVector Gstar = get_lower_triangle(G);
NumericMatrix actual = getOR(Astar, Gstar);
IntegerVector Astar = get_triangle(A);
IntegerVector Gstar = get_triangle(G);
NumericVector actual = getOR(Astar, Gstar);

// get distribution of odds ratios on permuted data
NumericMatrix permuted = permuteOR(A, Gstar, p);

// calculate the FDR
actual(0, 6) = getFDR_under(actual(0, 4), permuted(_, 4));
actual(0, 7) = getFDR_over(actual(0, 4), permuted(_, 4));
List fdr = getFDR(actual(4), permuted(_, 4));
actual(6) = as<double>(fdr["under"]);
actual(7) = as<double>(fdr["over"]);

return actual;
}
Loading
0