MODULE = Math::SimpleHisto::XS PACKAGE = Math::SimpleHisto::XS
double
integral(self, from, to, type = 0)
simple_histo_1d* self
double from
double to
int type
PREINIT:
double* data;
unsigned int i, n;
double binsize;
bool invert = 0;
CODE:
/* TODO nonconstant bins */
if (from > to) {
binsize = from; /* abuse as temp var */
from = to;
to = binsize;
invert = 1;
}
data = self->data;
binsize = self->binsize;
/* FIXME handle both to/from being off limits on the same side*/
if (to >= self->max)
to = self->max;
if (from < self->min)
from = self->min;
/*for (i = 1; i < self->nbins; ++i)
printf("%u: %f ", i, data[i]);
printf("\n");
*/
switch(type) {
case INTEGRAL_CONSTANT:
if (self->bins == NULL) {
/* first (fractional) bin */
from = (from - self->min) / binsize;
i = (int)from;
from -= (double)i;
/* last (fractional) bin */
to = (to - self->min) / binsize;
n = (int)to;
to -= (double)n;
if (i == n) {
RETVAL = (to-from) * data[i];
}
else {
RETVAL = data[i] * (1.-from)
+ data[n] * to;
++i;
for (; i < n; ++i)
RETVAL += data[i];
}
}
else { /* variable bin size */
/* TODO optimize */
double* bins = self->bins;
unsigned int nbins = self->nbins;
i = find_bin_nonconstant(from, nbins, bins);
binsize = (bins[i+1]-bins[i]);
RETVAL = (bins[i+1]-from)/binsize * data[i]; /* distance from 'from' to upper boundary of bin times data in bin */
n = find_bin_nonconstant(to, nbins, bins);
if (i == n) {
RETVAL -= (bins[i+1]-to)/binsize * data[i];
}
else {
++i;
for (; i < n; ++i) {
RETVAL += data[i];
}
binsize = bins[n+1]-bins[n];
RETVAL += data[n] * (to-bins[n])/binsize;
}
}
break;
default:
croak("Invalid integration type");
};
if (invert)
RETVAL *= -1.;
OUTPUT: RETVAL
double
mean(self)
simple_histo_1d* self
CODE:
RETVAL = histo_mean(aTHX_ self);
OUTPUT: RETVAL
double
standard_deviation(self, ...)
simple_histo_1d* self
CODE:
if (items > 1)
RETVAL = histo_standard_deviation_with_mean(aTHX_ self, SvNV(ST(1)));
else
RETVAL = histo_standard_deviation(aTHX_ self);
OUTPUT: RETVAL
double
median(self)
simple_histo_1d* self
PREINIT:
CODE:
RETVAL = histo_median(aTHX_ self);
OUTPUT: RETVAL
double
median_absolute_deviation(self, ...)
simple_histo_1d* self
PREINIT:
double median;
double *x, *data;
unsigned int i, n;
simple_histo_1d* madhist;
CODE:
if (items == 2)
median = SvNV(ST(1));
else if (items == 1)
median = histo_median(aTHX_ self);
/* FIXME think hard about the optimal nbins here. Also wrt. variable bin size */
madhist = histo_alloc_new_fixed_bins(aTHX_ self->nbins, 0., self->width);
n = self->nbins;
data = self->data;
Newx(x, n, double);
if (self->bins == 0) {
double min = self->min;
double binsize = self->binsize;
for (i = 0; i < n; ++i) {
/* abs(median-bin_center) */
x[i] = fabs( median - (min + binsize * (i + 0.5)) );
}
}
else { /* variable bin size */
double* bins = self->bins;
for (i = 0; i < n; ++i) {
x[i] = fabs( median - 0.5*(bins[i]+bins[i+1]) );
}
}
histo_fill(madhist, n, x, data);
Safefree(x);
RETVAL = histo_median(aTHX_ madhist);
HS_DEALLOCATE(madhist);
OUTPUT: RETVAL