Fix assertion for proper Huffman merge pattern: 0 == 1 modulo 1.
[pspp] / src / descript.c
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3    Written by Ben Pfaff <blp@gnu.org>.
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
18    02111-1307, USA. */
19
20 /* FIXME: Many possible optimizations. */
21
22 #include <config.h>
23 #include "error.h"
24 #include <limits.h>
25 #include <math.h>
26 #include <stdlib.h>
27 #include "algorithm.h"
28 #include "alloc.h"
29 #include "case.h"
30 #include "casefile.h"
31 #include "command.h"
32 #include "dictionary.h"
33 #include "lexer.h"
34 #include "error.h"
35 #include "magic.h"
36 #include "moments.h"
37 #include "som.h"
38 #include "tab.h"
39 #include "var.h"
40 #include "vfm.h"
41
42 /* DESCRIPTIVES private data. */
43
44 struct dsc_proc;
45
46 /* Handling of missing values. */
47 enum dsc_missing_type
48   {
49     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
50     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
51   };
52
53 /* Describes properties of a distribution for the purpose of
54    calculating a Z-score. */
55 struct dsc_z_score
56   {
57     int src_idx;                /* Source index into case data. */
58     int dst_idx;                /* Destination index into case data. */
59     double mean;                /* Distribution mean. */
60     double std_dev;             /* Distribution standard deviation. */
61     struct variable *v;         /* Variable on which z-score is based. */
62   };
63
64 /* DESCRIPTIVES transformation (for calculating Z-scores). */
65 struct dsc_trns
66   {
67     struct trns_header h;
68     struct dsc_z_score *z_scores; /* Array of Z-scores. */
69     int z_score_cnt;            /* Number of Z-scores. */
70     struct variable **vars;     /* Variables for listwise missing checks. */
71     size_t var_cnt;             /* Number of variables. */
72     enum dsc_missing_type missing_type; /* Treatment of missing values. */
73     int include_user_missing;   /* Nonzero to include user-missing values. */
74   };
75
76 /* Statistics.  Used as bit indexes, so must be 32 or fewer. */
77 enum dsc_statistic
78   {
79     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
80     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
81     DSC_MAX, DSC_SUM, DSC_N_STATS,
82
83     /* Only valid as sort criteria. */
84     DSC_NAME = -2,              /* Sort by name. */
85     DSC_NONE = -1               /* Unsorted. */
86   };
87
88 /* Describes one statistic. */
89 struct dsc_statistic_info
90   {
91     const char *identifier;     /* Identifier. */
92     const char *name;           /* Full name. */
93     enum moment moment;         /* Highest moment needed to calculate. */
94   };
95
96 /* Table of statistics, indexed by DSC_*. */
97 static const struct dsc_statistic_info dsc_info[DSC_N_STATS] =
98   {
99     {"MEAN", N_("Mean"), MOMENT_MEAN},
100     {"SEMEAN", N_("S E Mean"), MOMENT_VARIANCE},
101     {"STDDEV", N_("Std Dev"), MOMENT_VARIANCE},
102     {"VARIANCE", N_("Variance"), MOMENT_VARIANCE},
103     {"KURTOSIS", N_("Kurtosis"), MOMENT_KURTOSIS},
104     {"SEKURTOSIS", N_("S E Kurt"), MOMENT_NONE},
105     {"SKEWNESS", N_("Skewness"), MOMENT_SKEWNESS},
106     {"SESKEWNESS", N_("S E Skew"), MOMENT_NONE},
107     {"RANGE", N_("Range"), MOMENT_NONE},
108     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
109     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
110     {"SUM", N_("Sum"), MOMENT_MEAN},
111   };
112
113 /* Statistics calculated by default if none are explicitly
114    requested. */
115 #define DEFAULT_STATS                                                   \
116         ((1ul << DSC_MEAN) | (1ul << DSC_STDDEV) | (1ul << DSC_MIN)     \
117          | (1ul << DSC_MAX))
118      
119 /* A variable specified on DESCRIPTIVES. */
120 struct dsc_var
121   {
122     struct variable *v;         /* Variable to calculate on. */
123     char z_name[9];             /* Name for z-score variable. */
124     double valid, missing;      /* Valid, missing counts. */
125     struct moments *moments;    /* Moments. */
126     double min, max;            /* Maximum and mimimum values. */
127     double stats[DSC_N_STATS];  /* All the stats' values. */
128   };
129
130 /* Output format. */
131 enum dsc_format 
132   {
133     DSC_LINE,           /* Abbreviated format. */
134     DSC_SERIAL          /* Long format. */
135   };
136
137 /* A DESCRIPTIVES procedure. */
138 struct dsc_proc 
139   {
140     /* Per-variable info. */
141     struct dsc_var *vars;       /* Variables. */
142     size_t var_cnt;             /* Number of variables. */
143
144     /* User options. */
145     enum dsc_missing_type missing_type; /* Treatment of missing values. */
146     int include_user_missing;   /* Nonzero to include user-missing values. */
147     int show_var_labels;        /* Nonzero to show variable labels. */
148     int show_index;             /* Nonzero to show variable index. */
149     enum dsc_format format;     /* Output format. */
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     int 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
162 /* Parsing. */
163 static enum dsc_statistic match_statistic (void);
164 static void free_dsc_proc (struct dsc_proc *);
165
166 /* Z-score functions. */
167 static int try_name (struct dsc_proc *dsc, char *name);
168 static int generate_z_varname (struct dsc_proc *dsc, char *z_name,
169                                const char *name, int *z_cnt);
170 static void dump_z_table (struct dsc_proc *);
171 static void setup_z_trns (struct dsc_proc *);
172
173 /* Procedure execution functions. */
174 static void calc_descriptives (const struct casefile *, void *dsc_);
175 static void display (struct dsc_proc *dsc);
176 \f
177 /* Parser and outline. */
178
179 /* Handles DESCRIPTIVES. */
180 int
181 cmd_descriptives (void)
182 {
183   struct dsc_proc *dsc;
184   struct variable **vars = NULL;
185   int var_cnt = 0;
186   int save_z_scores = 0;
187   int z_cnt = 0;
188   int i;
189
190   /* Create and initialize dsc. */
191   dsc = xmalloc (sizeof *dsc);
192   dsc->vars = NULL;
193   dsc->var_cnt = 0;
194   dsc->missing_type = DSC_VARIABLE;
195   dsc->include_user_missing = 0;
196   dsc->show_var_labels = 1;
197   dsc->show_index = 0;
198   dsc->format = DSC_LINE;
199   dsc->missing_listwise = 0.;
200   dsc->valid = 0.;
201   dsc->bad_warn = 1;
202   dsc->sort_by_stat = DSC_NONE;
203   dsc->sort_ascending = 1;
204   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
205
206   /* Parse DESCRIPTIVES. */
207   while (token != '.') 
208     {
209       if (lex_match_id ("MISSING"))
210         {
211           lex_match ('=');
212           while (token != '.' && token != '/') 
213             {
214               if (lex_match_id ("VARIABLE"))
215                 dsc->missing_type = DSC_VARIABLE;
216               else if (lex_match_id ("LISTWISE"))
217                 dsc->missing_type = DSC_LISTWISE;
218               else if (lex_match_id ("INCLUDE"))
219                 dsc->include_user_missing = 1;
220               else
221                 {
222                   lex_error (NULL);
223                   goto error;
224                 }
225               lex_match (',');
226             }
227         }
228       else if (lex_match_id ("SAVE"))
229         save_z_scores = 1;
230       else if (lex_match_id ("FORMAT")) 
231         {
232           lex_match ('=');
233           while (token != '.' && token != '/') 
234             {
235               if (lex_match_id ("LABELS"))
236                 dsc->show_var_labels = 1;
237               else if (lex_match_id ("NOLABELS"))
238                 dsc->show_var_labels = 0;
239               else if (lex_match_id ("INDEX"))
240                 dsc->show_index = 1;
241               else if (lex_match_id ("NOINDEX"))
242                 dsc->show_index = 0;
243               else if (lex_match_id ("LINE"))
244                 dsc->format = DSC_LINE;
245               else if (lex_match_id ("SERIAL"))
246                 dsc->format = DSC_SERIAL;
247               else
248                 {
249                   lex_error (NULL);
250                   goto error;
251                 }
252               lex_match (',');
253             }
254         }
255       else if (lex_match_id ("STATISTICS")) 
256         {
257           lex_match ('=');
258           dsc->show_stats = 0;
259           while (token != '.' && token != '/') 
260             {
261               if (lex_match (T_ALL)) 
262                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
263               else if (lex_match_id ("DEFAULT"))
264                 dsc->show_stats |= DEFAULT_STATS;
265               else
266                 dsc->show_stats |= 1ul << (match_statistic ());
267               lex_match (',');
268             }
269           if (dsc->show_stats == 0)
270             dsc->show_stats = DEFAULT_STATS;
271         }
272       else if (lex_match_id ("SORT")) 
273         {
274           lex_match ('=');
275           if (lex_match_id ("NAME"))
276             dsc->sort_by_stat = DSC_NAME;
277           else 
278             {
279               dsc->sort_by_stat = match_statistic ();
280               if (dsc->sort_by_stat == DSC_NONE )
281                 dsc->sort_by_stat = DSC_MEAN;
282             }
283           if (lex_match ('(')) 
284             {
285               if (lex_match_id ("A"))
286                 dsc->sort_ascending = 1;
287               else if (lex_match_id ("D"))
288                 dsc->sort_ascending = 0;
289               else
290                 lex_error (NULL);
291               lex_force_match (')');
292             }
293         }
294       else if (var_cnt == 0)
295         {
296           if (lex_look_ahead () == '=') 
297             {
298               lex_match_id ("VARIABLES");
299               lex_match ('=');
300             }
301
302           while (token != '.' && token != '/') 
303             {
304               int i;
305               
306               if (!parse_variables (default_dict, &vars, &var_cnt,
307                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
308                 goto error;
309
310               dsc->vars = xrealloc (dsc->vars, sizeof *dsc->vars * var_cnt);
311               for (i = dsc->var_cnt; i < var_cnt; i++)
312                 {
313                   struct dsc_var *dv = &dsc->vars[i];
314                   dv->v = vars[i];
315                   dv->z_name[0] = '\0';
316                   dv->moments = NULL;
317                 }
318               dsc->var_cnt = var_cnt;
319
320               if (lex_match ('(')) 
321                 {
322                   if (token != T_ID) 
323                     {
324                       lex_error (NULL);
325                       goto error;
326                     }
327                   if (try_name (dsc, tokid)) 
328                     {
329                       strcpy (dsc->vars[dsc->var_cnt - 1].z_name, tokid);
330                       z_cnt++;
331                     }
332                   else
333                     msg (SE, _("Z-score variable name %s would be"
334                                " a duplicate variable name."), tokid);
335                   lex_get ();
336                   if (!lex_force_match (')'))
337                     goto error;
338                 }
339             }
340         }
341       else 
342         {
343           lex_error (NULL);
344           goto error; 
345         }
346
347       lex_match ('/');
348     }
349   if (var_cnt == 0)
350     {
351       msg (SE, _("No variables specified."));
352       goto error;
353     }
354
355   /* Construct z-score varnames, show translation table. */
356   if (z_cnt || save_z_scores)
357     {
358       if (save_z_scores) 
359         {
360           int gen_cnt = 0;
361
362           for (i = 0; i < dsc->var_cnt; i++)
363             if (dsc->vars[i].z_name[0] == 0) 
364               {
365                 if (!generate_z_varname (dsc, dsc->vars[i].z_name,
366                                          dsc->vars[i].v->name, &gen_cnt))
367                   goto error;
368                 z_cnt++;
369               } 
370         }
371       dump_z_table (dsc);
372     }
373
374   /* Figure out statistics to display. */
375   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
376     dsc->show_stats |= 1ul << DSC_SESKEW;
377   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
378     dsc->show_stats |= 1ul << DSC_SEKURT;
379
380   /* Figure out which statistics to calculate. */
381   dsc->calc_stats = dsc->show_stats;
382   if (z_cnt > 0)
383     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
384   if (dsc->sort_by_stat >= 0)
385     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
386   if (dsc->show_stats & (1ul << DSC_SESKEW))
387     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
388   if (dsc->show_stats & (1ul << DSC_SEKURT))
389     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
390
391   /* Figure out maximum moment needed and allocate moments for
392      the variables. */
393   dsc->max_moment = MOMENT_NONE;
394   for (i = 0; i < DSC_N_STATS; i++) 
395     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
396       dsc->max_moment = dsc_info[i].moment;
397   if (dsc->max_moment != MOMENT_NONE)
398     for (i = 0; i < dsc->var_cnt; i++)
399       dsc->vars[i].moments = moments_create (dsc->max_moment);
400
401   /* Data pass. */
402   multipass_procedure_with_splits (calc_descriptives, dsc);
403
404   /* Z-scoring! */
405   if (z_cnt)
406     setup_z_trns (dsc);
407
408   /* Done. */
409   free (vars);
410   free_dsc_proc (dsc);
411   return CMD_SUCCESS;
412
413  error:
414   free (vars);
415   free_dsc_proc (dsc);
416   return CMD_FAILURE;
417 }
418
419 /* Returns the statistic named by the current token and skips past the token.
420    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
421    specifiers). Emits an error if the current token ID does not name a
422    statistic. */
423 static enum dsc_statistic
424 match_statistic (void) 
425 {
426   if (token == T_ID) 
427     {
428       enum dsc_statistic stat;
429
430       for (stat = 0; stat < DSC_N_STATS; stat++)
431         if (lex_match_id (dsc_info[stat].identifier)) 
432           return stat;
433
434       lex_get();
435       lex_error (_("expecting statistic name: reverting to default"));
436     }
437
438   return DSC_NONE;
439 }
440
441 /* Frees DSC. */
442 static void
443 free_dsc_proc (struct dsc_proc *dsc)
444 {
445   size_t i;
446
447   if (dsc == NULL)
448     return;
449   
450   for (i = 0; i < dsc->var_cnt; i++)
451     moments_destroy (dsc->vars[i].moments);
452   free (dsc->vars);
453   free (dsc);
454 }
455 \f
456 /* Z scores. */
457
458 /* Returns 0 if NAME is a duplicate of any existing variable name or
459    of any previously-declared z-var name; otherwise returns 1. */
460 static int
461 try_name (struct dsc_proc *dsc, char *name)
462 {
463   int i;
464
465   if (dict_lookup_var (default_dict, name) != NULL)
466     return 0;
467   for (i = 0; i < dsc->var_cnt; i++)
468     if (!strcmp (dsc->vars[i].z_name, name))
469       return 0;
470   return 1;
471 }
472
473 /* Generates a name for a Z-score variable based on a variable
474    named VAR_NAME, given that *Z_CNT generated variable names are
475    known to already exist.  If successful, returns nonzero and
476    copies the new name into Z_NAME.  On failure, returns zero. */
477 static int
478 generate_z_varname (struct dsc_proc *dsc, char *z_name,
479                     const char *var_name, int *z_cnt)
480 {
481   char name[10];
482
483   /* Try a name based on the original variable name. */
484   name[0] = 'Z';
485   strcpy (name + 1, var_name);
486   name[8] = '\0';
487   if (try_name (dsc, name))
488     {
489       strcpy (z_name, name);
490       return 1;
491     }
492
493   /* Generate a synthetic name. */
494   for (;;)
495     {
496       (*z_cnt)++;
497
498       if (*z_cnt <= 99)
499         sprintf (name, "ZSC%03d", *z_cnt);
500       else if (*z_cnt <= 108)
501         sprintf (name, "STDZ%02d", *z_cnt - 99);
502       else if (*z_cnt <= 117)
503         sprintf (name, "ZZZZ%02d", *z_cnt - 108);
504       else if (*z_cnt <= 126)
505         sprintf (name, "ZQZQ%02d", *z_cnt - 117);
506       else
507         {
508           msg (SE, _("Ran out of generic names for Z-score variables.  "
509                      "There are only 126 generic names: ZSC001-ZSC0999, "
510                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
511           return 0;
512         }
513       
514       if (try_name (dsc, name))
515         {
516           strcpy (z_name, name);
517           return 1;
518         }
519     }
520 }
521
522 /* Outputs a table describing the mapping between source
523    variables and Z-score variables. */
524 static void
525 dump_z_table (struct dsc_proc *dsc)
526 {
527   int cnt = 0;
528   struct tab_table *t;
529   
530   {
531     int i;
532     
533     for (i = 0; i < dsc->var_cnt; i++)
534       if (dsc->vars[i].z_name[0] != '\0')
535         cnt++;
536   }
537   
538   t = tab_create (2, cnt + 1, 0);
539   tab_title (t, 0, _("Mapping of variables to corresponding Z-scores."));
540   tab_columns (t, SOM_COL_DOWN, 1);
541   tab_headers (t, 0, 0, 1, 0);
542   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, cnt);
543   tab_hline (t, TAL_2, 0, 1, 1);
544   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
545   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
546   tab_dim (t, tab_natural_dimensions);
547
548   {
549     int i, y;
550     
551     for (i = 0, y = 1; i < dsc->var_cnt; i++)
552       if (dsc->vars[i].z_name[0] != '\0')
553         {
554           tab_text (t, 0, y, TAB_LEFT, dsc->vars[i].v->name);
555           tab_text (t, 1, y++, TAB_LEFT, dsc->vars[i].z_name);
556         }
557   }
558   
559   tab_submit (t);
560 }
561
562 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
563    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
564    SYSMIS 3) score is user missing and they were not included in the original
565    analyis. 4) any of the variables in the original analysis were missing
566    (either system or user-missing values that weren't included).
567 */
568 static int
569 descriptives_trns_proc (struct trns_header *trns, struct ccase * c,
570                         int case_idx UNUSED)
571 {
572   struct dsc_trns *t = (struct dsc_trns *) trns;
573   struct dsc_z_score *z;
574   struct variable **vars;
575   int all_sysmis = 0;
576
577   if (t->missing_type == DSC_LISTWISE)
578     {
579       assert(t->vars);
580       for (vars = t->vars; vars < t->vars + t->var_cnt; vars++)
581         {
582           double score = case_num (c, (*vars)->fv);
583           if ( score == SYSMIS || (!t->include_user_missing 
584                                    && is_num_user_missing(score, *vars)) )
585             {
586               all_sysmis = 1;
587               break;
588             }
589         }
590     }
591       
592   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
593     {
594       double input = case_num (c, z->src_idx);
595       double *output = &case_data_rw (c, z->dst_idx)->f;
596
597       if (z->mean == SYSMIS || z->std_dev == SYSMIS 
598           || all_sysmis || input == SYSMIS 
599           || (!t->include_user_missing && is_num_user_missing(input, z->v)))
600         *output = SYSMIS;
601       else
602         *output = (input - z->mean) / z->std_dev;
603     }
604   return -1;
605 }
606
607 /* Frees a descriptives_trns struct. */
608 static void
609 descriptives_trns_free (struct trns_header * trns)
610 {
611   struct dsc_trns *t = (struct dsc_trns *) trns;
612
613   free (t->z_scores);
614   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
615   free (t->vars);
616 }
617
618 /* Sets up a transformation to calculate Z scores. */
619 static void
620 setup_z_trns (struct dsc_proc *dsc)
621 {
622   struct dsc_trns *t;
623   int cnt, i;
624
625   for (cnt = i = 0; i < dsc->var_cnt; i++)
626     if (dsc->vars[i].z_name[0] != '\0')
627       cnt++;
628
629   t = xmalloc (sizeof *t);
630   t->h.proc = descriptives_trns_proc;
631   t->h.free = descriptives_trns_free;
632   t->z_scores = xmalloc (cnt * sizeof *t->z_scores);
633   t->z_score_cnt = cnt;
634   t->missing_type = dsc->missing_type;
635   t->include_user_missing = dsc->include_user_missing;
636   if ( t->missing_type == DSC_LISTWISE )
637     {
638       t->var_cnt = dsc->var_cnt;
639       t->vars = xmalloc(t->var_cnt * sizeof *t->vars);
640       for (i = 0; i < t->var_cnt; i++)
641         t->vars[i] = dsc->vars[i].v;
642     }
643   else
644     {
645       t->var_cnt = 0;
646       t->vars = NULL;
647     }
648   
649
650   for (cnt = i = 0; i < dsc->var_cnt; i++)
651     {
652       struct dsc_var *dv = &dsc->vars[i];
653       if (dv->z_name[0] != '\0')
654         {
655           struct dsc_z_score *z;
656           char *cp;
657           struct variable *dst_var;
658
659           dst_var = dict_create_var_assert (default_dict, dv->z_name, 0);
660           dst_var->init = 0;
661           if (dv->v->label)
662             {
663               dst_var->label = xmalloc (strlen (dv->v->label) + 12);
664               cp = stpcpy (dst_var->label, _("Z-score of "));
665               strcpy (cp, dv->v->label);
666             }
667           else
668             {
669               dst_var->label = xmalloc (strlen (dv->v->name) + 12);
670               cp = stpcpy (dst_var->label, _("Z-score of "));
671               strcpy (cp, dv->v->name);
672             }
673
674           z = &t->z_scores[cnt++];
675           z->src_idx = dv->v->fv;
676           z->dst_idx = dst_var->fv;
677           z->mean = dv->stats[DSC_MEAN];
678           z->std_dev = dv->stats[DSC_STDDEV];
679           z->v = dv->v;
680         }
681     }
682
683   add_transformation ((struct trns_header *) t);
684 }
685 \f
686 /* Statistical calculation. */
687
688 static int listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
689
690 /* Calculates and displays descriptive statistics for the cases
691    in CF. */
692 static void
693 calc_descriptives (const struct casefile *cf, void *dsc_) 
694 {
695   struct dsc_proc *dsc = dsc_;
696   struct casereader *reader;
697   struct ccase c;
698   int i;
699
700   for (i = 0; i < dsc->var_cnt; i++)
701     {
702       struct dsc_var *dv = &dsc->vars[i];
703       
704       dv->valid = dv->missing = 0.0;
705       if (dv->moments != NULL)
706         moments_clear (dv->moments);
707       dv->min = DBL_MAX;
708       dv->max = -DBL_MAX;
709     }
710   dsc->missing_listwise = 0.;
711   dsc->valid = 0.;
712
713   /* First pass to handle most of the work. */
714   for (reader = casefile_get_reader (cf);
715        casereader_read (reader, &c);
716        case_destroy (&c))
717     {
718       double weight = dict_get_case_weight (default_dict, &c, &dsc->bad_warn);
719       if (weight <= 0.0) 
720         continue;
721        
722       /* Check for missing values. */
723       if (listwise_missing (dsc, &c)) 
724         {
725           dsc->missing_listwise += weight;
726           if (dsc->missing_type == DSC_LISTWISE)
727             continue; 
728         }
729       dsc->valid += weight;
730
731       for (i = 0; i < dsc->var_cnt; i++) 
732         {
733           struct dsc_var *dv = &dsc->vars[i];
734           double x = case_num (&c, dv->v->fv);
735           
736           if (dsc->missing_type != DSC_LISTWISE
737               && (x == SYSMIS
738                   || (!dsc->include_user_missing
739                       && is_num_user_missing (x, dv->v))))
740             {
741               dv->missing += weight;
742               continue;
743             }
744
745           if (dv->moments != NULL) 
746             moments_pass_one (dv->moments, x, weight);
747
748           if (x < dv->min)
749             dv->min = x;
750           if (x > dv->max)
751             dv->max = x;
752         }
753     }
754   casereader_destroy (reader);
755
756   /* Second pass for higher-order moments. */
757   if (dsc->max_moment > MOMENT_MEAN) 
758     {
759       for (reader = casefile_get_reader (cf);
760            casereader_read (reader, &c);
761            case_destroy (&c))
762         {
763           double weight = dict_get_case_weight (default_dict, &c, 
764                                                 &dsc->bad_warn);
765           if (weight <= 0.0)
766             continue;
767       
768           /* Check for missing values. */
769           if (listwise_missing (dsc, &c) 
770               && dsc->missing_type == DSC_LISTWISE)
771             continue; 
772
773           for (i = 0; i < dsc->var_cnt; i++) 
774             {
775               struct dsc_var *dv = &dsc->vars[i];
776               double x = case_num (&c, dv->v->fv);
777           
778               if (dsc->missing_type != DSC_LISTWISE
779                   && (x == SYSMIS
780                       || (!dsc->include_user_missing
781                           && is_num_user_missing (x, dv->v))))
782                 continue;
783
784               if (dv->moments != NULL)
785                 moments_pass_two (dv->moments, x, weight);
786             }
787         }
788       casereader_destroy (reader);
789     }
790   
791   /* Calculate results. */
792   for (i = 0; i < dsc->var_cnt; i++)
793     {
794       struct dsc_var *dv = &dsc->vars[i];
795       double W;
796       int j;
797
798       for (j = 0; j < DSC_N_STATS; j++)
799         dv->stats[j] = SYSMIS;
800
801       dv->valid = W = dsc->valid - dv->missing;
802
803       if (dv->moments != NULL)
804         moments_calculate (dv->moments, NULL,
805                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
806                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
807       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
808           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
809         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
810       if (dsc->calc_stats & (1ul << DSC_STDDEV)
811           && dv->stats[DSC_VARIANCE] != SYSMIS)
812         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
813       if (dsc->calc_stats & (1ul << DSC_SEKURT)) 
814         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
815             dv->stats[DSC_SEKURT] = calc_sekurt (W);
816       if (dsc->calc_stats & (1ul << DSC_SESKEW)
817           && dv->stats[DSC_SKEWNESS] != SYSMIS)
818         dv->stats[DSC_SESKEW] = calc_seskew (W);
819       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
820                               ? SYSMIS : dv->max - dv->min);
821       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
822       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
823       if (dsc->calc_stats & (1ul << DSC_SUM))
824         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
825     }
826
827   /* Output results. */
828   display (dsc);
829 }
830
831 /* Returns nonzero if any of the descriptives variables in DSC's
832    variable list have missing values in case C, zero otherwise. */
833 static int
834 listwise_missing (struct dsc_proc *dsc, const struct ccase *c) 
835 {
836   int i;
837
838   for (i = 0; i < dsc->var_cnt; i++)
839     {
840       struct dsc_var *dv = &dsc->vars[i];
841       double x = case_num (c, dv->v->fv);
842
843       if (x == SYSMIS
844           || (!dsc->include_user_missing && is_num_user_missing (x, dv->v)))
845         return 1;
846     }
847   return 0;
848 }
849 \f
850 /* Statistical display. */
851
852 static algo_compare_func descriptives_compare_dsc_vars;
853
854 /* Displays a table of descriptive statistics for DSC. */
855 static void
856 display (struct dsc_proc *dsc)
857 {
858   int i, j;
859   int nc;
860   struct tab_table *t;
861
862   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
863   for (i = 0; i < DSC_N_STATS; i++)
864     if (dsc->show_stats & (1ul << i))
865       nc++;
866
867   if (dsc->sort_by_stat != DSC_NONE)
868     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
869           descriptives_compare_dsc_vars, dsc);
870
871   t = tab_create (nc, dsc->var_cnt + 1, 0);
872   tab_headers (t, 1, 0, 1, 0);
873   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
874   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
875   tab_hline (t, TAL_2, 0, nc - 1, 1);
876   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
877   tab_dim (t, tab_natural_dimensions);
878
879   nc = 0;
880   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
881   if (dsc->format == DSC_SERIAL)
882     {
883       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
884       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
885     }
886   else
887     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
888
889   for (i = 0; i < DSC_N_STATS; i++)
890     if (dsc->show_stats & (1ul << i))
891       {
892         const char *title = gettext (dsc_info[i].name);
893         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
894       }
895
896   for (i = 0; i < dsc->var_cnt; i++)
897     {
898       struct dsc_var *dv = &dsc->vars[i];
899
900       nc = 0;
901       tab_text (t, nc++, i + 1, TAB_LEFT, dv->v->name);
902       tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", dv->valid);
903       if (dsc->format == DSC_SERIAL)
904         tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", dv->missing);
905       for (j = 0; j < DSC_N_STATS; j++)
906         if (dsc->show_stats & (1ul << j))
907           tab_float (t, nc++, i + 1, TAB_NONE, dv->stats[j], 10, 3);
908     }
909
910   tab_title (t, 1, _("Valid cases = %g; cases with missing value(s) = %g."),
911              dsc->valid, dsc->missing_listwise);
912
913   tab_submit (t);
914 }
915
916 /* Compares `struct dsc_var's A and B according to the ordering
917    specified by CMD. */
918 static int
919 descriptives_compare_dsc_vars (const void *a_, const void *b_, void *dsc_)
920 {
921   const struct dsc_var *a = a_;
922   const struct dsc_var *b = b_;
923   struct dsc_proc *dsc = dsc_;
924
925   int result;
926
927   if (dsc->sort_by_stat == DSC_NAME)
928     result = strcmp (a->v->name, b->v->name);
929   else 
930     {
931       double as = a->stats[dsc->sort_by_stat];
932       double bs = b->stats[dsc->sort_by_stat];
933
934       result = as < bs ? -1 : as > bs;
935     }
936
937   if (!dsc->sort_ascending)
938     result = -result;
939
940   return result;
941 }