Skip to content

Commit

Permalink
int32 for .state
Browse files Browse the repository at this point in the history
  • Loading branch information
@ms609
ms609 committed Sep 27, 2023
1 parent 6b048b0 commit b5215c3
Showing 1 changed file with 44 additions and 42 deletions.
86 changes: 44 additions & 42 deletions src/tree_distances.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,26 +87,26 @@ List cpp_robinson_foulds_info (const RawMatrix x, const RawMatrix y,

/* Dynamic allocation 20% faster for 105 tips, but VLA not permitted in C11 */
splitbit b_complement[SL_MAX_SPLITS][SL_MAX_BINS];
for (int16 i = 0; i != b.n_splits; i++) {
for (int16 bin = 0; bin != last_bin; ++bin) {
for (int32 i = 0; i != b.n_splits; i++) {
for (int32 bin = 0; bin != last_bin; ++bin) {
b_complement[i][bin] = ~b.state[i][bin];
}
b_complement[i][last_bin] = b.state[i][last_bin] ^ unset_mask;
}

for (int16 ai = 0; ai != a.n_splits; ++ai) {
for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 ai = 0; ai != a.n_splits; ++ai) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {

bool all_match = true, all_complement = true;

for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
if ((a.state[ai][bin] != b.state[bi][bin])) {
all_match = false;
break;
}
}
if (!all_match) {
for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
if ((a.state[ai][bin] != b_complement[bi][bin])) {
all_complement = false;
break;
Expand All @@ -115,7 +115,7 @@ List cpp_robinson_foulds_info (const RawMatrix x, const RawMatrix y,
}
if (all_match || all_complement) {
int16 leaves_in_split = 0;
for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
leaves_in_split += count_bits(a.state[ai][bin]);
}

Expand Down Expand Up @@ -154,20 +154,20 @@ List cpp_matching_split_distance (const RawMatrix x, const RawMatrix y,
cost** score = new cost*[most_splits];
for (int16 i = most_splits; i--; ) score[i] = new cost[most_splits];

for (int16 ai = 0; ai != a.n_splits; ++ai) {
for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 ai = 0; ai != a.n_splits; ++ai) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {
score[ai][bi] = 0;
for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
score[ai][bi] += count_bits(a.state[ai][bin] ^ b.state[bi][bin]);
}
if (score[ai][bi] > half_tips) score[ai][bi] = n_tips - score[ai][bi];
}
for (int16 bi = b.n_splits; bi < most_splits; ++bi) {
for (int32 bi = b.n_splits; bi < most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
for (int16 ai = a.n_splits; ai < most_splits; ++ai) {
for (int16 bi = 0; bi != most_splits; ++bi) {
for (int32 ai = a.n_splits; ai < most_splits; ++ai) {
for (int32 bi = 0; bi != most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
Expand Down Expand Up @@ -217,18 +217,18 @@ List cpp_jaccard_similarity (const RawMatrix x, const RawMatrix y,
cost** score = new cost*[most_splits];
for (int16 i = most_splits; i--; ) score[i] = new cost[most_splits];

for (int16 ai = 0; ai != a.n_splits; ++ai) {
for (int32 ai = 0; ai != a.n_splits; ++ai) {

const int16
na = a.in_split[ai],
nA = n_tips - na
;

for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {

// x divides tips into a|A; y divides tips into b|B
int16 a_and_b = 0;
for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
a_and_b += count_bits(a.state[ai][bin] & b.state[bi][bin]);
}

Expand Down Expand Up @@ -282,12 +282,12 @@ List cpp_jaccard_similarity (const RawMatrix x, const RawMatrix y,
}
}
}
for (int16 bi = b.n_splits; bi < most_splits; ++bi) {
for (int32 bi = b.n_splits; bi < most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
for (int16 ai = a.n_splits; ai < most_splits; ++ai) {
for (int16 bi = 0; bi != most_splits; ++bi) {
for (int32 ai = a.n_splits; ai < most_splits; ++ai) {
for (int32 bi = 0; bi != most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
Expand Down Expand Up @@ -332,14 +332,14 @@ List cpp_msi_distance (const RawMatrix x, const RawMatrix y,

splitbit different[SL_MAX_BINS];

for (int16 ai = 0; ai != a.n_splits; ++ai) {
for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 ai = 0; ai != a.n_splits; ++ai) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {
int16
n_different = 0,
n_a_only = 0,
n_a_and_b = 0
;
for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
different[bin] = a.state[ai][bin] ^ b.state[bi][bin];
n_different += count_bits(different[bin]);
n_a_only += count_bits(a.state[ai][bin] & different[bin]);
Expand All @@ -351,12 +351,12 @@ List cpp_msi_distance (const RawMatrix x, const RawMatrix y,
((max_score / max_possible) *
mmsi_score(n_same, n_a_and_b, n_different, n_a_only));
}
for (int16 bi = b.n_splits; bi < most_splits; ++bi) {
for (int32 bi = b.n_splits; bi < most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
for (int16 ai = a.n_splits; ai < most_splits; ++ai) {
for (int16 bi = 0; bi < most_splits; ++bi) {
for (int32 ai = a.n_splits; ai < most_splits; ++ai) {
for (int32 bi = 0; bi < most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
Expand Down Expand Up @@ -388,10 +388,12 @@ List cpp_msi_distance (const RawMatrix x, const RawMatrix y,
// [[Rcpp::export]]
List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y,
const IntegerVector nTip) {
Rcpp::Rcout << "\n\n\n\nASOUIGHAUI SFH ASHF KLASJHF KAJSF \n\n\n\n";
if (x.cols() != y.cols()) {
Rcpp::stop("Input splits must address same number of tips.");
}
const SplitList a(x), b(y);
Rcpp::Rcout << "\n\n\n\nSPLITS LISTED \n\n\n\n";
const bool a_has_more_splits = (a.n_splits > b.n_splits);
const int16
most_splits = a_has_more_splits ? a.n_splits : b.n_splits,
Expand All @@ -414,18 +416,18 @@ List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y,
NumericVector a_match(a.n_splits);
std::unique_ptr<int16[]> b_match = std::make_unique<int16[]>(b.n_splits);

for (int16 ai = 0; ai != a.n_splits; ++ai) {
for (int32 ai = 0; ai != a.n_splits; ++ai) {
if (a_match[ai]) continue;
const int16
na = a.in_split[ai],
nA = n_tips - na
;

for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {

// x divides tips into a|A; y divides tips into b|B
int16 a_and_b = 0;
for (int16 bin = 0; bin != a.n_bins; ++bin) {
for (int32 bin = 0; bin != a.n_bins; ++bin) {
a_and_b += count_bits(a.state[ai][bin] & b.state[bi][bin]);
}

Expand Down Expand Up @@ -461,7 +463,7 @@ List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y,
);
}
}
for (int16 bi = b.n_splits; bi < most_splits; ++bi) {
for (int32 bi = b.n_splits; bi < most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
Expand All @@ -482,21 +484,21 @@ List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y,

if (exact_matches) {
int16 a_pos = 0;
for (int16 ai = 0; ai != a.n_splits; ++ai) {
for (int32 ai = 0; ai != a.n_splits; ++ai) {
if (a_match[ai]) continue;
int16 b_pos = 0;
for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {
if (b_match[bi]) continue;
score[a_pos][b_pos] = score[ai][bi];
b_pos++;
}
for (int16 bi = lap_dim - a_extra_splits; bi < lap_dim; ++bi) {
for (int32 bi = lap_dim - a_extra_splits; bi < lap_dim; ++bi) {
score[a_pos][bi] = max_score;
}
a_pos++;
}
for (int16 ai = lap_dim - b_extra_splits; ai < lap_dim; ++ai) {
for (int16 bi = 0; bi != lap_dim; ++bi) {
for (int32 ai = lap_dim - b_extra_splits; ai < lap_dim; ++ai) {
for (int32 bi = 0; bi != lap_dim; ++bi) {
score[ai][bi] = max_score;
}
}
Expand All @@ -512,7 +514,7 @@ List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y,

std::unique_ptr<int16[]> lap_decode = std::make_unique<int16[]>(lap_dim);
int16 fuzzy_match = 0;
for (int16 bi = 0; bi != b.n_splits; ++bi) {
for (int32 bi = 0; bi != b.n_splits; ++bi) {
if (!b_match[bi]) {
assert(fuzzy_match < lap_dim);
lap_decode[fuzzy_match++] = bi + 1;
Expand Down Expand Up @@ -547,8 +549,8 @@ List cpp_mutual_clustering (const RawMatrix x, const RawMatrix y,
return List::create(Named("score") = final_score,
_["matching"] = final_matching);
} else {
for (int16 ai = a.n_splits; ai < most_splits; ++ai) {
for (int16 bi = 0; bi != most_splits; ++bi) {
for (int32 ai = a.n_splits; ai < most_splits; ++ai) {
for (int32 bi = 0; bi != most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
Expand Down Expand Up @@ -597,8 +599,8 @@ List cpp_shared_phylo (const RawMatrix x, const RawMatrix y,
cost** score = new cost*[most_splits];
for (int16 i = most_splits; i--; ) score[i] = new cost[most_splits];

for (int16 ai = a.n_splits; ai--; ) {
for (int16 bi = b.n_splits; bi--; ) {
for (int32 ai = a.n_splits; ai--; ) {
for (int32 bi = b.n_splits; bi--; ) {
const double spi_over = spi_overlap(a.state[ai], b.state[bi], n_tips,
a.in_split[ai], b.in_split[bi],
a.n_bins);
Expand All @@ -608,12 +610,12 @@ List cpp_shared_phylo (const RawMatrix x, const RawMatrix y,
max_score;

}
for (int16 bi = b.n_splits; bi < most_splits; ++bi) {
for (int32 bi = b.n_splits; bi < most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
for (int16 ai = a.n_splits; ai < most_splits; ++ai) {
for (int16 bi = 0; bi != most_splits; ++bi) {
for (int32 ai = a.n_splits; ai < most_splits; ++ai) {
for (int32 bi = 0; bi != most_splits; ++bi) {
score[ai][bi] = max_score;
}
}
Expand Down

0 comments on commit b5215c3

Please sign in to comment.