Merge remote-tracking branch 'origin/master' into sheet
[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/tab.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 z_score_cnt;            /* Number of Z-scores. */
75     const struct variable **vars;     /* Variables for listwise missing checks. */
76     size_t var_cnt;             /* 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     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 /* Output format. */
140 enum dsc_format
141   {
142     DSC_LINE,           /* Abbreviated format. */
143     DSC_SERIAL          /* Long format. */
144   };
145
146 /* A DESCRIPTIVES procedure. */
147 struct dsc_proc
148   {
149     /* Per-variable info. */
150     struct dsc_var *vars;       /* Variables. */
151     size_t var_cnt;             /* Number of variables. */
152
153     /* User options. */
154     enum dsc_missing_type missing_type; /* Treatment of missing values. */
155     enum mv_class exclude;      /* Classes of missing values to exclude. */
156     int show_var_labels;        /* Nonzero to show variable labels. */
157     int show_index;             /* Nonzero to show variable index. */
158     enum dsc_format format;     /* Output format. */
159
160     /* Accumulated results. */
161     double missing_listwise;    /* Sum of weights of cases missing listwise. */
162     double valid;               /* Sum of weights of valid cases. */
163     bool bad_warn;               /* Warn if bad weight found. */
164     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
165     int sort_ascending;         /* !0: ascending order; 0: descending. */
166     unsigned long show_stats;   /* Statistics to display. */
167     unsigned long calc_stats;   /* Statistics to calculate. */
168     enum moment max_moment;     /* Highest moment needed for stats. */
169
170     /* Z scores. */
171     struct casewriter *z_writer; /* Mean and stddev per SPLIT FILE group. */
172   };
173
174 /* Parsing. */
175 static enum dsc_statistic match_statistic (struct lexer *);
176 static void free_dsc_proc (struct dsc_proc *);
177
178 /* Z-score functions. */
179 static bool try_name (const struct dictionary *dict,
180                       struct dsc_proc *dsc, const char *name);
181 static char *generate_z_varname (const struct dictionary *dict,
182                                  struct dsc_proc *dsc,
183                                  const char *name, int *z_cnt);
184 static void dump_z_table (struct dsc_proc *);
185 static void setup_z_trns (struct dsc_proc *, struct dataset *);
186
187 /* Procedure execution functions. */
188 static void calc_descriptives (struct dsc_proc *, struct casereader *,
189                                struct dataset *);
190 static void display (struct dsc_proc *dsc);
191 \f
192 /* Parser and outline. */
193
194 /* Handles DESCRIPTIVES. */
195 int
196 cmd_descriptives (struct lexer *lexer, struct dataset *ds)
197 {
198   struct dictionary *dict = dataset_dict (ds);
199   struct dsc_proc *dsc;
200   const struct variable **vars = NULL;
201   size_t var_cnt = 0;
202   int save_z_scores = 0;
203   int z_cnt = 0;
204   size_t i;
205   bool ok;
206
207   struct casegrouper *grouper;
208   struct casereader *group;
209
210   /* Create and initialize dsc. */
211   dsc = xmalloc (sizeof *dsc);
212   dsc->vars = NULL;
213   dsc->var_cnt = 0;
214   dsc->missing_type = DSC_VARIABLE;
215   dsc->exclude = MV_ANY;
216   dsc->show_var_labels = 1;
217   dsc->show_index = 0;
218   dsc->format = DSC_LINE;
219   dsc->missing_listwise = 0.;
220   dsc->valid = 0.;
221   dsc->bad_warn = 1;
222   dsc->sort_by_stat = DSC_NONE;
223   dsc->sort_ascending = 1;
224   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
225   dsc->z_writer = NULL;
226
227   /* Parse DESCRIPTIVES. */
228   while (lex_token (lexer) != T_ENDCMD)
229     {
230       if (lex_match_id (lexer, "MISSING"))
231         {
232           lex_match (lexer, T_EQUALS);
233           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
234             {
235               if (lex_match_id (lexer, "VARIABLE"))
236                 dsc->missing_type = DSC_VARIABLE;
237               else if (lex_match_id (lexer, "LISTWISE"))
238                 dsc->missing_type = DSC_LISTWISE;
239               else if (lex_match_id (lexer, "INCLUDE"))
240                 dsc->exclude = MV_SYSTEM;
241               else
242                 {
243                   lex_error (lexer, NULL);
244                   goto error;
245                 }
246               lex_match (lexer, T_COMMA);
247             }
248         }
249       else if (lex_match_id (lexer, "SAVE"))
250         save_z_scores = 1;
251       else if (lex_match_id (lexer, "FORMAT"))
252         {
253           lex_match (lexer, T_EQUALS);
254           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
255             {
256               if (lex_match_id (lexer, "LABELS"))
257                 dsc->show_var_labels = 1;
258               else if (lex_match_id (lexer, "NOLABELS"))
259                 dsc->show_var_labels = 0;
260               else if (lex_match_id (lexer, "INDEX"))
261                 dsc->show_index = 1;
262               else if (lex_match_id (lexer, "NOINDEX"))
263                 dsc->show_index = 0;
264               else if (lex_match_id (lexer, "LINE"))
265                 dsc->format = DSC_LINE;
266               else if (lex_match_id (lexer, "SERIAL"))
267                 dsc->format = DSC_SERIAL;
268               else
269                 {
270                   lex_error (lexer, NULL);
271                   goto error;
272                 }
273               lex_match (lexer, T_COMMA);
274             }
275         }
276       else if (lex_match_id (lexer, "STATISTICS"))
277         {
278           lex_match (lexer, T_EQUALS);
279           dsc->show_stats = 0;
280           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
281             {
282               if (lex_match (lexer, T_ALL))
283                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
284               else if (lex_match_id (lexer, "DEFAULT"))
285                 dsc->show_stats |= DEFAULT_STATS;
286               else
287                 dsc->show_stats |= 1ul << (match_statistic (lexer));
288               lex_match (lexer, T_COMMA);
289             }
290           if (dsc->show_stats == 0)
291             dsc->show_stats = DEFAULT_STATS;
292         }
293       else if (lex_match_id (lexer, "SORT"))
294         {
295           lex_match (lexer, T_EQUALS);
296           if (lex_match_id (lexer, "NAME"))
297             dsc->sort_by_stat = DSC_NAME;
298           else
299             {
300               dsc->sort_by_stat = match_statistic (lexer);
301               if (dsc->sort_by_stat == DSC_NONE )
302                 dsc->sort_by_stat = DSC_MEAN;
303             }
304           if (lex_match (lexer, T_LPAREN))
305             {
306               if (lex_match_id (lexer, "A"))
307                 dsc->sort_ascending = 1;
308               else if (lex_match_id (lexer, "D"))
309                 dsc->sort_ascending = 0;
310               else
311                 lex_error (lexer, NULL);
312               if (! lex_force_match (lexer, T_RPAREN))
313                 goto error;
314             }
315         }
316       else if (var_cnt == 0)
317         {
318           if (lex_next_token (lexer, 1) == T_EQUALS)
319             {
320               lex_match_id (lexer, "VARIABLES");
321               lex_match (lexer, T_EQUALS);
322             }
323
324           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
325             {
326               int i;
327
328               if (!parse_variables_const (lexer, dict, &vars, &var_cnt,
329                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
330                 goto error;
331
332               dsc->vars = xnrealloc ((void *)dsc->vars, var_cnt, sizeof *dsc->vars);
333               for (i = dsc->var_cnt; i < var_cnt; i++)
334                 {
335                   struct dsc_var *dv = &dsc->vars[i];
336                   dv->v = vars[i];
337                   dv->z_name = NULL;
338                   dv->moments = NULL;
339                 }
340               dsc->var_cnt = var_cnt;
341
342               if (lex_match (lexer, T_LPAREN))
343                 {
344                   if (lex_token (lexer) != T_ID)
345                     {
346                       lex_error (lexer, NULL);
347                       goto error;
348                     }
349                   if (try_name (dict, dsc, lex_tokcstr (lexer)))
350                     {
351                       struct dsc_var *dsc_var = &dsc->vars[dsc->var_cnt - 1];
352                       dsc_var->z_name = xstrdup (lex_tokcstr (lexer));
353                       z_cnt++;
354                     }
355                   else
356                     msg (SE, _("Z-score variable name %s would be"
357                                " a duplicate variable name."), lex_tokcstr (lexer));
358                   lex_get (lexer);
359                   if (!lex_force_match (lexer, T_RPAREN))
360                     goto error;
361                 }
362             }
363         }
364       else
365         {
366           lex_error (lexer, NULL);
367           goto error;
368         }
369
370       lex_match (lexer, T_SLASH);
371     }
372   if (var_cnt == 0)
373     {
374       msg (SE, _("No variables specified."));
375       goto error;
376     }
377
378   /* Construct z-score varnames, show translation table. */
379   if (z_cnt || save_z_scores)
380     {
381       struct caseproto *proto;
382
383       if (save_z_scores)
384         {
385           int gen_cnt = 0;
386
387           for (i = 0; i < dsc->var_cnt; i++)
388             {
389               struct dsc_var *dsc_var = &dsc->vars[i];
390               if (dsc_var->z_name == NULL)
391                 {
392                   const char *name = var_get_name (dsc_var->v);
393                   dsc_var->z_name = generate_z_varname (dict, dsc, name,
394                                                         &gen_cnt);
395                   if (dsc_var->z_name == NULL)
396                     goto error;
397
398                   z_cnt++;
399                 }
400             }
401         }
402
403       /* It would be better to handle Z scores correctly (however we define
404          that) when TEMPORARY is in effect, but in the meantime this at least
405          prevents a use-after-free error.  See bug #38786.  */
406       if (proc_make_temporary_transformations_permanent (ds))
407         msg (SW, _("DESCRIPTIVES with Z scores ignores TEMPORARY.  "
408                    "Temporary transformations will be made permanent."));
409
410       proto = caseproto_create ();
411       for (i = 0; i < 1 + 2 * z_cnt; i++)
412         proto = caseproto_add_width (proto, 0);
413       dsc->z_writer = autopaging_writer_create (proto);
414       caseproto_unref (proto);
415
416       dump_z_table (dsc);
417     }
418
419   /* Figure out statistics to display. */
420   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
421     dsc->show_stats |= 1ul << DSC_SESKEW;
422   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
423     dsc->show_stats |= 1ul << DSC_SEKURT;
424
425   /* Figure out which statistics to calculate. */
426   dsc->calc_stats = dsc->show_stats;
427   if (z_cnt > 0)
428     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
429   if (dsc->sort_by_stat >= 0)
430     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
431   if (dsc->show_stats & (1ul << DSC_SESKEW))
432     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
433   if (dsc->show_stats & (1ul << DSC_SEKURT))
434     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
435
436   /* Figure out maximum moment needed and allocate moments for
437      the variables. */
438   dsc->max_moment = MOMENT_NONE;
439   for (i = 0; i < DSC_N_STATS; i++)
440     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
441       dsc->max_moment = dsc_info[i].moment;
442   if (dsc->max_moment != MOMENT_NONE)
443     for (i = 0; i < dsc->var_cnt; i++)
444       dsc->vars[i].moments = moments_create (dsc->max_moment);
445
446   /* Data pass. */
447   grouper = casegrouper_create_splits (proc_open_filtering (ds, false), dict);
448   while (casegrouper_get_next_group (grouper, &group))
449     calc_descriptives (dsc, group, ds);
450   ok = casegrouper_destroy (grouper);
451   ok = proc_commit (ds) && ok;
452
453   /* Z-scoring! */
454   if (ok && z_cnt)
455     setup_z_trns (dsc, ds);
456
457   /* Done. */
458   free (vars);
459   free_dsc_proc (dsc);
460   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
461
462  error:
463   free (vars);
464   free_dsc_proc (dsc);
465   return CMD_FAILURE;
466 }
467
468 /* Returns the statistic named by the current token and skips past the token.
469    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
470    specifiers). Emits an error if the current token ID does not name a
471    statistic. */
472 static enum dsc_statistic
473 match_statistic (struct lexer *lexer)
474 {
475   if (lex_token (lexer) == T_ID)
476     {
477       enum dsc_statistic stat;
478
479       for (stat = 0; stat < DSC_N_STATS; stat++)
480         if (lex_match_id (lexer, dsc_info[stat].identifier))
481           return stat;
482
483       lex_get (lexer);
484       lex_error (lexer, _("expecting statistic name: reverting to default"));
485     }
486
487   return DSC_NONE;
488 }
489
490 /* Frees DSC. */
491 static void
492 free_dsc_proc (struct dsc_proc *dsc)
493 {
494   size_t i;
495
496   if (dsc == NULL)
497     return;
498
499   for (i = 0; i < dsc->var_cnt; i++)
500     {
501       struct dsc_var *dsc_var = &dsc->vars[i];
502       free (dsc_var->z_name);
503       moments_destroy (dsc_var->moments);
504     }
505   casewriter_destroy (dsc->z_writer);
506   free (dsc->vars);
507   free (dsc);
508 }
509 \f
510 /* Z scores. */
511
512 /* Returns false if NAME is a duplicate of any existing variable name or
513    of any previously-declared z-var name; otherwise returns true. */
514 static bool
515 try_name (const struct dictionary *dict, struct dsc_proc *dsc,
516           const char *name)
517 {
518   size_t i;
519
520   if (dict_lookup_var (dict, name) != NULL)
521     return false;
522   for (i = 0; i < dsc->var_cnt; i++)
523     {
524       struct dsc_var *dsc_var = &dsc->vars[i];
525       if (dsc_var->z_name != NULL && !utf8_strcasecmp (dsc_var->z_name, name))
526         return false;
527     }
528   return true;
529 }
530
531 /* Generates a name for a Z-score variable based on a variable
532    named VAR_NAME, given that *Z_CNT generated variable names are
533    known to already exist.  If successful, returns the new name
534    as a dynamically allocated string.  On failure, returns NULL. */
535 static char *
536 generate_z_varname (const struct dictionary *dict, struct dsc_proc *dsc,
537                     const char *var_name, int *z_cnt)
538 {
539   char *z_name, *trunc_name;
540
541   /* Try a name based on the original variable name. */
542   z_name = xasprintf ("Z%s", var_name);
543   trunc_name = utf8_encoding_trunc (z_name, dict_get_encoding (dict),
544                                     ID_MAX_LEN);
545   free (z_name);
546   if (try_name (dict, dsc, trunc_name))
547     return trunc_name;
548   free (trunc_name);
549
550   /* Generate a synthetic name. */
551   for (;;)
552     {
553       char name[8];
554
555       (*z_cnt)++;
556
557       if (*z_cnt <= 99)
558         sprintf (name, "ZSC%03d", *z_cnt);
559       else if (*z_cnt <= 108)
560         sprintf (name, "STDZ%02d", *z_cnt - 99);
561       else if (*z_cnt <= 117)
562         sprintf (name, "ZZZZ%02d", *z_cnt - 108);
563       else if (*z_cnt <= 126)
564         sprintf (name, "ZQZQ%02d", *z_cnt - 117);
565       else
566         {
567           msg (SE, _("Ran out of generic names for Z-score variables.  "
568                      "There are only 126 generic names: ZSC001-ZSC0999, "
569                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
570           return NULL;
571         }
572
573       if (try_name (dict, dsc, name))
574         return xstrdup (name);
575     }
576   NOT_REACHED();
577 }
578
579 /* Outputs a table describing the mapping between source
580    variables and Z-score variables. */
581 static void
582 dump_z_table (struct dsc_proc *dsc)
583 {
584   size_t cnt = 0;
585   struct tab_table *t;
586
587   {
588     size_t i;
589
590     for (i = 0; i < dsc->var_cnt; i++)
591       if (dsc->vars[i].z_name != NULL)
592         cnt++;
593   }
594
595   t = tab_create (2, cnt + 1);
596   tab_title (t, _("Mapping of variables to corresponding Z-scores."));
597   tab_headers (t, 0, 0, 1, 0);
598   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, cnt);
599   tab_hline (t, TAL_2, 0, 1, 1);
600   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
601   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
602
603   {
604     size_t i, y;
605
606     for (i = 0, y = 1; i < dsc->var_cnt; i++)
607       if (dsc->vars[i].z_name != NULL)
608         {
609           tab_text (t, 0, y, TAB_LEFT, var_to_string (dsc->vars[i].v));
610           tab_text (t, 1, y++, TAB_LEFT, dsc->vars[i].z_name);
611         }
612   }
613
614   tab_submit (t);
615 }
616
617 static void
618 descriptives_set_all_sysmis_zscores (const struct dsc_trns *t, struct ccase *c)
619 {
620   const struct dsc_z_score *z;
621
622   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
623     case_data_rw (c, z->z_var)->f = SYSMIS;
624 }
625
626 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
627    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
628    SYSMIS 3) score is user missing and they were not included in the original
629    analyis. 4) any of the variables in the original analysis were missing
630    (either system or user-missing values that weren't included).
631 */
632 static int
633 descriptives_trns_proc (void *trns_, struct ccase **c,
634                         casenumber case_idx UNUSED)
635 {
636   struct dsc_trns *t = trns_;
637   struct dsc_z_score *z;
638   const struct variable **vars;
639
640   *c = case_unshare (*c);
641
642   if (t->filter)
643     {
644       double f = case_num (*c, t->filter);
645       if (f == 0.0 || var_is_num_missing (t->filter, f, MV_ANY))
646         {
647           descriptives_set_all_sysmis_zscores (t, *c);
648           return TRNS_CONTINUE;
649         }
650     }
651
652   if (t->count <= 0)
653     {
654       struct ccase *z_case;
655
656       z_case = casereader_read (t->z_reader);
657       if (z_case)
658         {
659           size_t z_idx = 0;
660
661           t->count = case_num_idx (z_case, z_idx++);
662           for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
663             {
664               z->mean = case_num_idx (z_case, z_idx++);
665               z->std_dev = case_num_idx (z_case, z_idx++);
666             }
667           case_unref (z_case);
668         }
669       else
670         {
671           if (t->ok)
672             {
673               msg (SE, _("Internal error processing Z scores"));
674               t->ok = false;
675             }
676           descriptives_set_all_sysmis_zscores (t, *c);
677           return TRNS_CONTINUE;
678         }
679     }
680   t->count--;
681
682   if (t->missing_type == DSC_LISTWISE)
683     {
684       assert(t->vars);
685       for (vars = t->vars; vars < t->vars + t->var_cnt; vars++)
686         {
687           double score = case_num (*c, *vars);
688           if (var_is_num_missing (*vars, score, t->exclude))
689             {
690               descriptives_set_all_sysmis_zscores (t, *c);
691               return TRNS_CONTINUE;
692             }
693         }
694     }
695
696   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
697     {
698       double input = case_num (*c, z->src_var);
699       double *output = &case_data_rw (*c, z->z_var)->f;
700
701       if (z->mean == SYSMIS || z->std_dev == SYSMIS
702           || var_is_num_missing (z->src_var, input, t->exclude))
703         *output = SYSMIS;
704       else
705         *output = (input - z->mean) / z->std_dev;
706     }
707   return TRNS_CONTINUE;
708 }
709
710 /* Frees a descriptives_trns struct. */
711 static bool
712 descriptives_trns_free (void *trns_)
713 {
714   struct dsc_trns *t = trns_;
715   bool ok = t->ok && !casereader_error (t->z_reader);
716
717   free (t->z_scores);
718   casereader_destroy (t->z_reader);
719   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
720   free (t->vars);
721   free (t);
722
723   return ok;
724 }
725
726 /* Sets up a transformation to calculate Z scores. */
727 static void
728 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
729 {
730   struct dsc_trns *t;
731   size_t cnt, i;
732
733   for (cnt = i = 0; i < dsc->var_cnt; i++)
734     if (dsc->vars[i].z_name != NULL)
735       cnt++;
736
737   t = xmalloc (sizeof *t);
738   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
739   t->z_score_cnt = cnt;
740   t->missing_type = dsc->missing_type;
741   t->exclude = dsc->exclude;
742   if ( t->missing_type == DSC_LISTWISE )
743     {
744       t->var_cnt = dsc->var_cnt;
745       t->vars = xnmalloc (t->var_cnt, sizeof *t->vars);
746       for (i = 0; i < t->var_cnt; i++)
747         t->vars[i] = dsc->vars[i].v;
748     }
749   else
750     {
751       t->var_cnt = 0;
752       t->vars = NULL;
753     }
754   t->filter = dict_get_filter (dataset_dict (ds));
755   t->z_reader = casewriter_make_reader (dsc->z_writer);
756   t->count = 0;
757   t->ok = true;
758   dsc->z_writer = NULL;
759
760   for (cnt = i = 0; i < dsc->var_cnt; i++)
761     {
762       struct dsc_var *dv = &dsc->vars[i];
763       if (dv->z_name != NULL)
764         {
765           struct dsc_z_score *z;
766           struct variable *dst_var;
767           char *label;
768
769           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
770
771           label = xasprintf (_("Z-score of %s"),var_to_string (dv->v));
772           var_set_label (dst_var, label);
773           free (label);
774
775           z = &t->z_scores[cnt++];
776           z->src_var = dv->v;
777           z->z_var = dst_var;
778         }
779     }
780
781   add_transformation (ds,
782                       descriptives_trns_proc, descriptives_trns_free, t);
783 }
784 \f
785 /* Statistical calculation. */
786
787 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
788
789 /* Calculates and displays descriptive statistics for the cases
790    in CF. */
791 static void
792 calc_descriptives (struct dsc_proc *dsc, struct casereader *group,
793                    struct dataset *ds)
794 {
795   struct variable *filter = dict_get_filter (dataset_dict (ds));
796   struct casereader *pass1, *pass2;
797   casenumber count;
798   struct ccase *c;
799   size_t z_idx;
800   size_t i;
801
802   c = casereader_peek (group, 0);
803   if (c == NULL)
804     {
805       casereader_destroy (group);
806       return;
807     }
808   output_split_file_values (ds, c);
809   case_unref (c);
810
811   group = casereader_create_filter_weight (group, dataset_dict (ds),
812                                            NULL, NULL);
813
814   pass1 = group;
815   pass2 = dsc->max_moment <= MOMENT_MEAN ? NULL : casereader_clone (pass1);
816
817   for (i = 0; i < dsc->var_cnt; i++)
818     {
819       struct dsc_var *dv = &dsc->vars[i];
820
821       dv->valid = dv->missing = 0.0;
822       if (dv->moments != NULL)
823         moments_clear (dv->moments);
824       dv->min = DBL_MAX;
825       dv->max = -DBL_MAX;
826     }
827   dsc->missing_listwise = 0.;
828   dsc->valid = 0.;
829
830   /* First pass to handle most of the work. */
831   count = 0;
832   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
833     {
834       double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
835
836       if (filter)
837         {
838           double f = case_num (c, filter);
839           if (f == 0.0 || var_is_num_missing (filter, f, MV_ANY))
840             continue;
841         }
842
843       /* Check for missing values. */
844       if (listwise_missing (dsc, c))
845         {
846           dsc->missing_listwise += weight;
847           if (dsc->missing_type == DSC_LISTWISE)
848             continue;
849         }
850       dsc->valid += weight;
851
852       for (i = 0; i < dsc->var_cnt; i++)
853         {
854           struct dsc_var *dv = &dsc->vars[i];
855           double x = case_num (c, dv->v);
856
857           if (var_is_num_missing (dv->v, x, dsc->exclude))
858             {
859               dv->missing += weight;
860               continue;
861             }
862
863           if (dv->moments != NULL)
864             moments_pass_one (dv->moments, x, weight);
865
866           if (x < dv->min)
867             dv->min = x;
868           if (x > dv->max)
869             dv->max = x;
870         }
871
872       count++;
873     }
874   if (!casereader_destroy (pass1))
875     {
876       casereader_destroy (pass2);
877       return;
878     }
879
880   /* Second pass for higher-order moments. */
881   if (dsc->max_moment > MOMENT_MEAN)
882     {
883       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
884         {
885           double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
886
887           if (filter)
888             {
889               double f = case_num (c, filter);
890               if (f == 0.0 || var_is_num_missing (filter, f, MV_ANY))
891                 continue;
892             }
893
894           /* Check for missing values. */
895           if (dsc->missing_type == DSC_LISTWISE && listwise_missing (dsc, c))
896             continue;
897
898           for (i = 0; i < dsc->var_cnt; i++)
899             {
900               struct dsc_var *dv = &dsc->vars[i];
901               double x = case_num (c, dv->v);
902
903               if (var_is_num_missing (dv->v, x, dsc->exclude))
904                 continue;
905
906               if (dv->moments != NULL)
907                 moments_pass_two (dv->moments, x, weight);
908             }
909         }
910       if (!casereader_destroy (pass2))
911         return;
912     }
913
914   /* Calculate results. */
915   if (dsc->z_writer && count > 0)
916     {
917       c = case_create (casewriter_get_proto (dsc->z_writer));
918       z_idx = 0;
919       case_data_rw_idx (c, z_idx++)->f = count;
920     }
921   else
922     c = NULL;
923
924   for (i = 0; i < dsc->var_cnt; i++)
925     {
926       struct dsc_var *dv = &dsc->vars[i];
927       double W;
928       int j;
929
930       for (j = 0; j < DSC_N_STATS; j++)
931         dv->stats[j] = SYSMIS;
932
933       dv->valid = W = dsc->valid - dv->missing;
934
935       if (dv->moments != NULL)
936         moments_calculate (dv->moments, NULL,
937                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
938                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
939       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
940           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
941         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
942       if (dsc->calc_stats & (1ul << DSC_STDDEV)
943           && dv->stats[DSC_VARIANCE] != SYSMIS)
944         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
945       if (dsc->calc_stats & (1ul << DSC_SEKURT))
946         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
947             dv->stats[DSC_SEKURT] = calc_sekurt (W);
948       if (dsc->calc_stats & (1ul << DSC_SESKEW)
949           && dv->stats[DSC_SKEWNESS] != SYSMIS)
950         dv->stats[DSC_SESKEW] = calc_seskew (W);
951       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
952                               ? SYSMIS : dv->max - dv->min);
953       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
954       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
955       if (dsc->calc_stats & (1ul << DSC_SUM))
956         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
957
958       if (dv->z_name && c != NULL)
959         {
960           case_data_rw_idx (c, z_idx++)->f = dv->stats[DSC_MEAN];
961           case_data_rw_idx (c, z_idx++)->f = dv->stats[DSC_STDDEV];
962         }
963     }
964
965   if (c != NULL)
966     casewriter_write (dsc->z_writer, c);
967
968   /* Output results. */
969   display (dsc);
970 }
971
972 /* Returns true if any of the descriptives variables in DSC's
973    variable list have missing values in case C, false otherwise. */
974 static bool
975 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
976 {
977   size_t i;
978
979   for (i = 0; i < dsc->var_cnt; i++)
980     {
981       struct dsc_var *dv = &dsc->vars[i];
982       double x = case_num (c, dv->v);
983
984       if (var_is_num_missing (dv->v, x, dsc->exclude))
985         return true;
986     }
987   return false;
988 }
989 \f
990 /* Statistical display. */
991
992 static algo_compare_func descriptives_compare_dsc_vars;
993
994 /* Displays a table of descriptive statistics for DSC. */
995 static void
996 display (struct dsc_proc *dsc)
997 {
998   size_t i;
999   int nc;
1000   struct tab_table *t;
1001
1002   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
1003   for (i = 0; i < DSC_N_STATS; i++)
1004     if (dsc->show_stats & (1ul << i))
1005       nc++;
1006
1007   if (dsc->sort_by_stat != DSC_NONE)
1008     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
1009           descriptives_compare_dsc_vars, dsc);
1010
1011   t = tab_create (nc, dsc->var_cnt + 1);
1012   tab_headers (t, 1, 0, 1, 0);
1013   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
1014   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
1015   tab_hline (t, TAL_2, 0, nc - 1, 1);
1016   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
1017
1018   nc = 0;
1019   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
1020   if (dsc->format == DSC_SERIAL)
1021     {
1022       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
1023       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
1024     }
1025   else
1026     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
1027
1028   for (i = 0; i < DSC_N_STATS; i++)
1029     if (dsc->show_stats & (1ul << i))
1030       {
1031         const char *title = gettext (dsc_info[i].name);
1032         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
1033       }
1034
1035   for (i = 0; i < dsc->var_cnt; i++)
1036     {
1037       struct dsc_var *dv = &dsc->vars[i];
1038       size_t j;
1039
1040       nc = 0;
1041       tab_text (t, nc++, i + 1, TAB_LEFT, var_to_string (dv->v));
1042       tab_text_format (t, nc++, i + 1, 0, "%.*g", DBL_DIG + 1, dv->valid);
1043       if (dsc->format == DSC_SERIAL)
1044         tab_text_format (t, nc++, i + 1, 0, "%.*g", DBL_DIG + 1, dv->missing);
1045
1046       for (j = 0; j < DSC_N_STATS; j++)
1047         if (dsc->show_stats & (1ul << j))
1048           tab_double (t, nc++, i + 1, TAB_NONE, dv->stats[j], NULL, RC_OTHER);
1049     }
1050
1051   tab_title (t, _("Valid cases = %.*g; cases with missing value(s) = %.*g."),
1052              DBL_DIG + 1, dsc->valid,
1053              DBL_DIG + 1, dsc->missing_listwise);
1054
1055   tab_submit (t);
1056 }
1057
1058 /* Compares `struct dsc_var's A and B according to the ordering
1059    specified by CMD. */
1060 static int
1061 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
1062 {
1063   const struct dsc_var *a = a_;
1064   const struct dsc_var *b = b_;
1065   const struct dsc_proc *dsc = dsc_;
1066
1067   int result;
1068
1069   if (dsc->sort_by_stat == DSC_NAME)
1070     result = utf8_strcasecmp (var_get_name (a->v), var_get_name (b->v));
1071   else
1072     {
1073       double as = a->stats[dsc->sort_by_stat];
1074       double bs = b->stats[dsc->sort_by_stat];
1075
1076       result = as < bs ? -1 : as > bs;
1077     }
1078
1079   if (!dsc->sort_ascending)
1080     result = -result;
1081
1082   return result;
1083 }