treewide: Replace <name>_cnt by n_<name>s and <name>_cap by allocated_<name>.
[pspp] / src / language / stats / descriptives.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 1997-2000, 2009-2014 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 #include <float.h>
20 #include <limits.h>
21 #include <math.h>
22 #include <stdlib.h>
23
24 #include "data/casegrouper.h"
25 #include "data/casereader.h"
26 #include "data/casewriter.h"
27 #include "data/dataset.h"
28 #include "data/dictionary.h"
29 #include "data/transformations.h"
30 #include "data/variable.h"
31 #include "language/command.h"
32 #include "language/dictionary/split-file.h"
33 #include "language/lexer/lexer.h"
34 #include "language/lexer/variable-parser.h"
35 #include "libpspp/array.h"
36 #include "libpspp/assertion.h"
37 #include "libpspp/compiler.h"
38 #include "libpspp/i18n.h"
39 #include "libpspp/message.h"
40 #include "math/moments.h"
41 #include "output/pivot-table.h"
42
43 #include "gl/xalloc.h"
44
45 #include "gettext.h"
46 #define _(msgid) gettext (msgid)
47 #define N_(msgid) msgid
48
49 /* DESCRIPTIVES private data. */
50
51 struct dsc_proc;
52
53 /* Handling of missing values. */
54 enum dsc_missing_type
55   {
56     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
57     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
58   };
59
60 /* Describes properties of a distribution for the purpose of
61    calculating a Z-score. */
62 struct dsc_z_score
63   {
64     const struct variable *src_var;   /* Variable on which z-score is based. */
65     struct variable *z_var;     /* New z-score variable. */
66     double mean;                /* Distribution mean. */
67     double std_dev;             /* Distribution standard deviation. */
68   };
69
70 /* DESCRIPTIVES transformation (for calculating Z-scores). */
71 struct dsc_trns
72   {
73     struct dsc_z_score *z_scores; /* Array of Z-scores. */
74     int n_z_scores;               /* Number of Z-scores. */
75     const struct variable **vars;     /* Variables for listwise missing checks. */
76     size_t n_vars;              /* Number of variables. */
77     enum dsc_missing_type missing_type; /* Treatment of missing values. */
78     enum mv_class exclude;      /* Classes of missing values to exclude. */
79     const struct variable *filter;    /* Dictionary FILTER BY variable. */
80     struct casereader *z_reader; /* Reader for count, mean, stddev. */
81     casenumber count;            /* Number left in this SPLIT FILE group.*/
82     bool ok;
83   };
84
85 /* Statistics.  Used as bit indexes, so must be 32 or fewer. */
86 enum dsc_statistic
87   {
88     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
89     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
90     DSC_MAX, DSC_SUM, DSC_N_STATS,
91
92     /* Only valid as sort criteria. */
93     DSC_NAME = -2,              /* Sort by name. */
94     DSC_NONE = -1               /* Unsorted. */
95   };
96
97 /* Describes one statistic. */
98 struct dsc_statistic_info
99   {
100     const char *identifier;     /* Identifier. */
101     const char *name;           /* Full name. */
102     enum moment moment;         /* Highest moment needed to calculate. */
103   };
104
105 /* Table of statistics, indexed by DSC_*. */
106 static const struct dsc_statistic_info dsc_info[DSC_N_STATS] =
107   {
108     {"MEAN", N_("Mean"), MOMENT_MEAN},
109     {"SEMEAN", N_("S.E. Mean"), MOMENT_VARIANCE},
110     {"STDDEV", N_("Std Dev"), MOMENT_VARIANCE},
111     {"VARIANCE", N_("Variance"), MOMENT_VARIANCE},
112     {"KURTOSIS", N_("Kurtosis"), MOMENT_KURTOSIS},
113     {"SEKURTOSIS", N_("S.E. Kurt"), MOMENT_NONE},
114     {"SKEWNESS", N_("Skewness"), MOMENT_SKEWNESS},
115     {"SESKEWNESS", N_("S.E. Skew"), MOMENT_NONE},
116     {"RANGE", N_("Range"), MOMENT_NONE},
117     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
118     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
119     {"SUM", N_("Sum"), MOMENT_MEAN},
120   };
121
122 /* Statistics calculated by default if none are explicitly
123    requested. */
124 #define DEFAULT_STATS                                                   \
125         ((1ul << DSC_MEAN) | (1ul << DSC_STDDEV) | (1ul << DSC_MIN)     \
126          | (1ul << DSC_MAX))
127
128 /* A variable specified on DESCRIPTIVES. */
129 struct dsc_var
130   {
131     const struct variable *v;         /* Variable to calculate on. */
132     char *z_name;                     /* Name for z-score variable. */
133     double valid, missing;      /* Valid, missing counts. */
134     struct moments *moments;    /* Moments. */
135     double min, max;            /* Maximum and mimimum values. */
136     double stats[DSC_N_STATS];  /* All the stats' values. */
137   };
138
139 /* A DESCRIPTIVES procedure. */
140 struct dsc_proc
141   {
142     /* Per-variable info. */
143     struct dictionary *dict;    /* Dictionary. */
144     struct dsc_var *vars;       /* Variables. */
145     size_t n_vars;              /* Number of variables. */
146
147     /* User options. */
148     enum dsc_missing_type missing_type; /* Treatment of missing values. */
149     enum mv_class exclude;      /* Classes of missing values to exclude. */
150
151     /* Accumulated results. */
152     double missing_listwise;    /* Sum of weights of cases missing listwise. */
153     double valid;               /* Sum of weights of valid cases. */
154     bool bad_warn;               /* Warn if bad weight found. */
155     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
156     int sort_ascending;         /* !0: ascending order; 0: descending. */
157     unsigned long show_stats;   /* Statistics to display. */
158     unsigned long calc_stats;   /* Statistics to calculate. */
159     enum moment max_moment;     /* Highest moment needed for stats. */
160
161     /* Z scores. */
162     struct casewriter *z_writer; /* Mean and stddev per SPLIT FILE group. */
163   };
164
165 /* Parsing. */
166 static enum dsc_statistic match_statistic (struct lexer *);
167 static void free_dsc_proc (struct dsc_proc *);
168
169 /* Z-score functions. */
170 static bool try_name (const struct dictionary *dict,
171                       struct dsc_proc *dsc, const char *name);
172 static char *generate_z_varname (const struct dictionary *dict,
173                                  struct dsc_proc *dsc,
174                                  const char *name, int *n_zs);
175 static void dump_z_table (struct dsc_proc *);
176 static void setup_z_trns (struct dsc_proc *, struct dataset *);
177
178 /* Procedure execution functions. */
179 static void calc_descriptives (struct dsc_proc *, struct casereader *,
180                                struct dataset *);
181 static void display (struct dsc_proc *dsc);
182 \f
183 /* Parser and outline. */
184
185 /* Handles DESCRIPTIVES. */
186 int
187 cmd_descriptives (struct lexer *lexer, struct dataset *ds)
188 {
189   struct dictionary *dict = dataset_dict (ds);
190   struct dsc_proc *dsc;
191   const struct variable **vars = NULL;
192   size_t n_vars = 0;
193   int save_z_scores = 0;
194   int n_zs = 0;
195   size_t i;
196   bool ok;
197
198   struct casegrouper *grouper;
199   struct casereader *group;
200
201   /* Create and initialize dsc. */
202   dsc = xmalloc (sizeof *dsc);
203   dsc->dict = dict;
204   dsc->vars = NULL;
205   dsc->n_vars = 0;
206   dsc->missing_type = DSC_VARIABLE;
207   dsc->exclude = MV_ANY;
208   dsc->missing_listwise = 0.;
209   dsc->valid = 0.;
210   dsc->bad_warn = 1;
211   dsc->sort_by_stat = DSC_NONE;
212   dsc->sort_ascending = 1;
213   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
214   dsc->z_writer = NULL;
215
216   /* Parse DESCRIPTIVES. */
217   while (lex_token (lexer) != T_ENDCMD)
218     {
219       if (lex_match_id (lexer, "MISSING"))
220         {
221           lex_match (lexer, T_EQUALS);
222           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
223             {
224               if (lex_match_id (lexer, "VARIABLE"))
225                 dsc->missing_type = DSC_VARIABLE;
226               else if (lex_match_id (lexer, "LISTWISE"))
227                 dsc->missing_type = DSC_LISTWISE;
228               else if (lex_match_id (lexer, "INCLUDE"))
229                 dsc->exclude = MV_SYSTEM;
230               else
231                 {
232                   lex_error (lexer, NULL);
233                   goto error;
234                 }
235               lex_match (lexer, T_COMMA);
236             }
237         }
238       else if (lex_match_id (lexer, "SAVE"))
239         save_z_scores = 1;
240       else if (lex_match_id (lexer, "FORMAT"))
241         {
242           lex_match (lexer, T_EQUALS);
243           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
244             {
245               if (lex_match_id (lexer, "LABELS")
246                   || lex_match_id (lexer, "NOLABELS")
247                   || lex_match_id (lexer, "INDEX")
248                   || lex_match_id (lexer, "NOINDEX")
249                   || lex_match_id (lexer, "LINE")
250                   || lex_match_id (lexer, "SERIAL"))
251                 {
252                   /* Ignore. */
253                 }
254               else
255                 {
256                   lex_error (lexer, NULL);
257                   goto error;
258                 }
259               lex_match (lexer, T_COMMA);
260             }
261         }
262       else if (lex_match_id (lexer, "STATISTICS"))
263         {
264           lex_match (lexer, T_EQUALS);
265           dsc->show_stats = 0;
266           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
267             {
268               if (lex_match (lexer, T_ALL))
269                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
270               else if (lex_match_id (lexer, "DEFAULT"))
271                 dsc->show_stats |= DEFAULT_STATS;
272               else
273                 {
274                   enum dsc_statistic s = match_statistic (lexer);
275                   if (s == DSC_NONE)
276                     {
277                       lex_error (lexer, NULL);
278                       goto error;
279                     }
280                   dsc->show_stats |= 1ul << s;
281                 }
282               lex_match (lexer, T_COMMA);
283             }
284           if (dsc->show_stats == 0)
285             dsc->show_stats = DEFAULT_STATS;
286         }
287       else if (lex_match_id (lexer, "SORT"))
288         {
289           lex_match (lexer, T_EQUALS);
290           if (lex_match_id (lexer, "NAME"))
291             dsc->sort_by_stat = DSC_NAME;
292           else
293             {
294               dsc->sort_by_stat = match_statistic (lexer);
295               if (dsc->sort_by_stat == DSC_NONE)
296                 dsc->sort_by_stat = DSC_MEAN;
297             }
298           if (lex_match (lexer, T_LPAREN))
299             {
300               if (lex_match_id (lexer, "A"))
301                 dsc->sort_ascending = 1;
302               else if (lex_match_id (lexer, "D"))
303                 dsc->sort_ascending = 0;
304               else
305                 lex_error (lexer, NULL);
306               if (! lex_force_match (lexer, T_RPAREN))
307                 goto error;
308             }
309         }
310       else if (n_vars == 0)
311         {
312           if (lex_next_token (lexer, 1) == T_EQUALS)
313             {
314               lex_match_id (lexer, "VARIABLES");
315               lex_match (lexer, T_EQUALS);
316             }
317
318           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
319             {
320               int i;
321
322               if (!parse_variables_const (lexer, dict, &vars, &n_vars,
323                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
324                 goto error;
325
326               dsc->vars = xnrealloc ((void *)dsc->vars, n_vars, sizeof *dsc->vars);
327               for (i = dsc->n_vars; i < n_vars; i++)
328                 {
329                   struct dsc_var *dv = &dsc->vars[i];
330                   dv->v = vars[i];
331                   dv->z_name = NULL;
332                   dv->moments = NULL;
333                 }
334               dsc->n_vars = n_vars;
335
336               if (lex_match (lexer, T_LPAREN))
337                 {
338                   if (lex_token (lexer) != T_ID)
339                     {
340                       lex_error (lexer, NULL);
341                       goto error;
342                     }
343                   if (try_name (dict, dsc, lex_tokcstr (lexer)))
344                     {
345                       struct dsc_var *dsc_var = &dsc->vars[dsc->n_vars - 1];
346                       dsc_var->z_name = xstrdup (lex_tokcstr (lexer));
347                       n_zs++;
348                     }
349                   else
350                     msg (SE, _("Z-score variable name %s would be"
351                                " a duplicate variable name."), lex_tokcstr (lexer));
352                   lex_get (lexer);
353                   if (!lex_force_match (lexer, T_RPAREN))
354                     goto error;
355                 }
356             }
357         }
358       else
359         {
360           lex_error (lexer, NULL);
361           goto error;
362         }
363
364       lex_match (lexer, T_SLASH);
365     }
366   if (n_vars == 0)
367     {
368       msg (SE, _("No variables specified."));
369       goto error;
370     }
371
372   /* Construct z-score varnames, show translation table. */
373   if (n_zs || save_z_scores)
374     {
375       struct caseproto *proto;
376
377       if (save_z_scores)
378         {
379           int n_gens = 0;
380
381           for (i = 0; i < dsc->n_vars; i++)
382             {
383               struct dsc_var *dsc_var = &dsc->vars[i];
384               if (dsc_var->z_name == NULL)
385                 {
386                   const char *name = var_get_name (dsc_var->v);
387                   dsc_var->z_name = generate_z_varname (dict, dsc, name,
388                                                         &n_gens);
389                   if (dsc_var->z_name == NULL)
390                     goto error;
391
392                   n_zs++;
393                 }
394             }
395         }
396
397       /* It would be better to handle Z scores correctly (however we define
398          that) when TEMPORARY is in effect, but in the meantime this at least
399          prevents a use-after-free error.  See bug #38786.  */
400       if (proc_make_temporary_transformations_permanent (ds))
401         msg (SW, _("DESCRIPTIVES with Z scores ignores TEMPORARY.  "
402                    "Temporary transformations will be made permanent."));
403
404       proto = caseproto_create ();
405       for (i = 0; i < 1 + 2 * n_zs; i++)
406         proto = caseproto_add_width (proto, 0);
407       dsc->z_writer = autopaging_writer_create (proto);
408       caseproto_unref (proto);
409
410       dump_z_table (dsc);
411     }
412
413   /* Figure out statistics to display. */
414   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
415     dsc->show_stats |= 1ul << DSC_SESKEW;
416   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
417     dsc->show_stats |= 1ul << DSC_SEKURT;
418
419   /* Figure out which statistics to calculate. */
420   dsc->calc_stats = dsc->show_stats;
421   if (n_zs > 0)
422     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
423   if (dsc->sort_by_stat >= 0)
424     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
425   if (dsc->show_stats & (1ul << DSC_SESKEW))
426     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
427   if (dsc->show_stats & (1ul << DSC_SEKURT))
428     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
429
430   /* Figure out maximum moment needed and allocate moments for
431      the variables. */
432   dsc->max_moment = MOMENT_NONE;
433   for (i = 0; i < DSC_N_STATS; i++)
434     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
435       dsc->max_moment = dsc_info[i].moment;
436   if (dsc->max_moment != MOMENT_NONE)
437     for (i = 0; i < dsc->n_vars; i++)
438       dsc->vars[i].moments = moments_create (dsc->max_moment);
439
440   /* Data pass. */
441   grouper = casegrouper_create_splits (proc_open_filtering (ds, false), dict);
442   while (casegrouper_get_next_group (grouper, &group))
443     calc_descriptives (dsc, group, ds);
444   ok = casegrouper_destroy (grouper);
445   ok = proc_commit (ds) && ok;
446
447   /* Z-scoring! */
448   if (ok && n_zs)
449     setup_z_trns (dsc, ds);
450
451   /* Done. */
452   free (vars);
453   free_dsc_proc (dsc);
454   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
455
456  error:
457   free (vars);
458   free_dsc_proc (dsc);
459   return CMD_FAILURE;
460 }
461
462 /* Returns the statistic named by the current token and skips past the token.
463    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
464    specifiers). Emits an error if the current token ID does not name a
465    statistic. */
466 static enum dsc_statistic
467 match_statistic (struct lexer *lexer)
468 {
469   if (lex_token (lexer) == T_ID)
470     {
471       enum dsc_statistic stat;
472
473       for (stat = 0; stat < DSC_N_STATS; stat++)
474         if (lex_match_id (lexer, dsc_info[stat].identifier))
475           return stat;
476
477       lex_get (lexer);
478       lex_error (lexer, _("expecting statistic name: reverting to default"));
479     }
480
481   return DSC_NONE;
482 }
483
484 /* Frees DSC. */
485 static void
486 free_dsc_proc (struct dsc_proc *dsc)
487 {
488   size_t i;
489
490   if (dsc == NULL)
491     return;
492
493   for (i = 0; i < dsc->n_vars; i++)
494     {
495       struct dsc_var *dsc_var = &dsc->vars[i];
496       free (dsc_var->z_name);
497       moments_destroy (dsc_var->moments);
498     }
499   casewriter_destroy (dsc->z_writer);
500   free (dsc->vars);
501   free (dsc);
502 }
503 \f
504 /* Z scores. */
505
506 /* Returns false if NAME is a duplicate of any existing variable name or
507    of any previously-declared z-var name; otherwise returns true. */
508 static bool
509 try_name (const struct dictionary *dict, struct dsc_proc *dsc,
510           const char *name)
511 {
512   size_t i;
513
514   if (dict_lookup_var (dict, name) != NULL)
515     return false;
516   for (i = 0; i < dsc->n_vars; i++)
517     {
518       struct dsc_var *dsc_var = &dsc->vars[i];
519       if (dsc_var->z_name != NULL && !utf8_strcasecmp (dsc_var->z_name, name))
520         return false;
521     }
522   return true;
523 }
524
525 /* Generates a name for a Z-score variable based on a variable
526    named VAR_NAME, given that *Z_CNT generated variable names are
527    known to already exist.  If successful, returns the new name
528    as a dynamically allocated string.  On failure, returns NULL. */
529 static char *
530 generate_z_varname (const struct dictionary *dict, struct dsc_proc *dsc,
531                     const char *var_name, int *n_zs)
532 {
533   char *z_name, *trunc_name;
534
535   /* Try a name based on the original variable name. */
536   z_name = xasprintf ("Z%s", var_name);
537   trunc_name = utf8_encoding_trunc (z_name, dict_get_encoding (dict),
538                                     ID_MAX_LEN);
539   free (z_name);
540   if (try_name (dict, dsc, trunc_name))
541     return trunc_name;
542   free (trunc_name);
543
544   /* Generate a synthetic name. */
545   for (;;)
546     {
547       char name[16];
548
549       (*n_zs)++;
550
551       if (*n_zs <= 99)
552         sprintf (name, "ZSC%03d", *n_zs);
553       else if (*n_zs <= 108)
554         sprintf (name, "STDZ%02d", *n_zs - 99);
555       else if (*n_zs <= 117)
556         sprintf (name, "ZZZZ%02d", *n_zs - 108);
557       else if (*n_zs <= 126)
558         sprintf (name, "ZQZQ%02d", *n_zs - 117);
559       else
560         {
561           msg (SE, _("Ran out of generic names for Z-score variables.  "
562                      "There are only 126 generic names: ZSC001-ZSC0999, "
563                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
564           return NULL;
565         }
566
567       if (try_name (dict, dsc, name))
568         return xstrdup (name);
569     }
570   NOT_REACHED();
571 }
572
573 /* Outputs a table describing the mapping between source
574    variables and Z-score variables. */
575 static void
576 dump_z_table (struct dsc_proc *dsc)
577 {
578   struct pivot_table *table = pivot_table_create (
579     N_("Mapping of Variables to Z-scores"));
580
581   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Names"),
582                           N_("Source"), N_("Target"));
583
584   struct pivot_dimension *names = pivot_dimension_create (
585     table, PIVOT_AXIS_ROW, N_("Variables"));
586   names->hide_all_labels = true;
587
588   for (size_t i = 0; i < dsc->n_vars; i++)
589     if (dsc->vars[i].z_name != NULL)
590       {
591         int row = pivot_category_create_leaf (names->root,
592                                               pivot_value_new_number (i));
593
594         pivot_table_put2 (table, 0, row,
595                           pivot_value_new_variable (dsc->vars[i].v));
596         pivot_table_put2 (table, 1, row,
597                           pivot_value_new_user_text (dsc->vars[i].z_name, -1));
598       }
599
600   pivot_table_submit (table);
601 }
602
603 static void
604 descriptives_set_all_sysmis_zscores (const struct dsc_trns *t, struct ccase *c)
605 {
606   const struct dsc_z_score *z;
607
608   for (z = t->z_scores; z < t->z_scores + t->n_z_scores; z++)
609     *case_num_rw (c, z->z_var) = SYSMIS;
610 }
611
612 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
613    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
614    SYSMIS 3) score is user missing and they were not included in the original
615    analyis. 4) any of the variables in the original analysis were missing
616    (either system or user-missing values that weren't included).
617 */
618 static int
619 descriptives_trns_proc (void *trns_, struct ccase **c,
620                         casenumber case_idx UNUSED)
621 {
622   struct dsc_trns *t = trns_;
623   struct dsc_z_score *z;
624   const struct variable **vars;
625
626   *c = case_unshare (*c);
627
628   if (t->filter)
629     {
630       double f = case_num (*c, t->filter);
631       if (f == 0.0 || var_is_num_missing (t->filter, f, MV_ANY))
632         {
633           descriptives_set_all_sysmis_zscores (t, *c);
634           return TRNS_CONTINUE;
635         }
636     }
637
638   if (t->count <= 0)
639     {
640       struct ccase *z_case;
641
642       z_case = casereader_read (t->z_reader);
643       if (z_case)
644         {
645           size_t z_idx = 0;
646
647           t->count = case_num_idx (z_case, z_idx++);
648           for (z = t->z_scores; z < t->z_scores + t->n_z_scores; z++)
649             {
650               z->mean = case_num_idx (z_case, z_idx++);
651               z->std_dev = case_num_idx (z_case, z_idx++);
652             }
653           case_unref (z_case);
654         }
655       else
656         {
657           if (t->ok)
658             {
659               msg (SE,  _("Internal error processing Z scores.  "
660                           "Please report this to %s."),
661                    PACKAGE_BUGREPORT);
662               t->ok = false;
663             }
664           descriptives_set_all_sysmis_zscores (t, *c);
665           return TRNS_CONTINUE;
666         }
667     }
668   t->count--;
669
670   if (t->missing_type == DSC_LISTWISE)
671     {
672       assert(t->vars);
673       for (vars = t->vars; vars < t->vars + t->n_vars; vars++)
674         {
675           double score = case_num (*c, *vars);
676           if (var_is_num_missing (*vars, score, t->exclude))
677             {
678               descriptives_set_all_sysmis_zscores (t, *c);
679               return TRNS_CONTINUE;
680             }
681         }
682     }
683
684   for (z = t->z_scores; z < t->z_scores + t->n_z_scores; z++)
685     {
686       double input = case_num (*c, z->src_var);
687       double *output = case_num_rw (*c, z->z_var);
688
689       if (z->mean == SYSMIS || z->std_dev == SYSMIS
690           || var_is_num_missing (z->src_var, input, t->exclude))
691         *output = SYSMIS;
692       else
693         *output = (input - z->mean) / z->std_dev;
694     }
695   return TRNS_CONTINUE;
696 }
697
698 /* Frees a descriptives_trns struct. */
699 static bool
700 descriptives_trns_free (void *trns_)
701 {
702   struct dsc_trns *t = trns_;
703   bool ok = t->ok && !casereader_error (t->z_reader);
704
705   free (t->z_scores);
706   casereader_destroy (t->z_reader);
707   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
708   free (t->vars);
709   free (t);
710
711   return ok;
712 }
713
714 /* Sets up a transformation to calculate Z scores. */
715 static void
716 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
717 {
718   struct dsc_trns *t;
719   size_t cnt, i;
720
721   for (cnt = i = 0; i < dsc->n_vars; i++)
722     if (dsc->vars[i].z_name != NULL)
723       cnt++;
724
725   t = xmalloc (sizeof *t);
726   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
727   t->n_z_scores = cnt;
728   t->missing_type = dsc->missing_type;
729   t->exclude = dsc->exclude;
730   if (t->missing_type == DSC_LISTWISE)
731     {
732       t->n_vars = dsc->n_vars;
733       t->vars = xnmalloc (t->n_vars, sizeof *t->vars);
734       for (i = 0; i < t->n_vars; i++)
735         t->vars[i] = dsc->vars[i].v;
736     }
737   else
738     {
739       t->n_vars = 0;
740       t->vars = NULL;
741     }
742   t->filter = dict_get_filter (dataset_dict (ds));
743   t->z_reader = casewriter_make_reader (dsc->z_writer);
744   t->count = 0;
745   t->ok = true;
746   dsc->z_writer = NULL;
747
748   for (cnt = i = 0; i < dsc->n_vars; i++)
749     {
750       struct dsc_var *dv = &dsc->vars[i];
751       if (dv->z_name != NULL)
752         {
753           struct dsc_z_score *z;
754           struct variable *dst_var;
755           char *label;
756
757           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
758
759           label = xasprintf (_("Z-score of %s"),var_to_string (dv->v));
760           var_set_label (dst_var, label);
761           free (label);
762
763           z = &t->z_scores[cnt++];
764           z->src_var = dv->v;
765           z->z_var = dst_var;
766         }
767     }
768
769   add_transformation (ds,
770                       descriptives_trns_proc, descriptives_trns_free, t);
771 }
772 \f
773 /* Statistical calculation. */
774
775 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
776
777 /* Calculates and displays descriptive statistics for the cases
778    in CF. */
779 static void
780 calc_descriptives (struct dsc_proc *dsc, struct casereader *group,
781                    struct dataset *ds)
782 {
783   const struct variable *filter = dict_get_filter (dataset_dict (ds));
784   struct casereader *pass1, *pass2;
785   casenumber count;
786   struct ccase *c;
787   size_t z_idx;
788   size_t i;
789
790   c = casereader_peek (group, 0);
791   if (c == NULL)
792     {
793       casereader_destroy (group);
794       return;
795     }
796   output_split_file_values (ds, c);
797   case_unref (c);
798
799   group = casereader_create_filter_weight (group, dataset_dict (ds),
800                                            NULL, NULL);
801
802   pass1 = group;
803   pass2 = dsc->max_moment <= MOMENT_MEAN ? NULL : casereader_clone (pass1);
804
805   for (i = 0; i < dsc->n_vars; i++)
806     {
807       struct dsc_var *dv = &dsc->vars[i];
808
809       dv->valid = dv->missing = 0.0;
810       if (dv->moments != NULL)
811         moments_clear (dv->moments);
812       dv->min = DBL_MAX;
813       dv->max = -DBL_MAX;
814     }
815   dsc->missing_listwise = 0.;
816   dsc->valid = 0.;
817
818   /* First pass to handle most of the work. */
819   count = 0;
820   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
821     {
822       double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
823
824       if (filter)
825         {
826           double f = case_num (c, filter);
827           if (f == 0.0 || var_is_num_missing (filter, f, MV_ANY))
828             continue;
829         }
830
831       /* Check for missing values. */
832       if (listwise_missing (dsc, c))
833         {
834           dsc->missing_listwise += weight;
835           if (dsc->missing_type == DSC_LISTWISE)
836             continue;
837         }
838       dsc->valid += weight;
839
840       for (i = 0; i < dsc->n_vars; i++)
841         {
842           struct dsc_var *dv = &dsc->vars[i];
843           double x = case_num (c, dv->v);
844
845           if (var_is_num_missing (dv->v, x, dsc->exclude))
846             {
847               dv->missing += weight;
848               continue;
849             }
850
851           if (dv->moments != NULL)
852             moments_pass_one (dv->moments, x, weight);
853
854           if (x < dv->min)
855             dv->min = x;
856           if (x > dv->max)
857             dv->max = x;
858         }
859
860       count++;
861     }
862   if (!casereader_destroy (pass1))
863     {
864       casereader_destroy (pass2);
865       return;
866     }
867
868   /* Second pass for higher-order moments. */
869   if (dsc->max_moment > MOMENT_MEAN)
870     {
871       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
872         {
873           double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
874
875           if (filter)
876             {
877               double f = case_num (c, filter);
878               if (f == 0.0 || var_is_num_missing (filter, f, MV_ANY))
879                 continue;
880             }
881
882           /* Check for missing values. */
883           if (dsc->missing_type == DSC_LISTWISE && listwise_missing (dsc, c))
884             continue;
885
886           for (i = 0; i < dsc->n_vars; i++)
887             {
888               struct dsc_var *dv = &dsc->vars[i];
889               double x = case_num (c, dv->v);
890
891               if (var_is_num_missing (dv->v, x, dsc->exclude))
892                 continue;
893
894               if (dv->moments != NULL)
895                 moments_pass_two (dv->moments, x, weight);
896             }
897         }
898       if (!casereader_destroy (pass2))
899         return;
900     }
901
902   /* Calculate results. */
903   if (dsc->z_writer && count > 0)
904     {
905       c = case_create (casewriter_get_proto (dsc->z_writer));
906       z_idx = 0;
907       *case_num_rw_idx (c, z_idx++) = count;
908     }
909   else
910     c = NULL;
911
912   for (i = 0; i < dsc->n_vars; i++)
913     {
914       struct dsc_var *dv = &dsc->vars[i];
915       double W;
916       int j;
917
918       for (j = 0; j < DSC_N_STATS; j++)
919         dv->stats[j] = SYSMIS;
920
921       dv->valid = W = dsc->valid - dv->missing;
922
923       if (dv->moments != NULL)
924         moments_calculate (dv->moments, NULL,
925                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
926                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
927       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
928           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
929         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
930       if (dsc->calc_stats & (1ul << DSC_STDDEV)
931           && dv->stats[DSC_VARIANCE] != SYSMIS)
932         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
933       if (dsc->calc_stats & (1ul << DSC_SEKURT))
934         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
935             dv->stats[DSC_SEKURT] = calc_sekurt (W);
936       if (dsc->calc_stats & (1ul << DSC_SESKEW)
937           && dv->stats[DSC_SKEWNESS] != SYSMIS)
938         dv->stats[DSC_SESKEW] = calc_seskew (W);
939       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
940                               ? SYSMIS : dv->max - dv->min);
941       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
942       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
943       if (dsc->calc_stats & (1ul << DSC_SUM))
944         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
945
946       if (dv->z_name && c != NULL)
947         {
948           *case_num_rw_idx (c, z_idx++) = dv->stats[DSC_MEAN];
949           *case_num_rw_idx (c, z_idx++) = dv->stats[DSC_STDDEV];
950         }
951     }
952
953   if (c != NULL)
954     casewriter_write (dsc->z_writer, c);
955
956   /* Output results. */
957   display (dsc);
958 }
959
960 /* Returns true if any of the descriptives variables in DSC's
961    variable list have missing values in case C, false otherwise. */
962 static bool
963 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
964 {
965   size_t i;
966
967   for (i = 0; i < dsc->n_vars; i++)
968     {
969       struct dsc_var *dv = &dsc->vars[i];
970       double x = case_num (c, dv->v);
971
972       if (var_is_num_missing (dv->v, x, dsc->exclude))
973         return true;
974     }
975   return false;
976 }
977 \f
978 /* Statistical display. */
979
980 static algo_compare_func descriptives_compare_dsc_vars;
981
982 /* Displays a table of descriptive statistics for DSC. */
983 static void
984 display (struct dsc_proc *dsc)
985 {
986   struct pivot_table *table = pivot_table_create (
987     N_("Descriptive Statistics"));
988   pivot_table_set_weight_var (table, dict_get_weight (dsc->dict));
989
990   struct pivot_dimension *statistics = pivot_dimension_create (
991     table, PIVOT_AXIS_COLUMN, N_("Statistics"));
992   pivot_category_create_leaf_rc (
993     statistics->root, pivot_value_new_text (N_("N")), PIVOT_RC_COUNT);
994   for (int i = 0; i < DSC_N_STATS; i++)
995     if (dsc->show_stats & (1ul << i))
996       pivot_category_create_leaf (statistics->root,
997                                   pivot_value_new_text (dsc_info[i].name));
998
999   if (dsc->sort_by_stat != DSC_NONE)
1000     sort (dsc->vars, dsc->n_vars, sizeof *dsc->vars,
1001           descriptives_compare_dsc_vars, dsc);
1002
1003   struct pivot_dimension *variables = pivot_dimension_create (
1004     table, PIVOT_AXIS_ROW, N_("Variable"));
1005   for (size_t i = 0; i < dsc->n_vars; i++)
1006     {
1007       const struct dsc_var *dv = &dsc->vars[i];
1008
1009       int row = pivot_category_create_leaf (variables->root,
1010                                             pivot_value_new_variable (dv->v));
1011
1012       int column = 0;
1013       pivot_table_put2 (table, column++, row,
1014                         pivot_value_new_number (dv->valid));
1015
1016       for (int j = 0; j < DSC_N_STATS; j++)
1017         if (dsc->show_stats & (1ul << j))
1018           {
1019             union value v = { .f = dv->stats[j] };
1020             struct pivot_value *pv = (j == DSC_MIN || j == DSC_MAX
1021                                       ? pivot_value_new_var_value (dv->v, &v)
1022                                       : pivot_value_new_number (dv->stats[j]));
1023             pivot_table_put2 (table, column++, row, pv);
1024           }
1025     }
1026
1027   int row = pivot_category_create_leaves (
1028     variables->root, N_("Valid N (listwise)"), N_("Missing N (listwise)"));
1029   pivot_table_put2 (table, 0, row, pivot_value_new_number (dsc->valid));
1030   pivot_table_put2 (table, 0, row + 1,
1031                     pivot_value_new_number (dsc->missing_listwise));
1032   pivot_table_submit (table);
1033 }
1034
1035 /* Compares `struct dsc_var's A and B according to the ordering
1036    specified by CMD. */
1037 static int
1038 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
1039 {
1040   const struct dsc_var *a = a_;
1041   const struct dsc_var *b = b_;
1042   const struct dsc_proc *dsc = dsc_;
1043
1044   int result;
1045
1046   if (dsc->sort_by_stat == DSC_NAME)
1047     result = utf8_strcasecmp (var_get_name (a->v), var_get_name (b->v));
1048   else
1049     {
1050       double as = a->stats[dsc->sort_by_stat];
1051       double bs = b->stats[dsc->sort_by_stat];
1052
1053       result = as < bs ? -1 : as > bs;
1054     }
1055
1056   if (!dsc->sort_ascending)
1057     result = -result;
1058
1059   return result;
1060 }