src/output/charts/plot-hist.c

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