<commit_before>// nonpareil - Calculation of nonpareil curves
// @author <NAME> <<EMAIL> at gmail dot com>
// @version 3.0
// @license artistic license 2.0

/*
  ROADMAP

  o To validate "logarithmic sampling".

*/

#include <math.h>
#include <getopt.h>
#include "enveomics/universal.h"
#include "enveomics/multinode.h"
#include "enveomics/sequence.h"
#include "enveomics/nonpareil_mating.h"
#include "enveomics/nonpareil_sampling.h"
#include "enveomics/References.h"
#include "enveomics/KmerCounter.h"
#include <string>

#define LARGEST_PATH 4096
#define NP_VERSION 3.000

using namespace std;
int processID;
int processes;

void help(const char *msg){
   if(processID==0 && msg!=NULL && strlen(msg) != 0) cerr << endl << msg << endl << endl;
   if(processID==0) cerr << "DESCRIPTION"<< endl
   	<< "   Nonpareil uses the redundancy of the reads in metagenomic datasets to estimate the average coverage and predict the" << endl
	<< "   amount of sequences that will be required to achieve 'nearly complete coverage'." << endl
	<< endl
   	<< "USAGE" << endl
  << "   nonpareil -s sequences.fa -T alignment -b output [options]" << endl
	<< "   nonpareil -s sequences.fa -T kmer -f fastq -b output [options]" << endl
	<< "   nonpareil -h" << endl
	<< "   nonpareil -V" << endl
	<< endl
   	<< "MANDATORY ARGUMENTS" << endl
	<< "   -s <str> : Path to the (input) file containing the sequences.  This is lowercase S." << endl
	<< endl
	<< "COMMON OPTIONS" << endl
	<< "   -f <str> : The format of the sequences.  Can be 'fasta' or 'fastq'.  By default: 'fasta'." << endl
	<< "   -b <str> : Path to the prefix for all the output files.  Replaces the options: -a, -C, -l, and -o; generating files" << endl
	<< "              with the suffixes .npa, npc, .npl, and .npo, respectively, unless explicitly set." << endl
	<< "   -d <num> : Subsample iteratively applying this factor to the number of reads, resulting in logarithmic subsampling." << endl
	<< "              Use -d 0 to fall back to linear sampling, controlled by -m, -M, & -i (this was the default before v2.4)." << endl
	<< "              By default: 0.7." << endl
	<< "   -n <int> : Number of sub-samples to generate per point.  If it is not a multiple of the number of threads (see -t)," << endl
	<< "              it is rounded to the next (upper) multiple.  By default: 1024." << endl
	<< "   -L <num> : Minimum overlapping percentage of the aligned region on the largest sequence. The similarity (see -S) is" << endl
	<< "              evaluated for the aligned region only.  By default: 50." << endl
	<< "   -X <int> : Maximum number of reads to use as query.  This is capital X.  By default, 1,000 reads." << endl
	<< "   -q <str> : Path to the (input) file containing a second dataset to be used as query, for dataset comparisons.  This" << endl
	<< "              option is currently experimental." << endl
	<< "   -R <int> : Maximum RAM usage in Mib.  Ideally this value should be larger than the sequences to analyze (discarding" << endl
	<< "              non-sequence elements like headers or quality).  This is particularly important when running in multiple" << endl
	<< "              cores (see -t).  This value is approximated.  By default 1024." << endl
	<< "              Maximum value in this version: " << (UINT_MAX/1024) << endl
	<< "   -t <int> : Number of threads.  Highest efficiency when the number of sub-samples (see -n) is multiple of the number" << endl
	<< "              of threads.  By default: 2." << endl
	<< "   -v <int> : Verbosity level, for debugging purposes.  By default 7.  This is lowercase V." << endl
	<< "   -V       : Show version information and exit.  This is uppercase V." << endl
	<< "   -h       : Display this message and exit." << endl
  << "   -k <int> : kmer length. By default: 24" << endl
  << "   -T <str> : Nonpareil algorithm kmer or alignment accepted"
	<< endl
	<< "See all supported arguments and additional documentation at http://nonpareil.readthedocs.org or execute man nonpareil." << endl
	<< endl;
   finalize_multinode();
   if(processID==0){ exit(1); }else{ exit(0); }
}

int main(int argc, char *argv[]) {
   init_multinode(argc, argv, processID, processes);
   if(processID==0) cout << "Nonpareil v" << NP_VERSION << endl;
   if(argc<=1) help("");

   // Vars
   char		*file, *format=(char *)"fasta", *nonpareiltype, *alldata, *cntfile, *outfile, *namFile,
   		*seqFile, *baseout, *qfile, *qNamFile, *qSeqFile;
   double	min=0.0, max=1.0, itv=0.01, qry_portion=0, min_sim=0.95, ovl=0.50, *sample_result,
   		avg_seq_len, adj_avg_seq_len, divide=0.7, q_avg_seq_len;
   int		v=7, largest_seq, rseed=time(NULL), n=1024, k=24, thr=2, ram=1024, *mates, samples_no,
   		sample_i, sample_after_20, sampling_points, q_largest_seq;
   unsigned int	q_total_seqs, lines_in_ram, hX=1000, qry_seqs_no, ram_Kb, required_ram_Kb;
   unsigned long long int total_seqs;
   bool		n_as_mismatch=false, portion_label=false, revcom=true, ok, autoadjust=false,
   		alt_query=false;
   matepar_t	matepar;
   samplepar_t 	samplepar;

   alldata = (char *)"";
   cntfile = (char *)"";
   outfile = (char *)"-";


   // GetOpt
   int		optchr;
   while ((optchr = getopt (argc, argv, "Aa:b:cC:d:f:Fhi:k:l:L:m:M:n:No:q:r:R:s:S:t:T:v:Vx:X:")) != EOF)
      switch(optchr) {
	 case 'a': alldata = optarg;		break;
	 case 'A': autoadjust = true;		break;
	 case 'b': baseout = optarg;		break;
	 case 'c': revcom = false;		break;
	 case 'C': cntfile = optarg;		break;
	 case 'd': divide = atof(optarg);	break;
	 case 'f': format = optarg;		break;
	 case 'F': portion_label = true;	break;
	 case 'h': help("");			break;
	 case 'i': itv = atof(optarg);		break;
   case 'k': k = atoi(optarg);  break;
	 case 'l': open_log(optarg);		break;
	 case 'L': ovl=atof(optarg)/100.0;	break;
	 case 'm': min = atof(optarg);		break;
	 case 'M': max = atof(optarg);		break;
	 case 'n': n = atoi(optarg);		break;
	 case 'N': n_as_mismatch=true;		break;
	 case 'o': outfile = optarg;		break;
	 case 'q': qfile = optarg; alt_query = true;	break;
	 case 'r': rseed=atoi(optarg);		break;
	 case 'R': ram = atoi(optarg);		break;
   case 's': file = optarg;		break;
	 case 'S': min_sim=atof(optarg);	break;
	 case 't': thr = atoi(optarg);		break;
	 case 'T': nonpareiltype = optarg; break; // For backwards compatibility
	 case 'v': v = atoi(optarg);		break;
	 case 'V': finalize_multinode(); return 0;
	 case 'x': qry_portion = atof(optarg);	break;
	 case 'X': hX = atoi(optarg);		break;
     }
   // Initialize
   set_verbosity(v);
   if(strlen(nonpareiltype)==0) help("");
   if(strcmp(nonpareiltype,"kmer")!=0 && strcmp(nonpareiltype,"alignment")!=0)
        help("Bad argument for -t option, accepted values are kmer or alignment");
   if(strlen(file)==0) help("");
   if(strcmp(format, "fasta")!=0 & strcmp(format, "fastq")!=0)
   				help("Unsupported value for -f option");
   if(min<0 | min>1)		help("Bad argument for -m option, accepted values are numbers in the range [0, 1]");
   if(max<0 | max>1)		help("Bad argument for -M option, accepted values are numbers in the range [0, 1]");
   if(itv<=0 | itv>1)		help("Bad argument for -i option, accepted values are numbers in the range (0, 1]");
   if(ovl<=0.0 | ovl>1.0)	help("Bad argument for -L option, accepted values are numbers in the range (0, 100]");
   if(thr<=0)			help("Bad argumement for -t option, accepted values are positive non-zero integers");
   if(n<=0)			help("Bad argument for -n option, accepted values are positive non-zero integers");
   if(ram<=0)			help("Bad argument for -R option, accepted values are positive non-zero integers");
   if(min_sim<=0 | min_sim>1)	help("Bad argument for -S option, accepted values are numbers in the range (0, 1]");
   if(thr<=0)			help("Bad argument for -t option, accepted values are positive non-zero integers");
   if(qry_portion<0 | qry_portion>1)
   				help("Bad argument for -x option, accepted values are numbers in the range (0, 1]");
   if(divide<0 | divide>=1)	help("Bad argument for -d option, accepted values are numbers in the range (0, 1)");
   char alldataTmp[LARGEST_PATH], outfileTmp[LARGEST_PATH], cntfileTmp[LARGEST_PATH];
   if(baseout && (strlen(baseout)>0)){
      if(!alldata || (strlen(alldata)<=0)){ sprintf(alldataTmp, "%s.npa", baseout); alldata=alldataTmp; }
      if(!cntfile || (strlen(cntfile)<=0)){ sprintf(cntfileTmp, "%s.npc", baseout); cntfile=cntfileTmp; }
      if(!outfile || (strcmp(outfile, "-")==0)){ sprintf(outfileTmp, "%s.npo", baseout); outfile=outfileTmp; }
      if(processID==0 && !log_is_open()) {
         char logfile[LARGEST_PATH];
	 sprintf(logfile, "%s.npl", baseout);
	 open_log(logfile);
      }
   }

   if(alldata && (strlen(alldata)>0)) remove(alldata);
   if(cntfile && (strlen(cntfile)>0)) remove(cntfile);
   if(outfile && (strlen(outfile)>0) & (strcmp(outfile, "-")!=0)) remove(outfile);

   if(strcmp(nonpareiltype,"kmer")==0) {
     if(alt_query) {
       if(strcmp(format,"fasta")==0) {
         cout << "testing" << endl;
         ifstream qifs((string(qfile)));
         say("1ss$","reading query", file);
         FastaReader qfastaReader(qifs);
         References references = References(qfastaReader, k, alt_query);
         say("1s$","Started counting");
         ifstream ifs((string(file)));
         FastaReader fastaReader(ifs);
         KmerCounter counter = KmerCounter(references, fastaReader, string(cntfile));
         mates = new int[hX];
         counter.getCounts(mates);
         avg_seq_len = counter.getAvgLen();
         total_seqs = counter.getTotalSeqs();
         qry_seqs_no = counter.getTotalQSeqs();
         adj_avg_seq_len = avg_seq_len - k + 1;
         say("1sus$", "Read file with ", total_seqs, " sequences");
         say("1sfs$", "Average read length is ", avg_seq_len, "bp");
         goto restart_samples;
       }
     }
     else {
       if(strcmp(format,"fastq")==0) {
         ifstream ifs((string(file)));
         say("1ss$","reading ", file);
         FastqReader fastqReader(ifs);
         say("1sus$","Picking ", hX, " random sequences");
         References references = References(fastqReader, hX, k);
         say("1s$","Started counting");
         KmerCounter counter = KmerCounter(references, fastqReader, string(cntfile));
         mates = new int[hX];
         counter.getCounts(mates);<commit_msg>fixed segfault<commit_after>// nonpareil - Calculation of nonpareil curves
// @author <NAME> <<EMAIL> at g<EMAIL> dot com>
// @version 3.0
// @license artistic license 2.0

/*
  ROADMAP

  o To validate "logarithmic sampling".

*/

#include <math.h>
#include <getopt.h>
#include "enveomics/universal.h"
#include "enveomics/multinode.h"
#include "enveomics/sequence.h"
#include "enveomics/nonpareil_mating.h"
#include "enveomics/nonpareil_sampling.h"
#include "enveomics/References.h"
#include "enveomics/KmerCounter.h"
#include <string>

#define LARGEST_PATH 4096
#define NP_VERSION 3.000

using namespace std;
int processID;
int processes;

void help(const char *msg){
   if(processID==0 && msg!=NULL && strlen(msg) != 0) cerr << endl << msg << endl << endl;
   if(processID==0) cerr << "DESCRIPTION"<< endl
   	<< "   Nonpareil uses the redundancy of the reads in metagenomic datasets to estimate the average coverage and predict the" << endl
	<< "   amount of sequences that will be required to achieve 'nearly complete coverage'." << endl
	<< endl
   	<< "USAGE" << endl
  << "   nonpareil -s sequences.fa -T alignment -b output [options]" << endl
	<< "   nonpareil -s sequences.fa -T kmer -f fastq -b output [options]" << endl
	<< "   nonpareil -h" << endl
	<< "   nonpareil -V" << endl
	<< endl
   	<< "MANDATORY ARGUMENTS" << endl
	<< "   -s <str> : Path to the (input) file containing the sequences.  This is lowercase S." << endl
	<< endl
	<< "COMMON OPTIONS" << endl
	<< "   -f <str> : The format of the sequences.  Can be 'fasta' or 'fastq'.  By default: 'fasta'." << endl
	<< "   -b <str> : Path to the prefix for all the output files.  Replaces the options: -a, -C, -l, and -o; generating files" << endl
	<< "              with the suffixes .npa, npc, .npl, and .npo, respectively, unless explicitly set." << endl
	<< "   -d <num> : Subsample iteratively applying this factor to the number of reads, resulting in logarithmic subsampling." << endl
	<< "              Use -d 0 to fall back to linear sampling, controlled by -m, -M, & -i (this was the default before v2.4)." << endl
	<< "              By default: 0.7." << endl
	<< "   -n <int> : Number of sub-samples to generate per point.  If it is not a multiple of the number of threads (see -t)," << endl
	<< "              it is rounded to the next (upper) multiple.  By default: 1024." << endl
	<< "   -L <num> : Minimum overlapping percentage of the aligned region on the largest sequence. The similarity (see -S) is" << endl
	<< "              evaluated for the aligned region only.  By default: 50." << endl
	<< "   -X <int> : Maximum number of reads to use as query.  This is capital X.  By default, 1,000 reads." << endl
	<< "   -q <str> : Path to the (input) file containing a second dataset to be used as query, for dataset comparisons.  This" << endl
	<< "              option is currently experimental." << endl
	<< "   -R <int> : Maximum RAM usage in Mib.  Ideally this value should be larger than the sequences to analyze (discarding" << endl
	<< "              non-sequence elements like headers or quality).  This is particularly important when running in multiple" << endl
	<< "              cores (see -t).  This value is approximated.  By default 1024." << endl
	<< "              Maximum value in this version: " << (UINT_MAX/1024) << endl
	<< "   -t <int> : Number of threads.  Highest efficiency when the number of sub-samples (see -n) is multiple of the number" << endl
	<< "              of threads.  By default: 2." << endl
	<< "   -v <int> : Verbosity level, for debugging purposes.  By default 7.  This is lowercase V." << endl
	<< "   -V       : Show version information and exit.  This is uppercase V." << endl
	<< "   -h       : Display this message and exit." << endl
  << "   -k <int> : kmer length. By default: 24" << endl
  << "   -T <str> : Nonpareil algorithm kmer or alignment accepted"
	<< endl
	<< "See all supported arguments and additional documentation at http://nonpareil.readthedocs.org or execute man nonpareil." << endl
	<< endl;
   finalize_multinode();
   if(processID==0){ exit(1); }else{ exit(0); }
}

int main(int argc, char *argv[]) {
   init_multinode(argc, argv, processID, processes);
   if(processID==0) cout << "Nonpareil v" << NP_VERSION << endl;
   if(argc<=1) help("");

   // Vars
   char		*file, *format=(char *)"fasta", *nonpareiltype, *alldata, *cntfile, *outfile, *namFile,
   		*seqFile, *baseout, *qfile, *qNamFile, *qSeqFile;
   double	min=0.0, max=1.0, itv=0.01, qry_portion=0, min_sim=0.95, ovl=0.50, *sample_result,
   		avg_seq_len, adj_avg_seq_len, divide=0.7, q_avg_seq_len;
   int		v=7, largest_seq, rseed=time(NULL), n=1024, k=24, thr=2, ram=1024, *mates, samples_no,
   		sample_i, sample_after_20, sampling_points, q_largest_seq;
   unsigned int	q_total_seqs, lines_in_ram, hX=1000, qry_seqs_no, ram_Kb, required_ram_Kb;
   unsigned long long int total_seqs;
   bool		n_as_mismatch=false, portion_label=false, revcom=true, ok, autoadjust=false,
   		alt_query=false;
   matepar_t	matepar;
   samplepar_t 	samplepar;

   alldata = (char *)"";
   cntfile = (char *)"";
   outfile = (char *)"-";


   // GetOpt
   int		optchr;
   while ((optchr = getopt (argc, argv, "Aa:b:cC:d:f:Fhi:k:l:L:m:M:n:No:q:r:R:s:S:t:T:v:Vx:X:")) != EOF)
      switch(optchr) {
	 case 'a': alldata = optarg;		break;
	 case 'A': autoadjust = true;		break;
	 case 'b': baseout = optarg;		break;
	 case 'c': revcom = false;		break;
	 case 'C': cntfile = optarg;		break;
	 case 'd': divide = atof(optarg);	break;
	 case 'f': format = optarg;		break;
	 case 'F': portion_label = true;	break;
	 case 'h': help("");			break;
	 case 'i': itv = atof(optarg);		break;
   case 'k': k = atoi(optarg);  break;
	 case 'l': open_log(optarg);		break;
	 case 'L': ovl=atof(optarg)/100.0;	break;
	 case 'm': min = atof(optarg);		break;
	 case 'M': max = atof(optarg);		break;
	 case 'n': n = atoi(optarg);		break;
	 case 'N': n_as_mismatch=true;		break;
	 case 'o': outfile = optarg;		break;
	 case 'q': qfile = optarg; alt_query = true;	break;
	 case 'r': rseed=atoi(optarg);		break;
	 case 'R': ram = atoi(optarg);		break;
   case 's': file = optarg;		break;
	 case 'S': min_sim=atof(optarg);	break;
	 case 't': thr = atoi(optarg);		break;
	 case 'T': nonpareiltype = optarg; break; // For backwards compatibility
	 case 'v': v = atoi(optarg);		break;
	 case 'V': finalize_multinode(); return 0;
	 case 'x': qry_portion = atof(optarg);	break;
	 case 'X': hX = atoi(optarg);		break;
     }
   // Initialize
   set_verbosity(v);
   if(strlen(nonpareiltype)==0) help("");
   if(strcmp(nonpareiltype,"kmer")!=0 && strcmp(nonpareiltype,"alignment")!=0)
        help("Bad argument for -t option, accepted values are kmer or alignment");
   if(strlen(file)==0) help("");
   if(strcmp(format, "fasta")!=0 & strcmp(format, "fastq")!=0)
   				help("Unsupported value for -f option");
   if(min<0 | min>1)		help("Bad argument for -m option, accepted values are numbers in the range [0, 1]");
   if(max<0 | max>1)		help("Bad argument for -M option, accepted values are numbers in the range [0, 1]");
   if(itv<=0 | itv>1)		help("Bad argument for -i option, accepted values are numbers in the range (0, 1]");
   if(ovl<=0.0 | ovl>1.0)	help("Bad argument for -L option, accepted values are numbers in the range (0, 100]");
   if(thr<=0)			help("Bad argumement for -t option, accepted values are positive non-zero integers");
   if(n<=0)			help("Bad argument for -n option, accepted values are positive non-zero integers");
   if(ram<=0)			help("Bad argument for -R option, accepted values are positive non-zero integers");
   if(min_sim<=0 | min_sim>1)	help("Bad argument for -S option, accepted values are numbers in the range (0, 1]");
   if(thr<=0)			help("Bad argument for -t option, accepted values are positive non-zero integers");
   if(qry_portion<0 | qry_portion>1)
   				help("Bad argument for -x option, accepted values are numbers in the range (0, 1]");
   if(divide<0 | divide>=1)	help("Bad argument for -d option, accepted values are numbers in the range (0, 1)");
   char alldataTmp[LARGEST_PATH], outfileTmp[LARGEST_PATH], cntfileTmp[LARGEST_PATH];
   if(baseout && (strlen(baseout)>0)){
      if(!alldata || (strlen(alldata)<=0)){ sprintf(alldataTmp, "%s.npa", baseout); alldata=alldataTmp; }
      if(!cntfile || (strlen(cntfile)<=0)){ sprintf(cntfileTmp, "%s.npc", baseout); cntfile=cntfileTmp; }
      if(!outfile || (strcmp(outfile, "-")==0)){ sprintf(outfileTmp, "%s.npo", baseout); outfile=outfileTmp; }
      if(processID==0 && !log_is_open()) {
         char logfile[LARGEST_PATH];
	 sprintf(logfile, "%s.npl", baseout);
	 open_log(logfile);
      }
   }

   if(alldata && (strlen(alldata)>0)) remove(alldata);
   if(cntfile && (strlen(cntfile)>0)) remove(cntfile);
   if(outfile && (strlen(outfile)>0) & (strcmp(outfile, "-")!=0)) remove(outfile);

   if(strcmp(nonpareiltype,"kmer")==0) {
     if(alt_query) {
       if(strcmp(format,"fasta")==0) {
         cout << "testing" << endl;
         ifstream qifs((string(qfile)));
         say("1ss$","reading query", file);
         FastaReader qfastaReader(qifs);
         References references = References(qfastaReader, k, alt_query);
         say("1s$","Started counting");
         ifstream ifs((string(file)));
         FastaReader fastaReader(ifs);
         KmerCounter counter = KmerCounter(references, fastaReader, string(cntfile));
         mates = new int[counter.getTotalSeqs()];
         counter.getCounts(mates);
         avg_seq_len = counter.getAvgLen();
         total_seqs = counter.getTotalSeqs();
         qry_seqs_no = counter.getTotalQSeqs();
         adj_avg_seq_len = avg_seq_len - k + 1;
         say("1sus$", "Read file with ", total_seqs, " sequences");
         say("1sfs$", "Average read length is ", avg_seq_len, "bp");
         goto restart_samples;
       }
     }
     else {
       if(strcmp(format,"fastq")==0) {
         ifstream ifs((string(file)));
         say("1ss$","reading ", file);
         FastqReader fastqReader(ifs);
         say("1sus$","Picking ", hX, " random sequences");
         References references = References(fastqReader, hX, k);
         say("1s$","Started counting");
         KmerCounter counter = KmerCounter(references, fastqReader, string(cntfile));
         mates = new int[hX];
         counter.getCounts(mates);
{'fraction_non_alphanumeric': 0.15301578081130796, 'fraction_numerical': 0.014328589408461613, 'mean_word_length': 3.3871310256954765, 'pattern_counts': {'":': 0, '<': 414, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1720959', 'n_tokens_mistral': 7462, 'n_tokens_neox': 7240, 'n_words': 2142}