Merge "master" into "output".
[pspp-builds.git] / src / language / stats / examine.q
index 08077942340234031939cded549e87886d5119ee..1d399b5a1888908a4b08c8cfe2874da0af039bed 100644 (file)
@@ -49,7 +49,7 @@
 #include <libpspp/str.h>
 #include <math/moments.h>
 #include <output/charts/box-whisker.h>
-#include <output/charts/cartesian.h>
+#include <output/charts/np-plot.h>
 #include <output/manager.h>
 #include <output/table.h>
 
@@ -104,11 +104,11 @@ struct factor_metrics
   struct percentile **ptl;
   size_t n_ptiles;
 
-  struct statistic *tukey_hinges;
-  struct statistic *box_whisker;
-  struct statistic *trimmed_mean;
-  struct statistic *histogram;
-  struct order_stats *np;
+  struct tukey_hinges *tukey_hinges;
+  struct box_whisker *box_whisker;
+  struct trimmed_mean *trimmed_mean;
+  struct histogram *histogram;
+  struct np *np;
 
   /* Three quartiles indexing into PTL */
   struct percentile **quartiles;
@@ -179,12 +179,12 @@ factor_destroy (struct xfactor *fctr)
          moments1_destroy (result->metrics[v].moments);
          extrema_destroy (result->metrics[v].minima);
          extrema_destroy (result->metrics[v].maxima);
-         statistic_destroy (result->metrics[v].trimmed_mean);
-         statistic_destroy (result->metrics[v].tukey_hinges);
-         statistic_destroy (result->metrics[v].box_whisker);
-         statistic_destroy (result->metrics[v].histogram);
+         statistic_destroy (&result->metrics[v].trimmed_mean->parent.parent);
+         statistic_destroy (&result->metrics[v].tukey_hinges->parent.parent);
+         statistic_destroy (&result->metrics[v].box_whisker->parent.parent);
+         statistic_destroy (&result->metrics[v].histogram->parent);
          for (i = 0 ; i < result->metrics[v].n_ptiles; ++i)
-           statistic_destroy ((struct statistic *) result->metrics[v].ptl[i]);
+           statistic_destroy (&result->metrics[v].ptl[i]->parent.parent);
          free (result->metrics[v].ptl);
          free (result->metrics[v].quartiles);
          casereader_destroy (result->metrics[v].up_reader);
@@ -320,82 +320,6 @@ cmd_examine (struct lexer *lexer, struct dataset *ds)
 };
 
 
-/* Plot the normal and detrended normal plots for RESULT.
-   Label the plots with LABEL */
-static void
-np_plot (struct np *np, const char *label)
-{
-  double yfirst = 0, ylast = 0;
-
-  double x_lower;
-  double x_upper;
-  double slack;
-
-  /* Normal Plot */
-  struct chart *np_chart;
-
-  /* Detrended Normal Plot */
-  struct chart *dnp_chart;
-
-  /* The slope and intercept of the ideal normal probability line */
-  const double slope = 1.0 / np->stddev;
-  const double intercept = -np->mean / np->stddev;
-
-  if ( np->n < 1.0 )
-    {
-      msg (MW, _("Not creating plot because data set is empty."));
-      return ;
-    }
-
-  np_chart = chart_create ();
-  dnp_chart = chart_create ();
-
-  if ( !np_chart || ! dnp_chart )
-    return ;
-
-  chart_write_title (np_chart, _("Normal Q-Q Plot of %s"), label);
-  chart_write_xlabel (np_chart, _("Observed Value"));
-  chart_write_ylabel (np_chart, _("Expected Normal"));
-
-  chart_write_title (dnp_chart, _("Detrended Normal Q-Q Plot of %s"),
-                    label);
-  chart_write_xlabel (dnp_chart, _("Observed Value"));
-  chart_write_ylabel (dnp_chart, _("Dev from Normal"));
-
-  yfirst = gsl_cdf_ugaussian_Pinv (1 / (np->n + 1));
-  ylast = gsl_cdf_ugaussian_Pinv (np->n / (np->n + 1));
-
-  /* Need to make sure that both the scatter plot and the ideal fit into the
-     plot */
-  x_lower = MIN (np->y_min, (yfirst - intercept) / slope) ;
-  x_upper = MAX (np->y_max, (ylast  - intercept) / slope) ;
-  slack = (x_upper - x_lower)  * 0.05 ;
-
-  chart_write_xscale (np_chart, x_lower - slack, x_upper + slack, 5);
-  chart_write_xscale (dnp_chart, np->y_min, np->y_max, 5);
-
-  chart_write_yscale (np_chart, yfirst, ylast, 5);
-  chart_write_yscale (dnp_chart, np->dns_min, np->dns_max, 5);
-
-  {
-    struct casereader *reader = casewriter_make_reader (np->writer);
-    struct ccase *c;
-    while ((c = casereader_read (reader)) != NULL)
-      {
-       chart_datum (np_chart, 0, case_data_idx (c, NP_IDX_Y)->f, case_data_idx (c, NP_IDX_NS)->f);
-       chart_datum (dnp_chart, 0, case_data_idx (c, NP_IDX_Y)->f, case_data_idx (c, NP_IDX_DNS)->f);
-
-       case_unref (c);
-      }
-    casereader_destroy (reader);
-  }
-
-  chart_line (dnp_chart, 0, 0, np->y_min, np->y_max , CHART_DIM_X);
-  chart_line (np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
-
-  chart_submit (np_chart);
-  chart_submit (dnp_chart);
-}
 
 
 static void
@@ -412,20 +336,37 @@ show_npplot (const struct variable **dependent_var,
           ll != ll_null (&fctr->result_list);
           ll = ll_next (ll))
        {
-         struct string str;
+         struct string label;
          const struct factor_result *result =
            ll_data (ll, struct factor_result, ll);
-
-         ds_init_empty (&str);
-         ds_put_format (&str, "%s ", var_get_name (dependent_var[v]));
-
-         factor_to_string (fctr, result, &str);
-
-         np_plot ((struct np*) result->metrics[v].np, ds_cstr(&str));
-
-         statistic_destroy ((struct statistic *)result->metrics[v].np);
-
-         ds_destroy (&str);
+          struct chart *npp, *dnpp;
+          struct casereader *reader;
+          struct np *np;
+
+         ds_init_empty (&label);
+         ds_put_format (&label, "%s ", var_get_name (dependent_var[v]));
+         factor_to_string (fctr, result, &label);
+
+          np = result->metrics[v].np;
+          reader = casewriter_make_reader (np->writer);
+          npp = np_plot_create (np, reader, ds_cstr (&label));
+          dnpp = dnp_plot_create (np, reader, ds_cstr (&label));
+
+         ds_destroy (&label);
+
+          if (npp == NULL || dnpp == NULL)
+            {
+              msg (MW, _("Not creating NP plot because data set is empty."));
+              chart_unref (npp);
+              chart_unref (dnpp);
+            }
+          else
+            {
+              chart_submit (npp);
+              chart_submit (dnpp);
+            }
+
+         statistic_destroy (&np->parent.parent);
        }
     }
 }
@@ -448,15 +389,25 @@ show_histogram (const struct variable **dependent_var,
          struct string str;
          const struct factor_result *result =
            ll_data (ll, struct factor_result, ll);
+          struct histogram *histogram;
+          double mean, var, n;
+
+          histogram = result->metrics[v].histogram;
+          if (histogram == NULL)
+            {
+              /* Probably all values are SYSMIS. */
+              continue;
+            }
 
          ds_init_empty (&str);
          ds_put_format (&str, "%s ", var_get_name (dependent_var[v]));
 
          factor_to_string (fctr, result, &str);
 
-         histogram_plot ((struct histogram *) result->metrics[v].histogram,
-                         ds_cstr (&str),
-                         (struct moments1 *) result->metrics[v].moments);
+          moments1_calculate (result->metrics[v].moments,
+                              &n, &mean, &var, NULL,  NULL);
+          chart_submit (histogram_chart_create (histogram, ds_cstr (&str),
+                                                n, mean, sqrt (var), false));
 
          ds_destroy (&str);
        }
@@ -474,25 +425,18 @@ show_boxplot_groups (const struct variable **dependent_var,
 
   for (v = 0; v < n_dep_var; ++v)
     {
-      struct ll *ll;
-      int f = 0;
-      struct chart *ch = chart_create ();
+      const struct factor_result *result;
+      struct boxplot *boxplot;
       double y_min = DBL_MAX;
       double y_max = -DBL_MAX;
+      char *title;
 
-      for (ll = ll_head (&fctr->result_list);
-          ll != ll_null (&fctr->result_list);
-          ll = ll_next (ll))
+      ll_for_each (result, struct factor_result, ll, &fctr->result_list)
        {
+          struct factor_metrics *metrics = &result->metrics[v];
+         const struct ll_list *max_list = extrema_list (metrics->maxima);
+         const struct ll_list *min_list = extrema_list (metrics->minima);
          const struct extremum  *max, *min;
-         const struct factor_result *result =
-           ll_data (ll, struct factor_result, ll);
-
-         const struct ll_list *max_list =
-           extrema_list (result->metrics[v].maxima);
-
-         const struct ll_list *min_list =
-           extrema_list (result->metrics[v].minima);
 
          if ( ll_is_empty (max_list))
            {
@@ -500,52 +444,34 @@ show_boxplot_groups (const struct variable **dependent_var,
              continue;
            }
 
-         max = (const struct extremum *)
-           ll_data (ll_head(max_list), struct extremum, ll);
-
-          min = (const struct extremum *)
-           ll_data (ll_head (min_list), struct extremum, ll);
+         max = ll_data (ll_head(max_list), struct extremum, ll);
+          min = ll_data (ll_head (min_list), struct extremum, ll);
 
          y_max = MAX (y_max, max->value);
          y_min = MIN (y_min, min->value);
        }
 
-      boxplot_draw_yscale (ch, y_max, y_min);
-
-      if ( fctr->indep_var[0])
-       chart_write_title (ch, _("Boxplot of %s vs. %s"),
+      if (fctr->indep_var[0])
+       title = xasprintf (_("Boxplot of %s vs. %s"),
                           var_to_string (dependent_var[v]),
-                          var_to_string (fctr->indep_var[0]) );
+                          var_to_string (fctr->indep_var[0]));
       else
-       chart_write_title (ch, _("Boxplot of %s"),
-                          var_to_string (dependent_var[v]));
+       title = xasprintf (_("Boxplot of %s"),
+                           var_to_string (dependent_var[v]));
+      boxplot = boxplot_create (y_min, y_max, title);
+      free (title);
 
-      for (ll = ll_head (&fctr->result_list);
-          ll != ll_null (&fctr->result_list);
-          ll = ll_next (ll))
+      ll_for_each (result, struct factor_result, ll, &fctr->result_list)
        {
-         const struct factor_result *result =
-           ll_data (ll, struct factor_result, ll);
-
-         struct string str;
-         const double box_width = (ch->data_right - ch->data_left)
-           / (ll_count (&fctr->result_list) * 2.0 ) ;
-
-         const double box_centre = (f++ * 2 + 1) * box_width + ch->data_left;
-
-         ds_init_empty (&str);
+          struct factor_metrics *metrics = &result->metrics[v];
+         struct string str = DS_EMPTY_INITIALIZER;
          factor_to_string_concise (fctr, result, &str);
-
-         boxplot_draw_boxplot (ch,
-                               box_centre, box_width,
-                               (const struct box_whisker *)
-                                result->metrics[v].box_whisker,
-                               ds_cstr (&str));
-
+          boxplot_add_box (boxplot, metrics->box_whisker, ds_cstr (&str));
+          metrics->box_whisker = NULL;
          ds_destroy (&str);
        }
 
-      chart_submit (ch);
+      chart_submit (boxplot_get_chart (boxplot));
     }
 }
 
@@ -558,74 +484,42 @@ show_boxplot_variables (const struct variable **dependent_var,
                        )
 
 {
+  const struct factor_result *result;
   int v;
-  struct ll *ll;
-  const struct ll_list *result_list = &fctr->result_list;
-
-  for (ll = ll_head (result_list);
-       ll != ll_null (result_list);
-       ll = ll_next (ll))
 
+  ll_for_each (result, struct factor_result, ll, &fctr->result_list)
     {
       struct string title;
-      struct chart *ch = chart_create ();
       double y_min = DBL_MAX;
       double y_max = -DBL_MAX;
-
-      const struct factor_result *result =
-       ll_data (ll, struct factor_result, ll);
-
-      const double box_width = (ch->data_right - ch->data_left)
-       / (n_dep_var * 2.0 ) ;
+      struct boxplot *boxplot;
 
       for (v = 0; v < n_dep_var; ++v)
        {
-         const struct ll *max_ll =
-           ll_head (extrema_list (result->metrics[v].maxima));
-         const struct ll *min_ll =
-           ll_head (extrema_list (result->metrics[v].minima));
-
-         const struct extremum  *max =
-           (const struct extremum *) ll_data (max_ll, struct extremum, ll);
-
-          const struct extremum  *min =
-           (const struct extremum *) ll_data (min_ll, struct extremum, ll);
+          const struct factor_metrics *metrics = &result->metrics[v];
+         const struct ll *max_ll = ll_head (extrema_list (metrics->maxima));
+         const struct ll *min_ll = ll_head (extrema_list (metrics->minima));
+         const struct extremum *max = ll_data (max_ll, struct extremum, ll);
+          const struct extremum *min = ll_data (min_ll, struct extremum, ll);
 
          y_max = MAX (y_max, max->value);
          y_min = MIN (y_min, min->value);
        }
 
-
-      boxplot_draw_yscale (ch, y_max, y_min);
-
       ds_init_empty (&title);
       factor_to_string (fctr, result, &title);
-
-#if 0
-      ds_put_format (&title, "%s = ", var_get_name (fctr->indep_var[0]));
-      var_append_value_name (fctr->indep_var[0], &result->value[0], &title);
-#endif
-
-      chart_write_title (ch, ds_cstr (&title));
+      boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
       ds_destroy (&title);
 
       for (v = 0; v < n_dep_var; ++v)
        {
-         struct string str;
-         const double box_centre = (v * 2 + 1) * box_width + ch->data_left;
-
-         ds_init_empty (&str);
-         ds_init_cstr (&str, var_get_name (dependent_var[v]));
-
-         boxplot_draw_boxplot (ch,
-                               box_centre, box_width,
-                               (const struct box_whisker *) result->metrics[v].box_whisker,
-                               ds_cstr (&str));
-
-         ds_destroy (&str);
+          struct factor_metrics *metrics = &result->metrics[v];
+          boxplot_add_box (boxplot, metrics->box_whisker,
+                           var_get_name (dependent_var[v]));
+          metrics->box_whisker = NULL;
        }
 
-      chart_submit (ch);
+      chart_submit (boxplot_get_chart (boxplot));
     }
 }
 
@@ -674,16 +568,14 @@ output_examine (const struct dictionary *dict)
       if ( cmd.sbc_percentiles)
        show_percentiles (dependent_vars, n_dependent_vars, factor);
 
-      if (cmd.a_plot[XMN_PLT_BOXPLOT] &&
-         cmd.cmp == XMN_GROUPS)
-       show_boxplot_groups (dependent_vars, n_dependent_vars, factor);
-
-
-      if (cmd.a_plot[XMN_PLT_BOXPLOT] &&
-         cmd.cmp == XMN_VARIABLES)
-       show_boxplot_variables (dependent_vars, n_dependent_vars,
-                               factor);
-
+      if (cmd.a_plot[XMN_PLT_BOXPLOT])
+        {
+          if (cmd.cmp == XMN_GROUPS)
+            show_boxplot_groups (dependent_vars, n_dependent_vars, factor);
+          else if (cmd.cmp == XMN_VARIABLES)
+            show_boxplot_variables (dependent_vars, n_dependent_vars, factor);
+        }
+      
       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
        show_histogram (dependent_vars, n_dependent_vars, factor);
 
@@ -989,15 +881,13 @@ examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
 
          metric->n_ptiles = percentile_list.n_data;
 
-         metric->ptl = xcalloc (metric->n_ptiles,
-                                sizeof (struct percentile *));
+         metric->ptl = xcalloc (metric->n_ptiles, sizeof *metric->ptl);
 
          metric->quartiles = xcalloc (3, sizeof (*metric->quartiles));
 
          for (i = 0 ; i < metric->n_ptiles; ++i)
            {
-             metric->ptl[i] = (struct percentile *)
-               percentile_create (percentile_list.data[i] / 100.0, metric->n_valid);
+             metric->ptl[i] = percentile_create (percentile_list.data[i] / 100.0, metric->n_valid);
 
              if ( percentile_list.data[i] == 25)
                metric->quartiles[0] = metric->ptl[i];
@@ -1018,18 +908,18 @@ examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
              n_os ++;
            }
 
-         os = xcalloc (sizeof (struct order_stats *), n_os);
+          os = xcalloc (n_os, sizeof *os);
 
          for (i = 0 ; i < metric->n_ptiles ; ++i )
            {
-             os[i] = (struct order_stats *) metric->ptl[i];
+             os[i] = &metric->ptl[i]->parent;
            }
 
-         os[i] = (struct order_stats *) metric->tukey_hinges;
-         os[i+1] = (struct order_stats *) metric->trimmed_mean;
+         os[i] = &metric->tukey_hinges->parent;
+         os[i+1] = &metric->trimmed_mean->parent;
 
          if (cmd->a_plot[XMN_PLT_NPPLOT])
-           os[i+2] = metric->np;
+           os[i+2] = &metric->np->parent;
 
          order_stats_accumulate (os, n_os,
                                  casereader_clone (metric->up_reader),
@@ -1080,7 +970,7 @@ examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
            {
              struct factor_metrics *metric = &result->metrics[v];
              if ( metric->histogram)
-               histogram_add ((struct histogram *) metric->histogram,
+               histogram_add (metric->histogram,
                               case_data (c, dependent_vars[v])->f, weight);
            }
          case_unref (c);
@@ -1096,13 +986,13 @@ examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
          struct factor_metrics *metric = &result->metrics[v];
           int n_vals = caseproto_get_n_widths (casereader_get_proto (
                                                  metric->up_reader));
+          struct order_stats *os;
 
          metric->box_whisker =
-           box_whisker_create ((struct tukey_hinges *) metric->tukey_hinges,
-                               cmd->v_id, n_vals - 1);
+           box_whisker_create ( metric->tukey_hinges, cmd->v_id, n_vals - 1);
 
-         order_stats_accumulate ((struct order_stats **) &metric->box_whisker,
-                                 1,
+          os = &metric->box_whisker->parent;
+         order_stats_accumulate ( &os, 1,
                                  casereader_clone (metric->up_reader),
                                  wv, dependent_vars[v], MV_ANY);
        }
@@ -1242,10 +1132,10 @@ show_summary (const struct variable **dependent_var, int n_dep_var,
 
   n_cols = heading_columns + 6;
 
-  tbl = tab_create (n_cols, n_rows, 0);
+  tbl = tab_create (n_cols, n_rows);
   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 
-  tab_dim (tbl, tab_natural_dimensions, NULL);
+  tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
 
   /* Outline the box */
   tab_box (tbl,
@@ -1479,10 +1369,10 @@ show_descriptives (const struct variable **dependent_var,
 
   n_cols = heading_columns + 2;
 
-  tbl = tab_create (n_cols, n_rows, 0);
+  tbl = tab_create (n_cols, n_rows);
   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 
-  tab_dim (tbl, tab_natural_dimensions, NULL);
+  tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
 
   /* Outline the box */
   tab_box (tbl,
@@ -1658,7 +1548,7 @@ show_descriptives (const struct variable **dependent_var,
          tab_double (tbl, n_cols - 2,
                     heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
                     TAB_CENTER,
-                    trimmed_mean_calculate ((struct trimmed_mean *) result->metrics[v].trimmed_mean),
+                    trimmed_mean_calculate (result->metrics[v].trimmed_mean),
                     NULL);
 
 
@@ -1791,10 +1681,10 @@ show_extremes (const struct variable **dependent_var,
 
   n_cols = heading_columns + 2;
 
-  tbl = tab_create (n_cols, n_rows, 0);
+  tbl = tab_create (n_cols, n_rows);
   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 
-  tab_dim (tbl, tab_natural_dimensions, NULL);
+  tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
 
   /* Outline the box */
   tab_box (tbl,
@@ -1995,10 +1885,10 @@ show_percentiles (const struct variable **dependent_var,
 
   n_cols = heading_columns + n_percentiles;
 
-  tbl = tab_create (n_cols, n_rows, 0);
+  tbl = tab_create (n_cols, n_rows);
   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 
-  tab_dim (tbl, tab_natural_dimensions, NULL);
+  tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
 
   /* Outline the box */
   tab_box (tbl,
@@ -2078,8 +1968,7 @@ show_percentiles (const struct variable **dependent_var,
 
          tab_vline (tbl, TAL_1, n_cols - n_percentiles -1, heading_rows, n_rows - 1);
 
-         tukey_hinges_calculate ((struct tukey_hinges *) result->metrics[v].tukey_hinges,
-                                 hinges);
+         tukey_hinges_calculate (result->metrics[v].tukey_hinges, hinges);
 
          for (j = 0; j < n_percentiles; ++j)
            {