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

   1 /*
   2   PSPP - a program for statistical analysis.
   3   Copyright (C) 2012, 2013  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 struct exploratory_stats
  85 {
  86   double missing;
  87   double non_missing;
  88
  89   struct moments *mom;
  90
  91   double minimum;
  92   double maximum;
  93
  94   /* Total weight */
  95   double cc;
  96
  97   /* The minimum weight */
  98   double cmin;
  99 };
 100
 101
 102 struct graph
 103 {
 104   struct pool *pool;
 105
 106   size_t n_dep_vars;
 107   const struct variable **dep_vars;
 108   struct exploratory_stats *es;
 109
 110   enum mv_class dep_excl;
 111   enum mv_class fctr_excl;
 112
 113   const struct dictionary *dict;
 114
 115   bool missing_pw;
 116
 117   /* ------------ Graph ---------------- */
 118   enum chart_type chart_type;
 119   enum scatter_type scatter_type;
 120   const struct variable *byvar;
 121 };
 122
 123
 124 static void
 125 show_scatterplot (const struct graph *cmd, const struct casereader *input)
 126 {
 127   struct string title;
 128   struct scatterplot_chart *scatterplot;
 129   bool byvar_overflow = false;
 130
 131   ds_init_empty (&title);
 132
 133   if (cmd->byvar)
 134     {
 135       ds_put_format (&title, _("%s vs. %s by %s"),
 136                            var_to_string (cmd->dep_vars[0]),
 137                            var_to_string (cmd->dep_vars[1]),
 138                            var_to_string (cmd->byvar));
 139     }
 140   else
 141     {
 142       ds_put_format (&title, _("%s vs. %s"),
 143                      var_to_string (cmd->dep_vars[0]),
 144                      var_to_string (cmd->dep_vars[1]));
 145     }
 146
 147
 148   scatterplot = scatterplot_create(input,
 149                                    cmd->dep_vars[0],
 150                                    cmd->dep_vars[1],
 151                                    cmd->byvar,
 152                                    &byvar_overflow,
 153                                    ds_cstr (&title),
 154                                    cmd->es[0].minimum, cmd->es[0].maximum,
 155                                    cmd->es[1].minimum, cmd->es[1].maximum);
 156   scatterplot_chart_submit(scatterplot);
 157   ds_destroy(&title);
 158
 159   if (byvar_overflow)
 160     {
 161       msg (MW, _("Maximum number of scatterplot categories reached."
 162                  "Your BY variable has too many distinct values."
 163                  "The colouring of the plot will not be correct"));
 164     }
 165
 166
 167 }
 168
 169 static void
 170 show_histogr (const struct graph *cmd, const struct casereader *input)
 171 {
 172   struct histogram *histogram;
 173   struct ccase *c;
 174   struct casereader *reader;
 175
 176   {
 177     /* Sturges Rule */
 178     double bin_width = fabs (cmd->es[0].minimum - cmd->es[0].maximum)
 179       / (1 + log2 (cmd->es[0].cc))
 180       ;
 181
 182     histogram =
 183       histogram_create (bin_width, cmd->es[0].minimum, cmd->es[0].maximum);
 184   }
 185
 186
 187   for (reader=casereader_clone(input);(c = casereader_read (reader)) != NULL; case_unref (c))
 188     {
 189       const struct variable *var = cmd->dep_vars[0];
 190       const double x = case_data (c, var)->f;
 191       const double weight = dict_get_case_weight(cmd->dict,c,NULL);
 192       moments_pass_two (cmd->es[0].mom, x, weight);
 193       histogram_add (histogram, x, weight);
 194     }
 195   casereader_destroy(reader);
 196
 197
 198   {
 199     double n, mean, var;
 200
 201     struct string label;
 202
 203     ds_init_cstr (&label,
 204                   var_to_string (cmd->dep_vars[0]));
 205
 206     moments_calculate (cmd->es[0].mom, &n, &mean, &var, NULL, NULL);
 207
 208     chart_item_submit
 209       ( histogram_chart_create (histogram->gsl_hist,
 210                                 ds_cstr (&label), n, mean,
 211                                 sqrt (var), false));
 212
 213     statistic_destroy(&histogram->parent);
 214     ds_destroy (&label);
 215   }
 216 }
 217
 218 static void
 219 cleanup_exploratory_stats (struct graph *cmd)
 220 {
 221   int v;
 222
 223   for (v = 0; v < cmd->n_dep_vars; ++v)
 224     {
 225       moments_destroy (cmd->es[v].mom);
 226     }
 227 }
 228
 229
 230 static void
 231 run_graph (struct graph *cmd, struct casereader *input)
 232 {
 233   struct ccase *c;
 234   struct casereader *reader;
 235
 236
 237   cmd->es = pool_calloc(cmd->pool,cmd->n_dep_vars,sizeof(struct exploratory_stats));
 238   for(int v=0;v<cmd->n_dep_vars;v++)
 239     {
 240       cmd->es[v].mom = moments_create (MOMENT_KURTOSIS);
 241       cmd->es[v].cmin = DBL_MAX;
 242       cmd->es[v].maximum = -DBL_MAX;
 243       cmd->es[v].minimum =  DBL_MAX;
 244     }
 245   /* Always remove cases listwise. This is correct for */
 246   /* the histogram because there is only one variable  */
 247   /* and a simple bivariate scatterplot                */
 248   /* if ( cmd->missing_pw == false)                    */
 249     input = casereader_create_filter_missing (input,
 250                                               cmd->dep_vars,
 251                                               cmd->n_dep_vars,
 252                                               cmd->dep_excl,
 253                                               NULL,
 254                                               NULL);
 255
 256   for (reader = casereader_clone (input);
 257        (c = casereader_read (reader)) != NULL; case_unref (c))
 258     {
 259       const double weight = dict_get_case_weight(cmd->dict,c,NULL);
 260       for(int v=0;v<cmd->n_dep_vars;v++)
 261         {
 262           const struct variable *var = cmd->dep_vars[v];
 263           const double x = case_data (c, var)->f;
 264
 265           if (var_is_value_missing (var, case_data (c, var), cmd->dep_excl))
 266             {
 267               cmd->es[v].missing += weight;
 268               continue;
 269             }
 270
 271           if (x > cmd->es[v].maximum)
 272             cmd->es[v].maximum = x;
 273
 274           if (x < cmd->es[v].minimum)
 275             cmd->es[v].minimum =  x;
 276
 277           cmd->es[v].non_missing += weight;
 278
 279           moments_pass_one (cmd->es[v].mom, x, weight);
 280
 281           cmd->es[v].cc += weight;
 282
 283           if (cmd->es[v].cmin > weight)
 284             cmd->es[v].cmin = weight;
 285         }
 286     }
 287   casereader_destroy (reader);
 288
 289   switch (cmd->chart_type)
 290     {
 291     case CT_HISTOGRAM:
 292       reader = casereader_clone(input);
 293       show_histogr(cmd,reader);
 294       casereader_destroy(reader);
 295       break;
 296     case CT_SCATTERPLOT:
 297       reader = casereader_clone(input);
 298       show_scatterplot(cmd,reader);
 299       casereader_destroy(reader);
 300       break;
 301     default:
 302       NOT_REACHED ();
 303       break;
 304     };
 305
 306   casereader_destroy(input);
 307
 308   cleanup_exploratory_stats (cmd);
 309 }
 310
 311
 312 int
 313 cmd_graph (struct lexer *lexer, struct dataset *ds)
 314 {
 315   struct graph graph;
 316
 317   graph.missing_pw = false;
 318
 319   graph.pool = pool_create ();
 320
 321   graph.dep_excl = MV_ANY;
 322   graph.fctr_excl = MV_ANY;
 323
 324   graph.dict = dataset_dict (ds);
 325
 326
 327   /* ---------------- graph ------------------ */
 328   graph.dep_vars = NULL;
 329   graph.chart_type = CT_NONE;
 330   graph.scatter_type = ST_BIVARIATE;
 331   graph.byvar = NULL;
 332
 333   while (lex_token (lexer) != T_ENDCMD)
 334     {
 335       lex_match (lexer, T_SLASH);
 336
 337       if (lex_match_id(lexer, "HISTOGRAM"))
 338         {
 339           if (graph.chart_type != CT_NONE)
 340             {
 341               lex_error(lexer, _("Only one chart type is allowed."));
 342               goto error;
 343             }
 344           if (!lex_force_match (lexer, T_EQUALS))
 345             goto error;
 346           graph.chart_type = CT_HISTOGRAM;
 347           if (!parse_variables_const (lexer, graph.dict,
 348                                       &graph.dep_vars, &graph.n_dep_vars,
 349                                       PV_NO_DUPLICATE | PV_NUMERIC))
 350             goto error;
 351           if (graph.n_dep_vars > 1)
 352             {
 353               lex_error(lexer, _("Only one variable is allowed"));
 354               goto error;
 355             }
 356         }
 357       else if (lex_match_id (lexer, "SCATTERPLOT"))
 358         {
 359           if (graph.chart_type != CT_NONE)
 360             {
 361               lex_error(lexer, _("Only one chart type is allowed."));
 362               goto error;
 363             }
 364           graph.chart_type = CT_SCATTERPLOT;
 365           if (lex_match (lexer, T_LPAREN))
 366             {
 367               if (lex_match_id (lexer, "BIVARIATE"))
 368                 {
 369                   /* This is the default anyway */
 370                 }
 371               else if (lex_match_id (lexer, "OVERLAY"))
 372                 {
 373                   lex_error(lexer, _("%s is not yet implemented."),"OVERLAY");
 374                   goto error;
 375                 }
 376               else if (lex_match_id (lexer, "MATRIX"))
 377                 {
 378                   lex_error(lexer, _("%s is not yet implemented."),"MATRIX");
 379                   goto error;
 380                 }
 381               else if (lex_match_id (lexer, "XYZ"))
 382                 {
 383                   lex_error(lexer, _("%s is not yet implemented."),"XYZ");
 384                   goto error;
 385                 }
 386               else
 387                 {
 388                   lex_error_expecting(lexer, "BIVARIATE", NULL);
 389                   goto error;
 390                 }
 391               if (!lex_force_match (lexer, T_RPAREN))
 392                 goto error;
 393             }
 394           if (!lex_force_match (lexer, T_EQUALS))
 395             goto error;
 396
 397           if (!parse_variables_const (lexer, graph.dict,
 398                                       &graph.dep_vars, &graph.n_dep_vars,
 399                                       PV_NO_DUPLICATE | PV_NUMERIC))
 400             goto error;
 401
 402           if (graph.scatter_type == ST_BIVARIATE && graph.n_dep_vars != 1)
 403             {
 404               lex_error(lexer, _("Only one variable allowed"));
 405               goto error;
 406             }
 407
 408           if (!lex_force_match (lexer, T_WITH))
 409             goto error;
 410
 411           if (!parse_variables_const (lexer, graph.dict,
 412                                       &graph.dep_vars, &graph.n_dep_vars,
 413                                       PV_NO_DUPLICATE | PV_NUMERIC | PV_APPEND))
 414             goto error;
 415
 416           if (graph.scatter_type == ST_BIVARIATE && graph.n_dep_vars != 2)
 417             {
 418               lex_error(lexer, _("Only one variable allowed"));
 419               goto error;
 420             }
 421
 422           if (lex_match(lexer, T_BY))
 423             {
 424               const struct variable *v = NULL;
 425               if (!lex_match_variable (lexer,graph.dict,&v))
 426                 {
 427                   lex_error(lexer, _("Variable expected"));
 428                   goto error;
 429                 }
 430               graph.byvar = v;
 431             }
 432         }
 433       else if (lex_match_id (lexer, "BAR"))
 434         {
 435           lex_error (lexer, _("%s is not yet implemented."),"BAR");
 436           goto error;
 437         }
 438       else if (lex_match_id (lexer, "LINE"))
 439         {
 440           lex_error (lexer, _("%s is not yet implemented."),"LINE");
 441           goto error;
 442         }
 443       else if (lex_match_id (lexer, "PIE"))
 444         {
 445           lex_error (lexer, _("%s is not yet implemented."),"PIE");
 446           goto error;
 447         }
 448       else if (lex_match_id (lexer, "ERRORBAR"))
 449         {
 450           lex_error (lexer, _("%s is not yet implemented."),"ERRORBAR");
 451           goto error;
 452         }
 453       else if (lex_match_id (lexer, "PARETO"))
 454         {
 455           lex_error (lexer, _("%s is not yet implemented."),"PARETO");
 456           goto error;
 457         }
 458       else if (lex_match_id (lexer, "TITLE"))
 459         {
 460           lex_error (lexer, _("%s is not yet implemented."),"TITLE");
 461           goto error;
 462         }
 463       else if (lex_match_id (lexer, "SUBTITLE"))
 464         {
 465           lex_error (lexer, _("%s is not yet implemented."),"SUBTITLE");
 466           goto error;
 467         }
 468       else if (lex_match_id (lexer, "FOOTNOTE"))
 469         {
 470           lex_error (lexer, _("%s is not yet implemented."),"FOOTNOTE");
 471           lex_error (lexer, _("FOOTNOTE is not implemented yet for GRAPH"));
 472           goto error;
 473         }
 474       else if (lex_match_id (lexer, "MISSING"))
 475         {
 476           lex_match (lexer, T_EQUALS);
 477
 478           while (lex_token (lexer) != T_ENDCMD
 479                  && lex_token (lexer) != T_SLASH)
 480             {
 481               if (lex_match_id (lexer, "LISTWISE"))
 482                 {
 483                   graph.missing_pw = false;
 484                 }
 485               else if (lex_match_id (lexer, "VARIABLE"))
 486                 {
 487                   graph.missing_pw = true;
 488                 }
 489               else if (lex_match_id (lexer, "EXCLUDE"))
 490                 {
 491                   graph.dep_excl = MV_ANY;
 492                 }
 493               else if (lex_match_id (lexer, "INCLUDE"))
 494                 {
 495                   graph.dep_excl = MV_SYSTEM;
 496                 }
 497               else if (lex_match_id (lexer, "REPORT"))
 498                 {
 499                   graph.fctr_excl = MV_NEVER;
 500                 }
 501               else if (lex_match_id (lexer, "NOREPORT"))
 502                 {
 503                   graph.fctr_excl = MV_ANY;
 504                 }
 505               else
 506                 {
 507                   lex_error (lexer, NULL);
 508                   goto error;
 509                 }
 510             }
 511         }
 512       else
 513         {
 514           lex_error (lexer, NULL);
 515           goto error;
 516         }
 517     }
 518
 519   if (graph.chart_type == CT_NONE)
 520     {
 521       lex_error_expecting(lexer,"HISTOGRAM","SCATTERPLOT",NULL);
 522       goto error;
 523     }
 524
 525
 526   {
 527     struct casegrouper *grouper;
 528     struct casereader *group;
 529     bool ok;
 530
 531     grouper = casegrouper_create_splits (proc_open (ds), graph.dict);
 532     while (casegrouper_get_next_group (grouper, &group))
 533       run_graph (&graph, group);
 534     ok = casegrouper_destroy (grouper);
 535     ok = proc_commit (ds) && ok;
 536   }
 537
 538   free (graph.dep_vars);
 539   pool_destroy (graph.pool);
 540
 541   return CMD_SUCCESS;
 542
 543  error:
 544   free (graph.dep_vars);
 545   pool_destroy (graph.pool);
 546
 547   return CMD_FAILURE;
 548 }