pintos-os.org Git - pspp/blob - 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, X_mode, M2, M3, M4;
1061   struct freq *f;
1062   int most_often;
1063
1064   /* Calculate the mean and  mode */
1065   X_bar = 0.0;
1066   most_often = -1;
1067   X_mode = SYSMIS;
1068   for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1069     {
1070       /* mean */
1071       X_bar += f->v.f * f->c;
1072
1073       /* mode */
1074       if (most_often < f->c )
1075         {
1076           most_often=f->c;
1077           X_mode= f->v.f;
1078         }
1079       else if ( most_often == f->c )
1080         {
1081           /* if there are 2 values , then mode is undefined */
1082           X_mode=SYSMIS;
1083         }
1084     }
1085   X_bar /= W;
1086
1087   /* Calculate moments about the mean. */
1088   M2 = M3 = M4 = 0.0;
1089   for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1090     {
1091       double dev = f->v.f - X_bar;
1092       double tmp;
1093       tmp = dev * dev;
1094       M2 += f->c * tmp;
1095       tmp *= dev;
1096       M3 += f->c * tmp;
1097       tmp *= dev;
1098       M4 += f->c * tmp;
1099     }
1100
1101   /* Formulas below are taken from _SPSS Statistical Algorithms_. */
1102   d[frq_min] = v->p.frq.tab.valid[0].v.f;
1103   d[frq_max] = v->p.frq.tab.missing[-1].v.f;
1104   d[frq_mode] = X_mode;
1105   d[frq_range] = d[frq_max] - d[frq_min];
1106   d[frq_median] = SYSMIS;
1107   d[frq_mean] = X_bar;
1108   d[frq_sum] = X_bar * W;
1109   d[frq_variance] = M2 / (W - 1);
1110   d[frq_stddev] = sqrt (d[frq_variance]);
1111   d[frq_semean] = d[frq_stddev] / sqrt (W);
1112   if (W >= 3.0 && d[frq_variance] > 0)
1113     {
1114       double S = d[frq_stddev];
1115       d[frq_skew] = (W * M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
1116       d[frq_seskew] = sqrt (6.0 * W * (W - 1.0)
1117                             / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
1118     }
1119   else
1120     {
1121       d[frq_skew] = d[frq_seskew] = SYSMIS;
1122     }
1123   if (W >= 4.0 && d[frq_variance] > 0)
1124     {
1125       double S2 = d[frq_variance];
1126       double SE_g1 = d[frq_seskew];
1127
1128       d[frq_kurt] = ((W * (W + 1.0) * M4 - 3.0 * M2 * M2 * (W - 1.0))
1129                      / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2));
1130       d[frq_sekurt] = sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1)
1131                             / ((W - 3.0) * (W + 5.0)));
1132     }
1133   else
1134     {
1135       d[frq_kurt] = d[frq_sekurt] = SYSMIS;
1136     }
1137 }
1138
1139 /* Displays a table of all the statistics requested for variable V. */
1140 static void
1141 dump_statistics (struct variable * v, int show_varname)
1142 {
1143   double stat_value[frq_n_stats];
1144   struct tab_table *t;
1145   int i, r;
1146
1147   if (v->type == ALPHA)
1148     return;
1149   if (v->p.frq.tab.n_valid == 0)
1150     {
1151       msg (SW, _("No valid data for variable %s; statistics not displayed."),
1152            v->name);
1153       return;
1154     }
1155   calc_stats (v, stat_value);
1156
1157   t = tab_create (2, n_stats, 0);
1158   tab_dim (t, tab_natural_dimensions);
1159   tab_vline (t, TAL_1 | TAL_SPACING, 1, 0, n_stats - 1);
1160   for (i = r = 0; i < frq_n_stats; i++)
1161     if (stats & BIT_INDEX (i))
1162       {
1163         tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1164                       gettext (st_name[i].s10));
1165         tab_float (t, 1, r, TAB_NONE, stat_value[i], 11, 3);
1166         r++;
1167       }
1168
1169   tab_columns (t, SOM_COL_DOWN, 1);
1170   if (show_varname)
1171     {
1172       if (v->label)
1173         tab_title (t, 1, "%s: %s", v->name, v->label);
1174       else
1175         tab_title (t, 0, v->name);
1176     }
1177   else
1178     tab_flags (t, SOMF_NO_TITLE);
1179
1180   tab_submit (t);
1181 }
1182 \f
1183 #if 0
1184 /* Statistical calculation. */
1185
1186 static int degree[6];
1187 static int maxdegree, minmax;
1188
1189 static void stat_func (struct freq *, VISIT, int);
1190 static void calc_stats (int);
1191 static void display_stats (int);
1192
1193 /* mapping of data[]:
1194  * 0=>8
1195  * 1=>9
1196  * 2=>10
1197  * index 3: number of modes found (detects multiple modes)
1198  * index 4: number of nodes processed, for calculation of median
1199  * 5=>11
1200  *
1201  * mapping of dbl[]:
1202  * index 0-3: sum of X**i
1203  * index 4: minimum
1204  * index 5: maximum
1205  * index 6: mode
1206  * index 7: median
1207  * index 8: number of cases, valid and missing
1208  * index 9: number of valid cases
1209  * index 10: maximum frequency found, for calculation of mode
1210  * index 11: maximum frequency
1211  */
1212 static void
1213 out_stats (int i)
1214 {
1215   int j;
1216
1217   if (cur_var->type == ALPHA)
1218     return;
1219   for (j = 0; j < 8; j++)
1220     cur_var->dbl[j] = 0.;
1221   cur_var->dbl[10] = 0;
1222   cur_var->dbl[4] = DBL_MAX;
1223   cur_var->dbl[5] = -DBL_MAX;
1224   for (j = 2; j < 5; j++)
1225     cur_var->data[j] = 0;
1226   cur_var->p.frq.median_ncases = cur_var->p.frq.t.valid_cases / 2;
1227   avlwalk (cur_var->p.frq.t.f, stat_func, LEFT_TO_RIGHT);
1228   calc_stats (i);
1229   display_stats (i);
1230 }
1231
1232 static void
1233 calc_stats (int i)
1234 {
1235   struct variable *v;
1236   double n;
1237   double *d;
1238
1239   v = v_variables[i];
1240   n = v->p.frq.t.valid_cases;
1241   d = v->dbl;
1242
1243   if (n < 2 || (n < 3 && stat[FRQ_ST_7]))
1244     {
1245       warn (_("only %g case%s for variable %s, statistics not "
1246             "computed"), n, n == 1 ? "" : "s", v->name);
1247       return;
1248     }
1249   if (stat[FRQ_ST_9])
1250     v->res[FRQ_ST_9] = d[5] - d[4];
1251   if (stat[FRQ_ST_10])
1252     v->res[FRQ_ST_10] = d[4];
1253   if (stat[FRQ_ST_11])
1254     v->res[FRQ_ST_11] = d[5];
1255   if (stat[FRQ_ST_12])
1256     v->res[FRQ_ST_12] = d[0];
1257   if (stat[FRQ_ST_1] || stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_6] || stat[FRQ_ST_7])
1258     {
1259       v->res[FRQ_ST_1] = calc_mean (d, n);
1260       v->res[FRQ_ST_6] = calc_variance (d, n);
1261     }
1262   if (stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_7])
1263     v->res[FRQ_ST_5] = calc_stddev (v->res[FRQ_ST_6]);
1264   if (stat[FRQ_ST_2])
1265     v->res[FRQ_ST_2] = calc_semean (v->res[FRQ_ST_5], n);
1266   if (stat[FRQ_ST_7])
1267     {
1268       v->res[FRQ_ST_7] = calc_kurt (d, n, v->res[FRQ_ST_6]);
1269       v->res[FRQ_ST_14] = calc_sekurt (n);
1270     }
1271   if (stat[FRQ_ST_8])
1272     {
1273       v->res[FRQ_ST_8] = calc_skew (d, n, v->res[FRQ_ST_5]);
1274       v->res[FRQ_ST_15] = calc_seskew (n);
1275     }
1276   if (stat[FRQ_ST_MODE])
1277     {
1278       v->res[FRQ_ST_MODE] = v->dbl[6];
1279       if (v->data[3] > 1)
1280         warn (_("The variable %s has %d modes.  The lowest of these "
1281               "is the one given in the table."), v->name, v->data[3]);
1282     }
1283   if (stat[FRQ_ST_MEDIAN])
1284     v->res[FRQ_ST_MEDIAN] = v->dbl[7];
1285 }
1286
1287 static void
1288 stat_func (struct freq * x, VISIT order, int param)
1289 {
1290   double d, f;
1291
1292   if (order != INORDER)
1293     return;
1294   f = d = x->v.f;
1295   cur_var->dbl[0] += (d * x->c);
1296   switch (maxdegree)
1297     {
1298     case 1:
1299       f *= d;
1300       cur_var->dbl[1] += (f * x->c);
1301       break;
1302     case 2:
1303       f *= d;
1304       cur_var->dbl[1] += (f * x->c);
1305       f *= d;
1306       cur_var->dbl[2] += (f * x->c);
1307       break;
1308     case 3:
1309       f *= d;
1310       cur_var->dbl[1] += (f * x->c);
1311       f *= d;
1312       cur_var->dbl[2] += (f * x->c);
1313       f *= d;
1314       cur_var->dbl[3] += (f * x->c);
1315       break;
1316     }
1317   if (minmax)
1318     {
1319       if (d < cur_var->dbl[4])
1320         cur_var->dbl[4] = d;
1321       if (d > cur_var->dbl[5])
1322         cur_var->dbl[5] = d;
1323     }
1324   if (x->c > cur_var->dbl[10])
1325     {
1326       cur_var->data[3] = 1;
1327       cur_var->dbl[10] = x->c;
1328       cur_var->dbl[6] = x->v.f;
1329     }
1330   else if (x->c == cur_var->dbl[10])
1331     cur_var->data[3]++;
1332   if (cur_var->data[4] < cur_var->p.frq.median_ncases
1333       && cur_var->data[4] + x->c >= cur_var->p.frq.median_ncases)
1334     cur_var->dbl[7] = x->v.f;
1335   cur_var->data[4] += x->c;
1336 }
1337 \f
1338 /* Statistical display. */
1339 static int column, ncolumns;
1340
1341 static void outstat (char *, double);
1342
1343 static void
1344 display_stats (int i)
1345 {
1346   statname *sp;
1347   struct variable *v;
1348   int nlines;
1349
1350   v = v_variables[i];
1351   ncolumns = (margin_width + 3) / 26;
1352   if (ncolumns < 1)
1353     ncolumns = 1;
1354   nlines = sc / ncolumns + (sc % ncolumns > 0);
1355   if (nlines == 2 && sc == 4)
1356     ncolumns = 2;
1357   if (nlines == 3 && sc == 9)
1358     ncolumns = 3;
1359   if (nlines == 4 && sc == 12)
1360     ncolumns = 3;
1361   column = 0;
1362   for (sp = st_name; sp->s != -1; sp++)
1363     if (stat[sp->s] == 1)
1364       outstat (gettext (sp->s10), v->res[sp->s]);
1365   if (column)
1366     out_eol ();
1367   blank_line ();
1368 }
1369
1370 static void
1371 outstat (char *label, double value)
1372 {
1373   char buf[128], *cp;
1374   int dw, n;
1375
1376   cp = &buf[0];
1377   if (!column)
1378     out_header ();
1379   else
1380     {
1381       memset (buf, ' ', 3);
1382       cp = &buf[3];
1383     }
1384   dw = 4;
1385   n = nsprintf (cp, "%-10s %12.4f", label, value);
1386   while (n > 23 && dw > 0)
1387     n = nsprintf (cp, "%-10s %12.*f", label, --dw, value);
1388   outs (buf);
1389   column++;
1390   if (column == ncolumns)
1391     {
1392       column = 0;
1393       out_eol ();
1394     }
1395 }
1396 \f
1397 /* Graphs. */
1398
1399 static rect pb, gb;             /* Page border, graph border. */
1400 static int px, py;              /* Page width, height. */
1401 static int ix, iy;              /* Inch width, height. */
1402
1403 static void draw_bar_chart (int);
1404 static void draw_histogram (int);
1405 static int scale_dep_axis (int);
1406
1407 static void
1408 out_graphs (int i)
1409 {
1410   struct variable *v;
1411
1412   v = v_variables[i];
1413   if (avlcount (cur_var->p.frq.t.f) < 2
1414       || (chart == HIST && v_variables[i]->type == ALPHA))
1415     return;
1416   if (driver_id && set_highres == 1)
1417     {
1418       char *text;
1419
1420       graf_page_size (&px, &py, &ix, &iy);
1421       graf_feed_page ();
1422
1423       /* Calculate borders. */
1424       pb.x1 = ix;
1425       pb.y1 = iy;
1426       pb.x2 = px - ix;
1427       pb.y2 = py - iy;
1428       gb.x1 = pb.x1 + ix;
1429       gb.y1 = pb.y1 + iy;
1430       gb.x2 = pb.x2 - ix / 2;
1431       gb.y2 = pb.y2 - iy;
1432
1433       /* Draw borders. */
1434       graf_frame_rect (COMPONENTS (pb));
1435       graf_frame_rect (COMPONENTS (gb));
1436
1437       /* Draw axis labels. */
1438       graf_font_size (iy / 4);  /* 18-point text */
1439       text = format == PERCENT ? _("Percentage") : _("Frequency");
1440       graf_text (pb.x1 + max (ix, iy) / 4 + max (ix, iy) / 16, gb.y2, text,
1441                  SIDEWAYS);
1442       text = v->label ? v->label : v->name;
1443       graf_text (gb.x1, pb.y2 - iy / 4, text, UPRIGHT);
1444
1445       /* Draw axes, chart proper. */
1446       if (chart == BAR ||
1447           (chart == HBAR
1448        && (avlcount (cur_var->p.frq.t.f) || v_variables[i]->type == ALPHA)))
1449         draw_bar_chart (i);
1450       else
1451         draw_histogram (i);
1452
1453       graf_eject_page ();
1454     }
1455   if (set_lowres == 1 || (set_lowres == 2 && (!driver_id || !set_highres)))
1456     {
1457       static warned;
1458
1459       /* Do character-based graphs. */
1460       if (!warned)
1461         {
1462           warn (_("low-res graphs not implemented"));
1463           warned = 1;
1464         }
1465     }
1466 }
1467
1468 #if __GNUC__ && !__CHECKER__
1469 #define BIG_TYPE long long
1470 #else /* !__GNUC__ */
1471 #define BIG_TYPE double
1472 #endif /* !__GNUC__ */
1473
1474 static void
1475 draw_bar_chart (int i)
1476 {
1477   int bar_width, bar_spacing;
1478   int w, max, row;
1479   double val;
1480   struct freq *f;
1481   rect r;
1482   AVLtraverser *t = NULL;
1483
1484   w = (px - ix * 7 / 2) / avlcount (cur_var->p.frq.t.f);
1485   bar_width = w * 2 / 3;
1486   bar_spacing = w - bar_width;
1487
1488 #if !ALLOW_HUGE_BARS
1489   if (bar_width > ix / 2)
1490     bar_width = ix / 2;
1491 #endif /* !ALLOW_HUGE_BARS */
1492
1493   max = scale_dep_axis (cur_var->p.frq.t.max_freq);
1494
1495   row = 0;
1496   r.x1 = gb.x1 + bar_spacing / 2;
1497   r.x2 = r.x1 + bar_width;
1498   r.y2 = gb.y2;
1499   graf_fill_color (255, 0, 0);
1500   for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1501        f = avltrav (cur_var->p.frq.t.f, &t))
1502     {
1503       char buf2[64];
1504       char *buf;
1505
1506       val = f->c;
1507       if (format == PERCENT)
1508         val = val * 100 / cur_var->p.frq.t.valid_cases;
1509       r.y1 = r.y2 - val * (height (gb) - 1) / max;
1510       graf_fill_rect (COMPONENTS (r));
1511       graf_frame_rect (COMPONENTS (r));
1512       buf = get_val_lab (cur_var, f->v, 0);
1513       if (!buf)
1514         if (cur_var->type == ALPHA)
1515           buf = f->v.s;
1516         else
1517           {
1518             sprintf (buf2, "%g", f->v.f);
1519             buf = buf2;
1520           }
1521       graf_text (r.x1 + bar_width / 2,
1522                  gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1523       row ^= 1;
1524       r.x1 += bar_width + bar_spacing;
1525       r.x2 += bar_width + bar_spacing;
1526     }
1527   graf_fill_color (0, 0, 0);
1528 }
1529
1530 #define round_down(X, V)                        \
1531         (floor ((X) / (V)) * (V))
1532 #define round_up(X, V)                          \
1533         (ceil ((X) / (V)) * (V))
1534
1535 static void
1536 draw_histogram (int i)
1537 {
1538   double lower, upper, interval;
1539   int bars[MAX_HIST_BARS + 1], top, j;
1540   int err, addend, rem, nbars, row, max_freq;
1541   char buf[25];
1542   rect r;
1543   struct freq *f;
1544   AVLtraverser *t = NULL;
1545
1546   lower = min == SYSMIS ? cur_var->dbl[4] : min;
1547   upper = max == SYSMIS ? cur_var->dbl[5] : max;
1548   if (upper - lower >= 10)
1549     {
1550       double l, u;
1551
1552       u = round_up (upper, 5);
1553       l = round_down (lower, 5);
1554       nbars = (u - l) / 5;
1555       if (nbars * 2 + 1 <= MAX_HIST_BARS)
1556         {
1557           nbars *= 2;
1558           u = round_up (upper, 2.5);
1559           l = round_down (lower, 2.5);
1560           if (l + 1.25 <= lower && u - 1.25 >= upper)
1561             nbars--, lower = l + 1.25, upper = u - 1.25;
1562           else if (l + 1.25 <= lower)
1563             lower = l + 1.25, upper = u + 1.25;
1564           else if (u - 1.25 >= upper)
1565             lower = l - 1.25, upper = u - 1.25;
1566           else
1567             nbars++, lower = l - 1.25, upper = u + 1.25;
1568         }
1569       else if (nbars < MAX_HIST_BARS)
1570         {
1571           if (l + 2.5 <= lower && u - 2.5 >= upper)
1572             nbars--, lower = l + 2.5, upper = u - 2.5;
1573           else if (l + 2.5 <= lower)
1574             lower = l + 2.5, upper = u + 2.5;
1575           else if (u - 2.5 >= upper)
1576             lower = l - 2.5, upper = u - 2.5;
1577           else
1578             nbars++, lower = l - 2.5, upper = u + 2.5;
1579         }
1580       else
1581         nbars = MAX_HIST_BARS;
1582     }
1583   else
1584     {
1585       nbars = avlcount (cur_var->p.frq.t.f);
1586       if (nbars > MAX_HIST_BARS)
1587         nbars = MAX_HIST_BARS;
1588     }
1589   if (nbars < MIN_HIST_BARS)
1590     nbars = MIN_HIST_BARS;
1591   interval = (upper - lower) / nbars;
1592
1593   memset (bars, 0, sizeof (int[nbars + 1]));
1594   if (lower >= upper)
1595     {
1596       msg (SE, _("Could not make histogram for %s for specified "
1597            "minimum %g and maximum %g; please discard graph."), cur_var->name,
1598            lower, upper);
1599       return;
1600     }
1601   for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1602        f = avltrav (cur_var->p.frq.t.f, &t))
1603     if (f->v.f == upper)
1604       bars[nbars - 1] += f->c;
1605     else if (f->v.f >= lower && f->v.f < upper)
1606       bars[(int) ((f->v.f - lower) / interval)] += f->c;
1607   bars[nbars - 1] += bars[nbars];
1608   for (j = top = 0; j < nbars; j++)
1609     if (bars[j] > top)
1610       top = bars[j];
1611   max_freq = top;
1612   top = scale_dep_axis (top);
1613
1614   err = row = 0;
1615   addend = width (gb) / nbars;
1616   rem = width (gb) % nbars;
1617   r.x1 = gb.x1;
1618   r.x2 = r.x1 + addend;
1619   r.y2 = gb.y2;
1620   err += rem;
1621   graf_fill_color (255, 0, 0);
1622   for (j = 0; j < nbars; j++)
1623     {
1624       int w;
1625
1626       r.y1 = r.y2 - (BIG_TYPE) bars[j] * (height (gb) - 1) / top;
1627       graf_fill_rect (COMPONENTS (r));
1628       graf_frame_rect (COMPONENTS (r));
1629       sprintf (buf, "%g", lower + interval / 2 + interval * j);
1630       graf_text (r.x1 + addend / 2,
1631                  gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1632       row ^= 1;
1633       w = addend;
1634       err += rem;
1635       while (err >= addend)
1636         {
1637           w++;
1638           err -= addend;
1639         }
1640       r.x1 = r.x2;
1641       r.x2 = r.x1 + w;
1642     }
1643   if (normal)
1644     {
1645       double x, y, variance, mean, step, factor;
1646
1647       variance = cur_var->res[FRQ_ST_VARIANCE];
1648       mean = cur_var->res[FRQ_ST_MEAN];
1649       factor = (1. / (sqrt (2. * PI * variance))
1650                 * cur_var->p.frq.t.valid_cases * interval);
1651       graf_polyline_begin ();
1652       for (x = lower, step = (upper - lower) / (POLYLINE_DENSITY);
1653            x <= upper; x += step)
1654         {
1655           y = factor * exp (-square (x - mean) / (2. * variance));
1656           debug_printf (("(%20.10f, %20.10f)\n", x, y));
1657           graf_polyline_point (gb.x1 + (x - lower) / (upper - lower) * width (gb),
1658                                gb.y2 - y * (height (gb) - 1) / top);
1659         }
1660       graf_polyline_end ();
1661     }
1662   graf_fill_color (0, 0, 0);
1663 }
1664
1665 static int
1666 scale_dep_axis (int max)
1667 {
1668   int j, s, x, y, ty, by;
1669   char buf[10];
1670
1671   x = 10, s = 2;
1672   if (scale != SYSMIS && max < scale)
1673     x = scale, s = scale / 5;
1674   else if (format == PERCENT)
1675     {
1676       max = ((BIG_TYPE) 100 * cur_var->p.frq.t.max_freq
1677              / cur_var->p.frq.t.valid_cases + 1);
1678       if (max < 5)
1679         x = 5, s = 1;
1680       else if (max < 10)
1681         x = 10, s = 2;
1682       else if (max < 25)
1683         x = 25, s = 5;
1684       else if (max < 50)
1685         x = 50, s = 10;
1686       else
1687         max = 100, s = 20;
1688     }
1689   else                          /* format==FREQ */
1690     /* Uses a progression of 10, 20, 50, 100, 200, 500, ... */
1691     for (;;)
1692       {
1693         if (x > max)
1694           break;
1695         x *= 2;
1696         s *= 2;
1697         if (x > max)
1698           break;
1699         x = x / 2 * 5;
1700         s = s / 2 * 5;
1701         if (x > max)
1702           break;
1703         x *= 2;
1704         s *= 2;
1705       }
1706   graf_font_size (iy / 9);      /* 8-pt text */
1707   for (j = 0; j <= x; j += s)
1708     {
1709       y = gb.y2 - (BIG_TYPE) j *(height (gb) - 1) / x;
1710       ty = y - iy / 64;
1711       by = y + iy / 64;
1712       if (ty < gb.y1)
1713         ty += iy / 64, by += iy / 64;
1714       else if (by > gb.y2)
1715         ty -= iy / 64, by -= iy / 64;
1716       graf_fill_rect (gb.x1 - ix / 16, ty, gb.x1, by);
1717       sprintf (buf, "%d", j);
1718       graf_text (gb.x1 - ix / 8, (ty + by) / 2, buf, CRJUST);
1719     }
1720   return x;
1721 }
1722 \f
1723 /* Percentiles. */
1724
1725 static void ungrouped_pcnt (int i);
1726 static int grouped_interval_pcnt (int i);
1727 static void out_pcnt (double, double);
1728
1729 static void
1730 out_percentiles (int i)
1731 {
1732   if (cur_var->type == ALPHA || !n_percentiles)
1733     return;
1734
1735   outs_line (_("Percentile    Value     "
1736              "Percentile    Value     "
1737              "Percentile    Value"));
1738   blank_line ();
1739
1740   column = 0;
1741   if (!g_var[i])
1742     ungrouped_pcnt (i);
1743   else if (g_var[i] == 1)
1744     grouped_interval_pcnt (i);
1745 #if 0
1746   else if (g_var[i] == -1)
1747     grouped_pcnt (i);
1748   else
1749     grouped_boundaries_pcnt (i);
1750 #else /* !0 */
1751   else
1752     warn (_("this form of percentiles not supported"));
1753 #endif
1754   if (column)
1755     out_eol ();
1756 }
1757
1758 static void
1759 out_pcnt (double pcnt, double value)
1760 {
1761   if (!column)
1762     out_header ();
1763   else
1764     outs ("     ");
1765   out ("%7.2f%13.3f", pcnt * 100., value);
1766   column++;
1767   if (column == 3)
1768     {
1769       out_eol ();
1770       column = 0;
1771     }
1772 }
1773
1774 static void
1775 ungrouped_pcnt (int i)
1776 {
1777   AVLtraverser *t = NULL;
1778   struct freq *f;
1779   double *p, *e;
1780   int sum;
1781
1782   p = percentiles;
1783   e = &percentiles[n_percentiles];
1784   sum = 0;
1785   for (f = avltrav (cur_var->p.frq.t.f, &t);
1786        f && p < e; f = avltrav (cur_var->p.frq.t.f, &t))
1787     {
1788       sum += f->c;
1789       while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1790         out_pcnt (*p++, f->v.f);
1791     }
1792 }
1793
1794
1795 static int
1796 grouped_interval_pcnt (int i)
1797 {
1798   AVLtraverser * t = NULL;
1799   struct freq * f, *fp;
1800   double *p, *e, w;
1801   int sum, psum;
1802
1803   p = percentiles;
1804   e = &percentiles[n_percentiles];
1805   w = gl_var[i][0];
1806   sum = psum = 0;
1807   for (fp = 0, f = avltrav (cur_var->p.frq.t.f, &t);
1808        f && p < e;
1809        fp = f, f = avltrav (cur_var->p.frq.t.f, &t))
1810     {
1811       if (fp)
1812         if (fabs (f->v.f - fp->v.f) < w)
1813           {
1814             out_eol ();
1815             column = 0;
1816             return msg (SE, _("Difference between %g and %g is "
1817                               "too small for grouping interval %g."), f->v.f,
1818                         fp->v.f, w);
1819           }
1820       psum = sum;
1821       sum += f->c;
1822       while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1823         {
1824           out_pcnt (p[0], (((p[0] * cur_var->p.frq.t.valid_cases) - psum) * w / f->c
1825                            + (f->v.f - w / 2)));
1826           p++;
1827         }
1828     }
1829   return 1;
1830 }
1831 #endif
1832
1833 /*
1834    Local Variables:
1835    mode: c
1836    End:
1837 */