X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;ds=sidebyside;f=src%2Fmath%2Fhistogram.c;h=3c324eff2485d583c3cff0c72e49c137c98aecca;hb=3380cd42c3831b98b71709a0c1fcb175e3e4a5d9;hp=499b8f4dacb846ba42036a48bcba81992e9a63bb;hpb=821d3990dbb9e13b58c6a334a0b1cf563daf4740;p=pspp diff --git a/src/math/histogram.c b/src/math/histogram.c index 499b8f4dac..3c324eff24 100644 --- a/src/math/histogram.c +++ b/src/math/histogram.c @@ -17,10 +17,12 @@ #include #include "math/histogram.h" +#include "math/decimal.h" #include #include +#include "data/settings.h" #include "libpspp/message.h" #include "libpspp/assertion.h" #include "libpspp/cast.h" @@ -57,68 +59,85 @@ destroy (struct statistic *s) } -struct histogram * -histogram_create (double bin_width, double min, double max) +static +double get_slack (double limit, double half_bin_width, int *n_half_bins) { - int bins; - struct histogram *h = xmalloc (sizeof *h); - struct statistic *stat = &h->parent; + double ipart, remainder; + + assert (half_bin_width > 0); + + remainder = modf (limit / half_bin_width, &ipart); + + /* In C modf and % behave in an unexpected (to me at any rate) manner + when presented with a negative value + + For example, modf (-7.0 / 3.0) returns -2.0 R -0.3333 + */ + + *n_half_bins = ipart; + + return remainder * half_bin_width; +} + + +/* This functions adjusts the upper and lower range of the histogram to make them fit BIN_WIDTH + MIN and MAX are the lowest and highest data to be plotted in the histogram. + ADJ_MIN and ADJ_MAX are locations of the adjusted values of MIN and MAX (the range will + always be equal or slightly larger). + Returns the number of bins. + + The "testing_assert" expressions in this function should be algebraically correct. + However, due to floating point rounding they could fail, especially when small numbers + are involved. In normal use, therefore, testing_assert does nothing. + */ +static int +adjust_bin_ranges (double bin_width, double min, double max, double *adj_min, double *adj_max) +{ const double half_bin_width = bin_width / 2.0; /* The lower and upper limits of the histogram, in units of half bin widths */ int lower_limit, upper_limit; - /* -1 if the lower end of the range contains more unused space - than the upper end. - +1 otherwise. */ - short sparse_end = 0; + double lower_slack = get_slack (min, half_bin_width, &lower_limit); + double upper_slack = -get_slack (max, half_bin_width, &upper_limit); - double ul, ll; - double lower_remainder = fabs (modf (min / half_bin_width, &ll)); - double upper_remainder = fabs (modf (max / half_bin_width, &ul)); + testing_assert (max > min); - - if (max == min) + /* If min is negative, then lower_slack may be less than zero. + In this case, the lower bound must be extended in the negative direction + so that it is less than OR EQUAL to min. + */ + if (lower_slack < 0) { - msg (MW, _("Not creating histogram because the data contains less than 2 distinct values")); - free (h); - return NULL; + lower_limit--; + lower_slack += half_bin_width; } - - assert (max > min); - - lower_limit = ll; - - /* If the minimum value is not aligned on a half bin width, - then the lower bound must be extended so that the histogram range includes it. */ - if (lower_remainder > 0) - lower_limit--; + testing_assert (lower_limit * half_bin_width <= min); /* However, the upper bound must be extended regardless, because histogram bins - span the range [lower, upper) */ - upper_limit = ul + 1; - - /* So, in the case of the maximum value coinciding with a half bin width, - the upper end will be the sparse end (because is got extended by a complete - half bin width). In other cases, it depends which got the bigger extension. */ - if (upper_remainder == 0) - sparse_end = +1; - else - sparse_end = lower_remainder < upper_remainder ? -1 : +1; + span the range [lower, upper). In other words, the upper bound must be + greater than max. + */ + upper_limit++;; + upper_slack += half_bin_width; + testing_assert (upper_limit * half_bin_width > max); /* The range must be an EVEN number of half bin_widths */ if ( (upper_limit - lower_limit) % 2) { /* Extend the range at the end which gives the least unused space */ - if (sparse_end == +1) - lower_limit--; + if (upper_slack > lower_slack) + { + lower_limit--; + lower_slack += half_bin_width; + } else - upper_limit++; - - /* Now the other end has more space */ - sparse_end *= -1; + { + upper_limit++; + upper_slack += half_bin_width; + } } /* But the range should be aligned to an ODD number of @@ -127,35 +146,112 @@ histogram_create (double bin_width, double min, double max) */ if ( lower_limit % 2 == 0) { - /* Shift the range away from the sparse end, EXCEPT if that is the upper end, - and it was extended to prevent the maximum value from getting lost */ - if (sparse_end == +1 && upper_remainder > 0) + /* If there is not enough slack at either end to perform a shift, + then we must extend the range so that there is. We must extend + by two half bin widths in order to preserve the EVEN condition + established above. Also, we extend on the end with the least + slack, in order to keep things as balanced as possible. */ + if ( upper_slack > lower_slack && upper_slack <= half_bin_width) { + lower_limit -= 2; + lower_slack += 2 * half_bin_width; + } + + if (lower_slack > upper_slack && lower_slack < half_bin_width) + { + upper_limit += 2; + upper_slack += 2 * half_bin_width; + } + + if (upper_slack > lower_slack) + { + testing_assert (upper_slack > half_bin_width); + + /* Adjust the range to the left */ lower_limit --; upper_limit --; + upper_slack -= half_bin_width; + lower_slack += half_bin_width; } else { + testing_assert (lower_slack >= half_bin_width); + + /* Adjust the range to the right */ lower_limit ++; upper_limit ++; + lower_slack -= half_bin_width; + upper_slack += half_bin_width; } } - bins = (upper_limit - lower_limit) / 2.0; - /* Force the number of bins to lie in a sensible range. - FIXME: We should redo the above calculations if this happens! */ - if (bins > 25) - bins = 25; + /* If there are any completely empty bins, then remove them, + since empty bins don't really add much information to the histogram. + */ + if (upper_slack > 2 * half_bin_width) + { + upper_slack -= 2 * half_bin_width; + upper_limit -=2; + } + + if (lower_slack >= 2 * half_bin_width) + { + lower_slack -= 2 * half_bin_width; + lower_limit +=2; + } + + *adj_min = lower_limit * half_bin_width; + *adj_max = upper_limit * half_bin_width; + + testing_assert (*adj_max > max); + testing_assert (*adj_min <= min); + + return (upper_limit - lower_limit) / 2.0; +} + + + +struct histogram * +histogram_create (double bin_width_in, double min, double max) +{ + struct decimal bin_width; + const int MAX_BINS = 25; + struct histogram *h; + struct statistic *stat; + int bins; + double adjusted_min, adjusted_max; + + if (max == min) + { + msg (MW, _("Not creating histogram because the data contains less than 2 distinct values")); + return NULL; + } + + assert (bin_width_in > 0); + chart_rounded_tick (bin_width_in, &bin_width); + bins = adjust_bin_ranges (decimal_to_double (&bin_width), + min, max, &adjusted_min, &adjusted_max); + + /* Force the number of bins to lie in a sensible range. */ + if (bins > MAX_BINS) + { + chart_rounded_tick ((max - min) / (double) (MAX_BINS - 1), &bin_width); + bins = adjust_bin_ranges (decimal_to_double (&bin_width), + min, max, &adjusted_min, &adjusted_max); + } + + /* Can this ever happen? */ if (bins < 1) bins = 1; + h = xmalloc (sizeof *h); + h->gsl_hist = gsl_histogram_alloc (bins); - gsl_histogram_set_ranges_uniform (h->gsl_hist, - lower_limit * half_bin_width, - upper_limit * half_bin_width); + gsl_histogram_set_ranges_uniform (h->gsl_hist, adjusted_min, adjusted_max); + stat = &h->parent; stat->accumulate = acc; stat->destroy = destroy;