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