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