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