/* See {dm_sample.h}. */
/* Last edited on 2006-12-25 19:59:50 by stolfi */

#define dm_sample_C_COPYRIGHT \
  "Copyright © 2006  by the State University of Campinas (UNICAMP)"

#include <dm_basic.h>
#include <dm_nucleic.h>
#include <dm_seq.h>

#include <vec.h>
#include <affirm.h>

#include <math.h>
#include <stdlib.h>
#include <ctype.h>
#include <assert.h>

#include <math.h>
#include <ctype.h>
#include <assert.h>

#define dm_sample_VERY_MIN (-128)
#define dm_sample_VERY_MAX (+127)
  /* True range of any {dm_sample_t} variable, including 
    invalid values. */

#define dm_sample_BIAS (-dm_sample_VERY_MIN)
  /* An integer such that {v + dm_sample_BIAS >= 0} for any 
    {dm_sample_t} {v}, valid or not. */

/* INTERNAL PROTOTYPES: */

static double *dm_decode_table = NULL;
  /* Sample decoding table.  The decoded value of signed byte {v}
    is {dm_decode_table[dm_sample_BIAS+v]}. */

void dm_decode_table_setup(void); 
  /* Allocates and initializes the table {dm_decode_table}. */

/* IMPLEMENTATIONS: */

void dm_decode_table_setup(void)
  { 
    auto double cumpr(double x);
      /* Cumulative probability distribution for one-sided
        (0,1)-normal variable. */
      
    double cumpr(double x) { return erf(x/M_SQRT2); }
    
    /* Paranoia: */
    assert(dm_sample_VERY_MIN <= dm_sample_VALID_MIN);
    assert(dm_sample_VERY_MAX >= dm_sample_VALID_MAX);
    assert(dm_sample_VALID_MIN + dm_sample_VALID_MAX == 0);
    
    int N = dm_sample_BIAS + dm_sample_VERY_MAX + 1;
    double *tb = (double *)notnull(malloc(N*sizeof(double)), "no mem");
    int v; /* Must be {int} and not {dm_sample_t} to avoid overflow. */
    /* Fill the table with {NaN}s: */
    for (v = dm_sample_VERY_MIN; v <= dm_sample_VERY_MAX; v++)
      { tb[dm_sample_BIAS + v] = NAN; }
    /* Now set valid entries: */
    tb[dm_sample_BIAS + 0] = 0.0;
    double xprev = 0.0;
    for (v = 1; v <= dm_sample_VALID_MAX; v++)
      { /* Map {v} linearly to a probability between 0 and 1: */
        double fv = ((double)v)/(dm_sample_VALID_MAX + 0.5);
        /* Find {x} such that {erf(x/sqrt(2)) = fv}: */
        double xinf = xprev;
        double xsup = 10.0;
        /* fprintf(stderr, "\n"); */
        /* fprintf(stderr, "looking for fv = %24.16e\n", fv); */
        /* fprintf(stderr, "xinf = %24.16e cumpr(xinf) = %24.16e\n", xinf, cumpr(xinf)); */
        /* fprintf(stderr, "xsup = %24.16e cumpr(xsup) = %24.16e\n", xsup, cumpr(xsup)); */
        double x;
        assert(cumpr(xinf) <= fv);
        assert(cumpr(xsup) >= fv);
        while(TRUE)
          { double xrad = (xsup - xinf)/2;
            x = xinf + xrad;
            x = fabs(x); /* Hack to ensure that {x} is truncated to double instead of extended. */
            /* fprintf(stderr, "%24.16e %24.16e %24.16e\n", xinf, x, xsup); */
            if ((x <= xinf) || (x >= xsup)) { break; }
            double ex = cumpr(x);
            if (fv < ex) 
              { xsup = x; }
            else if (fv > ex)
              { xinf = x; }
            else
              { break; }
          }
        /* fprintf(stderr, "%03d = %24.16e\n", v, x); */
        /* Store in table: */
        tb[dm_sample_BIAS + v] = +x;
        tb[dm_sample_BIAS - v] = -x;
        xprev = x;
      }
    dm_decode_table = tb;
  }       

double dm_sample_decode(dm_sample_t ev)
  { int iev = ev;
    demand((iev >= dm_sample_VALID_MIN) && (iev <= dm_sample_VALID_MAX), "invalid encoded value"); 
    if (dm_decode_table == NULL) { dm_decode_table_setup(); }
    int k = ev + dm_sample_BIAS; /* Should be in {1..255}. */
    return dm_decode_table[k];
  }

dm_sample_t dm_sample_encode(double dv)
  { /* Must check whether this is the best rounding: */
    double fv = erf(fabs(dv)/M_SQRT2);
    int v = (int)round(fv*(dm_sample_VALID_MAX + 0.5));
    if (v > dm_sample_VALID_MAX) { v = dm_sample_VALID_MAX; }
    if (dv < 0) { v = -v; }
    return v;
  }

double dm_sample_diffsq(dm_sample_t xs, dm_sample_t ys)
  { double D = dm_sample_decode(xs) - dm_sample_decode(ys);
    return D*D;
  }

vec_typeimpl(dm_sample_vec_t,dm_sample_vec,dm_sample_t);