prokrustean/prokrustean.unitig.count.cpp at main · KoslickiLab/prokrustean · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
// tbd: license comment
/*
 * prokrustean.cpp
 *  Created on: Jul 29, 2023
 *      Author: Adam Park
 */

#include <iostream>
#include <unistd.h>
#include "src/fm_index/index.hpp"
#include "src/fm_index/string.sdsl.hpp"
#include "src/construction/algorithms.hpp"
#include "src/prokrustean.hpp"
#include "src/prokrustean.support.hpp"
#include "src/application/dbg.count.hpp"

using namespace std;
using namespace sdsl;

string input_prokrustean;
string output_file;
int num_threads=8;
int from=-1;
int to=-1;
char TERM = '$';

void help(){

	cout << "unitig counting [options]" << endl <<
	"Input: prokrustean graph file." << endl <<
	"Output: maximal unitig counts count." << endl <<
	"Options:" << endl <<
	"-h          help" << endl <<
	"-p <arg>    (REQUIRED) prokrustean file name" << endl <<
	"-l <arg>    k range left. default: kmin." << endl <<
	"-r <arg>    k range right. default: largest sequence size." << endl <<
	"-o <arg>    output file name. Default: input file name + .unitig.count.txt" << endl <<
	"-t <arg>    thread count Default:" << num_threads << endl;
	exit(0);
}

int main(int argc, char** argv){

	if(argc < 2) help();
	int opt;
	while ((opt = getopt(argc, argv, "hp:l:r:o:t:")) != -1){
		switch (opt){
			case 'h':
				help();
			break;
			case 'p':
				input_prokrustean = string(optarg);
			break;
			case 'o':
				output_file = string(optarg);
			break;
			case 'l':
				from = stoi(string(optarg));
			break;
			case 'r':
				to = stoi(string(optarg));
			break;
			case 't':
				num_threads = stoi(string(optarg));
			break;
			default:
				help();
			return -1;
		}
	}

	if(input_prokrustean.size()==0){
		cout << "input empty" << endl;
		help();
	}
	if(output_file.size()==0) {
		output_file=input_prokrustean+".unitig.count.txt";
	};
	auto start_total = std::chrono::steady_clock::now();
	Prokrustean prokrustean;
	bool success=load_prokrustean(input_prokrustean, prokrustean);
	if(!success){
		cout << "loading failed " << input_prokrustean << endl;
		exit(0);
	}
	if(from==-1){
		from=prokrustean.kmin;
	} else if(from<prokrustean.kmin){
		cout << "from(l) value should be at least the kmin of prokrustean: " << prokrustean.kmin << endl;
		exit(0);
	}
	if(to==-1){
		to=prokrustean.sequences__size[0];
		for(auto &size: prokrustean.sequences__size){
			if(to<size) to=size;
		}
	} else if(to<from){
		cout << "to(r) value should be at least from(l) value: " << from << endl;
		exit(0);
	}

	cout << "k range " << from << ":" << to << endl;
	cout << "threads: " << num_threads << endl;

	prokrustean.print_abstract();

	auto start = std::chrono::steady_clock::now();

	ProkrusteanExtension ext(prokrustean);
	// count_left_right_character_extensions_parallel(ext, num_threads);

	// cout << "count left right" << (std::chrono::steady_clock::now()-start).count()/1000000 << "ms" << endl;
	// start = std::chrono::steady_clock::now();

	vector<uint64_t> output;
	count_maximal_unitigs_range_of_k_parallel(from, to, ext, output, num_threads);

	cout << (std::chrono::steady_clock::now()-start).count()/1000000 << "ms" << endl;
	cout << "save unitig counts... ";

	std::ofstream outputFile(output_file);
	for(int i=0; i<output.size(); i++){
		if(from<=i && i<=to){
//			outputFile << std::left << std::setw(20) << i;
//			outputFile << output[i] << endl;
			outputFile << i << "," << output[i] << endl;
		}
	}

	outputFile.close();
	cout << "total " << (std::chrono::steady_clock::now()-start_total).count()/1000000 << "ms" << endl;
}