#define PERL_NO_GET_CONTEXT
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"

#include "ppport.h"

#include "histogram.h"
#include "histo_perl_interf.h"

#include "mt.h"
#include "const-c.inc"

/* More HS_* macros to be found in histogram.h. Those here
 * are more chummy with perl than those in histogram.h, which only
 * currently use the memory allocation macros of perl. */

/* Ideally, HS_ASSERT_BIN_RANGE would be part of the histogram.h
 * API, but given that we know the unsigned data type here AND have
 * access to perl's croak conveniently, that seems premature cleanup. */
#define HS_ASSERT_BIN_RANGE(self, i) STMT_START {                                     \
  if (/* i < 0 || */ i >= self->nbins) {                                              \
    croak("Bin %u outside histogram range (highest bin index is %u", i, self->nbins); \
  } } STMT_END

#define HS_CLONE_GET_CLASS(classname, src, where) STMT_START {                        \
  if (!sv_isobject(src))                                                              \
    croak("Cannot call " #where "() on non-object");                                  \
  classname = sv_reftype(SvRV(src), TRUE);                                            \
  if ( !sv_isobject(src) || (SvTYPE(SvRV(src)) != SVt_PVMG) )                         \
    croak( "%s::" #where "() -- self is not a blessed SV reference", classname);      \
  } STMT_END


/* The following couple of lines are for the RNG, taken from Math::Random::MT */
typedef struct mt* Math__SimpleHisto__XS__RNG;

void*
U32ArrayPtr (pTHX_ int n)
{
  SV * sv = sv_2mortal( NEWSV( 0, n*sizeof(U32) ) );
  return SvPVX(sv);
}


MODULE = Math::SimpleHisto::XS    PACKAGE = Math::SimpleHisto::XS

REQUIRE: 2.2201

INCLUDE: const-xs.inc

INCLUDE: XS/rdgen.xs

INCLUDE: XS/construction.xs

INCLUDE: XS/bins.xs

INCLUDE: XS/aggregate.xs

void
DESTROY(self)
    simple_histo_1d* self
  CODE:
    HS_DEALLOCATE(self);

void
multiply_constant(self, factor = 1.)
    simple_histo_1d* self
    double factor
  CODE:
    if (factor < 0.)
      croak("Cannot multiply histogram with negative value %f", factor);
    histo_multiply_constant(self, factor);

void
add_histogram(self, operand)
    simple_histo_1d* self
    simple_histo_1d* operand
  ALIAS:
    subtract_histogram = 1
    multiply_histogram = 2
    divide_histogram = 3
  PREINIT:
    bool ok;
  CODE:
    if (ix == 0)
      ok = histo_add_histogram(self, operand);
    else if (ix == 1)
      ok = histo_subtract_histogram(self, operand);
    else if (ix == 2)
      ok = histo_multiply_histogram(self, operand);
    else
      ok = histo_divide_histogram(self, operand);

    if (!ok) {
      char *reason;
      if (ix == 0)
        reason = "add";
      else if (ix == 1)
        reason = "subtract";
      else if (ix == 2)
        reason = "multiply";
      else
        reason = "divide";
      croak("Failed to %s incompatible histogram. Binning not the same?", reason);
    }

void
normalize(self, normalization = 1.)
    simple_histo_1d* self
    double normalization
  CODE:
    if (normalization <= 0.)
      croak("Cannot normalize to %f", normalization);
    if (self->total == 0.)
      croak("Cannot normalize histogram without data");
    histo_multiply_constant(self, normalization / self->total);

void
fill(self, ...)
    simple_histo_1d* self
  CODE:
    if (items == 2) {
      SV* const x_tmp = ST(1);
      SvGETMAGIC(x_tmp);
      if (SvROK(x_tmp) && SvTYPE(SvRV(x_tmp)) == SVt_PVAV) {
        int i, n;
        SV** sv;
        double* x;
        AV* av = (AV*)SvRV(x_tmp);
        n = av_len(av);
        Newx(x, n+1, double);
        for (i = 0; i <= n; ++i) {
          sv = av_fetch(av, i, 0);
          if (UNLIKELY( sv == NULL )) {
            Safefree(x);
            croak("Shouldn't happen");
          }
          x[i] = SvNV(*sv);
        }
        histo_fill(self, n+1, x, NULL);
        Safefree(x);
      }
      else {
        double x = SvNV(ST(1));
        histo_fill(self, 1, &x, NULL);
      }
    }
    else if (items == 3) {
      SV* const x_tmp = ST(1);
      SV* const w_tmp = ST(2);
      SvGETMAGIC(x_tmp);
      SvGETMAGIC(w_tmp);
      if (SvROK(x_tmp) && SvTYPE(SvRV(x_tmp)) == SVt_PVAV) {
        int i, n;
        SV** sv;
        double *x, *w;
        AV *xav, *wav;
        if (UNLIKELY( !SvROK(w_tmp) || SvTYPE(SvRV(x_tmp)) != SVt_PVAV )) {
          croak("Need array of weights if using array of x values");
        }
        xav = (AV*)SvRV(x_tmp);
        wav = (AV*)SvRV(w_tmp);
        n = av_len(xav);
        if (UNLIKELY( av_len(wav) != n )) {
          croak("x and w array lengths differ");
        }

        Newx(x, n+1, double);
        Newx(w, n+1, double);
        for (i = 0; i <= n; ++i) {
          sv = av_fetch(xav, i, 0);
          if (UNLIKELY( sv == NULL )) {
            Safefree(x);
            Safefree(w);
            croak("Shouldn't happen");
          }
          x[i] = SvNV(*sv);

          sv = av_fetch(wav, i, 0);
          if (UNLIKELY( sv == NULL )) {
            Safefree(x);
            Safefree(w);
            croak("Shouldn't happen");
          }
          w[i] = SvNV(*sv);
        }
        histo_fill(self, n+1, x, w);
        Safefree(x);
        Safefree(w);
      }
      else {
        double x = SvNV(ST(1));
        double w = SvNV(ST(2));
        histo_fill(self, 1, &x, &w);
      }
    }
    else {
      croak("Invalid number of arguments to fill(self, ...)");
    }

void
fill_by_bin(self, ...)
    simple_histo_1d* self
  CODE:
    HS_INVALIDATE_CUMULATIVE(self);
    if (items == 2) {
      SV* const x_tmp = ST(1);
      SvGETMAGIC(x_tmp);
      if (SvROK(x_tmp) && SvTYPE(SvRV(x_tmp)) == SVt_PVAV) {
        int i, n;
        SV** sv;
        int* ibin;
        AV* av = (AV*)SvRV(x_tmp);
        n = av_len(av);
        Newx(ibin, n+1, int);
        for (i = 0; i <= n; ++i) {
          sv = av_fetch(av, i, 0);
          if (sv == NULL) {
            Safefree(ibin);
            croak("Shouldn't happen");
          }
          ibin[i] = SvIV(*sv);
        }
        histo_fill_by_bin(self, n+1, ibin, NULL);
        Safefree(ibin);
      }
      else {
        const int ibin = SvUV(ST(1));
        histo_fill_by_bin(self, 1, &ibin, NULL);
      }
    }
    else if (items == 3) {
      SV* const x_tmp = ST(1);
      SV* const w_tmp = ST(2);
      SvGETMAGIC(x_tmp);
      SvGETMAGIC(w_tmp);
      if (SvROK(x_tmp) && SvTYPE(SvRV(x_tmp)) == SVt_PVAV) {
        int i, n;
        SV** sv;
        int *ibin;
        double *w;
        AV *xav, *wav;
        if (!SvROK(w_tmp) || SvTYPE(SvRV(x_tmp)) != SVt_PVAV) {
          croak("Need array of weights if using array of bin numbers");
        }
        xav = (AV*)SvRV(x_tmp);
        wav = (AV*)SvRV(w_tmp);
        n = av_len(xav);
        if (av_len(wav) != n) {
          croak("ibin and w array lengths differ");
        }

        Newx(ibin, n+1, int);
        Newx(w, n+1, double);
        for (i = 0; i <= n; ++i) {
          sv = av_fetch(xav, i, 0);
          if (sv == NULL) {
            Safefree(ibin);
            Safefree(w);
            croak("Shouldn't happen");
          }
          ibin[i] = SvIV(*sv);

          sv = av_fetch(wav, i, 0);
          if (sv == NULL) {
            Safefree(ibin);
            Safefree(w);
            croak("Shouldn't happen");
          }
          w[i] = SvNV(*sv);
        }
        histo_fill_by_bin(self, n+1, ibin, w);
        Safefree(ibin);
        Safefree(w);
      }
      else {
        int ibin = SvIV(ST(1));
        double w = SvNV(ST(2));
        histo_fill_by_bin(self, 1, &ibin, &w);
      }
    }
    else {
      croak("Invalid number of arguments to fill(self, ...)");
    }


double
min(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->min;
  OUTPUT: RETVAL

double
max(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->max;
  OUTPUT: RETVAL

double
width(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->width;
  OUTPUT: RETVAL

double
overflow(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->overflow;
  OUTPUT: RETVAL

double
underflow(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->underflow;
  OUTPUT: RETVAL

double
total(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->total;
  OUTPUT: RETVAL

unsigned int
nbins(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->nbins;
  OUTPUT: RETVAL

unsigned int
highest_bin(self)
    simple_histo_1d* self
  CODE:
    /* I know. Trivial, but convienient! */
    RETVAL = self->nbins-1;
  OUTPUT: RETVAL

double
binsize(self, ibin = 0)
    simple_histo_1d* self
    unsigned int ibin
  CODE:
    HS_ASSERT_BIN_RANGE(self, ibin);
    if (self->bins == NULL)
      RETVAL = self->binsize;
    else
      RETVAL = self->bins[ibin+1] - self->bins[ibin];
  OUTPUT: RETVAL

unsigned int
nfills(self)
    simple_histo_1d* self
  CODE:
    RETVAL = self->nfills;
  OUTPUT: RETVAL


void
all_bin_contents(self)
    simple_histo_1d* self
  PREINIT:
    SV* rv;
  PPCODE:
    rv = histo_data_av(aTHX_ self);
    XPUSHs(sv_2mortal(rv));

void
set_all_bin_contents(self, new_data)
    simple_histo_1d* self
    AV* new_data
  PREINIT:
    unsigned int n, i;
    double* data;
    SV** elem;
  CODE:
    /* While this would be nicer in the histogram API, it will be much faster
     * to access the AV* on the fly instead of doing blanket conversion to remove
     * dependence on perl data structures, so this stays here for the time being. */
    HS_INVALIDATE_CUMULATIVE(self);
    n = self->nbins;
    if ((unsigned int)(av_len(new_data)+1) != n) {
      croak("Length of new data is %u, size of histogram is %u. That doesn't work.", (unsigned int)(av_len(new_data)+1), n);
    }
    data = self->data;
    for (i = 0; i < n; ++i) {
      elem = av_fetch(new_data, i, 0);
      if (elem == NULL) {
        croak("Shouldn't happen");
      }
      self->total -= data[i];
      data[i] = SvNV(*elem);
      self->total += data[i];
    }

double
bin_content(self, ibin)
    simple_histo_1d* self
    unsigned int ibin
  CODE:
    HS_ASSERT_BIN_RANGE(self, ibin);
    RETVAL = self->data[ibin];
  OUTPUT: RETVAL

void
set_bin_content(self, ibin, content)
    simple_histo_1d* self
    unsigned int ibin
    double content
  PPCODE:
    /* Would be nicer in the API, but again, this is faster. */
    HS_ASSERT_BIN_RANGE(self, ibin);
    HS_INVALIDATE_CUMULATIVE(self);
    self->total += content - self->data[ibin];
    self->data[ibin] = content;

void
set_underflow(self, content)
    simple_histo_1d* self
    double content
  PPCODE:
    /* This doesn't invalidate the INTERNAL cumulative histo */
    self->underflow = content;

void
set_overflow(self, content)
    simple_histo_1d* self
    double content
  PPCODE:
    /* This doesn't invalidate the INTERNAL cumulative histo */
    self->overflow = content;


void
set_nfills(self, nfills)
    simple_histo_1d* self
    unsigned int nfills
  PPCODE:
    /* This doesn't invalidate the INTERNAL cumulative histo */
    self->nfills = nfills;


#void
#binary_dump(self)
#    simple_histo_1d* self
#  PREINIT:
#    char* out;
#    SV* outSv;
#    double* tmp;
#    unsigned int size;
#  PPCODE:
#    size = sizeof(simple_histo_1d) + sizeof(double)*self->nbins;
#    outSv = newSVpvs("");
#    SvGROW(outSv, size+1);
#    printf("   %u\n", SvLEN(outSv));
#    out = SvPVX(outSv);
#    SvLEN_set(outSv, size);
#    printf("%u\n", SvLEN(outSv));
#    /*Newx(out, size+1, char);*/
#    tmp = self->data;
#    self->data = NULL;
#    Copy(self, out, sizeof(simple_histo_1d), char);
#    Copy(tmp, out+sizeof(simple_histo_1d), sizeof(double)*self->nbins, char);
#    out[size] = '\0';
#    printf("%u\n", SvLEN(outSv));
#    self->data = tmp;
#    XPUSHs(sv_2mortal(outSv));


void
rand(self, ...)
    simple_histo_1d* self
  PREINIT:
    double rndval;
    double retval;
    SV* rngsv;
    Math__SimpleHisto__XS__RNG rng;
    unsigned int ibin;
  PREINIT:
    simple_histo_1d* cum_hist;
  PPCODE:
    if (items > 1) {
      rngsv = ST(1);
    }
    else {
      rngsv = get_sv("Math::SimpleHisto::XS::RNG::Gen", 0);
      if (rngsv == 0) {
        croak("Cannot find default random number generator!");
      }
    }
    if (sv_derived_from(rngsv, "Math::SimpleHisto::XS::RNG")) {
      IV tmp = SvIV((SV*)SvRV(rngsv));
      rng = INT2PTR(Math__SimpleHisto__XS__RNG, tmp);
    }
    else
      Perl_croak(aTHX_ "%s: %s is not of type %s",
                  "Math::SimpleHisto::XS::rand",
                  "rng", "Math::SimpleHisto::XS::RNG");
    rndval = mt_genrand(rng);

    /* Get the properly normalized internal cumulative */
    HS_ASSERT_CUMULATIVE(self);
    cum_hist = self->cumulative_hist;

    /* This all operates on the cumulative histogram */
    ibin = rndval < cum_hist->data[0] ? 0 : find_bin_nonconstant(rndval, cum_hist->nbins, cum_hist->data);
    if (cum_hist->bins == 0) { /* constant bin size */
      retval = cum_hist->min + cum_hist->binsize * (double)(ibin+1);
      if (rndval > cum_hist->data[ibin]) {
        retval += cum_hist->binsize * (rndval - cum_hist->data[ibin])
                                      / (cum_hist->data[ibin+1] - cum_hist->data[ibin]);
      }
    }
    else { /* variable bin size */
      retval = cum_hist->bins[ibin+1];
      if (rndval > cum_hist->data[ibin]) {
        retval += (cum_hist->bins[ibin+1] - cum_hist->bins[ibin])
                  * (rndval - cum_hist->data[ibin])
                  / (cum_hist->data[ibin+1] - cum_hist->data[ibin]);
      }
    }
    XPUSHs(sv_2mortal(newSVnv(retval)));


void
_get_info(self)
    simple_histo_1d* self
  PREINIT:
    SV* data_ary;
    SV* bins_ary;
  PPCODE:
    /* min, max, nbins, nfills, overflow, underflow, dataref, binsref*/
    EXTEND(SP, 8);
    mPUSHn(self->min);
    mPUSHn(self->max);
    mPUSHu(self->nbins);
    mPUSHu(self->nfills);
    mPUSHn(self->overflow);
    mPUSHn(self->underflow);
    data_ary = histo_data_av(aTHX_ self);
    XPUSHs(sv_2mortal(data_ary));
    if (self->bins == NULL)
      bins_ary = &PL_sv_undef;
    else
      bins_ary = sv_2mortal(histo_bins_av(aTHX_ self));
    XPUSHs(bins_ary);