pintos-os.org Git - pspp/blob - src/language/stats/graph.c

   1 /*
   2   PSPP - a program for statistical analysis.
   3   Copyright (C) 2012, 2013, 2015 Free Software Foundation, Inc.
   4
   5   This program is free software: you can redistribute it and/or modify
   6   it under the terms of the GNU General Public License as published by
   7   the Free Software Foundation, either version 3 of the License, or
   8   (at your option) any later version.
   9
  10   This program is distributed in the hope that it will be useful,
  11   but WITHOUT ANY WARRANTY; without even the implied warranty of
  12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13   GNU 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, see <http://www.gnu.org/licenses/>.
  17 */
  18
  19 /*
  20  * This module implements the graph command
  21  */
  22
  23 #include <config.h>
  24
  25 #include <math.h>
  26 #include <gsl/gsl_cdf.h>
  27
  28 #include "libpspp/assertion.h"
  29 #include "libpspp/message.h"
  30 #include "libpspp/pool.h"
  31
  32
  33 #include "data/dataset.h"
  34 #include "data/dictionary.h"
  35 #include "data/casegrouper.h"
  36 #include "data/casereader.h"
  37 #include "data/casewriter.h"
  38 #include "data/caseproto.h"
  39 #include "data/subcase.h"
  40
  41
  42 #include "data/format.h"
  43
  44 #include "math/chart-geometry.h"
  45 #include "math/histogram.h"
  46 #include "math/moments.h"
  47 #include "math/sort.h"
  48 #include "math/order-stats.h"
  49 #include "output/charts/plot-hist.h"
  50 #include "output/charts/scatterplot.h"
  51
  52 #include "language/command.h"
  53 #include "language/lexer/lexer.h"
  54 #include "language/lexer/value-parser.h"
  55 #include "language/lexer/variable-parser.h"
  56
  57 #include "output/tab.h"
  58
  59 #include "gettext.h"
  60 #define _(msgid) gettext (msgid)
  61 #define N_(msgid) msgid
  62
  63 enum chart_type
  64   {
  65     CT_NONE,
  66     CT_BAR,
  67     CT_LINE,
  68     CT_PIE,
  69     CT_ERRORBAR,
  70     CT_HILO,
  71     CT_HISTOGRAM,
  72     CT_SCATTERPLOT,
  73     CT_PARETO
  74   };
  75
  76 enum scatter_type
  77   {
  78     ST_BIVARIATE,
  79     ST_OVERLAY,
  80     ST_MATRIX,
  81     ST_XYZ
  82   };
  83
  84 /* Variable index for histogram case */
  85 enum
  86   {
  87     HG_IDX_X,
  88     HG_IDX_WT
  89   };
  90
  91 struct exploratory_stats
  92 {
  93   double missing;
  94   double non_missing;
  95
  96   struct moments *mom;
  97
  98   double minimum;
  99   double maximum;
 100
 101   /* Total weight */
 102   double cc;
 103
 104   /* The minimum weight */
 105   double cmin;
 106 };
 107
 108
 109 struct graph
 110 {
 111   struct pool *pool;
 112
 113   size_t n_dep_vars;
 114   const struct variable **dep_vars;
 115   struct exploratory_stats *es;
 116
 117   enum mv_class dep_excl;
 118   enum mv_class fctr_excl;
 119
 120   const struct dictionary *dict;
 121
 122   bool missing_pw;
 123
 124   /* ------------ Graph ---------------- */
 125   enum chart_type chart_type;
 126   enum scatter_type scatter_type;
 127   const struct variable *byvar;
 128   /* A caseproto that contains the plot data */
 129   struct caseproto *gr_proto;
 130 };
 131
 132
 133 static void
 134 show_scatterplot (const struct graph *cmd, struct casereader *input)
 135 {
 136   struct string title;
 137   struct scatterplot_chart *scatterplot;
 138   bool byvar_overflow = false;
 139
 140   ds_init_empty (&title);
 141
 142   if (cmd->byvar)
 143     {
 144       ds_put_format (&title, _("%s vs. %s by %s"),
 145                            var_to_string (cmd->dep_vars[1]),
 146                            var_to_string (cmd->dep_vars[0]),
 147                            var_to_string (cmd->byvar));
 148     }
 149   else
 150     {
 151       ds_put_format (&title, _("%s vs. %s"),
 152                      var_to_string (cmd->dep_vars[1]),
 153                      var_to_string (cmd->dep_vars[0]));
 154     }
 155
 156
 157   scatterplot = scatterplot_create (input,
 158                                     var_to_string(cmd->dep_vars[0]),
 159                                     var_to_string(cmd->dep_vars[1]),
 160                                     cmd->byvar,
 161                                     &byvar_overflow,
 162                                     ds_cstr (&title),
 163                                     cmd->es[0].minimum, cmd->es[0].maximum,
 164                                     cmd->es[1].minimum, cmd->es[1].maximum);
 165   scatterplot_chart_submit (scatterplot);
 166   ds_destroy (&title);
 167
 168   if (byvar_overflow)
 169     {
 170       msg (MW, _("Maximum number of scatterplot categories reached."
 171                  "Your BY variable has too many distinct values."
 172                  "The coloring of the plot will not be correct"));
 173     }
 174 }
 175
 176 static void
 177 show_histogr (const struct graph *cmd, struct casereader *input)
 178 {
 179   struct histogram *histogram;
 180   struct ccase *c;
 181
 182   {
 183     /* Sturges Rule */
 184     double bin_width = fabs (cmd->es[0].minimum - cmd->es[0].maximum)
 185       / (1 + log2 (cmd->es[0].cc))
 186       ;
 187
 188     histogram =
 189       histogram_create (bin_width, cmd->es[0].minimum, cmd->es[0].maximum);
 190   }
 191
 192
 193   for (;(c = casereader_read (input)) != NULL; case_unref (c))
 194     {
 195       const double x      = case_data_idx (c, HG_IDX_X)->f;
 196       const double weight = case_data_idx (c, HG_IDX_WT)->f;
 197       moments_pass_two (cmd->es[0].mom, x, weight);
 198       histogram_add (histogram, x, weight);
 199     }
 200   casereader_destroy (input);
 201
 202
 203   {
 204     double n, mean, var;
 205
 206     struct string label;
 207
 208     ds_init_cstr (&label,
 209                   var_to_string (cmd->dep_vars[0]));
 210
 211     moments_calculate (cmd->es[0].mom, &n, &mean, &var, NULL, NULL);
 212
 213     chart_item_submit
 214       ( histogram_chart_create (histogram->gsl_hist,
 215                                 ds_cstr (&label), n, mean,
 216                                 sqrt (var), false));
 217
 218     statistic_destroy (&histogram->parent);
 219     ds_destroy (&label);
 220   }
 221 }
 222
 223 static void
 224 cleanup_exploratory_stats (struct graph *cmd)
 225 {
 226   int v;
 227
 228   for (v = 0; v < cmd->n_dep_vars; ++v)
 229     {
 230       moments_destroy (cmd->es[v].mom);
 231     }
 232 }
 233
 234
 235 static void
 236 run_graph (struct graph *cmd, struct casereader *input)
 237 {
 238   struct ccase *c;
 239   struct casereader *reader;
 240   struct casewriter *writer;
 241
 242   cmd->es = pool_calloc (cmd->pool,cmd->n_dep_vars,sizeof(struct exploratory_stats));
 243   for(int v=0;v<cmd->n_dep_vars;v++)
 244     {
 245       cmd->es[v].mom = moments_create (MOMENT_KURTOSIS);
 246       cmd->es[v].cmin = DBL_MAX;
 247       cmd->es[v].maximum = -DBL_MAX;
 248       cmd->es[v].minimum =  DBL_MAX;
 249     }
 250   /* Always remove cases listwise. This is correct for */
 251   /* the histogram because there is only one variable  */
 252   /* and a simple bivariate scatterplot                */
 253   /* if ( cmd->missing_pw == false)                    */
 254     input = casereader_create_filter_missing (input,
 255                                               cmd->dep_vars,
 256                                               cmd->n_dep_vars,
 257                                               cmd->dep_excl,
 258                                               NULL,
 259                                               NULL);
 260
 261   writer = autopaging_writer_create (cmd->gr_proto);
 262
 263   /* The case data is copied to a new writer        */
 264   /* The setup of the case depends on the Charttype */
 265   /* For Scatterplot x is assumed in dep_vars[0]    */
 266   /*                 y is assumed in dep_vars[1]    */
 267   /* For Histogram   x is assumed in dep_vars[0]    */
 268   assert(SP_IDX_X == 0 && SP_IDX_Y == 1 && HG_IDX_X == 0);
 269
 270   for (;(c = casereader_read (input)) != NULL; case_unref (c))
 271     {
 272       struct ccase *outcase = case_create (cmd->gr_proto);
 273       const double weight = dict_get_case_weight (cmd->dict,c,NULL);
 274       if (cmd->chart_type == CT_HISTOGRAM)
 275         case_data_rw_idx (outcase, HG_IDX_WT)->f = weight;
 276       if (cmd->chart_type == CT_SCATTERPLOT && cmd->byvar)
 277         value_copy (case_data_rw_idx (outcase, SP_IDX_BY),
 278                     case_data (c, cmd->byvar),
 279                     var_get_width (cmd->byvar));
 280       for(int v=0;v<cmd->n_dep_vars;v++)
 281         {
 282           const struct variable *var = cmd->dep_vars[v];
 283           const double x = case_data (c, var)->f;
 284
 285           if (var_is_value_missing (var, case_data (c, var), cmd->dep_excl))
 286             {
 287               cmd->es[v].missing += weight;
 288               continue;
 289             }
 290           /* Magically v value fits to SP_IDX_X, SP_IDX_Y, HG_IDX_X */
 291           case_data_rw_idx (outcase, v)->f = x;
 292
 293           if (x > cmd->es[v].maximum)
 294             cmd->es[v].maximum = x;
 295
 296           if (x < cmd->es[v].minimum)
 297             cmd->es[v].minimum =  x;
 298
 299           cmd->es[v].non_missing += weight;
 300
 301           moments_pass_one (cmd->es[v].mom, x, weight);
 302
 303           cmd->es[v].cc += weight;
 304
 305           if (cmd->es[v].cmin > weight)
 306             cmd->es[v].cmin = weight;
 307         }
 308       casewriter_write (writer,outcase);
 309     }
 310
 311   reader = casewriter_make_reader (writer);
 312
 313   switch (cmd->chart_type)
 314     {
 315     case CT_HISTOGRAM:
 316       show_histogr (cmd,reader);
 317       break;
 318     case CT_SCATTERPLOT:
 319       show_scatterplot (cmd,reader);
 320       break;
 321     default:
 322       NOT_REACHED ();
 323       break;
 324     };
 325
 326   casereader_destroy (input);
 327   cleanup_exploratory_stats (cmd);
 328 }
 329
 330
 331 int
 332 cmd_graph (struct lexer *lexer, struct dataset *ds)
 333 {
 334   struct graph graph;
 335
 336   graph.missing_pw = false;
 337
 338   graph.pool = pool_create ();
 339
 340   graph.dep_excl = MV_ANY;
 341   graph.fctr_excl = MV_ANY;
 342
 343   graph.dict = dataset_dict (ds);
 344
 345
 346   /* ---------------- graph ------------------ */
 347   graph.dep_vars = NULL;
 348   graph.chart_type = CT_NONE;
 349   graph.scatter_type = ST_BIVARIATE;
 350   graph.byvar = NULL;
 351   graph.gr_proto = caseproto_create ();
 352
 353   while (lex_token (lexer) != T_ENDCMD)
 354     {
 355       lex_match (lexer, T_SLASH);
 356
 357       if (lex_match_id (lexer, "HISTOGRAM"))
 358         {
 359           if (graph.chart_type != CT_NONE)
 360             {
 361               lex_error (lexer, _("Only one chart type is allowed."));
 362               goto error;
 363             }
 364           if (!lex_force_match (lexer, T_EQUALS))
 365             goto error;
 366           graph.chart_type = CT_HISTOGRAM;
 367           if (!parse_variables_const (lexer, graph.dict,
 368                                       &graph.dep_vars, &graph.n_dep_vars,
 369                                       PV_NO_DUPLICATE | PV_NUMERIC))
 370             goto error;
 371           if (graph.n_dep_vars > 1)
 372             {
 373               lex_error (lexer, _("Only one variable is allowed."));
 374               goto error;
 375             }
 376         }
 377       else if (lex_match_id (lexer, "SCATTERPLOT"))
 378         {
 379           if (graph.chart_type != CT_NONE)
 380             {
 381               lex_error (lexer, _("Only one chart type is allowed."));
 382               goto error;
 383             }
 384           graph.chart_type = CT_SCATTERPLOT;
 385           if (lex_match (lexer, T_LPAREN))
 386             {
 387               if (lex_match_id (lexer, "BIVARIATE"))
 388                 {
 389                   /* This is the default anyway */
 390                 }
 391               else if (lex_match_id (lexer, "OVERLAY"))
 392                 {
 393                   lex_error (lexer, _("%s is not yet implemented."),"OVERLAY");
 394                   goto error;
 395                 }
 396               else if (lex_match_id (lexer, "MATRIX"))
 397                 {
 398                   lex_error (lexer, _("%s is not yet implemented."),"MATRIX");
 399                   goto error;
 400                 }
 401               else if (lex_match_id (lexer, "XYZ"))
 402                 {
 403                   lex_error(lexer, _("%s is not yet implemented."),"XYZ");
 404                   goto error;
 405                 }
 406               else
 407                 {
 408                   lex_error_expecting (lexer, "BIVARIATE", NULL);
 409                   goto error;
 410                 }
 411               if (!lex_force_match (lexer, T_RPAREN))
 412                 goto error;
 413             }
 414           if (!lex_force_match (lexer, T_EQUALS))
 415             goto error;
 416
 417           if (!parse_variables_const (lexer, graph.dict,
 418                                       &graph.dep_vars, &graph.n_dep_vars,
 419                                       PV_NO_DUPLICATE | PV_NUMERIC))
 420             goto error;
 421
 422           if (graph.scatter_type == ST_BIVARIATE && graph.n_dep_vars != 1)
 423             {
 424               lex_error(lexer, _("Only one variable is allowed."));
 425               goto error;
 426             }
 427
 428           if (!lex_force_match (lexer, T_WITH))
 429             goto error;
 430
 431           if (!parse_variables_const (lexer, graph.dict,
 432                                       &graph.dep_vars, &graph.n_dep_vars,
 433                                       PV_NO_DUPLICATE | PV_NUMERIC | PV_APPEND))
 434             goto error;
 435
 436           if (graph.scatter_type == ST_BIVARIATE && graph.n_dep_vars != 2)
 437             {
 438               lex_error (lexer, _("Only one variable is allowed."));
 439               goto error;
 440             }
 441
 442           if (lex_match (lexer, T_BY))
 443             {
 444               const struct variable *v = NULL;
 445               if (!lex_match_variable (lexer,graph.dict,&v))
 446                 {
 447                   lex_error (lexer, _("Variable expected"));
 448                   goto error;
 449                 }
 450               graph.byvar = v;
 451             }
 452         }
 453       else if (lex_match_id (lexer, "BAR"))
 454         {
 455           lex_error (lexer, _("%s is not yet implemented."),"BAR");
 456           goto error;
 457         }
 458       else if (lex_match_id (lexer, "LINE"))
 459         {
 460           lex_error (lexer, _("%s is not yet implemented."),"LINE");
 461           goto error;
 462         }
 463       else if (lex_match_id (lexer, "PIE"))
 464         {
 465           lex_error (lexer, _("%s is not yet implemented."),"PIE");
 466           goto error;
 467         }
 468       else if (lex_match_id (lexer, "ERRORBAR"))
 469         {
 470           lex_error (lexer, _("%s is not yet implemented."),"ERRORBAR");
 471           goto error;
 472         }
 473       else if (lex_match_id (lexer, "PARETO"))
 474         {
 475           lex_error (lexer, _("%s is not yet implemented."),"PARETO");
 476           goto error;
 477         }
 478       else if (lex_match_id (lexer, "TITLE"))
 479         {
 480           lex_error (lexer, _("%s is not yet implemented."),"TITLE");
 481           goto error;
 482         }
 483       else if (lex_match_id (lexer, "SUBTITLE"))
 484         {
 485           lex_error (lexer, _("%s is not yet implemented."),"SUBTITLE");
 486           goto error;
 487         }
 488       else if (lex_match_id (lexer, "FOOTNOTE"))
 489         {
 490           lex_error (lexer, _("%s is not yet implemented."),"FOOTNOTE");
 491           lex_error (lexer, _("FOOTNOTE is not implemented yet for GRAPH"));
 492           goto error;
 493         }
 494       else if (lex_match_id (lexer, "MISSING"))
 495         {
 496           lex_match (lexer, T_EQUALS);
 497
 498           while (lex_token (lexer) != T_ENDCMD
 499                  && lex_token (lexer) != T_SLASH)
 500             {
 501               if (lex_match_id (lexer, "LISTWISE"))
 502                 {
 503                   graph.missing_pw = false;
 504                 }
 505               else if (lex_match_id (lexer, "VARIABLE"))
 506                 {
 507                   graph.missing_pw = true;
 508                 }
 509               else if (lex_match_id (lexer, "EXCLUDE"))
 510                 {
 511                   graph.dep_excl = MV_ANY;
 512                 }
 513               else if (lex_match_id (lexer, "INCLUDE"))
 514                 {
 515                   graph.dep_excl = MV_SYSTEM;
 516                 }
 517               else if (lex_match_id (lexer, "REPORT"))
 518                 {
 519                   graph.fctr_excl = MV_NEVER;
 520                 }
 521               else if (lex_match_id (lexer, "NOREPORT"))
 522                 {
 523                   graph.fctr_excl = MV_ANY;
 524                 }
 525               else
 526                 {
 527                   lex_error (lexer, NULL);
 528                   goto error;
 529                 }
 530             }
 531         }
 532       else
 533         {
 534           lex_error (lexer, NULL);
 535           goto error;
 536         }
 537     }
 538
 539   switch (graph.chart_type)
 540     {
 541     case CT_SCATTERPLOT:
 542       /* See scatterplot.h for the setup of the case prototype */
 543       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* x value - SP_IDX_X*/
 544       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* y value - SP_IDX_Y*/
 545       /* The byvar contains the plot categories for the different xy plot colors */
 546       if (graph.byvar) /* SP_IDX_BY */
 547         graph.gr_proto = caseproto_add_width (graph.gr_proto, var_get_width(graph.byvar));
 548       break;
 549     case CT_HISTOGRAM:
 550       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* x value      */
 551       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* weight value */
 552       break;
 553     case CT_NONE:
 554       lex_error_expecting (lexer,"HISTOGRAM","SCATTERPLOT",NULL);
 555       goto error;
 556     default:
 557       NOT_REACHED ();
 558       break;
 559     };
 560
 561   {
 562     struct casegrouper *grouper;
 563     struct casereader *group;
 564     bool ok;
 565
 566     grouper = casegrouper_create_splits (proc_open (ds), graph.dict);
 567     while (casegrouper_get_next_group (grouper, &group))
 568       run_graph (&graph, group);
 569     ok = casegrouper_destroy (grouper);
 570     ok = proc_commit (ds) && ok;
 571   }
 572
 573   free (graph.dep_vars);
 574   pool_destroy (graph.pool);
 575   caseproto_unref (graph.gr_proto);
 576
 577   return CMD_SUCCESS;
 578
 579  error:
 580   caseproto_unref (graph.gr_proto);
 581   free (graph.dep_vars);
 582   pool_destroy (graph.pool);
 583
 584   return CMD_FAILURE;
 585 }