Skip to content

Commit 272374b

Browse files
sumit-waliasumit-walia
authored
Merge pull request #12 from TurakhiaLab/placement_accuracy
Added support to handle polytomies Added functionality to re-estimate branch lengths without changing the topology
2 parents c5b98e1 + 14ba729 commit 272374b

19 files changed

Lines changed: 2817 additions & 425 deletions

‎cpu/mash_placement.hpp‎

Lines changed: 61 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,8 @@
1111
#define MASHPL_HPP
1212

1313
#include <stdint.h>
14+
#include <cmath>
15+
#include <functional>
1416
#include <iostream>
1517
#include <vector>
1618
#include <cstdio>
@@ -21,7 +23,8 @@
2123

2224
namespace MashPlacement
2325
{
24-
/** CLI and pipeline parameters; threads for TBB. */
26+
inline bool g_printBinaryNewick = false;
27+
2528
struct Param
2629
{
2730
uint64_t kmerSize;
@@ -327,6 +330,63 @@ namespace MashPlacement
327330
};
328331
static KPlacementDeviceArraysDC kplacementDeviceArraysDC;
329332

333+
inline void printNewickFromHostAdjacency(std::ostream& out, const std::vector<std::string>& name, int* h_head,
334+
int* h_e, int* h_nxt, double* h_len, int rootNode, bool binaryNewick) {
335+
const double kLenEps = 1e-12;
336+
auto hasKids = [h_head, h_e, h_nxt](int node, int from) -> bool {
337+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
338+
if (h_e[i] != from) {
339+
return true;
340+
}
341+
}
342+
return false;
343+
};
344+
std::function<void(int, int)> dfs;
345+
dfs = [&](int node, int from) {
346+
if (!hasKids(node, from)) {
347+
out << name[node];
348+
return;
349+
}
350+
if (binaryNewick) {
351+
out << "(";
352+
std::vector<int> pos;
353+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
354+
if (h_e[i] != from) {
355+
pos.push_back(i);
356+
}
357+
}
358+
for (size_t i = 0; i < pos.size(); ++i) {
359+
dfs(h_e[pos[i]], node);
360+
out << ":" << h_len[pos[i]] << (i + 1 == pos.size() ? ")" : ",");
361+
}
362+
return;
363+
}
364+
std::vector<std::pair<int, double>> kids;
365+
std::function<void(int, int)> addCollapsed;
366+
addCollapsed = [&](int n, int fr) {
367+
for (int i = h_head[n]; i != -1; i = h_nxt[i]) {
368+
int ch = h_e[i];
369+
if (ch == fr) {
370+
continue;
371+
}
372+
double L = h_len[i];
373+
if (std::fabs(L) <= kLenEps && hasKids(ch, n)) {
374+
addCollapsed(ch, n);
375+
} else {
376+
kids.push_back({ch, L});
377+
}
378+
}
379+
};
380+
addCollapsed(node, from);
381+
out << "(";
382+
for (size_t i = 0; i < kids.size(); ++i) {
383+
dfs(kids[i].first, node);
384+
out << ":" << kids[i].second << (i + 1 == kids.size() ? ")" : ",");
385+
}
386+
};
387+
dfs(rootNode, -1);
388+
}
389+
330390
};
331391

332392
#endif

‎src/divide_and_conquer/mash_placement.cuh‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,13 +11,17 @@
1111
#define MASHDC_CUH
1212

1313
#include <stdint.h>
14+
#include <cmath>
15+
#include <functional>
1416
#include <iostream>
1517
#include <vector>
1618
#include <cstdio>
1719
#include <string>
1820

1921
namespace MashPlacement
2022
{
23+
inline bool g_printBinaryNewick = false;
24+
2125
struct Param
2226
{
2327
uint64_t kmerSize;
@@ -225,6 +229,63 @@ namespace MashPlacement
225229
};
226230
static KPlacementDeviceArraysDC kplacementDeviceArraysDC;
227231

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

230291
#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
@@ -675,27 +675,6 @@ void MashPlacement::KPlacementDeviceArraysDC::printTreeDC(std::vector <std::stri
675675
double * h_len = new double[totalNumSequences*8];
676676
double * h_closest_dis = new double[totalNumSequences*20];
677677
int * h_closest_id = new int[totalNumSequences*20];
678-
std::function<void(int,int)> print=[&](int node, int from){
679-
if(h_nxt[h_head[node]]!=-1){
680-
// printf("(");
681-
output_ << "(";
682-
std::vector <int> pos;
683-
for(int i=h_head[node];i!=-1;i=h_nxt[i])
684-
if(h_e[i]!=from)
685-
pos.push_back(i);
686-
for(size_t i=0;i<pos.size();i++){
687-
print(h_e[pos[i]],node);
688-
// printf(":");
689-
// printf("%.5g%c",h_len[pos[i]],i+1==pos.size()?')':',');
690-
output_ << ":";
691-
output_ << h_len[pos[i]] << (i+1==pos.size()?')':',');
692-
}
693-
}
694-
// else std::cout<<name[node];
695-
else {
696-
output_<<name[node];
697-
}
698-
};
699678
auto err = cudaMemcpy(h_head, d_head, totalNumSequences*2*sizeof(int),cudaMemcpyDeviceToHost);
700679
if (err != cudaSuccess)
701680
{
@@ -728,8 +707,9 @@ void MashPlacement::KPlacementDeviceArraysDC::printTreeDC(std::vector <std::stri
728707
// }
729708
// std::cout << std::endl;
730709

731-
print(totalNumSequences+bd-2,-1);
732-
output_<<";\n";
710+
printNewickFromHostAdjacency(output_, name, h_head, h_e, h_nxt, h_len, totalNumSequences + bd - 2,
711+
g_printBinaryNewick);
712+
output_ << ";\n";
733713
}
734714

735715
__global__

‎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
/*

‎src/hip/mash_placement.cuh.hip‎

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;
@@ -305,6 +309,63 @@ namespace MashPlacement
305309
};
306310
static KPlacementDeviceArraysDC kplacementDeviceArraysDC;
307311

312+
inline void printNewickFromHostAdjacency(std::ostream& out, const std::vector<std::string>& name, int* h_head,
313+
int* h_e, int* h_nxt, double* h_len, int rootNode, bool binaryNewick) {
314+
const double kLenEps = 1e-12;
315+
auto hasKids = [h_head, h_e, h_nxt](int node, int from) -> bool {
316+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
317+
if (h_e[i] != from) {
318+
return true;
319+
}
320+
}
321+
return false;
322+
};
323+
std::function<void(int, int)> dfs;
324+
dfs = [&](int node, int from) {
325+
if (!hasKids(node, from)) {
326+
out << name[node];
327+
return;
328+
}
329+
if (binaryNewick) {
330+
out << "(";
331+
std::vector<int> pos;
332+
for (int i = h_head[node]; i != -1; i = h_nxt[i]) {
333+
if (h_e[i] != from) {
334+
pos.push_back(i);
335+
}
336+
}
337+
for (size_t i = 0; i < pos.size(); ++i) {
338+
dfs(h_e[pos[i]], node);
339+
out << ":" << h_len[pos[i]] << (i + 1 == pos.size() ? ")" : ",");
340+
}
341+
return;
342+
}
343+
std::vector<std::pair<int, double>> kids;
344+
std::function<void(int, int)> addCollapsed;
345+
addCollapsed = [&](int n, int fr) {
346+
for (int i = h_head[n]; i != -1; i = h_nxt[i]) {
347+
int ch = h_e[i];
348+
if (ch == fr) {
349+
continue;
350+
}
351+
double L = h_len[i];
352+
if (std::fabs(L) <= kLenEps && hasKids(ch, n)) {
353+
addCollapsed(ch, n);
354+
} else {
355+
kids.push_back({ch, L});
356+
}
357+
}
358+
};
359+
addCollapsed(node, from);
360+
out << "(";
361+
for (size_t i = 0; i < kids.size(); ++i) {
362+
dfs(kids[i].first, node);
363+
out << ":" << kids[i].second << (i + 1 == kids.size() ? ")" : ",");
364+
}
365+
};
366+
dfs(rootNode, -1);
367+
}
368+
308369
};
309370

310371
#endif

0 commit comments

Comments
 (0)