FREQUENCIES: Added new test for bug which causes a crash
[pspp] / src / language / stats / frequencies.q
index aab345b67ec7fa33dc351020fae9c406f97a3436..8a8ac6c4bd64b4633099c5db1b5f53e69774cb38 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 1997-9, 2000, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Copyright (C) 1997-9, 2000, 2007, 2009, 2010, 2011 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -23,9 +23,9 @@
 #include "data/case.h"
 #include "data/casegrouper.h"
 #include "data/casereader.h"
+#include "data/dataset.h"
 #include "data/dictionary.h"
 #include "data/format.h"
-#include "data/procedure.h"
 #include "data/settings.h"
 #include "data/value-labels.h"
 #include "data/variable.h"
@@ -43,6 +43,8 @@
 #include "libpspp/str.h"
 #include "math/histogram.h"
 #include "math/moments.h"
+#include "math/chart-geometry.h"
+
 #include "output/chart-item.h"
 #include "output/charts/piechart.h"
 #include "output/charts/plot-hist.h"
@@ -125,10 +127,6 @@ struct percentile
 {
   double p;        /* the %ile to be calculated */
   double value;    /* the %ile's value */
-  double x1;       /* The datum value <= the percentile */
-  double x2;       /* The datum value >= the percentile */
-  int flag;
-  int flag2;       /* Set to 1 if this percentile value has been found */
   bool show;       /* True to show this percentile in the statistics box. */
 };
 
@@ -192,7 +190,6 @@ struct var_freqs
 
     /* Variable attributes. */
     int width;
-    struct fmt_spec print;
   };
 
 struct frq_proc
@@ -222,8 +219,9 @@ struct frq_proc
 static void determine_charts (struct frq_proc *,
                               const struct cmd_frequencies *);
 
-static void calc_stats (const struct frq_proc *, const struct var_freqs *,
-                        double d[FRQ_N_STATS]);
+static void calc_stats (const struct var_freqs *, double d[FRQ_N_STATS]);
+static void calc_percentiles (const struct frq_proc *,
+                              const struct var_freqs *);
 
 static void precalc (struct frq_proc *, struct casereader *, struct dataset *);
 static void calc (struct frq_proc *, const struct ccase *,
@@ -399,10 +397,13 @@ determine_charts (struct frq_proc *frq, const struct cmd_frequencies *cmd)
       if (hist->x_min != SYSMIS && hist->x_max != SYSMIS
           && hist->x_min >= hist->x_max)
         {
-          msg (SE, _("MAX for histogram must be greater than or equal to MIN, "
-                     "but MIN was specified as %.15g and MAX as %.15g.  "
-                     "MIN and MAX will be ignored."),
-               hist->x_min, hist->x_max);
+          msg (SE, _("%s for histogram must be greater than or equal to %s, "
+                     "but %s was specified as %.15g and %s as %.15g.  "
+                     "%s and %s will be ignored."),
+              "MAX", "MIN", 
+              "MIN", hist->x_min, 
+              "MAX", hist->x_max,
+              "MIN", "MAX");
           hist->x_min = hist->x_max = SYSMIS;
         }
     }
@@ -420,9 +421,13 @@ determine_charts (struct frq_proc *frq, const struct cmd_frequencies *cmd)
       if (pie->x_min != SYSMIS && pie->x_max != SYSMIS
           && pie->x_min >= pie->x_max)
         {
-          msg (SE, _("MAX for pie chart must be greater than or equal to MIN, "
-                     "but MIN was specified as %.15g and MAX as %.15g.  "
-                     "MIN and MAX will be ignored."), pie->x_min, pie->x_max);
+          msg (SE, _("%s for pie chart must be greater than or equal to %s, "
+                     "but %s was specified as %.15g and %s as %.15g.  "
+                     "%s and %s will be ignored."), 
+              "MAX", "MIN", 
+              "MIN", pie->x_min,
+              "MAX", pie->x_max,
+              "MIN", "MAX");
           pie->x_min = pie->x_max = SYSMIS;
         }
     }
@@ -497,18 +502,21 @@ postcalc (struct frq_proc *frq, const struct dataset *ds)
          double d[FRQ_N_STATS];
          struct histogram *histogram;
 
-         calc_stats (frq, vf, d);
+         calc_stats (vf, d);
 
          histogram = freq_tab_to_hist (frq, &vf->tab, vf->var);
 
-          chart_item_submit (histogram_chart_create (
+         if ( histogram)
+           {
+             chart_item_submit (histogram_chart_create (
                                histogram->gsl_hist, var_to_string(vf->var),
                                vf->tab.valid_cases,
                                d[FRQ_MEAN],
                                d[FRQ_STDDEV],
                                frq->hist->draw_normal));
 
-         statistic_destroy (&histogram->parent);
+             statistic_destroy (&histogram->parent);
+           }
        }
 
       if (frq->pie)
@@ -621,10 +629,10 @@ frq_custom_variables (struct lexer *lexer, struct dataset *ds,
   size_t n_vars;
   size_t i;
 
-  lex_match (lexer, '=');
+  lex_match (lexer, T_EQUALS);
   if (lex_token (lexer) != T_ALL
       && (lex_token (lexer) != T_ID
-          || dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) == NULL))
+          || dict_lookup_var (dataset_dict (ds), lex_tokcstr (lexer)) == NULL))
     return 2;
 
   /* Get list of current variables, to avoid duplicates. */
@@ -649,7 +657,6 @@ frq_custom_variables (struct lexer *lexer, struct dataset *ds,
       vf->n_groups = 0;
       vf->groups = NULL;
       vf->width = var_get_width (var);
-      vf->print = *var_get_print_format (var);
     }
   frq->n_vars = n_vars;
 
@@ -665,8 +672,9 @@ frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequenc
 {
   struct frq_proc *frq = frq_;
 
-  lex_match (lexer, '=');
-  if ((lex_token (lexer) == T_ID && dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) != NULL)
+  lex_match (lexer, T_EQUALS);
+  if ((lex_token (lexer) == T_ID
+       && dict_lookup_var (dataset_dict (ds), lex_tokcstr (lexer)) != NULL)
       || lex_token (lexer) == T_ID)
     for (;;)
       {
@@ -683,7 +691,7 @@ frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequenc
        if (!parse_variables_const (lexer, dataset_dict (ds), &v, &n,
                               PV_NO_DUPLICATE | PV_NUMERIC))
          return 0;
-       if (lex_match (lexer, '('))
+       if (lex_match (lexer, T_LPAREN))
          {
            nl = ml = 0;
            dl = NULL;
@@ -696,14 +704,14 @@ frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequenc
                  }
                dl[nl++] = lex_tokval (lexer);
                lex_get (lexer);
-               lex_match (lexer, ',');
+               lex_match (lexer, T_COMMA);
              }
            /* Note that nl might still be 0 and dl might still be
               NULL.  That's okay. */
-           if (!lex_match (lexer, ')'))
+           if (!lex_match (lexer, T_RPAREN))
              {
                free (v);
-               msg (SE, _("`)' expected after GROUPED interval list."));
+                lex_error_expecting (lexer, "`)'", NULL_SENTINEL);
                return 0;
              }
          }
@@ -724,7 +732,7 @@ frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequenc
                   {
                     if (vf->groups != NULL)
                       msg (SE, _("Variables %s specified multiple times on "
-                                 "GROUPED subcommand."), var_get_name (v[i]));
+                                 "%s subcommand."), var_get_name (v[i]), "GROUPED");
                     else
                       {
                         vf->n_groups = nl;
@@ -733,21 +741,30 @@ frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequenc
                     goto found;
                   }
               }
-            msg (SE, _("Variables %s specified on GROUPED but not on "
-                       "VARIABLES."), var_get_name (v[i]));
+            msg (SE, _("Variables %s specified on %s but not on "
+                       "%s."), var_get_name (v[i]), "GROUPED", "VARIABLES");
 
           found:;
           }
 
        free (v);
-       if (!lex_match (lexer, '/'))
-         break;
-       if ((lex_token (lexer) != T_ID || dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) != NULL)
-            && lex_token (lexer) != T_ALL)
-         {
-           lex_put_back (lexer, '/');
-           break;
-         }
+        if (lex_token (lexer) != T_SLASH)
+          break;
+
+        if ((lex_next_token (lexer, 1) == T_ID
+             && dict_lookup_var (dataset_dict (ds),
+                                 lex_next_tokcstr (lexer, 1)))
+            || lex_next_token (lexer, 1) == T_ALL)
+          {
+            /* The token after the slash is a variable name.  Keep parsing. */
+            lex_get (lexer);
+          }
+        else
+          {
+            /* The token after the slash must be the start of a new
+               subcommand.  Let the caller see the slash. */
+            break;
+          }
       }
 
   return 1;
@@ -843,7 +860,7 @@ dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
       if (label != NULL)
         tab_text (t, 0, r, TAB_LEFT, label);
 
-      tab_value (t, 1, r, TAB_NONE, &f->value, ft->dict, &vf->print);
+      tab_value (t, 1, r, TAB_NONE, &f->value, vf->var, NULL);
       tab_double (t, 2, r, TAB_NONE, f->count, wfmt);
       tab_double (t, 3, r, TAB_NONE, percent, NULL);
       tab_double (t, 4, r, TAB_NONE, valid_percent, NULL);
@@ -860,7 +877,7 @@ dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
       if (label != NULL)
         tab_text (t, 0, r, TAB_LEFT, label);
 
-      tab_value (t, 1, r, TAB_NONE, &f->value, ft->dict, &vf->print);
+      tab_value (t, 1, r, TAB_NONE, &f->value, vf->var, NULL);
       tab_double (t, 2, r, TAB_NONE, f->count, wfmt);
       tab_double (t, 3, r, TAB_NONE,
                     f->count / ft->total_cases * 100.0, NULL);
@@ -883,114 +900,73 @@ dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
 \f
 /* Statistical display. */
 
-/* Calculates all the pertinent statistics for variable V, putting them in
-   array D[]. */
+static double
+calc_percentile (double p, double valid_cases, double x1, double x2)
+{
+  double s, dummy;
+
+  s = (settings_get_algorithm () != COMPATIBLE
+       ? modf ((valid_cases - 1) * p, &dummy)
+       : modf ((valid_cases + 1) * p - 1, &dummy));
+
+  return x1 + (x2 - x1) * s;
+}
+
+/* Calculates all of the percentiles for VF within FRQ. */
 static void
-calc_stats (const struct frq_proc *frq,
-            const struct var_freqs *vf, double d[FRQ_N_STATS])
+calc_percentiles (const struct frq_proc *frq, const struct var_freqs *vf)
 {
   const struct freq_tab *ft = &vf->tab;
   double W = ft->valid_cases;
-  struct moments *m;
-  struct freq *f=0;
-  double prev_value;
-  int most_often;
-  double X_mode;
-
+  const struct freq *f;
+  int percentile_idx;
   double rank;
-  int i = 0;
-  int idx;
-
-  /* Calculate percentiles. */
 
   assert (ft->n_valid > 0);
 
-  for (i = 0; i < frq->n_percentiles; i++)
-    {
-      struct percentile *pc = &frq->percentiles[i];
-
-      pc->flag = 0;
-      pc->flag2 = 0;
-    }
-
   rank = 0;
-  prev_value = SYSMIS;
-  for (idx = 0; idx < ft->n_valid; ++idx)
+  percentile_idx = 0;
+  for (f = ft->valid; f < ft->missing; f++)
     {
-      f = &ft->valid[idx];
-      rank += f->count ;
-      for (i = 0; i < frq->n_percentiles; i++)
+      rank += f->count;
+      for (; percentile_idx < frq->n_percentiles; percentile_idx++)
         {
-          struct percentile *pc = &frq->percentiles[i];
-         double tp;
+          struct percentile *pc = &frq->percentiles[percentile_idx];
+          double tp;
 
-         if ( pc->flag2  ) continue ;
+          tp = (settings_get_algorithm () == ENHANCED
+                ? (W - 1) * pc->p
+                : (W + 1) * pc->p - 1);
 
-         if ( settings_get_algorithm () != COMPATIBLE )
-           tp =
-             (ft->valid_cases - 1) *  pc->p;
-         else
-           tp =
-             (ft->valid_cases + 1) *  pc->p - 1;
+          if (rank <= tp)
+            break;
 
-         if ( pc->flag )
-           {
-             pc->x2 = f->value.f;
-             pc->x1 = prev_value;
-             pc->flag2 = 1;
-             continue;
-           }
-
-          if (rank >  tp )
-         {
-           if ( f->count > 1 && rank - (f->count - 1) > tp )
-             {
-               pc->x2 = pc->x1 = f->value.f;
-               pc->flag2 = 1;
-             }
-           else
-             {
-               pc->flag=1;
-             }
-
-           continue;
-         }
+          if (tp + 1 < rank || f + 1 >= ft->missing)
+            pc->value = f->value.f;
+          else
+            pc->value = calc_percentile (pc->p, W, f->value.f, f[1].value.f);
         }
-      prev_value = f->value.f;
     }
-
-  for (i = 0; i < frq->n_percentiles; i++)
+  for (; percentile_idx < frq->n_percentiles; percentile_idx++)
     {
-      struct percentile *pc = &frq->percentiles[i];
-
-      /* Catches the case when p == 100% */
-      if ( ! pc->flag2 )
-       pc->x1 = pc->x2 = f->value.f;
-
-      /*
-      printf("percentile %d (p==%.2f); X1 = %g; X2 = %g\n",
-            i,pc->p,pc->x1,pc->x2);
-      */
+      struct percentile *pc = &frq->percentiles[percentile_idx];
+      pc->value = ft->valid[ft->n_valid - 1].value.f;
     }
+}
 
-  for (i = 0; i < frq->n_percentiles; i++)
-    {
-      struct percentile *pc = &frq->percentiles[i];
-      double s;
-
-      double dummy;
-      if ( settings_get_algorithm () != COMPATIBLE )
-       {
-         s = modf((ft->valid_cases - 1) * pc->p , &dummy);
-       }
-      else
-       {
-         s = modf((ft->valid_cases + 1) * pc->p -1, &dummy);
-       }
-
-      pc->value = pc->x1 + (pc->x2 - pc->x1) * s ;
-    }
+/* Calculates all the pertinent statistics for VF, putting them in array
+   D[]. */
+static void
+calc_stats (const struct var_freqs *vf, double d[FRQ_N_STATS])
+{
+  const struct freq_tab *ft = &vf->tab;
+  double W = ft->valid_cases;
+  const struct freq *f;
+  struct moments *m;
+  int most_often;
+  double X_mode;
 
+  assert (ft->n_valid > 0);
 
   /* Calculate the mode. */
   most_often = -1;
@@ -1052,7 +1028,8 @@ dump_statistics (const struct frq_proc *frq, const struct var_freqs *vf,
           var_get_name (vf->var));
       return;
     }
-  calc_stats (frq, vf, stat_value);
+  calc_stats (vf, stat_value);
+  calc_percentiles (frq, vf);
 
   t = tab_create (3, frq->n_stats + frq->n_show_percentiles + 2);
 
@@ -1080,7 +1057,7 @@ dump_statistics (const struct frq_proc *frq, const struct var_freqs *vf,
   tab_double (t, 2, 0, TAB_NONE, ft->valid_cases, wfmt);
   tab_double (t, 2, 1, TAB_NONE, ft->total_cases - ft->valid_cases, wfmt);
 
-  for (i = 0; i < frq->n_percentiles; i++, r++)
+  for (i = 0; i < frq->n_percentiles; i++)
     {
       struct percentile *pc = &frq->percentiles[i];
 
@@ -1098,6 +1075,7 @@ dump_statistics (const struct frq_proc *frq, const struct var_freqs *vf,
         tab_fixed (t, 1, r, TAB_LEFT, pc->p * 100, 3, 0);
       tab_double (t, 2, r, TAB_NONE, pc->value,
                   var_get_print_format (vf->var));
+      r++;
     }
 
   tab_title (t, "%s", var_to_string (vf->var));
@@ -1148,10 +1126,9 @@ freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
 {
   double x_min, x_max, valid_freq;
   int i;
-
+  double bin_width;
   struct histogram *histogram;
   double iqr;
-  int bins;
 
   /* Find out the extremes of the x value, within the range to be included in
      the histogram, and sum the total frequency of those values. */
@@ -1171,19 +1148,13 @@ freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
 
   /* Freedman-Diaconis' choice of bin width. */
   iqr = calculate_iqr (frq);
-  if (iqr != SYSMIS)
-    {
-      double bin_width = 2 * iqr / pow (valid_freq, 1.0 / 3.0);
-      bins = (x_max - x_min) / bin_width;
-      if (bins < 5)
-        bins = 5;
-      else if (bins > 400)
-        bins = 400;
-    }
-  else
-    bins = 5;
+  bin_width = 2 * iqr / pow (valid_freq, 1.0 / 3.0);
+
+  histogram = histogram_create (bin_width, x_min, x_max);
+
+  if ( histogram == NULL)
+    return NULL;
 
-  histogram = histogram_create (bins, x_min, x_max);
   for (i = 0; i < ft->n_valid; i++)
     {
       const struct freq *f = &ft->valid[i];