#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h> // for mkdir
#include <inttypes.h>
#include <stdint.h>
#include <algorithm> // for max/min
#include <vector> // for sorting_kmers
#include <sys/time.h>

#define NNKS 4 // default minimal abundance for solidity

int max_memory; // the most memory one should alloc at any time, in MB
int order=0; // in kissnp, don't change it, it should be 0

#include "../minia/Bank.h"
#include "../minia/Hash16.h"
#include "../minia/Set.h"
#include "../minia/Pool.h"
#include "../minia/Bloom.h"
#include "../minia/Debloom.h"
#include "../minia/Utils.h"
#include "../minia/SortingCount.h"
#include "../minia/Terminator.h"
#include "../minia/Kmer.h"
#include "../minia/rvalues.h"
//#include "SNP.h"
//#include "Kmer_for_kissnp2.h"
#include "Starter_Bank.h"
#include "commons.h"

extern int size_seeds;
extern string prefix_trashable;
using namespace std;

/////////////////////// FOR MINIA //////////////////
extern int sizeKmer;
int64_t genome_size;
Bloom * bloo1;
int threshold;
BinaryBank * SolidKmers;
BranchingTerminator * terminator;


extern unsigned int nbits_nbseeds;
extern uint64_t  mask_nbseed ;
extern uint64_t  mask_offset_seed;
////////////////////////////////////////////////////


char * getVersion(){
    return (char *)"2.0.5.deb - Copyright INRIA - CeCILL License";
}

void print_usage_and_exit(char * name){
    fprintf (stderr, "NAME\nmapsembler, version %s\n", getVersion());
    fprintf (stderr, "\nSYNOPSIS\n%s <starters.fasta> <readsC1.fasta> [<readsC2.fasta> [<readsC3.fasta] ...] [-E] [-t extension_type] [-q value] [-k value] [-c value] [-d authorized_distance] [-e error_threshold] [-g value] [-i index_name] [-o name] [-m name] [-h]\n", name);
    fprintf (stderr, "\nDESCRIPTION\n");
    fprintf (stderr, "\t Mapsembler2 is a targeted assembly software. It takes as input a set of NGS raw reads (fasta or fastq, gzipped or not) and a set of input sequences (starters). It first determines if each starter is read-coherent, e.g. whether reads confirm the presence of each starter in the original sequence. Then for each read-coherent starter, Mapsembler2 outputs its sequence neighborhood as a linear sequence or as a graph, depending on the user choice.\n");
    
    fprintf (stderr, "\nOPTIONS\n");
    //fprintf (stderr, "\t -c Write the coverage of each position of the sequences. The fasta file will have 3 lines (instead of one) starting with \'>\' before the sequences themselves\n");
    fprintf (stderr, "\t -E Extend only: avoid the mapping+substarter generation phase\n");
    fprintf (stderr, "\t -t extension_type. Default: 1\n");
    fprintf (stderr, "\t    1: a strict sequence: any branching stops the extension\n");
    fprintf (stderr, "\t    2: a consensus sequence: contiging approach\n");
    fprintf (stderr, "\t    3: a strict graph: any branching is conserved in the graph\n");
    fprintf (stderr, "\t    4: a consensus graph: \"small\" polymorphism is merged, but \"large\" structures are represented\n");
    fprintf (stderr, "\t -q size_seed: will use seeds of length size_seed during the mapping process. Default: 25, value should be in [5-31]\n");
    fprintf (stderr, "\t -k size_kmers: Size of the k-mers used duriung the extension phase Default: 31. Accepted range, depends on the compilation (make k=42 for instance) \n");
    fprintf (stderr, "\t -c min_coverage: a sequence is covered by at least min_coverage coherent reads. Default: 2\n");
    fprintf (stderr, "\t -d authorized_distance: a substarter is distant by at most authorized_distance substitutions. Default: 2\n");
    fprintf (stderr, "\t -e error_threshold: a nucleotide is corrected if occurs less than error_threshold times. Default: 2\n");
    fprintf (stderr, "\t -g estimated_genome_size: estimation of the size of the genome whose reads come from. \n \t    It is in bp, does not need to be accurate, only controls memory usage. Default: 3 billion\n");
    fprintf (stderr, "\t -n node_len: limit max of nodes length. Default: 1000000\n");
    fprintf (stderr, "\t -l graph_max_depth: limit max of graph depth.Default: 10000\n");
    fprintf (stderr, "\t -i index_name: stores the index files in files starting with this prefix name. Can be re-used latter. Default: \"index\"\n");
    fprintf (stderr, "\t    IF THE FILE \"index_name.bloom\" EXISTS: the index is not re-created \n");
    fprintf (stderr, "\t -o file_name_prefix: where to write outputs. Default: \"res_mapsembler\" \n");
    fprintf (stderr, "\t -p search_mod: kind of prosses LARGEUR or PROFONDEUR. Default: LARGEUR \n"); 
    fprintf (stderr, "\t -m file_name: write in file \"file_name\" the reads mapped on starters\n");
    fprintf (stderr, "\t -h prints this message and exit\n");
    exit(0);
}



int main(int argc, char *argv[])
{
    
    setvbuf(stdout,(char*)NULL,_IONBF,0); // disable buffering for printf()
    size_seeds=25;
    int min_coverage=2; // minimal number of reads per positions extending the initial fragment
    int error_threshold=2; //comment todo
    int authorized_distance=2; // comment todo et comrendre diff
    char extension_type=1; // kind of extensions
    int max_graph_depth = 10000; //detph graph limit
    int max_nodes = 1000000; //node length limit
    int process_type = 1; //kind of process : 1 == LARGEUR, 2 == PROFONDEUR
    parcours_t search_mode = LARGEUR;
    bool export_mapped_reads=false;
    char * file_read_mapped;
    bool make_mapping = true;
    genome_size=3000000000;
    if(argc<2){
        print_usage_and_exit(argv[0]);
    }
    char * toCheck_file =strdup(argv[1]);
    FILE * out_results = stdout;
    char * index_file_prefix = new char [4096];
    char * res_prefix = new char [4096];
    strcpy(index_file_prefix,"index");
    strcpy(res_prefix,"res_mapsembler");
    
    ////////////////////////////////// GETTING THE OPTIONS /////////////////////////////////////
    sizeKmer=31;
    // GET ALL THE READ FILES
    // find the number of read sets
    int number_of_read_sets=0;
    while(number_of_read_sets+2<argc && argv[number_of_read_sets+2][0]!='-') number_of_read_sets++;
    char ** reads_file_names = (char **) malloc(sizeof(char *)*number_of_read_sets);
    
    // copy the read set file names
    number_of_read_sets=0;
    while(number_of_read_sets+2<argc && argv[number_of_read_sets+2][0]!='-'){
        reads_file_names[number_of_read_sets]=strdup(argv[number_of_read_sets+2]);
        number_of_read_sets++;
    }
    
    while (1)
    {
        int temoin = getopt (argc-number_of_read_sets-1, &argv[number_of_read_sets+1], "t:c:e:d:x:y:o:f:k:q:g:i:m:h:E");
        if (temoin == -1){
            break;
        }
        switch (temoin)
        {
            case 'E':
                make_mapping=false;
                printf(" Mapsembler will extend the starters, not checking their read-coherency nor generating substarters\n");
                break;
            case 't':
                extension_type=atoi(optarg);
                printf(" Extension_Type=%d\n", extension_type);
                break;
            case 'c':
                min_coverage=atoi(optarg);
                printf(" Will consider as read coherent if coverage is at least %d\n", min_coverage);
                break;
            case 'e':
                error_threshold=atoi(optarg);
                printf(" Will use error_threshold=%d\n", error_threshold);
                break;
            case 'd':
                authorized_distance=atoi(optarg);
                printf(" Will use authorized_distance=%d\n", authorized_distance);
                break;
            case 'g':
                genome_size=atoll(optarg);
                printf(" Will use genome_size=%llu\n", genome_size); //TODO format
                break;
	    case 'x':
               	max_nodes=atoi(optarg);
		printf(" Will use node_len=%d\n", max_nodes); //TODO format
                break;
	    case 'y':
		max_graph_depth=atoi(optarg);
                printf(" Will use max_graph_depth=%d\n", max_graph_depth); //TODO format
                break;
            case 'o':
                strcpy(res_prefix,optarg);
                printf(" Will output results in file prefixed with %s\n", optarg);
                break;
	    case 'f':
		process_type=atoi(optarg);
		if(process_type == 1){
		     	search_mode = LARGEUR; 
		}else if(process_type == 2){
			search_mode = PROFONDEUR;    
		}
		printf(" Will use process_type %d\n", process_type);
                break;
            case 'm':
                file_read_mapped=strdup(optarg);
                export_mapped_reads=true;
                printf(" Will output mapped reads in file: %s\n", file_read_mapped);
                break;
            case 'i':
                strcpy(index_file_prefix, optarg);
                printf(" Will store read index in files prefixed with: %s\n", index_file_prefix);
                break;
                
                
            case 'q':
                size_seeds=atoi(optarg);
                
                if(size_seeds<5){
                    fprintf(stderr,"%d too small to be used as seed length. kisSnpCheckReads continues with seeds of length 5\n",size_seeds);
                    size_seeds=5;
                }
                if(size_seeds>31){
                    fprintf(stderr,"%d too big to be used as seed length. kisSnpCheckReads continues with seeds of length 31\n",size_seeds);
                    size_seeds=31;
                }
                
                    printf(" Will use seeds of length %d for mapping\n", size_seeds);
                break;
                
            case 'k':
                sizeKmer=atoi(optarg);
                // let's make k even for now, because i havnt thought of how to handle palindromes (dont want to stop on them)
                if (sizeKmer%2==0) //TODO: verif comptaibilité minia
                {
                    sizeKmer-=1;
                    printf("Need odd kmer size to avoid palindromes. I've set kmer size to %d.\n",sizeKmer);
                }
                if (sizeKmer>(sizeof(kmer_type)*4))
                {
                    printf("Max kmer size on this compiled version is %d\n",(int)sizeof(kmer_type)*4);
                    exit(1);
                }
                
                printf(" Will use kmers of length %d for extensions\n", sizeKmer);
                break;
            case 'h':
                print_usage_and_exit(argv[0]);
                //	case 'v':
                //	  printf(" VERBOSE mode\n");
                //	  verbose=1;
                //	  break;
            case '-':
                if(strcmp("version", optarg)==0){
                    printf("kissReads version %s\n", getVersion());
                    exit(0);
                }
                printf(" what next ? %s\n",optarg);
                break;
            default:
                //			printf ("Unknown option %c\n", temoin);
                print_usage_and_exit(argv[0]);
        }
    }

    if ( argc  - optind <2)
    {
        print_usage_and_exit(argv[0]);
    }
    
    if(!make_mapping && export_mapped_reads){
        fprintf(stderr,"Warning options -E and -m are incompatible. I won't output mapped reads (disabling -m option)");
        export_mapped_reads=false;
    }
    
    
    sprintf(prefix,"%s_k_%d_c_%d", index_file_prefix, sizeKmer, min_coverage);
    
    char final_res_prefix [4096];
    sprintf(final_res_prefix,"%s_k_%d_q_%d_c_%d_t_%d", res_prefix, sizeKmer, size_seeds, min_coverage, (int)extension_type);
    
    init_static_variables();
    Bank *reads = new Bank(reads_file_names,number_of_read_sets);
    Starter_Bank *starter_bank;
    Extension_Bank * extensions;
    
    
    ///// USED BOTH FOR mapping and extensions
    nbits_nbseeds = 8*sizeof(uint64_t)- NBITS_OFFSET_SEED ;
    mask_nbseed  = ( 1ULL <<  (uint64_t) nbits_nbseeds  ) -1 ;
    mask_offset_seed = (1ULL << (NBITS_OFFSET_SEED)) -1 ;
    
    
   ////////////////////////////////// READ COHERENCE /////////////////////////////////////////
    if(make_mapping){
        printf("indexing starters\n");
        Fragment_Bank *fragment_bank= new Fragment_Bank(); // loading the fragment sequences
        fragment_bank->load_fragments(argv[1]);
        fragment_bank->index_fragments_and_init_vector();
        
        
        printf("map reads, checking the read coherency and discovering subfragments\n");
        fragment_bank->map_reads(reads, authorized_distance);
//        fragment_bank->show_alignments();
        if(export_mapped_reads) fragment_bank->output_mapped_reads(file_read_mapped);
        
        
        starter_bank = new Starter_Bank(fragment_bank);
        delete(fragment_bank);
        printf("generate substarters\n");
        extensions = starter_bank->generate_all_substarters(min_coverage, authorized_distance, error_threshold, sizeKmer);
        
        if(extensions->all_extensions.size()==0){
            printf("No starter was read coherent, we stop here\n");
            exit(0);
        }
    }
    else{
        Fragment_Bank *fragment_bank= new Fragment_Bank(); // loading the fragment sequences
        fragment_bank->load_fragments(argv[1]);
        starter_bank = new Starter_Bank(fragment_bank);
        extensions = starter_bank->generate_all_substarter_nomapping(sizeKmer);
    }

  
    ////////////////////////////////// EXTENSIONS //////////////////////////////////////////
    extensions->set_extension_type(extension_type);
    printf("Indexing reads, using minia approach, generating file prefixed with \"%s\"\n", prefix);
    
  
    
    kmerMask=(((kmer_type)1)<<(sizeKmer*2))-1;
    double lg2 = log(2);
      NBITS_PER_KMER = rvalues[sizeKmer][1];
    //NBITS_PER_KMER = log(16*sizeKmer*(lg2*lg2))/(lg2*lg2); // needed to process argv[5]
    // solidity
    nks =NNKS;
    nks = min_coverage;
    if (nks<=0) nks=1; // min abundance can't be 0
    int estimated_BL1 = max( (int)ceilf(log2f(genome_size * NBITS_PER_KMER )), 1);
    uint64_t estimated_nb_FP =  (uint64_t)(genome_size * 4 * powf(0.6,11)); // just indicative
    max_memory = max( (1LL << estimated_BL1)/8LL /1024LL/1024LL, 1LL );
    //  printf("estimated values: nbits Bloom %i, nb FP %lld, max memory %i MB\n",estimated_BL1,estimated_nb_FP,max_memory);
    
    
    
    
    
    // shortcuts to go directly to assembly using serialized bloom and serialized hash
    int START_FROM_SOLID_KMERS=0; // if = 0, construct the fasta file of solid kmers, if = 1, start directly from that file
    int LOAD_FALSE_POSITIVE_KMERS=0; // if = 0, construct the fasta file of false positive kmers (debloom), if = 1, load that file into the hashtable
    int NO_FALSE_POSITIVES_AT_ALL=0; // if = 0, normal behavior, if = 1, don't load false positives (will be a probabilistic de bruijn graph)
    
    //tester si le fichier prefix.bloom exist et faire:
    if( access( (string(prefix)+string(".debloom")).c_str(), F_OK ) != -1 &&
       access( (string(prefix)+string(".debloom2")).c_str(), F_OK ) != -1 &&
       access( (string(prefix)+string(".false_positive_kmers")).c_str(), F_OK ) != -1 &&
       access( (string(prefix)+string(".solid_kmers_binary")).c_str(), F_OK ) != -1){
        printf("THE INDEX %s.... were already created we use them\n", prefix);
        START_FROM_SOLID_KMERS=1;
        LOAD_FALSE_POSITIVE_KMERS=1;
    }
    else
        printf("CREATING THE index %s.... \n", prefix);
    
    
    //  fprintf (stderr,"taille cell %lu \n", sizeof(cell<kmer_type>));
    
    STARTWALL(0);
    // counter kmers, write solid kmers to disk, insert them into bloo1
    if (!START_FROM_SOLID_KMERS)
    {
        
        int max_disk_space = 0; // let dsk decide
        int verbose = 0;
        bool write_count = false;
        
        sorting_count(reads,prefix, max_memory, max_disk_space, write_count, verbose);
    }
    
    
    // debloom, write false positives to disk, insert them into false_positives
    if (! LOAD_FALSE_POSITIVE_KMERS)
    {
        debloom(order, max_memory);
    }
    
    bloo1 = bloom_create_bloo1((BloomCpt *)NULL);
    // load false positives from disk into false_positives
     false_positives = load_false_positives_cascading4();
    
    
    
    
    // load branching kmers
    //BinaryBank *SolidKmers = new BinaryBank(return_file_name(solid_kmers_file),sizeof(kmer_type),0);
    int LOAD_BRANCHING_KMERS=0;
    if (LOAD_BRANCHING_KMERS)
    {
        BinaryBank *BranchingKmers = new BinaryBank(return_file_name(branching_kmers_file),sizeof(kmer_type),false);
        terminator = new BranchingTerminator(BranchingKmers,SolidKmers, bloo1,false_positives);
        BranchingKmers->close();
 	delete BranchingKmers;
    }
    else{
        BinaryBank *SolidKmers = new BinaryBank(return_file_name(solid_kmers_file),sizeof(kmer_type),0);
        terminator = new BranchingTerminator(SolidKmers,genome_size, bloo1,false_positives);
        SolidKmers->close();
 	delete SolidKmers; 
    }
    
    
    printf("Generating %d extensions\n", (int) extensions->all_extensions.size());
    extensions->perform_extensions(prefix, max_graph_depth, max_nodes, search_mode);
    
    printf("Format results\n");
    starter_bank->format_results(extension_type, extensions, final_res_prefix);
    
    // CLEANING THE TEMP FILES
    string command;
//    command = string("rm -f ")+prefix_trashable.c_str()+string("*"); system((char *)command.c_str());
    //command = string("rm -f *.json_edges"); system((char *)command.c_str());
    //command = string("rm -f *.json_nodes"); system((char *)command.c_str());
    //command = string("rm -f *.json_substarters"); system((char *)command.c_str());
    
    STOPWALL(0,"Total");

    delete bloo1;
    delete terminator;
    delete starter_bank;
    delete extensions;
    delete reads;
    free(toCheck_file);
    for(int i=0;i<number_of_read_sets;i++) free(reads_file_names[i]); free(reads_file_names);
    
    printf("all done, results are in %s.%s\n", final_res_prefix, extension_type<3?"fasta":"json");
    
    return 0;
}