You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
- Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
- 02111-1307, USA. */
+ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ 02110-1301, USA. */
-#include <config.h>
-
-#include <stdio.h>
-#include <plot.h>
#include <math.h>
#include <gsl/gsl_histogram.h>
#include <assert.h>
-
-#include "hash.h"
-
-#include "var.h"
+#include "histogram.h"
#include "chart.h"
-/* Number of bins in which to divide data */
-#define BINS 7
-
-/* The approximate no of ticks on the y axis */
-#define YTICKS 10
-
-#ifndef M_PI
-#define M_PI ( 22.0 / 7.0 )
-#endif
-
-
-static double gaussian(double x, double mu, double sigma ) ;
-
-
-static double
-gaussian(double x, double mu, double sigma )
-{
- return (exp( - (( x - mu )* (x - mu) / (2.0 * sigma * sigma) ))
- / ( sigma * sqrt( M_PI * 2.0) )) ;
-}
-
-
-/* Write the legend of the chart */
-void
-histogram_write_legend(struct chart *ch, const struct normal_curve *norm)
-{
- char buf[100];
- pl_savestate_r(ch->lp);
-
- sprintf(buf,"N = %.2f",norm->N);
- pl_move_r(ch->lp, ch->legend_left, ch->data_bottom);
- pl_alabel_r(ch->lp,0,'b',buf);
-
- sprintf(buf,"Mean = %.1f",norm->mean);
- pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
- pl_alabel_r(ch->lp,0,'b',buf);
-
- sprintf(buf,"Std. Dev = %.2f",norm->stddev);
- pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5 * 2);
- pl_alabel_r(ch->lp,0,'b',buf);
-
- pl_restorestate_r(ch->lp);
-}
-
-static void hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar);
-
-
-static void
-hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar)
+gsl_histogram *
+histogram_create(double bins, double x_min, double x_max)
{
- double upper;
- double lower;
- double height;
-
- const size_t bins = gsl_histogram_bins(hist);
- const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
- const double width = (ch->data_right - ch->data_left) / (double) bins ;
-
+ int n;
+ double bin_width ;
+ double bin_width_2 ;
+ double upper_limit, lower_limit;
- assert ( 0 == gsl_histogram_get_range(hist, bar, &lower, &upper));
+ gsl_histogram *hist = gsl_histogram_alloc(bins);
- assert( upper >= lower);
- height = gsl_histogram_get(hist, bar) *
- (ch->data_top - ch->data_bottom) / gsl_histogram_max_val(hist);
-
- pl_savestate_r(ch->lp);
- pl_move_r(ch->lp,ch->data_left, ch->data_bottom);
- pl_fillcolorname_r(ch->lp, ch->fill_colour);
- pl_filltype_r(ch->lp,1);
-
-
- pl_fboxrel_r(ch->lp,
- x_pos, 0,
- x_pos + width, height);
-
- pl_restorestate_r(ch->lp);
-
- {
- char buf[5];
- snprintf(buf,5,"%g",(upper + lower) / 2.0);
- draw_tick(ch, TICK_ABSCISSA,
- x_pos + width / 2.0, buf);
- }
-
-
-}
-
-
-
-
-void
-histogram_plot(const gsl_histogram *hist,
- const char *factorname,
- const struct normal_curve *norm, short show_normal)
-{
- int i;
- int bins;
-
- struct chart ch;
-
- bins = gsl_histogram_bins(hist);
-
- chart_initialise(&ch);
- chart_write_title(&ch, _("HISTOGRAM"));
-
- chart_write_ylabel(&ch, _("Frequency"));
- chart_write_xlabel(&ch, factorname);
-
- chart_write_yscale(&ch, 0, gsl_histogram_max_val(hist), 5);
-
- for ( i = 0 ; i < bins ; ++i )
- hist_draw_bar(&ch, hist, i);
-
- histogram_write_legend(&ch, norm);
-
- if ( show_normal )
- {
- /* Draw the normal curve */
-
- double d ;
- double x_min, x_max, not_used ;
- double abscissa_scale ;
- double ordinate_scale ;
+ bin_width = chart_rounded_tick((x_max - x_min)/ bins);
+ bin_width_2 = bin_width / 2.0;
+
+ n = ceil( x_max / (bin_width_2) ) ;
+ if ( ! (n % 2 ) ) n++;
+ upper_limit = n * bin_width_2;
+ n = floor( x_min / (bin_width_2) ) ;
+ if ( ! (n % 2 ) ) n--;
+ lower_limit = n * bin_width_2;
- gsl_histogram_get_range(hist, 0, &x_min, ¬_used);
- gsl_histogram_get_range(hist, bins - 1, &x_max, ¬_used);
-
- abscissa_scale = (ch.data_right - ch.data_left) / (x_max - x_min);
- ordinate_scale = (ch.data_top - ch.data_bottom) /
- gsl_histogram_max_val(hist) ;
+ gsl_histogram_set_ranges_uniform(hist, lower_limit, upper_limit);
- pl_move_r(ch.lp, ch.data_left, ch.data_bottom);
- for( d = ch.data_left;
- d <= ch.data_right ;
- d += (ch.data_right - ch.data_left) / 100.0)
- {
- const double x = (d - ch.data_left) / abscissa_scale + x_min ;
-
- pl_fcont_r(ch.lp, d,
- ch.data_bottom + norm->N * ordinate_scale *
- gaussian(x, norm->mean, norm->stddev)
- );
- }
- pl_endpath_r(ch.lp);
- }
- chart_finalise(&ch);
+ return hist;
}