Skip to content

Commit

Permalink
Merge pull request #50 from suzannejin/master
Browse files Browse the repository at this point in the history 
speed up graflex again
  • Loading branch information
@suzannejin
suzannejin authored Sep 13, 2024
2 parents 9a8e679 + 1a06c77 commit b557940
Show file tree
Hide file tree
Showing 6 changed files with 245 additions and 154 deletions.
5 changes: 5 additions & 0 deletions NEWS.md
View file
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
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
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
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 comments on commit b557940

Please sign in to comment.