#define PROG_NAME "dm_compute_weights"
#define PROG_DESC "analysis of true/false candidate discrimination"
#define PROG_VERS "1.0"

/* Last edited on 2009-02-23 22:24:20 by stolfi */

#define dm_cand_filter_C_COPYRIGHT \
  "Copyright � 2006  by the State University of Campinas (UNICAMP)"

#define PROG_HELP \
  PROG_NAME " \\\n" \
  "  -seqA {ID_A} {NAME_A} {FILENAME_A} \\\n" \
  "  -seqB {ID_B} {NAME_B} {FILENAME_B} \\\n" \
  "  -candF {CAND_FILE_F} \\\n" \
  "  -candT {CAND_FILE_T} \\\n" \
  "  [-minCandSize {MINCANDSIZE} -maxCandSize {MAXCANDSIZE} ] \\\n" \
  "  [-adjustWeight {MAXITERATIONS} ] \\\n" \
  "  -outPrefix {OUT_PREFIX}"

#define PROG_INFO \
  "NAME\n" \
  "  " PROG_NAME " - " PROG_DESC "\n" \
  "\n" \
  "SYNOPSIS\n" \
  "  " PROG_HELP "\n" \
  "\n" \
  "DESCRIPTION\n" \
  "  This program  reads a candidate file {CAND_FILE_T} containing correct pairings,  " \
  " another one {CAND_FILE_F} containing random parings, " \
  " and the corresponding sequences {FILENAME_A} and {FILENAME_B}, and computes " \
  " a estimate for the weights for dm_match for such candidates.\n " \
  "\n" \
  "OUTPUT FILES\n" \
  "  All output files will have names starting with {OUT_PREFIX}.\n" \
  "\n" \
  "OPTIONS\n" \
  "\n" \
  argparser_help_info_HELP_INFO "\n" \
  "SEE ALSO\n" \
  "  dm_seq_filter(1)\n" \
  "\n" \
  "AUTHOR\n" \
  "  This program was created on 12/nov/2012 by J. Stolfi and R. Saracchini.\n" \
  "WARRANTY\n" \
  argparser_help_info_NO_WARRANTY "\n" \
  "\n" \
  "RIGHTS\n" \
  "  " dm_cand_filter_C_COPYRIGHT ".\n" \
  "\n" \
  argparser_help_info_STANDARD_RIGHTS  

#include <dm_basic.h>
#include <dm_cand.h>
#include <dm_classify.h>
#include <dm_seq.h>
#include <dm_nucleic.h>
#include <dm_score.h>
#include <dm_filter.h>

#include <msm_rung.h>
#include <msm_pairing.h>
#include <msm_cand.h>
#include <msm_cand_vec.h>
#include <msm_cand_refine.h>
#include <msm_dyn.h>
#include <msm_multi.h>

#include <affirm.h>
#include <wt_table.h>
#include <jsfile.h>
#include <argparser.h>

#include <math.h>
#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

typedef struct options_t 
  { /*Sequence parameters*/
    msm_seq_id_t idA,idB;
    char* nameA, *nameB;
    char* filenameA, *filenameB;
    /*Candidate parameters*/
    char* candFileF;
    char* candFileT;
    char *outPrefix;       /* Output file name prefix (minus extensions). */
    int minCandSize,maxCandSize;
    bool_t adjustWeight;
    int maxIter;
    bool_t writeSVM;
  } options_t;
  
int main(int argc, char**argv);

options_t *parse_options(int argc, char**argv);
  /* Parses the command line options, packs 
    them into a {options_t} record. */

  
msm_cand_vec_t readCDV(char* filename);
msm_cand_vec_t readCDV(char* filename){
  FILE* arq_cdv = open_read(filename,TRUE);
  msm_cand_vec_t cdv =  msm_cand_vec_read(arq_cdv);
  fclose(arq_cdv);
  return cdv;
}

dm_seq_t readEQS(char* filename);
dm_seq_t readEQS(char* filename){
  FILE* arq_seq = open_read(filename,TRUE);
  dm_seq_t seq = dm_seq_read(arq_seq,1);
  fclose(arq_seq);
  return seq;
}

void trimCDV(msm_cand_vec_t* orig, msm_cand_vec_t* dest, int min, int max);
void trimCDV(msm_cand_vec_t* orig, msm_cand_vec_t* dest, int min, int max){
    int N = orig->ne;
//     msm_cand_vec_expand(dest, N);
    int count = 0;
    int i;
    for(i = 0; i < N; i++){
      msm_cand_t* cd = &(orig->e[i]);
      msm_pairing_t* pr = cd->pr;
      int nr = msm_pairing_num_rungs(pr);
      if( (nr >= min) && (nr <= max)){
	dest->e[count] = *cd;
	count++;
      }
    }
    msm_cand_vec_trim(dest,count);
}


int main(int argc, char**argv)
  { 
    options_t *o = parse_options(argc, argv);
    
    msm_cand_vec_t cdvForig = readCDV(o->candFileF);
    msm_cand_vec_t cdvTorig = readCDV(o->candFileT);
    
    msm_cand_vec_t cdvF = msm_cand_vec_new(cdvForig.ne);
    msm_cand_vec_t cdvT = msm_cand_vec_new(cdvTorig.ne);
    
    
    trimCDV(&cdvForig,&cdvF,o->minCandSize,o->maxCandSize);
    fprintf(stderr,"False candidates - selected %d of %d\n",cdvF.ne, cdvForig.ne);
    trimCDV(&cdvTorig,&cdvT,o->minCandSize,o->maxCandSize);
    fprintf(stderr,"True  candidates - selected %d of %d\n",cdvT.ne, cdvTorig.ne);
    
    
    dm_seq_t seq_a = readEQS(o->filenameA);
    dm_seq_t seq_b = readEQS(o->filenameB);
    
    double wEql,  wDif,  wBrk,  wSkp, wUnp, K;
    
    dm_classify_determine_weights(&cdvT, &cdvF,&seq_a, &seq_b, &wEql, &wDif, &wBrk,&wSkp,&wUnp,&K);
    
    
    double vec[5] = { wEql,wDif,wBrk,wSkp,wUnp };
    double max = wEql;
    int i;
    for(i = 0; i < 5; i++){ if(max < vec[i]){ max = vec[i];} }
    if(max > 100.0){
      double ratio = 100.0/max;
      wEql = wEql*ratio;
      wDif = wDif*ratio;
      wBrk = wBrk*ratio; 
      wSkp = wSkp*ratio;
      wUnp = wUnp*ratio;
    }
    fprintf(stderr,"Computed weights\nEQL: %9.6lf\nDIF: %9.6lf\nBRK: %9.6lf\nSKP: %9.6lf\nUNP: %9.6lf\nK: %9.6lf\n",wEql,wDif,wBrk,wSkp,wUnp,K);
//     dm_classify_adjust_weights(&cdvT, &cdvF,&seq_a, &seq_b, &wEql, &wDif, &wBrk,&wSkp,&wUnp,&K);
//     fprintf(stderr,"Adjusted weights\nEQL: %9.6lf\nDIF: %9.6lf\nBRK: %9.6lf\nSKP: %9.6lf\nUNP: %9.6lf\nK: %9.6lf\n",wEql,wDif,wBrk,wSkp,wUnp,K);
    
    if(o->writeSVM){
      char* arq_svm_filename = NULL;
      asprintf(&arq_svm_filename,"%s-svm.txt",o->outPrefix);
      FILE* arq_svm = open_write(arq_svm_filename,TRUE);
      dm_classify_write_svm_file(&cdvT,&cdvF,&seq_a, &seq_b, arq_svm);
      fclose(arq_svm);
    }
    
    
    if(o->adjustWeight){
      fprintf(stderr,"Adjusting weights - %d iterations\n",o->maxIter);
       double nEql,  nDif,  nBrk,  nSkp, nUnp, nK;
       dm_classify_adjust_weights(&cdvT,&cdvF,&seq_a,&seq_b,o->maxIter, wEql,wDif,wBrk,wSkp,wUnp,K, &nEql, &nDif, &nBrk,&nSkp,&nUnp,&nK);
       wEql = nEql;
       wDif = nDif;
       wBrk = nBrk;
       wSkp = nSkp;
       K = nK;

    }
      
    char* outfilename = NULL;
    asprintf(&outfilename,"%s-wgt.txt",o->outPrefix);
    FILE* arq_weight = open_write(outfilename,TRUE);
    fprintf(arq_weight,"data %9.6lf %9.6lf move %9.6lf %9.6lf %9.6lf\n",wEql,wDif,wBrk,wSkp,wUnp);
    fclose(arq_weight);

    fprintf(stderr,"Computing distances \n");
    double* distT = dm_classify_compute_distances_to_weight_plane(&cdvT, &seq_a, &seq_b,wEql,wDif,wBrk,wSkp,wUnp,K);
    double* distF = dm_classify_compute_distances_to_weight_plane(&cdvF, &seq_a, &seq_b,wEql,wDif,wBrk,wSkp,wUnp,K);
    
    char* outfilename_stats = NULL;
    asprintf(&outfilename_stats,"%s-sts.txt",o->outPrefix);
    FILE* arq_weight_stats = open_write(outfilename_stats,TRUE);
    
    fprintf(arq_weight_stats,"#TRUE CANDIDATES\n\n");
    int fn,tp;
    fn = 0; tp = 0;
    for(i = 0; i < cdvT.ne; i++){
      double eql,dif,skp,brk,unp;
      int span,nrungs;
      dm_cand_count_cand_stats(&(cdvT.e[i]),&seq_a, &seq_b, &eql, &dif,&brk,&skp,&unp,&span,&nrungs);
      fprintf(arq_weight_stats,"[%d] : %+9.6lf - %9.6lf %9.6lf %9.6lf %9.6lf %9.6lf\n",i,distT[i],eql,dif,brk,skp,unp);
      if(distT[i] >= 0) tp++;
      else{ fn++; }
    }
    int tn,fp;
    tn = fp = 0;
    fprintf(arq_weight_stats,"\n#FALSE CANDIDATES\n\n");
    for(i = 0; i < cdvF.ne; i++){
      double eql,dif,skp,brk,unp;
      int span,nrungs;
      dm_cand_count_cand_stats(&(cdvF.e[i]),&seq_a, &seq_b, &eql, &dif,&brk,&skp,&unp,&span,&nrungs);
      fprintf(arq_weight_stats,"[%d] : %+9.6lf - %9.6lf %9.6lf %9.6lf %9.6lf %9.6lf\n",i,distF[i],eql,dif,brk,skp,unp);
      if(distF[i] < 0){
	tn++;
      }else{
	fp++;
      }
    }
    fprintf(arq_weight_stats,"\n#SUMMARY\n");
    fprintf(arq_weight_stats,"True  Positives %04d (%3.2lf%%)\n",tp,100.0*tp/(double)cdvT.ne);
    fprintf(arq_weight_stats,"False Negatives %04d (%3.2lf%%)\n",fn,100.0*fn/(double)cdvT.ne);
    fprintf(arq_weight_stats,"True  Negatives %04d (%3.2lf%%)\n",tn,100.0*tn/(double)cdvF.ne);
    fprintf(arq_weight_stats,"False Positives %04d (%3.2lf%%)\n",fp,100.0*fp/(double)cdvF.ne);
    fclose(arq_weight_stats);
    
    return 0;
  }


options_t *parse_options(int argc, char**argv)
  { 
    options_t *o = (options_t *)notnull(malloc(sizeof(options_t)), "no mem");
    
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
     
       
    argparser_get_keyword(pp, "-seqA");
    o->idA = argparser_get_next_int(pp,0,INT_MAX);
    o->nameA = argparser_get_next_non_keyword(pp);
    o->filenameA = argparser_get_next_non_keyword(pp);
    
    argparser_get_keyword(pp, "-seqB");
    o->idB = argparser_get_next_int(pp,0,INT_MAX);
    o->nameB = argparser_get_next_non_keyword(pp);
    o->filenameB = argparser_get_next_non_keyword(pp);
    
    argparser_get_keyword(pp,"-candF");
    o->candFileF = argparser_get_next_non_keyword(pp);
    
    argparser_get_keyword(pp,"-candT");
    o->candFileT = argparser_get_next_non_keyword(pp);
    
    argparser_get_keyword(pp,"-outPrefix");
    o->outPrefix = argparser_get_next(pp);
    
    o->writeSVM =  argparser_keyword_present(pp,"-writeSVM");
    
    o->minCandSize = 0 ;
    if(argparser_keyword_present(pp,"-minCandSize")){
      o->minCandSize =  argparser_get_next_int(pp,0,INT_MAX);
    }
    
    o->maxCandSize = INT_MAX ;
    if(argparser_keyword_present(pp,"-maxCandSize")){
      o->maxCandSize =  argparser_get_next_int(pp,0,INT_MAX);
    }
    
    if(o->maxCandSize < o->minCandSize){
      fprintf(stderr,"minCandSize is bigger than maxCandSize !\n");
      assert(FALSE);
    }
    
    o->adjustWeight = argparser_keyword_present(pp,"-adjustWeight");
    if(o->adjustWeight){
      o->maxIter =  argparser_get_next_int(pp,1,INT_MAX);
    }
    
    argparser_finish(pp);
    
    return o;
  }