/* See {dm_test_tools.h} */
/* Last edited on 2008-04-19 23:15:39 by stolfi */

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

#define _GNU_SOURCE
#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <jsrandom.h>

#include <affirm.h>
#include <bool.h>

#include <msm_basic.h>
#include <msm_seq_desc.h>
#include <msm_rung.h>
#include <msm_double_vec.h>
#include <msm_cand_refine.h>
#include <msm_ps_tools.h>
#include <msm_test_tools.h>

#include <dm_test_tools.h>
#include <dm_nucleic.h>
#include <dm_score.h>
#include <dm_seq.h>

msm_cand_t dm_test_tools_get_optimum_pairing
  ( dm_score_args_t *sc, 
    dm_seq_t *xp, 
    dm_seq_t *yp, 
    int delta, 
    msm_dyn_tableau_t *tb,
    int maxIter
  ) 
  { 
    auto double step_score
      ( msm_seq_desc_t *xd, 
        msm_seq_desc_t *yd,
        msm_rung_t g0,
        msm_rung_t g1
      );
    /* Assumes that {xd} and {yd} are the descriptors of {xp} and {yp},
      and tthat {g0} and {g1} are rungs between the 
      two sequences.  Returns the score of the step {g0->g1}. */
    
    return  msm_test_tools_get_optimum_pairing
      ( &(xp->sd), &(yp->sd), delta, &step_score, tb, maxIter );

    double step_score
      ( msm_seq_desc_t *xd, 
        msm_seq_desc_t *yd,
        msm_rung_t g0,
        msm_rung_t g1
      )
      { (void)msm_seq_desc_same_seq(xd, &(xp->sd), TRUE); 
        (void)msm_seq_desc_same_seq(yd, &(yp->sd), TRUE); 
        return dm_score_step(sc, g0, g1, xp, yp);
      }
  }

void dm_test_tools_seq_write_and_plot_named
  ( dm_seq_t *seq,
    char *title,
    char *name,
    char *tag
  )
  {
    fprintf(stderr, "sequence = "); 
    msm_seq_desc_write(stderr, " ( ", &(seq->sd), 1, 1, 4, " )"); 
    fprintf(stderr, "  samples = %4d  nsub = %2d\n", seq->dv.ne, seq->nsub);

    /* Check whether the sample vector is compatible with the seq descriptor: */
    
    int nsmp = seq->dv.ne;
    int npos = (seq->nsub*(nsmp-1)+1);
    demand(dm_seq_num_positions(seq) == npos, "inconsistent sequence descriptor");
    
    /* Ignore the {seq} for now: */
    dm_seq_write_named(seq, name, tag);
    
    /* Plotting: */
    msm_ps_tools_t *mps = msm_ps_tools_new_graph
      ( NULL, name, tag, 
        /*hGraphSize:*/ 165.0, 
        /*vGraphSize:*/  25.0, 
        /*scaleL:*/ FALSE, /*titleL:*/ FALSE, 
        /*scaleR:*/ FALSE, /*titleR:*/ FALSE,
        /*scaleB:*/ FALSE, /*titleB:*/ FALSE,
        /*scaleT:*/ FALSE, /*titleT:*/ FALSE, 
        /*fontSize:*/ dm_test_tools_FONT_SIZE,
        /*maxLabChars:*/ 5,
        /*mrg:*/ 1.0
      );
      
    /* Repack the samples as a single vector {smp[]}, for {msm_ps_tools_draw_graphs}: */
    int nplt = ( nsmp); /* Number of samples for plotting. */
    double smp[dm_CHANNELS*nplt]; 
    int i, c;
    for (i = 0; i < nplt; i++)
      for (c = 0; c < dm_CHANNELS; c++)
        { smp[nplt*c + i] = dm_seq_get_sample(seq, i, c); }

    msm_ps_tools_draw_graphs
      ( mps, dm_CHANNELS, nsmp, 
        NULL, smp, +INF, -INF
      );
  }

void dm_test_tools_seq_multi_write_and_plot_named
  ( dm_seq_t seq[],
    int maxLevel,
    char *title,
    char *name,
    char *tag
  )
  { int level;
    for (level = 0; level <= maxLevel; level++)
      { char *titlei = NULL; 
        asprintf(&titlei, "%s level %02d", title, level);
        char *tagi = NULL; 
        asprintf(&tagi, "%s-%02d", tag, level);
        dm_test_tools_seq_write_and_plot_named(&(seq[level]),titlei,name,tagi);
        free(titlei); free(tagi);
      }
  }

#define dm_rung_huge (msm_rung_t){{ INT_MAX, INT_MAX }}

void dm_test_tools_make_seq_pair
  ( char *borg, 
    double mutProb, 
    double delProb, 
    bool_t circ,
    dm_seq_id_t xid, 
    char *xtag,
    dm_seq_t *xP,
    dm_seq_id_t yid, 
    char *ytag,
    dm_seq_t *yP,
    msm_pairing_t **prP,
    char *outName
  )
  { /* Generate two mutated copies of {b}, save the indices of copied chars: */
    char *xdna, *ydna; /* The mutated copies. */
    msm_rung_vec_t xgv, ygv;  /* Index pairs for letters that were copied. */
    dm_nucleic_string_mutate(borg, TRUE, TRUE, mutProb, delProb, &xdna, &xgv);
    dm_nucleic_string_mutate(borg, TRUE, TRUE, mutProb, delProb, &ydna, &ygv);
    /* Merge the rung vectors {xgv,ygv} by component 0 to create the pairing {*prP}: */
    int ngmax = (xgv.ne < ygv.ne ? xgv.ne : ygv.ne);
    msm_rung_vec_t gv = msm_rung_vec_new(ngmax);
    int ng = 0; /* Output rungs will be {gv[0..ng-1]}. */
    int kx = 0, ky = 0;
    while ((kx < xgv.ne) && (ky < ygv.ne))
      { msm_rung_t xg = (kx >= xgv.ne ? dm_rung_huge : xgv.e[kx]);
        msm_rung_t yg = (ky >= ygv.ne ? dm_rung_huge : ygv.e[ky]);
        if (xg.c[0] < yg.c[0]) 
          { kx++; }
        else if (xg.c[0] > yg.c[0])
          { ky++; }
        else
          { msm_rung_t g = (msm_rung_t){{ xg.c[1], yg.c[1] }};
            if (ng > 0) { (void)msm_rung_step_is_increasing(gv.e[ng-1], g, /*die*/ TRUE); }
            msm_rung_vec_expand(&gv, ng);
            gv.e[ng] = g; ng++; 
            kx++; ky++; 
          }
      }
    demand(ng > 0, "empty output pairing");
    if (circ)
      { /* Add the closing rung: */
        msm_rung_t gini = gv.e[0];
        int nx = strlen(xdna);
        int ny = strlen(ydna);
        msm_rung_t grep = (msm_rung_t){{ gini.c[0] + nx, gini.c[1] + ny }};
        (void)msm_rung_step_is_increasing(gv.e[ng-1], grep, /*die*/ TRUE);
        msm_rung_vec_expand(&gv, ng);
        gv.e[ng] = grep; ng++; 
      }
    msm_rung_vec_trim(&gv, ng);
    /* Generate pairing from rungs: */
    msm_pairing_t *pr = msm_pairing_from_rung_vec(&gv);
    /* Generate comment: */
    char *cmt = NULL;
    asprintf(&cmt, "Mutated mp = %8.6f dp = %8.6f\n", mutProb, delProb);
    (*xP) = dm_seq_from_nucleic_string(xid, xtag, cmt, xdna);
    (*yP) = dm_seq_from_nucleic_string(yid, ytag, cmt, ydna);
    (*prP) = pr;
    /* Write sequences to disk: */
    dm_test_tools_write_generated_sequence(xdna, xP, outName);
    dm_test_tools_write_generated_sequence(ydna, yP, outName);
    free(xdna); free(xgv.e);
    free(ydna); free(ygv.e);
  }

void dm_test_tools_write_generated_sequence(char *b, dm_seq_t *s, char *outName)
  { /* Write sequence to disk: */
    char *fileTag = NULL;
    asprintf(&fileTag, "-%s", s->sd.name);
    dm_nucleic_string_write_named(b, s->cmt, outName, fileTag);
    dm_seq_write_named(s, outName, fileTag);
    free(fileTag);
  }
  

msm_cand_vec_t dm_test_tools_generate_cands
  ( int nCands, 
    bool_t homo, 
    dm_seq_t *x, 
    dm_seq_t *y, 
    msm_pairing_t *p
  )  
  {
//     demand(! msm_pairing_is_circular(p), "pairing must be open");
//     demand(! x->sd.circ, "X seq must be open");
//     demand(! y->sd.circ, "Y seq must be open");
    
    /* Compute the number {nSegs} of pieces into which {x,y,p} will be split: */
    int nSegs = nCands + 1;
    
    /* Divide {p} into {nSegs} pieces {*pc[i]}, each with span {~cdSpan}: */
    msm_pairing_t *pc[nSegs];
    dm_test_tools_split_pairing(p, nSegs, pc);
    
    /* Generate the candidates: */
    msm_cand_vec_t cdv = msm_cand_vec_new(nCands);
    int icd;
    for (icd = 0; icd < nCands; icd++)
      { msm_pairing_t *pr; /* Raw pairing for candidate. */
        if (homo)
          { /* Create a candidate from segment {icd} of the pairing: */
            pr = msm_pairing_copy(pc[icd]);
          }
        else
          { /* Create a candidate that pairs two totally unrelated segments of {x,y}: */
            /* Choose a segment {xseg} in the X sequence: */
            int xseg = icd;
            assert((xseg >= 0) && (xseg < nSegs));
            /* Choose the Y segment {yseg != xseg} that is to be paired with it: */
            int yseg;
            yseg = abrandom(0,nSegs-2);
            if (yseg >= xseg) { yseg++; }
            assert((yseg >= 0) && (yseg < nSegs));
            assert(xseg != yseg);
            /* Create a candidate for these two segments, with trivial pairing: */
            pr = dm_test_tools_make_cross_pairing(x, pc[xseg], y, pc[yseg]);
          }
        /* Interpolate {pr} to obtain an atomic pairing {pi}: */
        msm_pairing_t *pi = msm_pairing_interpolate(pr);
        /* Turn pairing {pi} into a candidate: */
        cdv.e[icd] = msm_cand_from_pairing(&(x->sd), &(y->sd), pi, 0.0);
        msm_pairing_free(pr); 
      }
    return cdv;
  }
  
  

void dm_test_tools_split_pairing(msm_pairing_t *p, int nSegs, msm_pairing_t *pc[])
  {
    (void)msm_pairing_is_increasing(p, TRUE);
    
    /* Get the count of fundamental rungs: */
    int nRungs = msm_pairing_num_rungs(p);

    /* Get the total span of {p} on the X and Y sequences, not counting the loopback step: */
    int totSpan = msm_pairing_sub_span(p, 0, nRungs-1, -1);
    
    /* Choose the total span {cdSpan} of each piece (and hence of each candidate): */
    int cdSpan = totSpan/nSegs;
    
    fprintf(stderr, "splitting a pairing with %d rungs into %d segments\n", nRungs, nSegs);
    int k = 0; /* Index of current rung in pairing {p}. */
    int iseg;
    for (iseg = 0; iseg < nSegs; iseg++)
      { assert(k < nRungs);
        int kini = k;
        while ((k < nRungs) && (msm_pairing_sub_span(p, kini, k, -1) <= cdSpan)) { k++; }
        int kfin = k - 1;
        assert(kini <= kfin);

        /* Show the selected segment"*/
        fprintf(stderr, "  segment %d = p[%d..%d]", iseg, kini, kfin);
        msm_rung_t gini = msm_pairing_get_rung(p, kini);
        msm_rung_t gfin = msm_pairing_get_rung(p, kfin);
        fprintf(stderr, "  spans X[%d..%d] (%d)", gini.c[0], gfin.c[0], gfin.c[0]-gini.c[0]+1);
        fprintf(stderr, " + Y[%d..%d] (%d)", gini.c[1], gfin.c[1], gfin.c[1]-gini.c[1]+1);

        pc[iseg] = msm_pairing_sub_copy(p, kini, kfin);
        assert(msm_rung_equal(gini, msm_pairing_get_rung(pc[iseg], 0)));
        assert(msm_rung_equal(gfin, msm_pairing_get_rung(pc[iseg], kfin-kini)));
        fprintf(stderr, " = (%d)\n", msm_pairing_span(pc[iseg], -1));
//         assert(! msm_pairing_is_circular(pc[iseg]));
      }
  }

msm_pairing_t *dm_test_tools_make_cross_pairing
  ( dm_seq_t *x, 
    msm_pairing_t *xp, 
    dm_seq_t *y, 
    msm_pairing_t *yp
  )
  {
    /* Obtain the initial and final rungs {gini,gfin} of the desired pairing: */
    msm_rung_t gini, gfin;
    int j;
    for (j = 0; j < 2; j++)
      { /* Grab the whole sequence {s} on side {j} and its position count {ns}: */
        dm_seq_t *s = (j == 0 ? x : y);
//         demand(! s->sd.circ, "seq must be open");
        int ns = s->sd.npos;
        /* Grab the pairing {p} that defines the segm on {s}, and its rung count {ng}: */
        msm_pairing_t *p = (j == 0 ? xp : yp);
//         demand(! msm_pairing_is_circular(p), "pairing must be open");
        int ng = msm_pairing_num_rungs(p);
        /* Set the indices of first and last spanned positions on {s}: */
        gini.c[j] = msm_pairing_get_rung(p, 0).c[j];
        demand((gini.c[j] >= 0) && (gini.c[j] < ns), "initial rung out of bouds");
        gfin.c[j] = msm_pairing_get_rung(p, ng-1).c[j];
        demand((gfin.c[j] >= 0) && (gfin.c[j] < ns), "final rung out of bounds");
      }
    /* Create a trivial pairing from {gini,gfin}: */
    msm_rung_vec_t gvt = msm_rung_vec_new(2);
    gvt.e[0] = gini; 
    gvt.e[1] = gfin;
    return msm_pairing_from_rung_vec(&gvt);
  }