src/output/charts/plot-hist.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 2004 Free Software Foundation, Inc.
   3
   4    This program is free software; you can redistribute it and/or
   5    modify it under the terms of the GNU General Public License as
   6    published by the Free Software Foundation; either version 2 of the
   7    License, or (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful, but
  10    WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12    General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  17    02110-1301, USA. */
  18
  19
  20 #include <config.h>
  21
  22 #include <stdio.h>
  23 #include <plot.h>
  24 #include <math.h>
  25 #include <gsl/gsl_histogram.h>
  26 #include <gsl/gsl_randist.h>
  27 #include <assert.h>
  28
  29 #include <output/charts/plot-hist.h>
  30 #include <output/charts/plot-chart.h>
  31
  32 #include <data/variable.h>
  33 #include <libpspp/hash.h>
  34 #include <output/chart.h>
  35
  36 #include "gettext.h"
  37 #define _(msgid) gettext (msgid)
  38
  39 /* Write the legend of the chart */
  40 void
  41 histogram_write_legend(struct chart *ch, const struct normal_curve *norm)
  42 {
  43   char buf[100];
  44   if ( !ch )
  45     return ;
  46
  47   pl_savestate_r(ch->lp);
  48
  49   sprintf(buf,"N = %.2f",norm->N);
  50   pl_move_r(ch->lp, ch->legend_left, ch->data_bottom);
  51   pl_alabel_r(ch->lp,0,'b',buf);
  52
  53   sprintf(buf,"Mean = %.1f",norm->mean);
  54   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5);
  55   pl_alabel_r(ch->lp,0,'b',buf);
  56
  57   sprintf(buf,"Std. Dev = %.2f",norm->stddev);
  58   pl_fmove_r(ch->lp,ch->legend_left,ch->data_bottom + ch->font_size * 1.5 * 2);
  59   pl_alabel_r(ch->lp,0,'b',buf);
  60
  61   pl_restorestate_r(ch->lp);
  62 }
  63
  64 static void hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar);
  65
  66
  67 static void
  68 hist_draw_bar(struct chart *ch, const gsl_histogram *hist, int bar)
  69 {
  70   if ( !ch )
  71     return ;
  72
  73
  74   {
  75     double upper;
  76     double lower;
  77     double height;
  78
  79     const size_t bins = gsl_histogram_bins(hist);
  80     const double x_pos = (ch->data_right - ch->data_left) * bar / (double) bins ;
  81     const double width = (ch->data_right - ch->data_left) / (double) bins ;
  82
  83
  84     assert ( 0 == gsl_histogram_get_range(hist, bar, &lower, &upper));
  85
  86     assert( upper >= lower);
  87
  88     height = gsl_histogram_get(hist, bar) *
  89       (ch->data_top - ch->data_bottom) / gsl_histogram_max_val(hist);
  90
  91     pl_savestate_r(ch->lp);
  92     pl_move_r(ch->lp,ch->data_left, ch->data_bottom);
  93     pl_fillcolorname_r(ch->lp, ch->fill_colour);
  94     pl_filltype_r(ch->lp,1);
  95
  96
  97     pl_fboxrel_r(ch->lp,
  98                  x_pos, 0,
  99                  x_pos + width, height);
 100
 101     pl_restorestate_r(ch->lp);
 102
 103     {
 104       char buf[5];
 105       snprintf(buf,5,"%g",(upper + lower) / 2.0);
 106       draw_tick(ch, TICK_ABSCISSA,
 107                 x_pos + width / 2.0, buf);
 108     }
 109   }
 110 }
 111
 112
 113
 114
 115 void
 116 histogram_plot(const gsl_histogram *hist,
 117                const char *factorname,
 118                const struct normal_curve *norm, short show_normal)
 119 {
 120   int i;
 121   int bins;
 122
 123   struct chart *ch;
 124
 125   ch = chart_create();
 126   chart_write_title(ch, _("HISTOGRAM"));
 127
 128   chart_write_ylabel(ch, _("Frequency"));
 129   chart_write_xlabel(ch, factorname);
 130
 131   if ( ! hist ) /* If this happens, probably all values are SYSMIS */
 132     {
 133       chart_submit(ch);
 134       return ;
 135     }
 136   else
 137     {
 138       bins = gsl_histogram_bins(hist);
 139     }
 140
 141   chart_write_yscale(ch, 0, gsl_histogram_max_val(hist), 5);
 142
 143   for ( i = 0 ; i < bins ; ++i )
 144       hist_draw_bar(ch, hist, i);
 145
 146   histogram_write_legend(ch, norm);
 147
 148   if ( show_normal  )
 149   {
 150     /* Draw the normal curve */
 151
 152     double d ;
 153     double x_min, x_max, not_used ;
 154     double abscissa_scale ;
 155     double ordinate_scale ;
 156     double range ;
 157
 158     gsl_histogram_get_range(hist, 0, &x_min, &not_used);
 159     range = not_used - x_min;
 160     gsl_histogram_get_range(hist, bins - 1, &not_used, &x_max);
 161     assert(range == x_max - not_used);
 162
 163     abscissa_scale = (ch->data_right - ch->data_left) / (x_max - x_min);
 164     ordinate_scale = (ch->data_top - ch->data_bottom) /
 165       gsl_histogram_max_val(hist) ;
 166
 167     pl_move_r(ch->lp, ch->data_left, ch->data_bottom);
 168     for( d = ch->data_left;
 169          d <= ch->data_right ;
 170          d += (ch->data_right - ch->data_left) / 100.0)
 171       {
 172         const double x = (d - ch->data_left) / abscissa_scale + x_min ;
 173         const double y = norm->N * range *
 174           gsl_ran_gaussian_pdf(x - norm->mean, norm->stddev);
 175
 176         pl_fcont_r(ch->lp,  d,  ch->data_bottom  + y * ordinate_scale);
 177
 178       }
 179     pl_endpath_r(ch->lp);
 180
 181   }
 182   chart_submit(ch);
 183 }
 184