src/frequencies.q

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  18    02111-1307, USA. */
  19
  20 /*
  21   TODO:
  22
  23   * Remember that histograms, bar charts need mean, stddev.
  24 */
  25
  26 #include <config.h>
  27 #include <assert.h>
  28 #include <math.h>
  29 #include <stdlib.h>
  30 #include "alloc.h"
  31 #include "avl.h"
  32 #include "bitvector.h"
  33 #include "hash.h"
  34 #include "pool.h"
  35 #include "command.h"
  36 #include "lexer.h"
  37 #include "error.h"
  38 #include "approx.h"
  39 #include "magic.h"
  40 #include "misc.h"
  41 #include "stats.h"
  42 #include "output.h"
  43 #include "som.h"
  44 #include "tab.h"
  45 #include "var.h"
  46 #include "vfm.h"
  47
  48 #undef DEBUGGING
  49 /*#define DEBUGGING 1 */
  50 #include "debug-print.h"
  51
  52 /* (specification)
  53    FREQUENCIES (frq_):
  54      *variables=custom;
  55      format=cond:condense/onepage(*n:onepage_limit,"%s>=0")/!standard,
  56             table:limit(n:limit,"%s>0")/notable/!table,
  57             labels:!labels/nolabels,
  58             sort:!avalue/dvalue/afreq/dfreq,
  59             spaces:!single/double,
  60             paging:newpage/!oldpage;
  61      missing=miss:include/!exclude;
  62      barchart(ba_)=:minimum(d:min),
  63             :maximum(d:max),
  64             scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0");
  65      histogram(hi_)=:minimum(d:min),
  66             :maximum(d:max),
  67             scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
  68             norm:!nonormal/normal,
  69             incr:increment(d:inc,"%s>0");
  70      hbar(hb_)=:minimum(d:min),
  71             :maximum(d:max),
  72             scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
  73             norm:!nonormal/normal,
  74             incr:increment(d:inc,"%s>0");
  75      grouped=custom;
  76      ntiles=custom;
  77      percentiles=custom;
  78      statistics[st_]=1|mean,2|semean,3|median,4|mode,5|stddev,6|variance,
  79             7|kurtosis,8|skewness,9|range,10|minimum,11|maximum,12|sum,
  80             13|default,14|seskewness,15|sekurtosis,all,none.
  81 */
  82 /* (declarations) */
  83 /* (functions) */
  84
  85 /* Description of a statistic. */
  86 struct frq_info
  87   {
  88     int st_indx;                /* Index into a_statistics[]. */
  89     const char *s10;            /* Identifying string. */
  90   };
  91
  92 /* Table of statistics, indexed by dsc_*. */
  93 static struct frq_info st_name[frq_n_stats + 1] =
  94 {
  95   {FRQ_ST_MEAN, N_("Mean")},
  96   {FRQ_ST_SEMEAN, N_("S.E. Mean")},
  97   {FRQ_ST_MEDIAN, N_("Median")},
  98   {FRQ_ST_MODE, N_("Mode")},
  99   {FRQ_ST_STDDEV, N_("Std Dev")},
 100   {FRQ_ST_VARIANCE, N_("Variance")},
 101   {FRQ_ST_KURTOSIS, N_("Kurtosis")},
 102   {FRQ_ST_SEKURTOSIS, N_("S.E. Kurt")},
 103   {FRQ_ST_SKEWNESS, N_("Skewness")},
 104   {FRQ_ST_SESKEWNESS, N_("S.E. Skew")},
 105   {FRQ_ST_RANGE, N_("Range")},
 106   {FRQ_ST_MINIMUM, N_("Minimum")},
 107   {FRQ_ST_MAXIMUM, N_("Maximum")},
 108   {FRQ_ST_SUM, N_("Sum")},
 109   {-1, 0},
 110 };
 111
 112 /* Percentiles to calculate. */
 113 static double *percentiles;
 114 static int n_percentiles;
 115
 116 /* Groups of statistics. */
 117 #define BI          BIT_INDEX
 118 #define frq_default                                                     \
 119         (BI (frq_mean) | BI (frq_stddev) | BI (frq_min) | BI (frq_max))
 120 #define frq_all                                                 \
 121         (BI (frq_sum) | BI(frq_min) | BI(frq_max)               \
 122          | BI(frq_mean) | BI(frq_semean) | BI(frq_stddev)       \
 123          | BI(frq_variance) | BI(frq_kurt) | BI(frq_sekurt)     \
 124          | BI(frq_skew) | BI(frq_seskew) | BI(frq_range)        \
 125          | BI(frq_range) | BI(frq_mode) | BI(frq_median))
 126
 127 /* Statistics; number of statistics. */
 128 static unsigned long stats;
 129 static int n_stats;
 130
 131 /* Types of graphs. */
 132 enum
 133   {
 134     GFT_NONE,                   /* Don't draw graphs. */
 135     GFT_BAR,                    /* Draw bar charts. */
 136     GFT_HIST,                   /* Draw histograms. */
 137     GFT_HBAR                    /* Draw bar charts or histograms at our discretion. */
 138   };
 139
 140 /* Parsed command. */
 141 static struct cmd_frequencies cmd;
 142
 143 /* Summary of the barchart, histogram, and hbar subcommands. */
 144 static int chart;               /* NONE/BAR/HIST/HBAR. */
 145 static double min, max;         /* Minimum, maximum on y axis. */
 146 static int format;              /* FREQ/PERCENT: Scaling of y axis. */
 147 static double scale, incr;      /* FIXME */
 148 static int normal;              /* FIXME */
 149
 150 /* Variables for which to calculate statistics. */
 151 static int n_variables;
 152 static struct variable **v_variables;
 153
 154 /* Arenas used to store semi-permanent storage. */
 155 static struct pool *int_pool;   /* Integer mode. */
 156 static struct pool *gen_pool;   /* General mode. */
 157
 158 /* Easier access to a_statistics. */
 159 #define stat cmd.a_statistics
 160
 161 static void determine_charts (void);
 162
 163 static void precalc (void);
 164 static int calc_weighting (struct ccase *);
 165 static int calc_no_weight (struct ccase *);
 166 static void postcalc (void);
 167
 168 static void postprocess_freq_tab (struct variable *);
 169 static void dump_full (struct variable *);
 170 static void dump_condensed (struct variable *);
 171 static void dump_statistics (struct variable *, int show_varname);
 172 static void cleanup_freq_tab (struct variable *);
 173
 174 static int compare_value_numeric_a (const void *, const void *, void *);
 175 static int compare_value_alpha_a (const void *, const void *, void *);
 176 static int compare_value_numeric_d (const void *, const void *, void *);
 177 static int compare_value_alpha_d (const void *, const void *, void *);
 178 static int compare_freq_numeric_a (const void *, const void *, void *);
 179 static int compare_freq_alpha_a (const void *, const void *, void *);
 180 static int compare_freq_numeric_d (const void *, const void *, void *);
 181 static int compare_freq_alpha_d (const void *, const void *, void *);
 182 \f
 183 /* Parser and outline. */
 184
 185 static int internal_cmd_frequencies (void);
 186
 187 int
 188 cmd_frequencies (void)
 189 {
 190   int result;
 191
 192   int_pool = pool_create ();
 193   result = internal_cmd_frequencies ();
 194   pool_destroy (int_pool);
 195   int_pool=0;
 196   pool_destroy (gen_pool);
 197   gen_pool=0;
 198   free (v_variables);
 199   v_variables=0;
 200   return result;
 201 }
 202
 203 static int
 204 internal_cmd_frequencies (void)
 205 {
 206   int (*calc) (struct ccase *);
 207   int i;
 208
 209   n_percentiles = 0;
 210   percentiles = NULL;
 211
 212   n_variables = 0;
 213   v_variables = NULL;
 214
 215   for (i = 0; i < default_dict.nvar; i++)
 216     default_dict.var[i]->foo = 0;
 217
 218   lex_match_id ("FREQUENCIES");
 219   if (!parse_frequencies (&cmd))
 220     return CMD_FAILURE;
 221
 222   if (cmd.onepage_limit == NOT_LONG)
 223     cmd.onepage_limit = 50;
 224
 225   /* Figure out statistics to calculate. */
 226   stats = 0;
 227   if (stat[FRQ_ST_DEFAULT] || !cmd.sbc_statistics)
 228     stats |= frq_default;
 229   if (stat[FRQ_ST_ALL])
 230     stats |= frq_all;
 231   if (cmd.sort != FRQ_AVALUE && cmd.sort != FRQ_DVALUE)
 232     stats &= ~frq_median;
 233   for (i = 0; i < frq_n_stats; i++)
 234     if (stat[st_name[i].st_indx])
 235       stats |= BIT_INDEX (i);
 236   if (stats & frq_kurt)
 237     stats |= frq_sekurt;
 238   if (stats & frq_skew)
 239     stats |= frq_seskew;
 240
 241   /* Calculate n_stats. */
 242   n_stats = 0;
 243   for (i = 0; i < frq_n_stats; i++)
 244     if ((stats & BIT_INDEX (i)))
 245       n_stats++;
 246
 247   /* Charting. */
 248   determine_charts ();
 249   if (chart != GFT_NONE || cmd.sbc_ntiles)
 250     cmd.sort = FRQ_AVALUE;
 251
 252   /* Do it! */
 253   update_weighting (&default_dict);
 254   calc = default_dict.weight_index == -1 ? calc_no_weight : calc_weighting;
 255   procedure (precalc, calc, postcalc);
 256
 257   return CMD_SUCCESS;
 258 }
 259
 260 /* Figure out which charts the user requested.  */
 261 static void
 262 determine_charts (void)
 263 {
 264   int count = (!!cmd.sbc_histogram) + (!!cmd.sbc_barchart) + (!!cmd.sbc_hbar);
 265
 266   if (!count)
 267     {
 268       chart = GFT_NONE;
 269       return;
 270     }
 271   else if (count > 1)
 272     {
 273       chart = GFT_HBAR;
 274       msg (SW, _("At most one of BARCHART, HISTOGRAM, or HBAR should be "
 275            "given.  HBAR will be assumed.  Argument values will be "
 276            "given precedence increasing along the order given."));
 277     }
 278   else if (cmd.sbc_histogram)
 279     chart = GFT_HIST;
 280   else if (cmd.sbc_barchart)
 281     chart = GFT_BAR;
 282   else
 283     chart = GFT_HBAR;
 284
 285   min = max = SYSMIS;
 286   format = FRQ_FREQ;
 287   scale = SYSMIS;
 288   incr = SYSMIS;
 289   normal = 0;
 290
 291   if (cmd.sbc_barchart)
 292     {
 293       if (cmd.ba_min != SYSMIS)
 294         min = cmd.ba_min;
 295       if (cmd.ba_max != SYSMIS)
 296         max = cmd.ba_max;
 297       if (cmd.ba_scale == FRQ_FREQ)
 298         {
 299           format = FRQ_FREQ;
 300           scale = cmd.ba_freq;
 301         }
 302       else if (cmd.ba_scale == FRQ_PERCENT)
 303         {
 304           format = FRQ_PERCENT;
 305           scale = cmd.ba_pcnt;
 306         }
 307     }
 308
 309   if (cmd.sbc_histogram)
 310     {
 311       if (cmd.hi_min != SYSMIS)
 312         min = cmd.hi_min;
 313       if (cmd.hi_max != SYSMIS)
 314         max = cmd.hi_max;
 315       if (cmd.hi_scale == FRQ_FREQ)
 316         {
 317           format = FRQ_FREQ;
 318           scale = cmd.hi_freq;
 319         }
 320       else if (cmd.hi_scale == FRQ_PERCENT)
 321         {
 322           format = FRQ_PERCENT;
 323           scale = cmd.ba_pcnt;
 324         }
 325       if (cmd.hi_norm)
 326         normal = 1;
 327       if (cmd.hi_incr == FRQ_INCREMENT)
 328         incr = cmd.hi_inc;
 329     }
 330
 331   if (cmd.sbc_hbar)
 332     {
 333       if (cmd.hb_min != SYSMIS)
 334         min = cmd.hb_min;
 335       if (cmd.hb_max != SYSMIS)
 336         max = cmd.hb_max;
 337       if (cmd.hb_scale == FRQ_FREQ)
 338         {
 339           format = FRQ_FREQ;
 340           scale = cmd.hb_freq;
 341         }
 342       else if (cmd.hb_scale == FRQ_PERCENT)
 343         {
 344           format = FRQ_PERCENT;
 345           scale = cmd.ba_pcnt;
 346         }
 347       if (cmd.hb_norm)
 348         normal = 1;
 349       if (cmd.hb_incr == FRQ_INCREMENT)
 350         incr = cmd.hb_inc;
 351     }
 352
 353   if (min != SYSMIS && max != SYSMIS && min >= max)
 354     {
 355       msg (SE, _("MAX must be greater than or equal to MIN, if both are "
 356            "specified.  However, MIN was specified as %g and MAX as %g.  "
 357            "MIN and MAX will be ignored."), min, max);
 358       min = max = SYSMIS;
 359     }
 360 }
 361
 362 /* Generate each calc_*(). */
 363 #define WEIGHTING 0
 364 #include "frequencies.g"
 365
 366 #define WEIGHTING 1
 367 #include "frequencies.g"
 368
 369 /* Prepares each variable that is the target of FREQUENCIES by setting
 370    up its hash table. */
 371 static void
 372 precalc (void)
 373 {
 374   int i;
 375
 376   pool_destroy (gen_pool);
 377   gen_pool = pool_create ();
 378
 379   for (i = 0; i < n_variables; i++)
 380     {
 381       struct variable *v = v_variables[i];
 382
 383       if (v->p.frq.tab.mode == FRQM_GENERAL)
 384         {
 385           avl_comparison_func compare;
 386           if (v->type == NUMERIC)
 387             compare = compare_value_numeric_a;
 388           else
 389             compare = compare_value_alpha_a;
 390           v->p.frq.tab.tree = avl_create (gen_pool, compare,
 391                                           (void *) v->width);
 392           v->p.frq.tab.n_missing = 0;
 393         }
 394       else
 395         {
 396           int j;
 397
 398           for (j = (v->p.frq.tab.max - v->p.frq.tab.min); j >= 0; j--)
 399             v->p.frq.tab.vector[j] = 0.0;
 400           v->p.frq.tab.out_of_range = 0.0;
 401           v->p.frq.tab.sysmis = 0.0;
 402         }
 403     }
 404 }
 405
 406 /* Finishes up with the variables after frequencies have been
 407    calculated.  Displays statistics, percentiles, ... */
 408 static void
 409 postcalc (void)
 410 {
 411   int i;
 412
 413   for (i = 0; i < n_variables; i++)
 414     {
 415       struct variable *v = v_variables[i];
 416       int n_categories;
 417       int dumped_freq_tab = 1;
 418
 419       postprocess_freq_tab (v);
 420
 421       /* Frequencies tables. */
 422       n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
 423       if (cmd.table == FRQ_TABLE
 424           || (cmd.table == FRQ_LIMIT && n_categories <= cmd.limit))
 425         switch (cmd.cond)
 426           {
 427           case FRQ_CONDENSE:
 428             dump_condensed (v);
 429             break;
 430           case FRQ_STANDARD:
 431             dump_full (v);
 432             break;
 433           case FRQ_ONEPAGE:
 434             if (n_categories > cmd.onepage_limit)
 435               dump_condensed (v);
 436             else
 437               dump_full (v);
 438             break;
 439           default:
 440             assert (0);
 441           }
 442       else
 443         dumped_freq_tab = 0;
 444
 445       /* Statistics. */
 446       if (n_stats)
 447         dump_statistics (v, !dumped_freq_tab);
 448
 449       cleanup_freq_tab (v);
 450     }
 451 }
 452
 453 /* Comparison function called by comparison_helper(). */
 454 static avl_comparison_func comparison_func;
 455
 456 /* Passed to comparison function by comparison_helper(). */
 457 static void *comparison_param;
 458
 459 /* Used by postprocess_freq_tab to re-sort frequency tables. */
 460 static int
 461 comparison_helper (const void *a, const void *b)
 462 {
 463   return comparison_func (&((struct freq *) a)->v,
 464                           &((struct freq *) b)->v, comparison_param);
 465 }
 466
 467 /* Used by postprocess_freq_tab to construct the array members valid,
 468    missing of freq_tab. */
 469 static void
 470 add_freq (void *data, void *param)
 471 {
 472   struct freq *f = data;
 473   struct variable *v = param;
 474
 475   v->p.frq.tab.total_cases += f->c;
 476
 477   if ((v->type == NUMERIC && f->v.f == SYSMIS)
 478       || (cmd.miss == FRQ_EXCLUDE && is_user_missing (&f->v, v)))
 479     {
 480       *v->p.frq.tab.missing++ = *f;
 481       v->p.frq.tab.valid_cases -= f->c;
 482     }
 483   else
 484     *v->p.frq.tab.valid++ = *f;
 485 }
 486
 487 static void
 488 postprocess_freq_tab (struct variable * v)
 489 {
 490   avl_comparison_func compare;
 491
 492   switch (cmd.sort | (v->type << 16))
 493     {
 494       /* Note that q2c generates tags beginning with 1000. */
 495     case FRQ_AVALUE | (NUMERIC << 16):
 496       compare = NULL;
 497       break;
 498     case FRQ_AVALUE | (ALPHA << 16):
 499       compare = NULL;
 500       break;
 501     case FRQ_DVALUE | (NUMERIC << 16):
 502       comparison_func = compare_value_numeric_d;
 503       break;
 504     case FRQ_DVALUE | (ALPHA << 16):
 505       compare = compare_value_alpha_d;
 506       break;
 507     case FRQ_AFREQ | (NUMERIC << 16):
 508       compare = compare_freq_numeric_a;
 509       break;
 510     case FRQ_AFREQ | (ALPHA << 16):
 511       compare = compare_freq_alpha_a;
 512       break;
 513     case FRQ_DFREQ | (NUMERIC << 16):
 514       compare = compare_freq_numeric_d;
 515       break;
 516     case FRQ_DFREQ | (ALPHA << 16):
 517       compare = compare_freq_alpha_d;
 518       break;
 519     default:
 520       assert (0);
 521     }
 522   comparison_func = compare;
 523
 524   if (v->p.frq.tab.mode == FRQM_GENERAL)
 525     {
 526       int total;
 527       struct freq_tab *ft = &v->p.frq.tab;
 528
 529       total = avl_count (ft->tree);
 530       ft->n_valid = total - ft->n_missing;
 531       ft->valid = xmalloc (sizeof (struct freq) * total);
 532       ft->missing = &ft->valid[ft->n_valid];
 533       ft->valid_cases = ft->total_cases = 0.0;
 534
 535       avl_walk (ft->tree, add_freq, (void *) v);
 536
 537       ft->valid -= ft->n_valid;
 538       ft->missing -= ft->n_missing;
 539       ft->valid_cases += ft->total_cases;
 540
 541       if (compare)
 542         {
 543           qsort (ft->valid, ft->n_valid, sizeof (struct freq), comparison_helper);
 544           qsort (ft->missing, ft->n_missing, sizeof (struct freq), comparison_helper);
 545         }
 546     }
 547   else
 548     assert (0);
 549 }
 550
 551 static void
 552 cleanup_freq_tab (struct variable * v)
 553 {
 554   if (v->p.frq.tab.mode == FRQM_GENERAL)
 555     {
 556       struct freq_tab *ft = &v->p.frq.tab;
 557
 558       free (ft->valid);
 559     }
 560   else
 561     assert (0);
 562 }
 563
 564 /* Parses the VARIABLES subcommand, adding to
 565    {n_variables,v_variables}. */
 566 static int
 567 frq_custom_variables (struct cmd_frequencies *cmd unused)
 568 {
 569   int mode;
 570   int min, max;
 571
 572   int old_n_variables = n_variables;
 573   int i;
 574
 575   lex_match ('=');
 576   if (token != T_ALL && (token != T_ID || !is_varname (tokid)))
 577     return 2;
 578
 579   if (!parse_variables (NULL, &v_variables, &n_variables,
 580                         PV_APPEND | PV_NO_SCRATCH))
 581     return 0;
 582
 583   for (i = old_n_variables; i < n_variables; i++)
 584     v_variables[i]->p.frq.tab.mode = FRQM_GENERAL;
 585
 586   if (!lex_match ('('))
 587     mode = FRQM_GENERAL;
 588   else
 589     {
 590       mode = FRQM_INTEGER;
 591       if (!lex_force_int ())
 592         return 0;
 593       min = lex_integer ();
 594       lex_get ();
 595       if (!lex_force_match (','))
 596         return 0;
 597       if (!lex_force_int ())
 598         return 0;
 599       max = lex_integer ();
 600       lex_get ();
 601       if (!lex_force_match (')'))
 602         return 0;
 603       if (max < min)
 604         {
 605           msg (SE, _("Upper limit of integer mode value range must be "
 606                      "greater than lower limit."));
 607           return 0;
 608         }
 609     }
 610
 611   for (i = old_n_variables; i < n_variables; i++)
 612     {
 613       struct variable *v = v_variables[i];
 614
 615       if (v->foo != 0)
 616         {
 617           msg (SE, _("Variable %s specified multiple times on VARIABLES "
 618                      "subcommand."), v->name);
 619           return 0;
 620         }
 621
 622       v->foo = 1;               /* Used simply as a marker. */
 623
 624       v->p.frq.tab.valid = v->p.frq.tab.missing = NULL;
 625
 626       if (mode == FRQM_INTEGER)
 627         {
 628           if (v->type != NUMERIC)
 629             {
 630               msg (SE, _("Integer mode specified, but %s is not a numeric "
 631                          "variable."), v->name);
 632               return 0;
 633             }
 634
 635           v->p.frq.tab.min = min;
 636           v->p.frq.tab.max = max;
 637           v->p.frq.tab.vector = pool_alloc (int_pool,
 638                                             sizeof (struct freq) * (max - min + 1));
 639         }
 640       else
 641         v->p.frq.tab.vector = NULL;
 642
 643       v->p.frq.n_groups = 0;
 644       v->p.frq.groups = NULL;
 645     }
 646   return 1;
 647 }
 648
 649 /* Parses the GROUPED subcommand, setting the frq.{n_grouped,grouped}
 650    fields of specified variables. */
 651 static int
 652 frq_custom_grouped (struct cmd_frequencies *cmd unused)
 653 {
 654   lex_match ('=');
 655   if ((token == T_ID && is_varname (tokid)) || token == T_ID)
 656     for (;;)
 657       {
 658         int i;
 659
 660         /* Max, current size of list; list itself. */
 661         int nl, ml;
 662         double *dl;
 663
 664         /* Variable list. */
 665         int n;
 666         struct variable **v;
 667
 668         if (!parse_variables (NULL, &v, &n, PV_NO_DUPLICATE | PV_NUMERIC))
 669           return 0;
 670         if (lex_match ('('))
 671           {
 672             nl = ml = 0;
 673             dl = NULL;
 674             while (token == T_NUM)
 675               {
 676                 if (nl >= ml)
 677                   {
 678                     ml += 16;
 679                     dl = pool_realloc (int_pool, dl, ml * sizeof (double));
 680                   }
 681                 dl[nl++] = tokval;
 682                 lex_get ();
 683                 lex_match (',');
 684               }
 685             /* Note that nl might still be 0 and dl might still be
 686                NULL.  That's okay. */
 687             if (!lex_match (')'))
 688               {
 689                 free (v);
 690                 msg (SE, _("`)' expected after GROUPED interval list."));
 691                 return 0;
 692               }
 693           }
 694         else
 695           nl = 0;
 696
 697         for (i = 0; i < n; i++)
 698           {
 699             if (v[i]->foo == 0)
 700               msg (SE, _("Variables %s specified on GROUPED but not on "
 701                    "VARIABLES."), v[i]->name);
 702             if (v[i]->p.frq.groups != NULL)
 703               msg (SE, _("Variables %s specified multiple times on GROUPED "
 704                    "subcommand."), v[i]->name);
 705             else
 706               {
 707                 v[i]->p.frq.n_groups = nl;
 708                 v[i]->p.frq.groups = dl;
 709               }
 710           }
 711         free (v);
 712         if (!lex_match ('/'))
 713           break;
 714         if ((token != T_ID || !is_varname (tokid)) && token != T_ALL)
 715           {
 716             lex_put_back ('/');
 717             break;
 718           }
 719       }
 720
 721   return 1;
 722 }
 723
 724 /* Adds X to the list of percentiles, keeping the list in proper
 725    order. */
 726 static void
 727 add_percentile (double x)
 728 {
 729   int i;
 730
 731   for (i = 0; i < n_percentiles; i++)
 732     if (x <= percentiles[i])
 733       break;
 734   if (i >= n_percentiles || tokval != percentiles[i])
 735     {
 736       percentiles = pool_realloc (int_pool, percentiles,
 737                                   (n_percentiles + 1) * sizeof (double));
 738       if (i < n_percentiles)
 739         memmove (&percentiles[i + 1], &percentiles[i],
 740                  (n_percentiles - i) * sizeof (double));
 741       percentiles[i] = x;
 742       n_percentiles++;
 743     }
 744 }
 745
 746 /* Parses the PERCENTILES subcommand, adding user-specified
 747    percentiles to the list. */
 748 static int
 749 frq_custom_percentiles (struct cmd_frequencies *cmd unused)
 750 {
 751   lex_match ('=');
 752   if (token != T_NUM)
 753     {
 754       msg (SE, _("Percentile list expected after PERCENTILES."));
 755       return 0;
 756     }
 757
 758   do
 759     {
 760       if (tokval <= 0 || tokval >= 100)
 761         {
 762           msg (SE, _("Percentiles must be greater than "
 763                      "0 and less than 100."));
 764           return 0;
 765         }
 766
 767       add_percentile (tokval / 100.0);
 768       lex_get ();
 769       lex_match (',');
 770     }
 771   while (token == T_NUM);
 772   return 1;
 773 }
 774
 775 /* Parses the NTILES subcommand, adding the percentiles that
 776    correspond to the specified evenly-distributed ntiles. */
 777 static int
 778 frq_custom_ntiles (struct cmd_frequencies *cmd unused)
 779 {
 780   int i;
 781
 782   lex_match ('=');
 783   if (!lex_force_int ())
 784     return 0;
 785   for (i = 1; i < lex_integer (); i++)
 786     add_percentile (1.0 / lex_integer () * i);
 787   lex_get ();
 788   return 1;
 789 }
 790 \f
 791 /* Comparison functions. */
 792
 793 /* Ascending numeric compare of values. */
 794 static int
 795 compare_value_numeric_a (const void *a, const void *b, void *foo unused)
 796 {
 797   return approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
 798 }
 799
 800 /* Ascending string compare of values. */
 801 static int
 802 compare_value_alpha_a (const void *a, const void *b, void *len)
 803 {
 804   return memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
 805 }
 806
 807 /* Descending numeric compare of values. */
 808 static int
 809 compare_value_numeric_d (const void *a, const void *b, void *foo unused)
 810 {
 811   return approx_compare (((struct freq *) b)->v.f, ((struct freq *) a)->v.f);
 812 }
 813
 814 /* Descending string compare of values. */
 815 static int
 816 compare_value_alpha_d (const void *a, const void *b, void *len)
 817 {
 818   return memcmp (((struct freq *) b)->v.s, ((struct freq *) a)->v.s, (int) len);
 819 }
 820
 821 /* Ascending numeric compare of frequency;
 822    secondary key on ascending numeric value. */
 823 static int
 824 compare_freq_numeric_a (const void *a, const void *b, void *foo unused)
 825 {
 826   int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
 827   return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
 828 }
 829
 830 /* Ascending numeric compare of frequency;
 831    secondary key on ascending string value. */
 832 static int
 833 compare_freq_alpha_a (const void *a, const void *b, void *len)
 834 {
 835   int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
 836   return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
 837 }
 838
 839 /* Descending numeric compare of frequency;
 840    secondary key on ascending numeric value. */
 841 static int
 842 compare_freq_numeric_d (const void *a, const void *b, void *foo unused)
 843 {
 844   int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
 845   return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
 846 }
 847
 848 /* Descending numeric compare of frequency;
 849    secondary key on ascending string value. */
 850 static int
 851 compare_freq_alpha_d (const void *a, const void *b, void *len)
 852 {
 853   int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
 854   return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
 855 }
 856 \f
 857 /* Frequency table display. */
 858
 859 /* Sets the widths of all the columns and heights of all the rows in
 860    table T for driver D. */
 861 static void
 862 full_dim (struct tab_table *t, struct outp_driver *d)
 863 {
 864   int lab = cmd.labels == FRQ_LABELS;
 865   int i;
 866
 867   if (lab)
 868     t->w[0] = min (tab_natural_width (t, d, 0), d->prop_em_width * 15);
 869   for (i = lab; i < lab + 5; i++)
 870     t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
 871   for (i = 0; i < t->nr; i++)
 872     t->h[i] = d->font_height;
 873 }
 874
 875 /* Displays a full frequency table for variable V. */
 876 static void
 877 dump_full (struct variable * v)
 878 {
 879   int n_categories;
 880   struct freq *f;
 881   struct tab_table *t;
 882   int r;
 883   double cum_percent = 0.0;
 884   double cum_freq = 0.0;
 885
 886   struct init
 887     {
 888       int c, r;
 889       const char *s;
 890     };
 891
 892   struct init *p;
 893
 894   static struct init vec[] =
 895   {
 896     {4, 0, N_("Valid")},
 897     {5, 0, N_("Cum")},
 898     {1, 1, N_("Value")},
 899     {2, 1, N_("Frequency")},
 900     {3, 1, N_("Percent")},
 901     {4, 1, N_("Percent")},
 902     {5, 1, N_("Percent")},
 903     {0, 0, NULL},
 904     {1, 0, NULL},
 905     {2, 0, NULL},
 906     {3, 0, NULL},
 907     {-1, -1, NULL},
 908   };
 909
 910   int lab = cmd.labels == FRQ_LABELS;
 911
 912   n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
 913   t = tab_create (5 + lab, n_categories + 3, 0);
 914   tab_headers (t, 0, 0, 2, 0);
 915   tab_dim (t, full_dim);
 916
 917   if (lab)
 918     tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value Label"));
 919   for (p = vec; p->s; p++)
 920     tab_text (t, p->c - (p->r ? !lab : 0), p->r,
 921                   TAB_CENTER | TAT_TITLE, gettext (p->s));
 922
 923   r = 2;
 924   for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
 925     {
 926       double percent, valid_percent;
 927
 928       cum_freq += f->c;
 929
 930       percent = f->c / v->p.frq.tab.total_cases * 100.0;
 931       valid_percent = f->c / v->p.frq.tab.valid_cases * 100.0;
 932       cum_percent += valid_percent;
 933
 934       if (lab)
 935         {
 936           char *label = get_val_lab (v, f->v, 0);
 937           if (label != NULL)
 938             tab_text (t, 0, r, TAB_LEFT, label);
 939         }
 940
 941       tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
 942       tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
 943       tab_float (t, 2 + lab, r, TAB_NONE, percent, 5, 1);
 944       tab_float (t, 3 + lab, r, TAB_NONE, valid_percent, 5, 1);
 945       tab_float (t, 4 + lab, r, TAB_NONE, cum_percent, 5, 1);
 946       r++;
 947     }
 948   for (; f < &v->p.frq.tab.valid[n_categories]; f++)
 949     {
 950       cum_freq += f->c;
 951
 952       if (lab)
 953         {
 954           char *label = get_val_lab (v, f->v, 0);
 955           if (label != NULL)
 956             tab_text (t, 0, r, TAB_LEFT, label);
 957         }
 958
 959       tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
 960       tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
 961       tab_float (t, 2 + lab, r, TAB_NONE,
 962                      f->c / v->p.frq.tab.total_cases * 100.0, 5, 1);
 963       tab_text (t, 3 + lab, r, TAB_NONE, _("Missing"));
 964       r++;
 965     }
 966
 967   tab_box (t, TAL_1, TAL_1,
 968            cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
 969            0, 0, 4 + lab, r);
 970   tab_hline (t, TAL_2, 0, 4 + lab, 2);
 971   tab_hline (t, TAL_2, 0, 4 + lab, r);
 972   tab_joint_text (t, 0, r, 0 + lab, r, TAB_RIGHT | TAT_TITLE, _("Total"));
 973   tab_vline (t, TAL_0, 1, r, r);
 974   tab_float (t, 1 + lab, r, TAB_NONE, cum_freq, 8, 0);
 975   tab_float (t, 2 + lab, r, TAB_NONE, 100.0, 5, 1);
 976   tab_float (t, 3 + lab, r, TAB_NONE, 100.0, 5, 1);
 977
 978   tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
 979   tab_submit (t);
 980 }
 981
 982 /* Sets the widths of all the columns and heights of all the rows in
 983    table T for driver D. */
 984 static void
 985 condensed_dim (struct tab_table *t, struct outp_driver *d)
 986 {
 987   int cum_w = max (outp_string_width (d, _("Cum")),
 988                    max (outp_string_width (d, _("Cum")),
 989                         outp_string_width (d, "000")));
 990
 991   int i;
 992
 993   for (i = 0; i < 2; i++)
 994     t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
 995   for (i = 2; i < 4; i++)
 996     t->w[i] = cum_w;
 997   for (i = 0; i < t->nr; i++)
 998     t->h[i] = d->font_height;
 999 }
1000
1001 /* Display condensed frequency table for variable V. */
1002 static void
1003 dump_condensed (struct variable * v)
1004 {
1005   int n_categories;
1006   struct freq *f;
1007   struct tab_table *t;
1008   int r;
1009   double cum_percent = 0.0;
1010
1011   n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
1012   t = tab_create (4, n_categories + 2, 0);
1013
1014   tab_headers (t, 0, 0, 2, 0);
1015   tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value"));
1016   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Freq"));
1017   tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1018   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Cum"));
1019   tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1020   tab_dim (t, condensed_dim);
1021
1022   r = 2;
1023   for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1024     {
1025       double percent;
1026
1027       percent = f->c / v->p.frq.tab.total_cases * 100.0;
1028       cum_percent += f->c / v->p.frq.tab.valid_cases * 100.0;
1029
1030       tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1031       tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1032       tab_float (t, 2, r, TAB_NONE, percent, 3, 0);
1033       tab_float (t, 3, r, TAB_NONE, cum_percent, 3, 0);
1034       r++;
1035     }
1036   for (; f < &v->p.frq.tab.valid[n_categories]; f++)
1037     {
1038       tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1039       tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1040       tab_float (t, 2, r, TAB_NONE,
1041                  f->c / v->p.frq.tab.total_cases * 100.0, 3, 0);
1042       r++;
1043     }
1044
1045   tab_box (t, TAL_1, TAL_1,
1046            cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
1047            0, 0, 3, r - 1);
1048   tab_hline (t, TAL_2, 0, 3, 2);
1049   tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
1050   tab_columns (t, SOM_COL_DOWN, 1);
1051   tab_submit (t);
1052 }
1053 \f
1054 /* Statistical display. */
1055
1056 /* Calculates all the pertinent statistics for variable V, putting
1057    them in array D[].  FIXME: This could be made much more optimal. */
1058 static void
1059 calc_stats (struct variable * v, double d[frq_n_stats])
1060 {
1061   double W = v->p.frq.tab.valid_cases;
1062   double X_bar, M2, M3, M4;
1063   struct freq *f;
1064
1065   /* Calculate the mean. */
1066   X_bar = 0.0;
1067   for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1068     X_bar += f->v.f * f->c;
1069   X_bar /= W;
1070
1071   /* Calculate moments about the mean. */
1072   M2 = M3 = M4 = 0.0;
1073   for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1074     {
1075       double dev = f->v.f - X_bar;
1076       double tmp;
1077       tmp = dev * dev;
1078       M2 += f->c * tmp;
1079       tmp *= dev;
1080       M3 += f->c * tmp;
1081       tmp *= dev;
1082       M4 += f->c * tmp;
1083     }
1084
1085   /* Formulas below are taken from _SPSS Statistical Algorithms_. */
1086   d[frq_min] = v->p.frq.tab.valid[0].v.f;
1087   d[frq_max] = v->p.frq.tab.missing[-1].v.f;
1088   d[frq_mode] = 0.0;
1089   d[frq_range] = d[frq_max] - d[frq_min];
1090   d[frq_median] = 0.0;
1091   d[frq_mean] = X_bar;
1092   d[frq_sum] = X_bar * W;
1093   d[frq_variance] = M2 / (W - 1);
1094   d[frq_stddev] = sqrt (d[frq_variance]);
1095   d[frq_semean] = d[frq_stddev] / sqrt (W);
1096   if (W >= 3.0 && d[frq_variance] > 0)
1097     {
1098       double S = d[frq_stddev];
1099       d[frq_skew] = (W * M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
1100       d[frq_seskew] = sqrt (6.0 * W * (W - 1.0)
1101                             / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
1102     }
1103   else
1104     {
1105       d[frq_skew] = d[frq_seskew] = SYSMIS;
1106     }
1107   if (W >= 4.0 && d[frq_variance] > 0)
1108     {
1109       double S2 = d[frq_variance];
1110       double SE_g1 = d[frq_seskew];
1111
1112       d[frq_kurt] = ((W * (W + 1.0) * M4 - 3.0 * M2 * M2 * (W - 1.0))
1113                      / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2));
1114       d[frq_sekurt] = sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1)
1115                             / ((W - 3.0) * (W + 5.0)));
1116     }
1117   else
1118     {
1119       d[frq_kurt] = d[frq_sekurt] = SYSMIS;
1120     }
1121 }
1122
1123 /* Displays a table of all the statistics requested for variable V. */
1124 static void
1125 dump_statistics (struct variable * v, int show_varname)
1126 {
1127   double stat_value[frq_n_stats];
1128   struct tab_table *t;
1129   int i, r;
1130
1131   if (v->type == ALPHA)
1132     return;
1133   if (v->p.frq.tab.n_valid == 0)
1134     {
1135       msg (SW, _("No valid data for variable %s; statistics not displayed."),
1136            v->name);
1137       return;
1138     }
1139   calc_stats (v, stat_value);
1140
1141   t = tab_create (2, n_stats, 0);
1142   tab_dim (t, tab_natural_dimensions);
1143   tab_vline (t, TAL_1 | TAL_SPACING, 1, 0, n_stats - 1);
1144   for (i = r = 0; i < frq_n_stats; i++)
1145     if (stats & BIT_INDEX (i))
1146       {
1147         tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1148                       gettext (st_name[i].s10));
1149         tab_float (t, 1, r, TAB_NONE, stat_value[i], 11, 3);
1150         r++;
1151       }
1152
1153   tab_columns (t, SOM_COL_DOWN, 1);
1154   if (show_varname)
1155     {
1156       if (v->label)
1157         tab_title (t, 1, "%s: %s", v->name, v->label);
1158       else
1159         tab_title (t, 0, v->name);
1160     }
1161   else
1162     tab_flags (t, SOMF_NO_TITLE);
1163
1164   tab_submit (t);
1165 }
1166 \f
1167 #if 0
1168 /* Statistical calculation. */
1169
1170 static int degree[6];
1171 static int maxdegree, minmax;
1172
1173 static void stat_func (struct freq *, VISIT, int);
1174 static void calc_stats (int);
1175 static void display_stats (int);
1176
1177 /* mapping of data[]:
1178  * 0=>8
1179  * 1=>9
1180  * 2=>10
1181  * index 3: number of modes found (detects multiple modes)
1182  * index 4: number of nodes processed, for calculation of median
1183  * 5=>11
1184  *
1185  * mapping of dbl[]:
1186  * index 0-3: sum of X**i
1187  * index 4: minimum
1188  * index 5: maximum
1189  * index 6: mode
1190  * index 7: median
1191  * index 8: number of cases, valid and missing
1192  * index 9: number of valid cases
1193  * index 10: maximum frequency found, for calculation of mode
1194  * index 11: maximum frequency
1195  */
1196 static void
1197 out_stats (int i)
1198 {
1199   int j;
1200
1201   if (cur_var->type == ALPHA)
1202     return;
1203   for (j = 0; j < 8; j++)
1204     cur_var->dbl[j] = 0.;
1205   cur_var->dbl[10] = 0;
1206   cur_var->dbl[4] = DBL_MAX;
1207   cur_var->dbl[5] = -DBL_MAX;
1208   for (j = 2; j < 5; j++)
1209     cur_var->data[j] = 0;
1210   cur_var->p.frq.median_ncases = cur_var->p.frq.t.valid_cases / 2;
1211   avlwalk (cur_var->p.frq.t.f, stat_func, LEFT_TO_RIGHT);
1212   calc_stats (i);
1213   display_stats (i);
1214 }
1215
1216 static void
1217 calc_stats (int i)
1218 {
1219   struct variable *v;
1220   double n;
1221   double *d;
1222
1223   v = v_variables[i];
1224   n = v->p.frq.t.valid_cases;
1225   d = v->dbl;
1226
1227   if (n < 2 || (n < 3 && stat[FRQ_ST_7]))
1228     {
1229       warn (_("only %g case%s for variable %s, statistics not "
1230             "computed"), n, n == 1 ? "" : "s", v->name);
1231       return;
1232     }
1233   if (stat[FRQ_ST_9])
1234     v->res[FRQ_ST_9] = d[5] - d[4];
1235   if (stat[FRQ_ST_10])
1236     v->res[FRQ_ST_10] = d[4];
1237   if (stat[FRQ_ST_11])
1238     v->res[FRQ_ST_11] = d[5];
1239   if (stat[FRQ_ST_12])
1240     v->res[FRQ_ST_12] = d[0];
1241   if (stat[FRQ_ST_1] || stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_6] || stat[FRQ_ST_7])
1242     {
1243       v->res[FRQ_ST_1] = calc_mean (d, n);
1244       v->res[FRQ_ST_6] = calc_variance (d, n);
1245     }
1246   if (stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_7])
1247     v->res[FRQ_ST_5] = calc_stddev (v->res[FRQ_ST_6]);
1248   if (stat[FRQ_ST_2])
1249     v->res[FRQ_ST_2] = calc_semean (v->res[FRQ_ST_5], n);
1250   if (stat[FRQ_ST_7])
1251     {
1252       v->res[FRQ_ST_7] = calc_kurt (d, n, v->res[FRQ_ST_6]);
1253       v->res[FRQ_ST_14] = calc_sekurt (n);
1254     }
1255   if (stat[FRQ_ST_8])
1256     {
1257       v->res[FRQ_ST_8] = calc_skew (d, n, v->res[FRQ_ST_5]);
1258       v->res[FRQ_ST_15] = calc_seskew (n);
1259     }
1260   if (stat[FRQ_ST_MODE])
1261     {
1262       v->res[FRQ_ST_MODE] = v->dbl[6];
1263       if (v->data[3] > 1)
1264         warn (_("The variable %s has %d modes.  The lowest of these "
1265               "is the one given in the table."), v->name, v->data[3]);
1266     }
1267   if (stat[FRQ_ST_MEDIAN])
1268     v->res[FRQ_ST_MEDIAN] = v->dbl[7];
1269 }
1270
1271 static void
1272 stat_func (struct freq * x, VISIT order, int param)
1273 {
1274   double d, f;
1275
1276   if (order != INORDER)
1277     return;
1278   f = d = x->v.f;
1279   cur_var->dbl[0] += (d * x->c);
1280   switch (maxdegree)
1281     {
1282     case 1:
1283       f *= d;
1284       cur_var->dbl[1] += (f * x->c);
1285       break;
1286     case 2:
1287       f *= d;
1288       cur_var->dbl[1] += (f * x->c);
1289       f *= d;
1290       cur_var->dbl[2] += (f * x->c);
1291       break;
1292     case 3:
1293       f *= d;
1294       cur_var->dbl[1] += (f * x->c);
1295       f *= d;
1296       cur_var->dbl[2] += (f * x->c);
1297       f *= d;
1298       cur_var->dbl[3] += (f * x->c);
1299       break;
1300     }
1301   if (minmax)
1302     {
1303       if (d < cur_var->dbl[4])
1304         cur_var->dbl[4] = d;
1305       if (d > cur_var->dbl[5])
1306         cur_var->dbl[5] = d;
1307     }
1308   if (x->c > cur_var->dbl[10])
1309     {
1310       cur_var->data[3] = 1;
1311       cur_var->dbl[10] = x->c;
1312       cur_var->dbl[6] = x->v.f;
1313     }
1314   else if (x->c == cur_var->dbl[10])
1315     cur_var->data[3]++;
1316   if (cur_var->data[4] < cur_var->p.frq.median_ncases
1317       && cur_var->data[4] + x->c >= cur_var->p.frq.median_ncases)
1318     cur_var->dbl[7] = x->v.f;
1319   cur_var->data[4] += x->c;
1320 }
1321 \f
1322 /* Statistical display. */
1323 static int column, ncolumns;
1324
1325 static void outstat (char *, double);
1326
1327 static void
1328 display_stats (int i)
1329 {
1330   statname *sp;
1331   struct variable *v;
1332   int nlines;
1333
1334   v = v_variables[i];
1335   ncolumns = (margin_width + 3) / 26;
1336   if (ncolumns < 1)
1337     ncolumns = 1;
1338   nlines = sc / ncolumns + (sc % ncolumns > 0);
1339   if (nlines == 2 && sc == 4)
1340     ncolumns = 2;
1341   if (nlines == 3 && sc == 9)
1342     ncolumns = 3;
1343   if (nlines == 4 && sc == 12)
1344     ncolumns = 3;
1345   column = 0;
1346   for (sp = st_name; sp->s != -1; sp++)
1347     if (stat[sp->s] == 1)
1348       outstat (gettext (sp->s10), v->res[sp->s]);
1349   if (column)
1350     out_eol ();
1351   blank_line ();
1352 }
1353
1354 static void
1355 outstat (char *label, double value)
1356 {
1357   char buf[128], *cp;
1358   int dw, n;
1359
1360   cp = &buf[0];
1361   if (!column)
1362     out_header ();
1363   else
1364     {
1365       memset (buf, ' ', 3);
1366       cp = &buf[3];
1367     }
1368   dw = 4;
1369   n = nsprintf (cp, "%-10s %12.4f", label, value);
1370   while (n > 23 && dw > 0)
1371     n = nsprintf (cp, "%-10s %12.*f", label, --dw, value);
1372   outs (buf);
1373   column++;
1374   if (column == ncolumns)
1375     {
1376       column = 0;
1377       out_eol ();
1378     }
1379 }
1380 \f
1381 /* Graphs. */
1382
1383 static rect pb, gb;             /* Page border, graph border. */
1384 static int px, py;              /* Page width, height. */
1385 static int ix, iy;              /* Inch width, height. */
1386
1387 static void draw_bar_chart (int);
1388 static void draw_histogram (int);
1389 static int scale_dep_axis (int);
1390
1391 static void
1392 out_graphs (int i)
1393 {
1394   struct variable *v;
1395
1396   v = v_variables[i];
1397   if (avlcount (cur_var->p.frq.t.f) < 2
1398       || (chart == HIST && v_variables[i]->type == ALPHA))
1399     return;
1400   if (driver_id && set_highres == 1)
1401     {
1402       char *text;
1403
1404       graf_page_size (&px, &py, &ix, &iy);
1405       graf_feed_page ();
1406
1407       /* Calculate borders. */
1408       pb.x1 = ix;
1409       pb.y1 = iy;
1410       pb.x2 = px - ix;
1411       pb.y2 = py - iy;
1412       gb.x1 = pb.x1 + ix;
1413       gb.y1 = pb.y1 + iy;
1414       gb.x2 = pb.x2 - ix / 2;
1415       gb.y2 = pb.y2 - iy;
1416
1417       /* Draw borders. */
1418       graf_frame_rect (COMPONENTS (pb));
1419       graf_frame_rect (COMPONENTS (gb));
1420
1421       /* Draw axis labels. */
1422       graf_font_size (iy / 4);  /* 18-point text */
1423       text = format == PERCENT ? _("Percentage") : _("Frequency");
1424       graf_text (pb.x1 + max (ix, iy) / 4 + max (ix, iy) / 16, gb.y2, text,
1425                  SIDEWAYS);
1426       text = v->label ? v->label : v->name;
1427       graf_text (gb.x1, pb.y2 - iy / 4, text, UPRIGHT);
1428
1429       /* Draw axes, chart proper. */
1430       if (chart == BAR ||
1431           (chart == HBAR
1432        && (avlcount (cur_var->p.frq.t.f) || v_variables[i]->type == ALPHA)))
1433         draw_bar_chart (i);
1434       else
1435         draw_histogram (i);
1436
1437       graf_eject_page ();
1438     }
1439   if (set_lowres == 1 || (set_lowres == 2 && (!driver_id || !set_highres)))
1440     {
1441       static warned;
1442
1443       /* Do character-based graphs. */
1444       if (!warned)
1445         {
1446           warn (_("low-res graphs not implemented"));
1447           warned = 1;
1448         }
1449     }
1450 }
1451
1452 #if __GNUC__ && !__CHECKER__
1453 #define BIG_TYPE long long
1454 #else /* !__GNUC__ */
1455 #define BIG_TYPE double
1456 #endif /* !__GNUC__ */
1457
1458 static void
1459 draw_bar_chart (int i)
1460 {
1461   int bar_width, bar_spacing;
1462   int w, max, row;
1463   double val;
1464   struct freq *f;
1465   rect r;
1466   AVLtraverser *t = NULL;
1467
1468   w = (px - ix * 7 / 2) / avlcount (cur_var->p.frq.t.f);
1469   bar_width = w * 2 / 3;
1470   bar_spacing = w - bar_width;
1471
1472 #if !ALLOW_HUGE_BARS
1473   if (bar_width > ix / 2)
1474     bar_width = ix / 2;
1475 #endif /* !ALLOW_HUGE_BARS */
1476
1477   max = scale_dep_axis (cur_var->p.frq.t.max_freq);
1478
1479   row = 0;
1480   r.x1 = gb.x1 + bar_spacing / 2;
1481   r.x2 = r.x1 + bar_width;
1482   r.y2 = gb.y2;
1483   graf_fill_color (255, 0, 0);
1484   for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1485        f = avltrav (cur_var->p.frq.t.f, &t))
1486     {
1487       char buf2[64];
1488       char *buf;
1489
1490       val = f->c;
1491       if (format == PERCENT)
1492         val = val * 100 / cur_var->p.frq.t.valid_cases;
1493       r.y1 = r.y2 - val * (height (gb) - 1) / max;
1494       graf_fill_rect (COMPONENTS (r));
1495       graf_frame_rect (COMPONENTS (r));
1496       buf = get_val_lab (cur_var, f->v, 0);
1497       if (!buf)
1498         if (cur_var->type == ALPHA)
1499           buf = f->v.s;
1500         else
1501           {
1502             sprintf (buf2, "%g", f->v.f);
1503             buf = buf2;
1504           }
1505       graf_text (r.x1 + bar_width / 2,
1506                  gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1507       row ^= 1;
1508       r.x1 += bar_width + bar_spacing;
1509       r.x2 += bar_width + bar_spacing;
1510     }
1511   graf_fill_color (0, 0, 0);
1512 }
1513
1514 #define round_down(X, V)                        \
1515         (floor ((X) / (V)) * (V))
1516 #define round_up(X, V)                          \
1517         (ceil ((X) / (V)) * (V))
1518
1519 static void
1520 draw_histogram (int i)
1521 {
1522   double lower, upper, interval;
1523   int bars[MAX_HIST_BARS + 1], top, j;
1524   int err, addend, rem, nbars, row, max_freq;
1525   char buf[25];
1526   rect r;
1527   struct freq *f;
1528   AVLtraverser *t = NULL;
1529
1530   lower = min == SYSMIS ? cur_var->dbl[4] : min;
1531   upper = max == SYSMIS ? cur_var->dbl[5] : max;
1532   if (upper - lower >= 10)
1533     {
1534       double l, u;
1535
1536       u = round_up (upper, 5);
1537       l = round_down (lower, 5);
1538       nbars = (u - l) / 5;
1539       if (nbars * 2 + 1 <= MAX_HIST_BARS)
1540         {
1541           nbars *= 2;
1542           u = round_up (upper, 2.5);
1543           l = round_down (lower, 2.5);
1544           if (l + 1.25 <= lower && u - 1.25 >= upper)
1545             nbars--, lower = l + 1.25, upper = u - 1.25;
1546           else if (l + 1.25 <= lower)
1547             lower = l + 1.25, upper = u + 1.25;
1548           else if (u - 1.25 >= upper)
1549             lower = l - 1.25, upper = u - 1.25;
1550           else
1551             nbars++, lower = l - 1.25, upper = u + 1.25;
1552         }
1553       else if (nbars < MAX_HIST_BARS)
1554         {
1555           if (l + 2.5 <= lower && u - 2.5 >= upper)
1556             nbars--, lower = l + 2.5, upper = u - 2.5;
1557           else if (l + 2.5 <= lower)
1558             lower = l + 2.5, upper = u + 2.5;
1559           else if (u - 2.5 >= upper)
1560             lower = l - 2.5, upper = u - 2.5;
1561           else
1562             nbars++, lower = l - 2.5, upper = u + 2.5;
1563         }
1564       else
1565         nbars = MAX_HIST_BARS;
1566     }
1567   else
1568     {
1569       nbars = avlcount (cur_var->p.frq.t.f);
1570       if (nbars > MAX_HIST_BARS)
1571         nbars = MAX_HIST_BARS;
1572     }
1573   if (nbars < MIN_HIST_BARS)
1574     nbars = MIN_HIST_BARS;
1575   interval = (upper - lower) / nbars;
1576
1577   memset (bars, 0, sizeof (int[nbars + 1]));
1578   if (lower >= upper)
1579     {
1580       msg (SE, _("Could not make histogram for %s for specified "
1581            "minimum %g and maximum %g; please discard graph."), cur_var->name,
1582            lower, upper);
1583       return;
1584     }
1585   for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1586        f = avltrav (cur_var->p.frq.t.f, &t))
1587     if (f->v.f == upper)
1588       bars[nbars - 1] += f->c;
1589     else if (f->v.f >= lower && f->v.f < upper)
1590       bars[(int) ((f->v.f - lower) / interval)] += f->c;
1591   bars[nbars - 1] += bars[nbars];
1592   for (j = top = 0; j < nbars; j++)
1593     if (bars[j] > top)
1594       top = bars[j];
1595   max_freq = top;
1596   top = scale_dep_axis (top);
1597
1598   err = row = 0;
1599   addend = width (gb) / nbars;
1600   rem = width (gb) % nbars;
1601   r.x1 = gb.x1;
1602   r.x2 = r.x1 + addend;
1603   r.y2 = gb.y2;
1604   err += rem;
1605   graf_fill_color (255, 0, 0);
1606   for (j = 0; j < nbars; j++)
1607     {
1608       int w;
1609
1610       r.y1 = r.y2 - (BIG_TYPE) bars[j] * (height (gb) - 1) / top;
1611       graf_fill_rect (COMPONENTS (r));
1612       graf_frame_rect (COMPONENTS (r));
1613       sprintf (buf, "%g", lower + interval / 2 + interval * j);
1614       graf_text (r.x1 + addend / 2,
1615                  gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1616       row ^= 1;
1617       w = addend;
1618       err += rem;
1619       while (err >= addend)
1620         {
1621           w++;
1622           err -= addend;
1623         }
1624       r.x1 = r.x2;
1625       r.x2 = r.x1 + w;
1626     }
1627   if (normal)
1628     {
1629       double x, y, variance, mean, step, factor;
1630
1631       variance = cur_var->res[FRQ_ST_VARIANCE];
1632       mean = cur_var->res[FRQ_ST_MEAN];
1633       factor = (1. / (sqrt (2. * PI * variance))
1634                 * cur_var->p.frq.t.valid_cases * interval);
1635       graf_polyline_begin ();
1636       for (x = lower, step = (upper - lower) / (POLYLINE_DENSITY);
1637            x <= upper; x += step)
1638         {
1639           y = factor * exp (-square (x - mean) / (2. * variance));
1640           debug_printf (("(%20.10f, %20.10f)\n", x, y));
1641           graf_polyline_point (gb.x1 + (x - lower) / (upper - lower) * width (gb),
1642                                gb.y2 - y * (height (gb) - 1) / top);
1643         }
1644       graf_polyline_end ();
1645     }
1646   graf_fill_color (0, 0, 0);
1647 }
1648
1649 static int
1650 scale_dep_axis (int max)
1651 {
1652   int j, s, x, y, ty, by;
1653   char buf[10];
1654
1655   x = 10, s = 2;
1656   if (scale != SYSMIS && max < scale)
1657     x = scale, s = scale / 5;
1658   else if (format == PERCENT)
1659     {
1660       max = ((BIG_TYPE) 100 * cur_var->p.frq.t.max_freq
1661              / cur_var->p.frq.t.valid_cases + 1);
1662       if (max < 5)
1663         x = 5, s = 1;
1664       else if (max < 10)
1665         x = 10, s = 2;
1666       else if (max < 25)
1667         x = 25, s = 5;
1668       else if (max < 50)
1669         x = 50, s = 10;
1670       else
1671         max = 100, s = 20;
1672     }
1673   else                          /* format==FREQ */
1674     /* Uses a progression of 10, 20, 50, 100, 200, 500, ... */
1675     for (;;)
1676       {
1677         if (x > max)
1678           break;
1679         x *= 2;
1680         s *= 2;
1681         if (x > max)
1682           break;
1683         x = x / 2 * 5;
1684         s = s / 2 * 5;
1685         if (x > max)
1686           break;
1687         x *= 2;
1688         s *= 2;
1689       }
1690   graf_font_size (iy / 9);      /* 8-pt text */
1691   for (j = 0; j <= x; j += s)
1692     {
1693       y = gb.y2 - (BIG_TYPE) j *(height (gb) - 1) / x;
1694       ty = y - iy / 64;
1695       by = y + iy / 64;
1696       if (ty < gb.y1)
1697         ty += iy / 64, by += iy / 64;
1698       else if (by > gb.y2)
1699         ty -= iy / 64, by -= iy / 64;
1700       graf_fill_rect (gb.x1 - ix / 16, ty, gb.x1, by);
1701       sprintf (buf, "%d", j);
1702       graf_text (gb.x1 - ix / 8, (ty + by) / 2, buf, CRJUST);
1703     }
1704   return x;
1705 }
1706 \f
1707 /* Percentiles. */
1708
1709 static void ungrouped_pcnt (int i);
1710 static int grouped_interval_pcnt (int i);
1711 static void out_pcnt (double, double);
1712
1713 static void
1714 out_percentiles (int i)
1715 {
1716   if (cur_var->type == ALPHA || !n_percentiles)
1717     return;
1718
1719   outs_line (_("Percentile    Value     "
1720              "Percentile    Value     "
1721              "Percentile    Value"));
1722   blank_line ();
1723
1724   column = 0;
1725   if (!g_var[i])
1726     ungrouped_pcnt (i);
1727   else if (g_var[i] == 1)
1728     grouped_interval_pcnt (i);
1729 #if 0
1730   else if (g_var[i] == -1)
1731     grouped_pcnt (i);
1732   else
1733     grouped_boundaries_pcnt (i);
1734 #else /* !0 */
1735   else
1736     warn (_("this form of percentiles not supported"));
1737 #endif
1738   if (column)
1739     out_eol ();
1740 }
1741
1742 static void
1743 out_pcnt (double pcnt, double value)
1744 {
1745   if (!column)
1746     out_header ();
1747   else
1748     outs ("     ");
1749   out ("%7.2f%13.3f", pcnt * 100., value);
1750   column++;
1751   if (column == 3)
1752     {
1753       out_eol ();
1754       column = 0;
1755     }
1756 }
1757
1758 static void
1759 ungrouped_pcnt (int i)
1760 {
1761   AVLtraverser *t = NULL;
1762   struct freq *f;
1763   double *p, *e;
1764   int sum;
1765
1766   p = percentiles;
1767   e = &percentiles[n_percentiles];
1768   sum = 0;
1769   for (f = avltrav (cur_var->p.frq.t.f, &t);
1770        f && p < e; f = avltrav (cur_var->p.frq.t.f, &t))
1771     {
1772       sum += f->c;
1773       while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1774         out_pcnt (*p++, f->v.f);
1775     }
1776 }
1777
1778
1779 static int
1780 grouped_interval_pcnt (int i)
1781 {
1782   AVLtraverser * t = NULL;
1783   struct freq * f, *fp;
1784   double *p, *e, w;
1785   int sum, psum;
1786
1787   p = percentiles;
1788   e = &percentiles[n_percentiles];
1789   w = gl_var[i][0];
1790   sum = psum = 0;
1791   for (fp = 0, f = avltrav (cur_var->p.frq.t.f, &t);
1792        f && p < e;
1793        fp = f, f = avltrav (cur_var->p.frq.t.f, &t))
1794     {
1795       if (fp)
1796         if (fabs (f->v.f - fp->v.f) < w)
1797           {
1798             out_eol ();
1799             column = 0;
1800             return msg (SE, _("Difference between %g and %g is "
1801                               "too small for grouping interval %g."), f->v.f,
1802                         fp->v.f, w);
1803           }
1804       psum = sum;
1805       sum += f->c;
1806       while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1807         {
1808           out_pcnt (p[0], (((p[0] * cur_var->p.frq.t.valid_cases) - psum) * w / f->c
1809                            + (f->v.f - w / 2)));
1810           p++;
1811         }
1812     }
1813   return 1;
1814 }
1815 #endif
1816
1817 /*
1818    Local Variables:
1819    mode: c
1820    End:
1821 */