Skip to content

Commit 3153cf6

Browse files
sumit-waliasumit-walia
committed
Added support to handles polytomies and for converting given tree branch lengths to distance based (DIPPER) branch lenghts (fixed topology)
1 parent 6c78f96 commit 3153cf6

20 files changed

Lines changed: 2756 additions & 359 deletions

‎cpu/mash_placement.cuh‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@
22
#define MASHPL_CUH
33

44
#include <stdint.h>
5+
#include <cmath>
6+
#include <functional>
57
#include <iostream>
68
#include <vector>
79
#include <cstdio>
@@ -13,6 +15,8 @@
1315

1416
namespace MashPlacement
1517
{
18+
inline bool g_printBinaryNewick = false;
19+
1620
struct Param
1721
{
1822
uint64_t kmerSize;
@@ -288,6 +292,63 @@ namespace MashPlacement
288292
};
289293
static KPlacementDeviceArraysDC kplacementDeviceArraysDC;
290294

295+
inline void printNewickFromHostAdjacency(std::ostream& out, const std::vector<std::string>& name, int* h_head,
296+
int* h_e, int* h_nxt, double* h_len, int rootNode, bool binaryNewick) {
297+
const double kLenEps = 1e-12;
298+
auto hasKids = [h_head, h_e, h_nxt](int node, int from) -> bool {
299+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
300+
if (h_e[i] != from) {
301+
return true;
302+
}
303+
}
304+
return false;
305+
};
306+
std::function<void(int, int)> dfs;
307+
dfs = [&](int node, int from) {
308+
if (!hasKids(node, from)) {
309+
out << name[node];
310+
return;
311+
}
312+
if (binaryNewick) {
313+
out << "(";
314+
std::vector<int> pos;
315+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
316+
if (h_e[i] != from) {
317+
pos.push_back(i);
318+
}
319+
}
320+
for (size_t i = 0; i < pos.size(); ++i) {
321+
dfs(h_e[pos[i]], node);
322+
out << ":" << h_len[pos[i]] << (i + 1 == pos.size() ? ")" : ",");
323+
}
324+
return;
325+
}
326+
std::vector<std::pair<int, double>> kids;
327+
std::function<void(int, int)> addCollapsed;
328+
addCollapsed = [&](int n, int fr) {
329+
for (int i = h_head[n]; i != -1; i = h_nxt[i]) {
330+
int ch = h_e[i];
331+
if (ch == fr) {
332+
continue;
333+
}
334+
double L = h_len[i];
335+
if (std::fabs(L) <= kLenEps && hasKids(ch, n)) {
336+
addCollapsed(ch, n);
337+
} else {
338+
kids.push_back({ch, L});
339+
}
340+
}
341+
};
342+
addCollapsed(node, from);
343+
out << "(";
344+
for (size_t i = 0; i < kids.size(); ++i) {
345+
dfs(kids[i].first, node);
346+
out << ":" << kids[i].second << (i + 1 == kids.size() ? ")" : ",");
347+
}
348+
};
349+
dfs(rootNode, -1);
350+
}
351+
291352
};
292353

293354
#endif

‎cpu/placement_close_k.cu‎

Lines changed: 2 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -551,26 +551,6 @@ void MashPlacement::KPlacementDeviceArrays::printTree(std::vector <std::string>
551551
double * h_len = new double[numSequences*8];
552552
double * h_closest_dis = new double[numSequences*20];
553553
int * h_closest_id = new int[numSequences*20];
554-
std::function<void(int,int)> print=[&](int node, int from){
555-
if(h_nxt[h_head[node]]!=-1){
556-
// printf("(");
557-
output_ << "(";
558-
std::vector <int> pos;
559-
for(int i=h_head[node];i!=-1;i=h_nxt[i])
560-
if(h_e[i]!=from)
561-
pos.push_back(i);
562-
for(size_t i=0;i<pos.size();i++){
563-
print(h_e[pos[i]],node);
564-
// printf(":");
565-
// printf("%.5g%c",h_len[pos[i]],i+1==pos.size()?')':',');
566-
output_ << ":";
567-
// output_ << "%.5g%c",h_len[pos[i]],i+1==pos.size()?')':',';
568-
output_ << h_len[pos[i]] << (i+1==pos.size()?')':',');
569-
}
570-
}
571-
// else std::cout<<name[node];
572-
else output_<<name[node];
573-
};
574554
for (int i = 0; i < numSequences*2; ++i) {
575555
h_head[i] = d_head[i];
576556
}
@@ -625,7 +605,8 @@ void MashPlacement::KPlacementDeviceArrays::printTree(std::vector <std::string>
625605
// std::cerr<<"\n";
626606

627607

628-
print(numSequences+bd-2,-1);
608+
printNewickFromHostAdjacency(output_, name, h_head, h_e, h_nxt, h_len, numSequences + bd - 2,
609+
g_printBinaryNewick);
629610
// std::cout<<";\n";
630611
output_<<";\n";
631612
}

‎src/divide_and_conquer/mash_placement.cuh‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@
22
#define MASHDC_CUH
33

44
#include <stdint.h>
5+
#include <cmath>
6+
#include <functional>
57
#include <iostream>
68
#include <vector>
79
#include <cstdio>
@@ -12,6 +14,8 @@
1214

1315
namespace MashPlacement
1416
{
17+
inline bool g_printBinaryNewick = false;
18+
1519
struct Param
1620
{
1721
uint64_t kmerSize;
@@ -222,6 +226,63 @@ namespace MashPlacement
222226
};
223227
static KPlacementDeviceArraysDC kplacementDeviceArraysDC;
224228

229+
inline void printNewickFromHostAdjacency(std::ostream& out, const std::vector<std::string>& name, int* h_head,
230+
int* h_e, int* h_nxt, double* h_len, int rootNode, bool binaryNewick) {
231+
const double kLenEps = 1e-12;
232+
auto hasKids = [h_head, h_e, h_nxt](int node, int from) -> bool {
233+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
234+
if (h_e[i] != from) {
235+
return true;
236+
}
237+
}
238+
return false;
239+
};
240+
std::function<void(int, int)> dfs;
241+
dfs = [&](int node, int from) {
242+
if (!hasKids(node, from)) {
243+
out << name[node];
244+
return;
245+
}
246+
if (binaryNewick) {
247+
out << "(";
248+
std::vector<int> pos;
249+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
250+
if (h_e[i] != from) {
251+
pos.push_back(i);
252+
}
253+
}
254+
for (size_t i = 0; i < pos.size(); ++i) {
255+
dfs(h_e[pos[i]], node);
256+
out << ":" << h_len[pos[i]] << (i + 1 == pos.size() ? ")" : ",");
257+
}
258+
return;
259+
}
260+
std::vector<std::pair<int, double>> kids;
261+
std::function<void(int, int)> addCollapsed;
262+
addCollapsed = [&](int n, int fr) {
263+
for (int i = h_head[n]; i != -1; i = h_nxt[i]) {
264+
int ch = h_e[i];
265+
if (ch == fr) {
266+
continue;
267+
}
268+
double L = h_len[i];
269+
if (std::fabs(L) <= kLenEps && hasKids(ch, n)) {
270+
addCollapsed(ch, n);
271+
} else {
272+
kids.push_back({ch, L});
273+
}
274+
}
275+
};
276+
addCollapsed(node, from);
277+
out << "(";
278+
for (size_t i = 0; i < kids.size(); ++i) {
279+
dfs(kids[i].first, node);
280+
out << ":" << kids[i].second << (i + 1 == kids.size() ? ")" : ",");
281+
}
282+
};
283+
dfs(rootNode, -1);
284+
}
285+
225286
};
226287

227288
#endif

‎src/divide_and_conquer/placement_close_k.cpp‎

Lines changed: 5 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#include <queue>
55
#include <chrono>
66
#include <iostream>
7+
#include <sstream>
78
#include <tuple>
89
#include <cassert>
910
#include <algorithm>
@@ -593,25 +594,10 @@ void MashPlacement::KPlacementDeviceArraysHostDC::deallocateHostArraysDC(){
593594

594595

595596
void MashPlacement::KPlacementDeviceArraysHostDC::printTreeCpuDC(std::vector <std::string> name){
596-
597-
auto print=[&](int node, int from, auto&& print)->void {
598-
if(h_nxt[h_head[node]]!=-1){
599-
printf("(");
600-
std::vector <int> pos;
601-
for(int i=h_head[node];i!=-1;i=h_nxt[i])
602-
if(h_e[i]!=from)
603-
pos.push_back(i);
604-
for(size_t i=0;i<pos.size();i++){
605-
print(h_e[pos[i]],node, print);
606-
printf(":");
607-
printf("%.5g%c",h_len[pos[i]],i+1==pos.size()?')':',');
608-
}
609-
}
610-
else std::cout<<name[node];
611-
};
612-
613-
print(totalNumSequences+bd-2,-1, print);
614-
std::cout<<";\n";
597+
std::ostringstream oss;
598+
printNewickFromHostAdjacency(oss, name, h_head, h_e, h_nxt, h_len, totalNumSequences + bd - 2,
599+
g_printBinaryNewick);
600+
std::cout << oss.str() << ";\n";
615601
}
616602

617603
/* Clusterting function on CPU - > might need modification

‎src/divide_and_conquer/placement_close_k.cu‎

Lines changed: 3 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -659,27 +659,6 @@ void MashPlacement::KPlacementDeviceArraysDC::printTreeDC(std::vector <std::stri
659659
double * h_len = new double[totalNumSequences*8];
660660
double * h_closest_dis = new double[totalNumSequences*20];
661661
int * h_closest_id = new int[totalNumSequences*20];
662-
std::function<void(int,int)> print=[&](int node, int from){
663-
if(h_nxt[h_head[node]]!=-1){
664-
// printf("(");
665-
output_ << "(";
666-
std::vector <int> pos;
667-
for(int i=h_head[node];i!=-1;i=h_nxt[i])
668-
if(h_e[i]!=from)
669-
pos.push_back(i);
670-
for(size_t i=0;i<pos.size();i++){
671-
print(h_e[pos[i]],node);
672-
// printf(":");
673-
// printf("%.5g%c",h_len[pos[i]],i+1==pos.size()?')':',');
674-
output_ << ":";
675-
output_ << h_len[pos[i]] << (i+1==pos.size()?')':',');
676-
}
677-
}
678-
// else std::cout<<name[node];
679-
else {
680-
output_<<name[node];
681-
}
682-
};
683662
auto err = cudaMemcpy(h_head, d_head, totalNumSequences*2*sizeof(int),cudaMemcpyDeviceToHost);
684663
if (err != cudaSuccess)
685664
{
@@ -712,8 +691,9 @@ void MashPlacement::KPlacementDeviceArraysDC::printTreeDC(std::vector <std::stri
712691
// }
713692
// std::cout << std::endl;
714693

715-
print(totalNumSequences+bd-2,-1);
716-
output_<<";\n";
694+
printNewickFromHostAdjacency(output_, name, h_head, h_e, h_nxt, h_len, totalNumSequences + bd - 2,
695+
g_printBinaryNewick);
696+
output_ << ";\n";
717697
}
718698

719699
/*

‎src/divide_and_conquer/tree_generation_divide_and_conquer.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,8 @@ void parseArguments(int argc, char** argv)
4343
("algorithm,a", po::value<std::string>(), "Algorithm selection (0 - default mode, 1 - force placement, 2 - force conventional NJ)")
4444
("placement-mode,p", po::value<std::string>(), "Placement mode selection (0 - exact mode, 1 - k-closest mode), default is k-closest")
4545
// ("batch-size,b", po::value<std::string>(), "Batch size for GPU processing (Default = 100000)")
46+
("print-binary-newick",
47+
"Print fully resolved binary Newick with explicit :0 branches; default collapses near-zero internal edges")
4648
("help,h", "Print help messages");
4749

4850
}
@@ -194,6 +196,8 @@ int main(int argc, char** argv) {
194196
return 1;
195197
}
196198

199+
MashPlacement::g_printBinaryNewick = vm.count("print-binary-newick") > 0;
200+
197201
// Kmer Size
198202
uint64_t k = 15;
199203
try {k= (uint64_t)std::stoi(vm["kmer-size"].as<std::string>());}

‎src/hip/divide_and_conquer/placement_close_k.cu.hip‎

Lines changed: 3 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -657,27 +657,6 @@ void MashPlacement::KPlacementDeviceArraysDC::printTreeDC(std::vector <std::stri
657657
double * h_len = new double[totalNumSequences*8];
658658
double * h_closest_dis = new double[totalNumSequences*20];
659659
int * h_closest_id = new int[totalNumSequences*20];
660-
std::function<void(int,int)> print=[&](int node, int from){
661-
if(h_nxt[h_head[node]]!=-1){
662-
// printf("(");
663-
output_ << "(";
664-
std::vector <int> pos;
665-
for(int i=h_head[node];i!=-1;i=h_nxt[i])
666-
if(h_e[i]!=from)
667-
pos.push_back(i);
668-
for(size_t i=0;i<pos.size();i++){
669-
print(h_e[pos[i]],node);
670-
// printf(":");
671-
// printf("%.5g%c",h_len[pos[i]],i+1==pos.size()?')':',');
672-
output_ << ":";
673-
output_ << h_len[pos[i]] << (i+1==pos.size()?')':',');
674-
}
675-
}
676-
// else std::cout<<name[node];
677-
else {
678-
output_<<name[node];
679-
}
680-
};
681660
auto err = hipMemcpy(h_head, d_head, totalNumSequences*2*sizeof(int),hipMemcpyDeviceToHost);
682661
if (err != hipSuccess)
683662
{
@@ -710,8 +689,9 @@ void MashPlacement::KPlacementDeviceArraysDC::printTreeDC(std::vector <std::stri
710689
// }
711690
// std::cout << std::endl;
712691

713-
print(totalNumSequences+bd-2,-1);
714-
output_<<";\n";
692+
printNewickFromHostAdjacency(output_, name, h_head, h_e, h_nxt, h_len, totalNumSequences + bd - 2,
693+
g_printBinaryNewick);
694+
output_ << ";\n";
715695
}
716696

717697
/*

0 commit comments

Comments
 (0)