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   if (NULL == histogram)
 193     {
 194       casereader_destroy (input);
 195       return;
 196     }
 197
 198   for (;(c = casereader_read (input)) != NULL; case_unref (c))
 199     {
 200       const double x      = case_data_idx (c, HG_IDX_X)->f;
 201       const double weight = case_data_idx (c, HG_IDX_WT)->f;
 202       moments_pass_two (cmd->es[0].mom, x, weight);
 203       histogram_add (histogram, x, weight);
 204     }
 205   casereader_destroy (input);
 206
 207
 208   {
 209     double n, mean, var;
 210
 211     struct string label;
 212
 213     ds_init_cstr (&label,
 214                   var_to_string (cmd->dep_vars[0]));
 215
 216     moments_calculate (cmd->es[0].mom, &n, &mean, &var, NULL, NULL);
 217
 218     chart_item_submit
 219       ( histogram_chart_create (histogram->gsl_hist,
 220                                 ds_cstr (&label), n, mean,
 221                                 sqrt (var), false));
 222
 223     statistic_destroy (&histogram->parent);
 224     ds_destroy (&label);
 225   }
 226 }
 227
 228 static void
 229 cleanup_exploratory_stats (struct graph *cmd)
 230 {
 231   int v;
 232
 233   for (v = 0; v < cmd->n_dep_vars; ++v)
 234     {
 235       moments_destroy (cmd->es[v].mom);
 236     }
 237 }
 238
 239
 240 static void
 241 run_graph (struct graph *cmd, struct casereader *input)
 242 {
 243   struct ccase *c;
 244   struct casereader *reader;
 245   struct casewriter *writer;
 246
 247   cmd->es = pool_calloc (cmd->pool,cmd->n_dep_vars,sizeof(struct exploratory_stats));
 248   for(int v=0;v<cmd->n_dep_vars;v++)
 249     {
 250       cmd->es[v].mom = moments_create (MOMENT_KURTOSIS);
 251       cmd->es[v].cmin = DBL_MAX;
 252       cmd->es[v].maximum = -DBL_MAX;
 253       cmd->es[v].minimum =  DBL_MAX;
 254     }
 255   /* Always remove cases listwise. This is correct for */
 256   /* the histogram because there is only one variable  */
 257   /* and a simple bivariate scatterplot                */
 258   /* if ( cmd->missing_pw == false)                    */
 259     input = casereader_create_filter_missing (input,
 260                                               cmd->dep_vars,
 261                                               cmd->n_dep_vars,
 262                                               cmd->dep_excl,
 263                                               NULL,
 264                                               NULL);
 265
 266   writer = autopaging_writer_create (cmd->gr_proto);
 267
 268   /* The case data is copied to a new writer        */
 269   /* The setup of the case depends on the Charttype */
 270   /* For Scatterplot x is assumed in dep_vars[0]    */
 271   /*                 y is assumed in dep_vars[1]    */
 272   /* For Histogram   x is assumed in dep_vars[0]    */
 273   assert(SP_IDX_X == 0 && SP_IDX_Y == 1 && HG_IDX_X == 0);
 274
 275   for (;(c = casereader_read (input)) != NULL; case_unref (c))
 276     {
 277       struct ccase *outcase = case_create (cmd->gr_proto);
 278       const double weight = dict_get_case_weight (cmd->dict,c,NULL);
 279       if (cmd->chart_type == CT_HISTOGRAM)
 280         case_data_rw_idx (outcase, HG_IDX_WT)->f = weight;
 281       if (cmd->chart_type == CT_SCATTERPLOT && cmd->byvar)
 282         value_copy (case_data_rw_idx (outcase, SP_IDX_BY),
 283                     case_data (c, cmd->byvar),
 284                     var_get_width (cmd->byvar));
 285       for(int v=0;v<cmd->n_dep_vars;v++)
 286         {
 287           const struct variable *var = cmd->dep_vars[v];
 288           const double x = case_data (c, var)->f;
 289
 290           if (var_is_value_missing (var, case_data (c, var), cmd->dep_excl))
 291             {
 292               cmd->es[v].missing += weight;
 293               continue;
 294             }
 295           /* Magically v value fits to SP_IDX_X, SP_IDX_Y, HG_IDX_X */
 296           case_data_rw_idx (outcase, v)->f = x;
 297
 298           if (x > cmd->es[v].maximum)
 299             cmd->es[v].maximum = x;
 300
 301           if (x < cmd->es[v].minimum)
 302             cmd->es[v].minimum =  x;
 303
 304           cmd->es[v].non_missing += weight;
 305
 306           moments_pass_one (cmd->es[v].mom, x, weight);
 307
 308           cmd->es[v].cc += weight;
 309
 310           if (cmd->es[v].cmin > weight)
 311             cmd->es[v].cmin = weight;
 312         }
 313       casewriter_write (writer,outcase);
 314     }
 315
 316   reader = casewriter_make_reader (writer);
 317
 318   switch (cmd->chart_type)
 319     {
 320     case CT_HISTOGRAM:
 321       show_histogr (cmd,reader);
 322       break;
 323     case CT_SCATTERPLOT:
 324       show_scatterplot (cmd,reader);
 325       break;
 326     default:
 327       NOT_REACHED ();
 328       break;
 329     };
 330
 331   casereader_destroy (input);
 332   cleanup_exploratory_stats (cmd);
 333 }
 334
 335
 336 int
 337 cmd_graph (struct lexer *lexer, struct dataset *ds)
 338 {
 339   struct graph graph;
 340
 341   graph.missing_pw = false;
 342
 343   graph.pool = pool_create ();
 344
 345   graph.dep_excl = MV_ANY;
 346   graph.fctr_excl = MV_ANY;
 347
 348   graph.dict = dataset_dict (ds);
 349
 350
 351   /* ---------------- graph ------------------ */
 352   graph.dep_vars = NULL;
 353   graph.chart_type = CT_NONE;
 354   graph.scatter_type = ST_BIVARIATE;
 355   graph.byvar = NULL;
 356   graph.gr_proto = caseproto_create ();
 357
 358   while (lex_token (lexer) != T_ENDCMD)
 359     {
 360       lex_match (lexer, T_SLASH);
 361
 362       if (lex_match_id (lexer, "HISTOGRAM"))
 363         {
 364           if (graph.chart_type != CT_NONE)
 365             {
 366               lex_error (lexer, _("Only one chart type is allowed."));
 367               goto error;
 368             }
 369           if (!lex_force_match (lexer, T_EQUALS))
 370             goto error;
 371           graph.chart_type = CT_HISTOGRAM;
 372           if (!parse_variables_const (lexer, graph.dict,
 373                                       &graph.dep_vars, &graph.n_dep_vars,
 374                                       PV_NO_DUPLICATE | PV_NUMERIC))
 375             goto error;
 376           if (graph.n_dep_vars > 1)
 377             {
 378               lex_error (lexer, _("Only one variable is allowed."));
 379               goto error;
 380             }
 381         }
 382       else if (lex_match_id (lexer, "SCATTERPLOT"))
 383         {
 384           if (graph.chart_type != CT_NONE)
 385             {
 386               lex_error (lexer, _("Only one chart type is allowed."));
 387               goto error;
 388             }
 389           graph.chart_type = CT_SCATTERPLOT;
 390           if (lex_match (lexer, T_LPAREN))
 391             {
 392               if (lex_match_id (lexer, "BIVARIATE"))
 393                 {
 394                   /* This is the default anyway */
 395                 }
 396               else if (lex_match_id (lexer, "OVERLAY"))
 397                 {
 398                   lex_error (lexer, _("%s is not yet implemented."),"OVERLAY");
 399                   goto error;
 400                 }
 401               else if (lex_match_id (lexer, "MATRIX"))
 402                 {
 403                   lex_error (lexer, _("%s is not yet implemented."),"MATRIX");
 404                   goto error;
 405                 }
 406               else if (lex_match_id (lexer, "XYZ"))
 407                 {
 408                   lex_error(lexer, _("%s is not yet implemented."),"XYZ");
 409                   goto error;
 410                 }
 411               else
 412                 {
 413                   lex_error_expecting (lexer, "BIVARIATE", NULL);
 414                   goto error;
 415                 }
 416               if (!lex_force_match (lexer, T_RPAREN))
 417                 goto error;
 418             }
 419           if (!lex_force_match (lexer, T_EQUALS))
 420             goto error;
 421
 422           if (!parse_variables_const (lexer, graph.dict,
 423                                       &graph.dep_vars, &graph.n_dep_vars,
 424                                       PV_NO_DUPLICATE | PV_NUMERIC))
 425             goto error;
 426
 427           if (graph.scatter_type == ST_BIVARIATE && graph.n_dep_vars != 1)
 428             {
 429               lex_error(lexer, _("Only one variable is allowed."));
 430               goto error;
 431             }
 432
 433           if (!lex_force_match (lexer, T_WITH))
 434             goto error;
 435
 436           if (!parse_variables_const (lexer, graph.dict,
 437                                       &graph.dep_vars, &graph.n_dep_vars,
 438                                       PV_NO_DUPLICATE | PV_NUMERIC | PV_APPEND))
 439             goto error;
 440
 441           if (graph.scatter_type == ST_BIVARIATE && graph.n_dep_vars != 2)
 442             {
 443               lex_error (lexer, _("Only one variable is allowed."));
 444               goto error;
 445             }
 446
 447           if (lex_match (lexer, T_BY))
 448             {
 449               const struct variable *v = NULL;
 450               if (!lex_match_variable (lexer,graph.dict,&v))
 451                 {
 452                   lex_error (lexer, _("Variable expected"));
 453                   goto error;
 454                 }
 455               graph.byvar = v;
 456             }
 457         }
 458       else if (lex_match_id (lexer, "BAR"))
 459         {
 460           lex_error (lexer, _("%s is not yet implemented."),"BAR");
 461           goto error;
 462         }
 463       else if (lex_match_id (lexer, "LINE"))
 464         {
 465           lex_error (lexer, _("%s is not yet implemented."),"LINE");
 466           goto error;
 467         }
 468       else if (lex_match_id (lexer, "PIE"))
 469         {
 470           lex_error (lexer, _("%s is not yet implemented."),"PIE");
 471           goto error;
 472         }
 473       else if (lex_match_id (lexer, "ERRORBAR"))
 474         {
 475           lex_error (lexer, _("%s is not yet implemented."),"ERRORBAR");
 476           goto error;
 477         }
 478       else if (lex_match_id (lexer, "PARETO"))
 479         {
 480           lex_error (lexer, _("%s is not yet implemented."),"PARETO");
 481           goto error;
 482         }
 483       else if (lex_match_id (lexer, "TITLE"))
 484         {
 485           lex_error (lexer, _("%s is not yet implemented."),"TITLE");
 486           goto error;
 487         }
 488       else if (lex_match_id (lexer, "SUBTITLE"))
 489         {
 490           lex_error (lexer, _("%s is not yet implemented."),"SUBTITLE");
 491           goto error;
 492         }
 493       else if (lex_match_id (lexer, "FOOTNOTE"))
 494         {
 495           lex_error (lexer, _("%s is not yet implemented."),"FOOTNOTE");
 496           lex_error (lexer, _("FOOTNOTE is not implemented yet for GRAPH"));
 497           goto error;
 498         }
 499       else if (lex_match_id (lexer, "MISSING"))
 500         {
 501           lex_match (lexer, T_EQUALS);
 502
 503           while (lex_token (lexer) != T_ENDCMD
 504                  && lex_token (lexer) != T_SLASH)
 505             {
 506               if (lex_match_id (lexer, "LISTWISE"))
 507                 {
 508                   graph.missing_pw = false;
 509                 }
 510               else if (lex_match_id (lexer, "VARIABLE"))
 511                 {
 512                   graph.missing_pw = true;
 513                 }
 514               else if (lex_match_id (lexer, "EXCLUDE"))
 515                 {
 516                   graph.dep_excl = MV_ANY;
 517                 }
 518               else if (lex_match_id (lexer, "INCLUDE"))
 519                 {
 520                   graph.dep_excl = MV_SYSTEM;
 521                 }
 522               else if (lex_match_id (lexer, "REPORT"))
 523                 {
 524                   graph.fctr_excl = MV_NEVER;
 525                 }
 526               else if (lex_match_id (lexer, "NOREPORT"))
 527                 {
 528                   graph.fctr_excl = MV_ANY;
 529                 }
 530               else
 531                 {
 532                   lex_error (lexer, NULL);
 533                   goto error;
 534                 }
 535             }
 536         }
 537       else
 538         {
 539           lex_error (lexer, NULL);
 540           goto error;
 541         }
 542     }
 543
 544   switch (graph.chart_type)
 545     {
 546     case CT_SCATTERPLOT:
 547       /* See scatterplot.h for the setup of the case prototype */
 548       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* x value - SP_IDX_X*/
 549       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* y value - SP_IDX_Y*/
 550       /* The byvar contains the plot categories for the different xy plot colors */
 551       if (graph.byvar) /* SP_IDX_BY */
 552         graph.gr_proto = caseproto_add_width (graph.gr_proto, var_get_width(graph.byvar));
 553       break;
 554     case CT_HISTOGRAM:
 555       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* x value      */
 556       graph.gr_proto = caseproto_add_width (graph.gr_proto, 0); /* weight value */
 557       break;
 558     case CT_NONE:
 559       lex_error_expecting (lexer,"HISTOGRAM","SCATTERPLOT",NULL);
 560       goto error;
 561     default:
 562       NOT_REACHED ();
 563       break;
 564     };
 565
 566   {
 567     struct casegrouper *grouper;
 568     struct casereader *group;
 569     bool ok;
 570
 571     grouper = casegrouper_create_splits (proc_open (ds), graph.dict);
 572     while (casegrouper_get_next_group (grouper, &group))
 573       run_graph (&graph, group);
 574     ok = casegrouper_destroy (grouper);
 575     ok = proc_commit (ds) && ok;
 576   }
 577
 578   free (graph.dep_vars);
 579   pool_destroy (graph.pool);
 580   caseproto_unref (graph.gr_proto);
 581
 582   return CMD_SUCCESS;
 583
 584  error:
 585   caseproto_unref (graph.gr_proto);
 586   free (graph.dep_vars);
 587   pool_destroy (graph.pool);
 588
 589   return CMD_FAILURE;
 590 }