/* Last edited on 2023-02-22 12:29:30 by stolfi */

#include <pz_types.h>
#include <pz_proc.h>
#include <pz_match.h>
#include <pz_segment.h>
#include <pz_double_chain.h>

#include <nget.h>
#include <fget.h>
#include <filefmt.h>
#include <jsfile.h>

#include <math.h>
#include <stdio.h>
#include <libiberty.h>
#include <values.h>

pz_double_chain_t pz_double_chain_cut
  ( pz_double_chain_t *c, 
    unsigned start,
    unsigned length 
  )
  {
    pz_double_chain_t x = double_vec_new(length);
    pz_double_chain_do_cut(c, start, length, &x);
    return x;
  }

void pz_double_chain_do_cut
  ( pz_double_chain_t *c, 
    unsigned start,
    unsigned length, 
    pz_double_chain_t *x 
  )
  {
    int n = c->nel;
    int ini = (start % n);
    int i, j;
    for (i= 0, j = ini; i < length; i++, j++)
      { if (j > n) { j -= n; }
        x->el[i] = c->el[j];
      }
  }

void pz_double_chain_reverse ( pz_double_chain_t *c )
  {
    double aux;
    int n = c->nel, m = n/2;
    int i, j;
    for (i = 0, j = n-1; i < m; i++, j--)
      { aux = c->el[i]; c->el[i] = c->el[j]; c->el[j] = aux; }
  }

void pz_double_chain_complement ( pz_double_chain_t *c )
  {
    int n = c->nel, i;
    for (i = 0; i <= n-1 ; i++) { c->el[i] = - c->el[i]; }
  }

pz_double_chain_t pz_double_chain_extract_segment
  ( pz_double_chain_t *c, 
    pz_segment_t s, 
    bool_t curvature 
  )
  {
    affirm(c->nel == s.tot, "wrong tot size");
    double_vec_t t = double_vec_new(s.ns);
    pz_double_chain_do_extract_segment(c, s, curvature, &t);
    return t;
  }

void pz_double_chain_do_extract_segment
  ( pz_double_chain_t *c, 
    pz_segment_t s, 
    bool_t curvature, 
    pz_double_chain_t *t 
  )
  {
    affirm(c->nel == s.tot, "wrong tot size");
    pz_double_chain_do_cut(c, s.ini, s.ns, t);
    if (s.rev)
      { pz_double_chain_reverse(t);
        if (curvature) { pz_double_chain_complement(t); }
      }
  }

double pz_double_chain_adjust_unit ( double givenUnit, pz_double_chain_t *c )
  {
    int n = c->nel;
    double avg = 0, dev = 0, big = 0;
    /* Compute average and max-abs value: */
    if (n > 0 )
      { int i;
        for (i = 0; i < n; i++)
          { double ci = c->el[i];
            avg = avg + ci; 
            if (fabs(ci) > big) { big = fabs(ci); }
          }
        avg = avg/(double)n;
      }
    /* Compute standard deviation: */
    if (n > 1)
      { int i;
        for (i = 0; i < n; i++)
          { double d = c->el[i] - avg;
            dev += d*d;
          }
        dev = sqrt(dev/((double)n-1));
      }
    return pz_adjust_unit(givenUnit, dev, big);
  }

#define FileVersion "97-10-29"

void pz_double_chain_write
  ( FILE *wr, 
    char *cmt, 
    pz_double_chain_t *c, 
    double stride, 
    double unit 
  )
  {
    auto int pz_double_chain_to_int( double x );
    int pz_double_chain_to_int( double x )
      {
        affirm(x > ((double)MININT), "out of range");
        affirm(x < ((double)MAXINT), "out of range");
        return pz_round(x);
      }
      
    filefmt_write_header(wr, "pz_double_chain_t", FileVersion);
    int ind = 0; /* Comment indentation. */
    filefmt_write_comment(wr, cmt, ind, '|');
    fprintf(wr, "stride = %g\n", stride);
    fprintf(wr, "unit = %g\n", unit);
    fprintf(wr, "samples = %d\n", c->nel);
    int i;
    for (i = 0; i < c->nel; i++)
      { fprintf(wr, "%d\n", pz_double_chain_to_int(c->el[i]/unit)); }
    filefmt_write_footer(wr, "pz_double_chain_t");
    fflush(wr);
  }

pz_double_chain_read_data_t  pz_double_chain_read ( FILE *rd, bool_t header_only )
  {
    auto double pz_double_chain_to_double ( int n );
    double pz_double_chain_to_double ( int n )
      {
        affirm(n > MININT, "out of range");
        affirm(n < MAXINT, "out of range");
        return ((double)n);
      }

    pz_double_chain_read_data_t d;
    filefmt_read_header(rd, "pz_double_chain_t", FileVersion);
    d.cmt = filefmt_read_comment(rd, '|');
    d.stride = nget_double(rd, "stride");  fget_skip_formatting_chars(rd);
    d.unit = nget_double(rd, "unit");      fget_skip_formatting_chars(rd);
    d.samples = nget_int32(rd, "samples");   fget_eol(rd);

    int n = d.samples;
    double unit = d.unit;
    if (header_only)
      { d.c = double_vec_new(n); }
    else
      { d.c = double_vec_new(n);
        int j;
        for (j = 0; j < n; j++)
          { fget_skip_formatting_chars(rd);
            d.c.el[j] = pz_double_chain_to_double(fget_int32(rd)) * unit;
          }
        filefmt_read_footer(rd, "pz_double_chain_t");
      }
    return d;
  }

#define NoData (pz_double_chain_read_data_t){ \
    /*cmt*/ "NOT READ", \
    /*samples*/ 0, \
    /*c*/ (double_vec_t){ 0, NULL }, \
    /*stride*/ 0, \
    /*unit*/ 0 \
  }

pz_double_chain_read_all_data_t pz_double_chain_read_all
  ( char *prefix,
    unsigned band,
    char *extension,
    bool_vec_t *sel,
    bool_t header_only,
    char *dir            /* default was "." */
  )
  {
    pz_double_chain_read_data_t *dd = 
      malloc(sel->nel * sizeof(pz_double_chain_read_data_t));
    double unitMin = +INFINITY;
    double unitMax = -INFINITY;
    char *cmt = ""; 
    
    fprintf(stderr, "reading numeric chains:\n");

    int k;
    for (k = 0; k < sel->nel; k++)
      {
        pz_double_chain_read_data_t *ddk = &(dd[k]); 
        if (sel->el[k])
          {
            fprintf(stderr, "  %04d: ", k);
            char *full_name = jsprintf("%s/%04d/%s%03d%s", dir, k, prefix, band, extension);
            FILE *rd = open_read(full_name, TRUE);
            if (rd != NULL)
              { *ddk = pz_double_chain_read(rd, header_only);
                fprintf(stderr, "read from %s\n", full_name);
                affirm(ddk->unit > 0, "bad unit");
                if (unitMin > unitMax) { cmt = ddk->cmt; }
                unitMin = fmin(unitMin, ddk->unit);
                unitMax = fmax(unitMax, ddk->unit);
                fclose(rd);
              }
            else  
              { *ddk = NoData;
                fprintf(stderr, "file %s not found.\n", full_name);
              }
          }
        else
          { *ddk = NoData; }
      }
    char *buf = NULL; 
    char *buf = jsprintf(
      "Chains read from %s/NNNN/%sBBB%s\nLast entry:\n%s", 
      dir, prefix, extension, cmt
    );
    return (pz_double_chain_read_all_data_t) {
      /* chData */  dd,
      /* unitMin */ unitMin,
      /* unitMax */ unitMax,
      /* cmt */     buf
    };
  }

double  pz_double_chain_delta
  ( pz_double_chain_t *c, 
    double stride, 
    int ini,
    int fin 
  )
  { int n = c->nel;
    double ri = c->el[ini % n] + stride * ((double)ini / n);
    double rf = c->el[fin % n] + stride * ((double)fin / n);
    return rf - ri;
  }

double pz_double_chain_value_from_arg
  ( pz_double_chain_t *c, 
    double stride, 
    double x 
  )
  { int n = c->nel;
    int ka = (int)floor(x);
    int kb = ka + 1;
    double dx = x - ((double)ka);
    double ya = c->el[ka % n] + stride * ((double)ka / n);
    double yb = c->el[kb % n] + stride * ((double)kb / n);
    return  ya + dx * (yb - ya);
  }

void pz_double_chain_values_from_args
  ( pz_double_chain_t *c,
    double stride,
    double_vec_t *x,
    double_vec_t *y 
  )
  { int n = c->nel;
    int m = x->nel;
    affirm(y->nel == m, "wrong size");
    int i;
    for (i = 0; i < m; i++)
      { int ka = floor(x->el[i]);
        int kb = ka + 1;
        double dx = x->el[i] - ((double)ka);
        double ya = c->el[ka % n] + stride * ((double)ka / n);
        double yb = c->el[kb % n] + stride * ((double)kb / n);
        y->el[i] = ya + dx * (yb - ya);
      }
  }

double pz_double_chain_arg_from_value
  ( pz_double_chain_t *c, 
    double stride, 
    double y 
  )
  {
    int n = c->nel;
    int i, j;
    /* Reduce {y} to the range {[c[0] _ c[0] + stride)} */
    affirm(stride > 0.0, "bad stride");
    int q = (int)floor((y - c->el[0])/stride);
    y = y - ((double)q) * stride;
    while (y >= c->el[0] + stride){ q++; y = y - stride; }
    /* Binary search for the interval {[c[i] _ c{i+1])} that contains {y}: */
    affirm((c->el[0] <= y) && (y < c->el[0] + stride), "bug");
    i = 0; j = n;
    while (i < j-1)
      { int m = (i+j) / 2;
        affirm(m < n, "bug");
        if (y < c->el[m]) { j = m; } else { i = m; }
      }
    /* Interpolation: */
    double fi = ((double)i) + ((double)q*n);
    double ci = c->el[i];
    int k = (i+1)/n;
    double cj = c->el[(i+1) % n] + stride * ((double)k);
    affirm((y >= ci) && (y <= cj), "bug");
    if (ci == cj)
      { return fi + 0.5; }
    else
      { return fi + (y-ci)/(cj-ci); }
  }

void pz_double_chain_args_from_values
  ( pz_double_chain_t *c,
    double stride,
    double_vec_t *y,
    double_vec_t *x 
  )
  {
    unsigned ia, ib;
    int k;
    double yy;

    int nc = c->nel;
    int ny = y->nel;
    double cLo = c->el[0];
    double cHi = cLo + stride;
    affirm(stride > 0.0, "bad stride");
    affirm(x->nel == ny, "wrong size");
    ia = 0; ib = 1;
    int j;
    for (j = 0; j < ny; j++)
      { /* Invariant: {ia IN [0..np-1]}, {ib == ia+1} */
        /* Reduce {y[j]} to the range {[c[0] _ c[0] + stride)} */
        yy = y->el[j];
        pz_reduce_to_period(&yy, cLo, cHi, &k);
        /* Sequential search for the interval {[c[ia] _ c[ib])} that contains {yy}: */
        affirm((cLo <= yy) && (yy < cHi), "bug");
        while ((ib < nc) && (yy > c->el[ib])) { ia = ib; ib++; }
        while (yy < c->el[ia]){ ib = ia; ia--; }
        /* Interpolation: */
        double fi = ((double)ia) + ((double)k*nc);
        double ci = c->el[ia];
        double cj = c->el[ib % nc] + stride * ((double)ib / nc);
        /* do */
        affirm((yy >= ci) && (yy <= cj), "bug");
        if (ci == cj)
          { x->el[j] = fi + 0.5; }
        else
          { x->el[j] = fi + (yy-ci)/(cj-ci); }
      }
  }

double  pz_double_chain_value_from_value
  ( pz_double_chain_t *a, 
    double aStride,
    pz_double_chain_t *b, 
    double bStride,
    double ya 
  )
  {  double x = pz_double_chain_arg_from_value(a, aStride, ya);
     return pz_double_chain_value_from_arg(b, bStride, x);
  }

double  pz_double_chain_arg_from_arg
  ( pz_double_chain_t *a, 
    double aStride,
    pz_double_chain_t *b, 
    double bStride,
    double xa 
  )
  {
    double y = pz_double_chain_value_from_arg(a, aStride, xa);
    return pz_double_chain_arg_from_value(b, bStride, y);
  }

pz_segment_t  pz_double_chain_map_segment
  ( pz_segment_t *sOld,
    pz_double_chain_t *tOld,
    pz_double_chain_t *tNew,
    double stride 
  )
  {
    int cvx = sOld->cvx;

    int iniOld = sOld->ini;
    int finOld = iniOld + sOld->ns - 1;
    /* DON'T take MOD! */
    double labelIni = pz_double_chain_value_from_arg(tOld, stride, ((double)iniOld));
    double labelFin = pz_double_chain_value_from_arg(tOld, stride, ((double)finOld));

    int totNew = tNew->nel;

    int iniNew = pz_round(pz_double_chain_arg_from_value(tNew, stride, labelIni));
    int finNew = pz_round(pz_double_chain_arg_from_value(tNew, stride, labelFin));
    
    affirm(iniNew <= finNew, "bug");
    return (pz_segment_t) {
        /*cvx */ cvx,
        /*tot */ totNew,
        /*ini */ iniNew % totNew,
        /*ns */  fmin(finNew - iniNew + 1, totNew + 1),
        /*rev */ sOld->rev
      };
  }

/*
  Copyright © 2001 Universidade Estadual de Campinas (UNICAMP).
  Authors: Helena C. G. Leitão and Jorge Stolfi.
  
  This file can be freely distributed, used, and modified, provided
  that this copyright and authorship notice is preserved, and that any
  modified versions are clearly marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the authors nor their employers chall be held
  responsible for any losses or damages that may result from its use.
*/