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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005, 2006, 2007, 2009, 2010, 2011, 2012, 2013, 2014, 2016 Free Software Foundation, Inc
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <math.h>
  20 #include <gsl/gsl_cdf.h>
  21
  22 #include "data/case.h"
  23 #include "data/casegrouper.h"
  24 #include "data/casereader.h"
  25 #include "data/dataset.h"
  26 #include "data/dictionary.h"
  27 #include "data/format.h"
  28 #include "data/variable.h"
  29 #include "data/subcase.h"
  30 #include "data/casewriter.h"
  31 #include "data/short-names.h"
  32 #include "language/command.h"
  33 #include "language/lexer/lexer.h"
  34 #include "language/lexer/variable-parser.h"
  35 #include "language/stats/sort-criteria.h"
  36 #include "math/sort.h"
  37 #include "libpspp/assertion.h"
  38 #include "libpspp/i18n.h"
  39 #include "libpspp/message.h"
  40 #include "libpspp/misc.h"
  41 #include "libpspp/pool.h"
  42 #include "libpspp/string-set.h"
  43 #include "libpspp/taint.h"
  44 #include "output/pivot-table.h"
  45
  46 #include "gettext.h"
  47 #define _(msgid) gettext (msgid)
  48 #define N_(msgid) (msgid)
  49
  50 struct rank;
  51
  52 typedef double (*rank_function_t) (const struct rank*, double c, double cc, double cc_1,
  53                                    int i, double w);
  54
  55 static double rank_proportion (const struct rank *, double c, double cc, double cc_1,
  56                                int i, double w);
  57
  58 static double rank_normal (const struct rank *, double c, double cc, double cc_1,
  59                            int i, double w);
  60
  61 static double rank_percent (const struct rank *, double c, double cc, double cc_1,
  62                             int i, double w);
  63
  64 static double rank_rfraction (const struct rank *, double c, double cc, double cc_1,
  65                               int i, double w);
  66
  67 static double rank_rank (const struct rank *, double c, double cc, double cc_1,
  68                          int i, double w);
  69
  70 static double rank_n (const struct rank *, double c, double cc, double cc_1,
  71                       int i, double w);
  72
  73 static double rank_savage (const struct rank *, double c, double cc, double cc_1,
  74                            int i, double w);
  75
  76 static double rank_ntiles (const struct rank *, double c, double cc, double cc_1,
  77                            int i, double w);
  78
  79
  80 enum rank_func
  81   {
  82     RANK,
  83     NORMAL,
  84     PERCENT,
  85     RFRACTION,
  86     PROPORTION,
  87     N,
  88     NTILES,
  89     SAVAGE,
  90     n_RANK_FUNCS
  91   };
  92
  93 static const struct fmt_spec dest_format[n_RANK_FUNCS] = {
  94   [RANK]       = { .type = FMT_F, .w = 9, .d = 3 },
  95   [NORMAL]     = { .type = FMT_F, .w = 6, .d = 4 },
  96   [PERCENT]    = { .type = FMT_F, .w = 6, .d = 2 },
  97   [RFRACTION]  = { .type = FMT_F, .w = 6, .d = 4 },
  98   [PROPORTION] = { .type = FMT_F, .w = 6, .d = 4 },
  99   [N]          = { .type = FMT_F, .w = 6, .d = 0 },
 100   [NTILES]     = { .type = FMT_F, .w = 3, .d = 0 },
 101   [SAVAGE]     = { .type = FMT_F, .w = 8, .d = 4 }
 102 };
 103
 104 static const char * const function_name[n_RANK_FUNCS] = {
 105   "RANK",
 106   "NORMAL",
 107   "PERCENT",
 108   "RFRACTION",
 109   "PROPORTION",
 110   "N",
 111   "NTILES",
 112   "SAVAGE"
 113 };
 114
 115 static const rank_function_t rank_func[n_RANK_FUNCS] = {
 116   rank_rank,
 117   rank_normal,
 118   rank_percent,
 119   rank_rfraction,
 120   rank_proportion,
 121   rank_n,
 122   rank_ntiles,
 123   rank_savage
 124 };
 125
 126
 127 enum ties
 128   {
 129     TIES_LOW,
 130     TIES_HIGH,
 131     TIES_MEAN,
 132     TIES_CONDENSE
 133   };
 134
 135 enum fraction
 136   {
 137     FRAC_BLOM,
 138     FRAC_RANKIT,
 139     FRAC_TUKEY,
 140     FRAC_VW
 141   };
 142
 143 struct rank_spec
 144 {
 145   enum rank_func rfunc;
 146   const char **dest_names;
 147   const char **dest_labels;
 148 };
 149
 150 /* If NEW_NAME exists in DICT or NEW_NAMES, returns NULL without changing
 151    anything.  Otherwise, inserts NEW_NAME in NEW_NAMES and returns the copy of
 152    NEW_NAME now in NEW_NAMES. */
 153 static const char *
 154 try_new_name (const char *new_name,
 155               const struct dictionary *dict, struct string_set *new_names)
 156 {
 157   return (!dict_lookup_var (dict, new_name)
 158           && string_set_insert (new_names, new_name)
 159           ? string_set_find_node (new_names, new_name)->string
 160           : NULL);
 161 }
 162
 163 /* Returns a variable name for storing ranks of a variable named SRC_NAME
 164    according to the rank function F.  The name chosen will not be one already in
 165    DICT or NEW_NAMES.
 166
 167    If successful, adds the new name to NEW_NAMES and returns the name added.
 168    If no name can be generated, returns NULL. */
 169 static const char *
 170 rank_choose_dest_name (struct dictionary *dict, struct string_set *new_names,
 171                        enum rank_func f, const char *src_name)
 172 {
 173   char *src_name_7;
 174   char name[128];
 175   const char *s;
 176   int i;
 177
 178   /* Try the first character of the ranking function followed by the first 7
 179      bytes of the srcinal variable name. */
 180   src_name_7 = utf8_encoding_trunc (src_name, dict_get_encoding (dict), 7);
 181   snprintf (name, sizeof name, "%c%s", function_name[f][0], src_name_7);
 182   free (src_name_7);
 183   s = try_new_name (name, dict, new_names);
 184   if (s != NULL)
 185     return s;
 186
 187   /* Try "fun###". */
 188   for (i = 1; i <= 999; i++)
 189     {
 190       sprintf (name, "%.3s%03d", function_name[f], i);
 191       s = try_new_name (name, dict, new_names);
 192       if (s != NULL)
 193         return s;
 194     }
 195
 196   /* Try "RNKfn##". */
 197   for (i = 1; i <= 99; i++)
 198     {
 199       sprintf (name, "RNK%.2s%02d", function_name[f], i);
 200       s = try_new_name (name, dict, new_names);
 201       if (s != NULL)
 202         return s;
 203     }
 204
 205   msg (ME, _("Cannot generate variable name for ranking %s with %s.  "
 206              "All candidates in use."),
 207        src_name, function_name[f]);
 208   return NULL;
 209 }
 210
 211 struct rank
 212 {
 213   struct dictionary *dict;
 214
 215   struct subcase sc;
 216
 217   const struct variable **vars;
 218   size_t n_vars;
 219
 220   const struct variable **group_vars;
 221   size_t n_group_vars;
 222
 223
 224   enum mv_class exclude;
 225
 226   struct rank_spec *rs;
 227   size_t n_rs;
 228
 229   enum ties ties;
 230
 231   enum fraction fraction;
 232   int k_ntiles;
 233
 234   bool print;
 235
 236   /* Pool on which cell functions may allocate data */
 237   struct pool *pool;
 238 };
 239
 240
 241 static void
 242 destroy_rank (struct rank *rank)
 243 {
 244  free (rank->vars);
 245  free (rank->group_vars);
 246  subcase_uninit (&rank->sc);
 247  pool_destroy (rank->pool);
 248 }
 249
 250 static bool
 251 parse_into (struct lexer *lexer, struct rank *cmd,
 252             struct string_set *new_names)
 253 {
 254   int var_count = 0;
 255   struct rank_spec *rs = NULL;
 256
 257   cmd->rs = pool_realloc (cmd->pool, cmd->rs, sizeof (*cmd->rs) * (cmd->n_rs + 1));
 258   rs = &cmd->rs[cmd->n_rs];
 259
 260   if (lex_match_id (lexer, "RANK"))
 261     {
 262       rs->rfunc = RANK;
 263     }
 264   else if (lex_match_id (lexer, "NORMAL"))
 265     {
 266       rs->rfunc = NORMAL;
 267     }
 268   else if (lex_match_id (lexer, "RFRACTION"))
 269     {
 270       rs->rfunc = RFRACTION;
 271     }
 272   else if (lex_match_id (lexer, "N"))
 273     {
 274       rs->rfunc = N;
 275     }
 276   else if (lex_match_id (lexer, "SAVAGE"))
 277     {
 278       rs->rfunc = SAVAGE;
 279     }
 280   else if (lex_match_id (lexer, "PERCENT"))
 281     {
 282       rs->rfunc = PERCENT;
 283     }
 284   else if (lex_match_id (lexer, "PROPORTION"))
 285     {
 286       rs->rfunc = PROPORTION;
 287     }
 288   else if (lex_match_id (lexer, "NTILES"))
 289     {
 290       if (!lex_force_match (lexer, T_LPAREN))
 291         return false;
 292
 293       if (! lex_force_int_range (lexer, "NTILES", 1, INT_MAX))
 294         return false;
 295
 296       cmd->k_ntiles = lex_integer (lexer);
 297       lex_get (lexer);
 298
 299       if (!lex_force_match (lexer, T_RPAREN))
 300         return false;
 301
 302       rs->rfunc = NTILES;
 303     }
 304   else
 305     {
 306       lex_error (lexer, NULL);
 307       return false;
 308     }
 309
 310   cmd->n_rs++;
 311   rs->dest_names = pool_calloc (cmd->pool, cmd->n_vars,
 312                                 sizeof *rs->dest_names);
 313
 314   if (lex_match_id (lexer, "INTO"))
 315     {
 316       while(lex_token (lexer) == T_ID)
 317         {
 318           const char *name = lex_tokcstr (lexer);
 319
 320           if (var_count >= subcase_get_n_fields (&cmd->sc))
 321             msg (SE, _("Too many variables in %s clause."), "INTO");
 322           else if (dict_lookup_var (cmd->dict, name) != NULL)
 323             msg (SE, _("Variable %s already exists."), name);
 324           else if (string_set_contains (new_names, name))
 325             msg (SE, _("Duplicate variable name %s."), name);
 326           else
 327             {
 328               string_set_insert (new_names, name);
 329               rs->dest_names[var_count++] = pool_strdup (cmd->pool, name);
 330               lex_get (lexer);
 331               continue;
 332             }
 333
 334           /* Error path. */
 335           return false;
 336         }
 337     }
 338
 339   return true;
 340 }
 341
 342 /* Hardly a rank function !! */
 343 static double
 344 rank_n (const struct rank *cmd UNUSED, double c UNUSED, double cc UNUSED, double cc_1 UNUSED,
 345         int i UNUSED, double w)
 346 {
 347   return w;
 348 }
 349
 350
 351 static double
 352 rank_rank (const struct rank *cmd, double c, double cc, double cc_1,
 353            int i, double w UNUSED)
 354 {
 355   double rank;
 356
 357   if (c >= 1.0)
 358     {
 359       switch (cmd->ties)
 360         {
 361         case TIES_LOW:
 362           rank = cc_1 + 1;
 363           break;
 364         case TIES_HIGH:
 365           rank = cc;
 366           break;
 367         case TIES_MEAN:
 368           rank = cc_1 + (c + 1.0)/ 2.0;
 369           break;
 370         case TIES_CONDENSE:
 371           rank = i;
 372           break;
 373         default:
 374           NOT_REACHED ();
 375         }
 376     }
 377   else
 378     {
 379       switch (cmd->ties)
 380         {
 381         case TIES_LOW:
 382           rank = cc_1;
 383           break;
 384         case TIES_HIGH:
 385           rank = cc;
 386           break;
 387         case TIES_MEAN:
 388           rank = cc_1 + c / 2.0 ;
 389           break;
 390         case TIES_CONDENSE:
 391           rank = i;
 392           break;
 393         default:
 394           NOT_REACHED ();
 395         }
 396     }
 397
 398   return rank;
 399 }
 400
 401
 402 static double
 403 rank_rfraction (const struct rank *cmd, double c, double cc, double cc_1,
 404                 int i, double w)
 405 {
 406   return rank_rank (cmd, c, cc, cc_1, i, w) / w ;
 407 }
 408
 409
 410 static double
 411 rank_percent (const struct rank *cmd, double c, double cc, double cc_1,
 412               int i, double w)
 413 {
 414   return rank_rank (cmd, c, cc, cc_1, i, w) * 100.0 / w ;
 415 }
 416
 417
 418 static double
 419 rank_proportion (const struct rank *cmd, double c, double cc, double cc_1,
 420                  int i, double w)
 421 {
 422   const double r =  rank_rank (cmd, c, cc, cc_1, i, w) ;
 423
 424   double f;
 425
 426   switch (cmd->fraction)
 427     {
 428     case FRAC_BLOM:
 429       f =  (r - 3.0/8.0) / (w + 0.25);
 430       break;
 431     case FRAC_RANKIT:
 432       f = (r - 0.5) / w ;
 433       break;
 434     case FRAC_TUKEY:
 435       f = (r - 1.0/3.0) / (w + 1.0/3.0);
 436       break;
 437     case FRAC_VW:
 438       f = r / (w + 1.0);
 439       break;
 440     default:
 441       NOT_REACHED ();
 442     }
 443
 444
 445   return (f > 0) ? f : SYSMIS;
 446 }
 447
 448 static double
 449 rank_normal (const struct rank *cmd, double c, double cc, double cc_1,
 450              int i, double w)
 451 {
 452   double f = rank_proportion (cmd, c, cc, cc_1, i, w);
 453
 454   return gsl_cdf_ugaussian_Pinv (f);
 455 }
 456
 457 static double
 458 rank_ntiles (const struct rank *cmd, double c, double cc, double cc_1,
 459              int i, double w)
 460 {
 461   double r = rank_rank (cmd, c, cc, cc_1, i, w);
 462
 463
 464   return (floor ((r * cmd->k_ntiles) / (w + 1)) + 1);
 465 }
 466
 467 /* Expected value of the order statistics from an exponential distribution */
 468 static double
 469 ee (int j, double w_star)
 470 {
 471   int k;
 472   double sum = 0.0;
 473
 474   for (k = 1 ; k <= j; k++)
 475     sum += 1.0 / (w_star + 1 - k);
 476
 477   return sum;
 478 }
 479
 480
 481 static double
 482 rank_savage (const struct rank *cmd UNUSED, double c, double cc, double cc_1,
 483              int i UNUSED, double w)
 484 {
 485   double int_part;
 486   const int i_1 = floor (cc_1);
 487   const int i_2 = floor (cc);
 488
 489   const double w_star = (modf (w, &int_part) == 0) ? w : floor (w) + 1;
 490
 491   const double g_1 = cc_1 - i_1;
 492   const double g_2 = cc - i_2;
 493
 494   /* The second factor is infinite, when the first is zero.
 495      Therefore, evaluate the second, only when the first is non-zero */
 496   const double expr1 =  (1 - g_1) ? (1 - g_1) * ee(i_1+1, w_star) : (1 - g_1);
 497   const double expr2 =  g_2 ? g_2 * ee (i_2+1, w_star) : g_2 ;
 498
 499   if (i_1 == i_2)
 500     return ee (i_1 + 1, w_star) - 1;
 501
 502   if (i_1 + 1 == i_2)
 503     return ((expr1 + expr2)/c) - 1;
 504
 505   if (i_1 + 2 <= i_2)
 506     {
 507       int j;
 508       double sigma = 0.0;
 509       for (j = i_1 + 2 ; j <= i_2; ++j)
 510         sigma += ee (j, w_star);
 511       return ((expr1 + expr2 + sigma) / c) -1;
 512     }
 513
 514   NOT_REACHED();
 515 }
 516
 517 static double
 518 sum_weights (const struct casereader *input, int weight_idx)
 519 {
 520   if (weight_idx == -1)
 521     return casereader_count_cases (input);
 522   else
 523     {
 524       struct casereader *pass;
 525       struct ccase *c;
 526       double w;
 527
 528       w = 0.0;
 529       pass = casereader_clone (input);
 530       for (; (c = casereader_read (pass)) != NULL; case_unref (c))
 531         w += case_num_idx (c, weight_idx);
 532       casereader_destroy (pass);
 533
 534       return w;
 535     }
 536 }
 537
 538 static void
 539 rank_sorted_file (struct casereader *input,
 540                   struct casewriter *output,
 541                   int weight_idx,
 542                   const struct rank *cmd)
 543 {
 544   struct casegrouper *tie_grouper;
 545   struct casereader *tied_cases;
 546   struct subcase input_var;
 547   int tie_group = 1;
 548   struct ccase *c;
 549   double cc = 0.0;
 550   double w;
 551
 552   /* Get total group weight. */
 553   w = sum_weights (input, weight_idx);
 554
 555   /* Do ranking. */
 556   subcase_init (&input_var, 0, 0, SC_ASCEND);
 557   tie_grouper = casegrouper_create_subcase (input, &input_var);
 558   subcase_uninit (&input_var);
 559   for (; casegrouper_get_next_group (tie_grouper, &tied_cases);
 560        casereader_destroy (tied_cases))
 561     {
 562       double tw = sum_weights (tied_cases, weight_idx);
 563       double cc_1 = cc;
 564       cc += tw;
 565
 566       taint_propagate (casereader_get_taint (tied_cases),
 567                        casewriter_get_taint (output));
 568
 569       /* Rank tied cases. */
 570       for (; (c = casereader_read (tied_cases)) != NULL; case_unref (c))
 571         {
 572           struct ccase *out_case;
 573           size_t i;
 574
 575           out_case = case_create (casewriter_get_proto (output));
 576           *case_num_rw_idx (out_case, 0) = case_num_idx (c, 1);
 577           for (i = 0; i < cmd->n_rs; ++i)
 578             {
 579               rank_function_t func = rank_func[cmd->rs[i].rfunc];
 580               double rank = func (cmd, tw, cc, cc_1, tie_group, w);
 581               *case_num_rw_idx (out_case, i + 1) = rank;
 582             }
 583
 584           casewriter_write (output, out_case);
 585         }
 586       tie_group++;
 587     }
 588   casegrouper_destroy (tie_grouper);
 589 }
 590
 591
 592 static bool
 593 rank_cmd (struct dataset *ds,  const struct rank *cmd);
 594
 595 static const char *
 596 fraction_name (const struct rank *cmd)
 597 {
 598   switch (cmd->fraction)
 599     {
 600     case FRAC_BLOM:   return "BLOM";
 601     case FRAC_RANKIT: return "RANKIT";
 602     case FRAC_TUKEY:  return "TUKEY";
 603     case FRAC_VW:     return "VW";
 604     default:          NOT_REACHED ();
 605     }
 606 }
 607
 608 /* Returns a label for a variable derived from SRC_VAR with function F. */
 609 static const char *
 610 create_var_label (struct rank *cmd, const struct variable *src_var,
 611                   enum rank_func f)
 612 {
 613   struct string label;
 614   const char *pool_label;
 615
 616   ds_init_empty (&label);
 617
 618   if (cmd->n_group_vars > 0)
 619     {
 620       struct string group_var_str;
 621       int g;
 622
 623       ds_init_empty (&group_var_str);
 624
 625       for (g = 0 ; g < cmd->n_group_vars ; ++g)
 626         {
 627           if (g > 0) ds_put_cstr (&group_var_str, " ");
 628           ds_put_cstr (&group_var_str, var_get_name (cmd->group_vars[g]));
 629         }
 630
 631       ds_put_format (&label, _("%s of %s by %s"), function_name[f],
 632                      var_get_name (src_var), ds_cstr (&group_var_str));
 633       ds_destroy (&group_var_str);
 634     }
 635   else
 636     ds_put_format (&label, _("%s of %s"),
 637                    function_name[f], var_get_name (src_var));
 638
 639   pool_label = pool_strdup (cmd->pool, ds_cstr (&label));
 640
 641   ds_destroy (&label);
 642
 643   return pool_label;
 644 }
 645
 646 int
 647 cmd_rank (struct lexer *lexer, struct dataset *ds)
 648 {
 649   struct string_set new_names;
 650   struct rank rank;
 651   struct rank_spec *rs;
 652
 653   subcase_init_empty (&rank.sc);
 654
 655   rank.rs = NULL;
 656   rank.n_rs = 0;
 657   rank.exclude = MV_ANY;
 658   rank.n_group_vars = 0;
 659   rank.group_vars = NULL;
 660   rank.dict = dataset_dict (ds);
 661   rank.ties = TIES_MEAN;
 662   rank.fraction = FRAC_BLOM;
 663   rank.print = true;
 664   rank.vars = NULL;
 665   rank.pool = pool_create ();
 666
 667   string_set_init (&new_names);
 668
 669   if (lex_match_id (lexer, "VARIABLES"))
 670     if (! lex_force_match (lexer, T_EQUALS))
 671       goto error;
 672
 673   if (!parse_sort_criteria (lexer, rank.dict,
 674                             &rank.sc,
 675                             &rank.vars, NULL))
 676     goto error;
 677
 678   rank.n_vars = rank.sc.n_fields;
 679
 680   if (lex_match (lexer, T_BY))
 681     {
 682       if (! parse_variables_const (lexer, rank.dict,
 683                                     &rank.group_vars, &rank.n_group_vars,
 684                                     PV_NO_DUPLICATE | PV_NO_SCRATCH))
 685         goto error;
 686     }
 687
 688
 689   while (lex_token (lexer) != T_ENDCMD)
 690     {
 691       if (! lex_force_match (lexer, T_SLASH))
 692         goto error;
 693       if (lex_match_id (lexer, "TIES"))
 694         {
 695           if (! lex_force_match (lexer, T_EQUALS))
 696             goto error;
 697           if (lex_match_id (lexer, "MEAN"))
 698             {
 699               rank.ties = TIES_MEAN;
 700             }
 701           else if (lex_match_id (lexer, "LOW"))
 702             {
 703               rank.ties = TIES_LOW;
 704             }
 705           else if (lex_match_id (lexer, "HIGH"))
 706             {
 707               rank.ties = TIES_HIGH;
 708             }
 709           else if (lex_match_id (lexer, "CONDENSE"))
 710             {
 711               rank.ties = TIES_CONDENSE;
 712             }
 713           else
 714             {
 715               lex_error (lexer, NULL);
 716               goto error;
 717             }
 718         }
 719       else if (lex_match_id (lexer, "FRACTION"))
 720         {
 721           if (! lex_force_match (lexer, T_EQUALS))
 722             goto error;
 723           if (lex_match_id (lexer, "BLOM"))
 724             {
 725               rank.fraction = FRAC_BLOM;
 726             }
 727           else if (lex_match_id (lexer, "TUKEY"))
 728             {
 729               rank.fraction = FRAC_TUKEY;
 730             }
 731           else if (lex_match_id (lexer, "VW"))
 732             {
 733               rank.fraction = FRAC_VW;
 734             }
 735           else if (lex_match_id (lexer, "RANKIT"))
 736             {
 737               rank.fraction = FRAC_RANKIT;
 738             }
 739           else
 740             {
 741               lex_error (lexer, NULL);
 742               goto error;
 743             }
 744         }
 745       else if (lex_match_id (lexer, "PRINT"))
 746         {
 747           if (! lex_force_match (lexer, T_EQUALS))
 748             goto error;
 749           if (lex_match_id (lexer, "YES"))
 750             {
 751               rank.print = true;
 752             }
 753           else if (lex_match_id (lexer, "NO"))
 754             {
 755               rank.print = false;
 756             }
 757           else
 758             {
 759               lex_error (lexer, NULL);
 760               goto error;
 761             }
 762         }
 763       else if (lex_match_id (lexer, "MISSING"))
 764         {
 765           if (! lex_force_match (lexer, T_EQUALS))
 766             goto error;
 767           if (lex_match_id (lexer, "INCLUDE"))
 768             {
 769               rank.exclude = MV_SYSTEM;
 770             }
 771           else if (lex_match_id (lexer, "EXCLUDE"))
 772             {
 773               rank.exclude = MV_ANY;
 774             }
 775           else
 776             {
 777               lex_error (lexer, NULL);
 778               goto error;
 779             }
 780         }
 781       else if (! parse_into (lexer, &rank, &new_names))
 782         goto error;
 783     }
 784
 785
 786   /* If no rank specs are given, then apply a default */
 787   if (rank.n_rs == 0)
 788     {
 789       struct rank_spec *rs;
 790
 791       rs = pool_calloc (rank.pool, 1, sizeof *rs);
 792       rs->rfunc = RANK;
 793       rs->dest_names = pool_calloc (rank.pool, rank.n_vars,
 794                                     sizeof *rs->dest_names);
 795
 796       rank.rs = rs;
 797       rank.n_rs = 1;
 798     }
 799
 800   /* Choose variable names for all rank destinations which haven't already been
 801      created with INTO. */
 802   for (rs = rank.rs; rs < &rank.rs[rank.n_rs]; rs++)
 803     {
 804       rs->dest_labels = pool_calloc (rank.pool, rank.n_vars,
 805                                      sizeof *rs->dest_labels);
 806       for (int v = 0 ; v < rank.n_vars ;  v ++)
 807         {
 808           const char **dst_name = &rs->dest_names[v];
 809           if (*dst_name == NULL)
 810             {
 811               *dst_name = rank_choose_dest_name (rank.dict, &new_names,
 812                                                  rs->rfunc,
 813                                                  var_get_name (rank.vars[v]));
 814               if (*dst_name == NULL)
 815                 goto error;
 816             }
 817
 818           rs->dest_labels[v] = create_var_label (&rank, rank.vars[v],
 819                                                  rs->rfunc);
 820         }
 821     }
 822
 823   if (rank.print)
 824     {
 825       struct pivot_table *table = pivot_table_create (
 826         N_("Variables Created by RANK"));
 827
 828       pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("New Variable"),
 829                               N_("New Variable"), N_("Function"),
 830                               N_("Fraction"), N_("Grouping Variables"));
 831
 832       struct pivot_dimension *variables = pivot_dimension_create (
 833         table, PIVOT_AXIS_ROW, N_("Existing Variable"),
 834         N_("Existing Variable"));
 835       variables->root->show_label = true;
 836
 837       for (size_t i = 0 ; i <  rank.n_rs ; ++i)
 838         {
 839           for (size_t v = 0 ; v < rank.n_vars ;  v ++)
 840             {
 841               int row_idx = pivot_category_create_leaf (
 842                 variables->root, pivot_value_new_variable (rank.vars[v]));
 843
 844               struct string group_vars = DS_EMPTY_INITIALIZER;
 845               for (int g = 0 ; g < rank.n_group_vars ; ++g)
 846                 {
 847                   if (g)
 848                     ds_put_byte (&group_vars, ' ');
 849                   ds_put_cstr (&group_vars, var_get_name (rank.group_vars[g]));
 850                 }
 851
 852               enum rank_func rfunc = rank.rs[i].rfunc;
 853               bool has_fraction = rfunc == NORMAL || rfunc == PROPORTION;
 854               const char *entries[] =
 855                 {
 856                   rank.rs[i].dest_names[v],
 857                   function_name[rank.rs[i].rfunc],
 858                   has_fraction ? fraction_name (&rank) : NULL,
 859                   rank.n_group_vars ? ds_cstr (&group_vars) : NULL,
 860                 };
 861               for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 862                 {
 863                   const char *entry = entries[j];
 864                   if (entry)
 865                     pivot_table_put2 (table, j, row_idx,
 866                                       pivot_value_new_user_text (entry, -1));
 867                 }
 868               ds_destroy (&group_vars);
 869             }
 870         }
 871
 872       pivot_table_submit (table);
 873     }
 874
 875   /* Do the ranking */
 876   rank_cmd (ds, &rank);
 877
 878   destroy_rank (&rank);
 879   string_set_destroy (&new_names);
 880   return CMD_SUCCESS;
 881
 882  error:
 883
 884   destroy_rank (&rank);
 885   string_set_destroy (&new_names);
 886   return CMD_FAILURE;
 887 }
 888
 889 /* RANK transformation. */
 890 struct rank_trns
 891   {
 892     int order_case_idx;
 893
 894     struct rank_trns_input_var *input_vars;
 895     size_t n_input_vars;
 896
 897     size_t n_funcs;
 898   };
 899
 900 struct rank_trns_input_var
 901   {
 902     struct casereader *input;
 903     struct ccase *current;
 904
 905     struct variable **output_vars;
 906   };
 907
 908 static void
 909 advance_ranking (struct rank_trns_input_var *iv)
 910 {
 911   case_unref (iv->current);
 912   iv->current = casereader_read (iv->input);
 913 }
 914
 915 static enum trns_result
 916 rank_trns_proc (void *trns_, struct ccase **c, casenumber case_idx UNUSED)
 917 {
 918   struct rank_trns *trns = trns_;
 919   double order = case_num_idx (*c, trns->order_case_idx);
 920   struct rank_trns_input_var *iv;
 921
 922   *c = case_unshare (*c);
 923   for (iv = trns->input_vars; iv < &trns->input_vars[trns->n_input_vars]; iv++)
 924     while (iv->current != NULL)
 925       {
 926         double iv_order = case_num_idx (iv->current, 0);
 927         if (iv_order == order)
 928           {
 929             size_t i;
 930
 931             for (i = 0; i < trns->n_funcs; i++)
 932               *case_num_rw (*c, iv->output_vars[i])
 933                 = case_num_idx (iv->current, i + 1);
 934             advance_ranking (iv);
 935             break;
 936           }
 937         else if (iv_order > order)
 938           break;
 939         else
 940           advance_ranking (iv);
 941       }
 942   return TRNS_CONTINUE;
 943 }
 944
 945 static bool
 946 rank_trns_free (void *trns_)
 947 {
 948   struct rank_trns *trns = trns_;
 949   struct rank_trns_input_var *iv;
 950
 951   for (iv = trns->input_vars; iv < &trns->input_vars[trns->n_input_vars]; iv++)
 952     {
 953       casereader_destroy (iv->input);
 954       case_unref (iv->current);
 955
 956       free (iv->output_vars);
 957     }
 958   free (trns->input_vars);
 959   free (trns);
 960
 961   return true;
 962 }
 963
 964 static const struct trns_class rank_trns_class = {
 965   .name = "RANK",
 966   .execute = rank_trns_proc,
 967   .destroy = rank_trns_free,
 968 };
 969
 970 static bool
 971 rank_cmd (struct dataset *ds, const struct rank *cmd)
 972 {
 973   struct dictionary *d = dataset_dict (ds);
 974   struct variable *weight_var = dict_get_weight (d);
 975   struct casewriter **outputs;
 976   struct variable *order_var;
 977   struct casereader *input;
 978   struct rank_trns *trns;
 979   bool ok = true;
 980   int i;
 981
 982   order_var = add_permanent_ordering_transformation (ds);
 983
 984   /* Create output files. */
 985   {
 986     struct caseproto *output_proto;
 987     struct subcase by_order;
 988
 989     output_proto = caseproto_create ();
 990     for (i = 0; i < cmd->n_rs + 1; i++)
 991       output_proto = caseproto_add_width (output_proto, 0);
 992
 993     subcase_init (&by_order, 0, 0, SC_ASCEND);
 994
 995     outputs = xnmalloc (cmd->n_vars, sizeof *outputs);
 996     for (i = 0; i < cmd->n_vars; i++)
 997       outputs[i] = sort_create_writer (&by_order, output_proto);
 998
 999     subcase_uninit (&by_order);
1000     caseproto_unref (output_proto);
1001   }
1002
1003   /* Open the active file and make one pass per input variable. */
1004   input = proc_open (ds);
1005   input = casereader_create_filter_weight (input, d, NULL, NULL);
1006   for (i = 0 ; i < cmd->n_vars ; ++i)
1007     {
1008       const struct variable *input_var = cmd->vars[i];
1009       struct casereader *input_pass;
1010       struct casegrouper *split_grouper;
1011       struct casereader *split_group;
1012       struct subcase rank_ordering;
1013       struct subcase projection;
1014       struct subcase split_vars;
1015       struct subcase group_vars;
1016       int weight_idx;
1017       int j;
1018
1019       /* Discard cases that have missing values of input variable. */
1020       input_pass = i == cmd->n_vars - 1 ? input : casereader_clone (input);
1021       input_pass = casereader_create_filter_missing (input_pass, &input_var, 1,
1022                                                      cmd->exclude, NULL, NULL);
1023
1024       /* Keep only the columns we really need, to save time and space when we
1025          sort them just below.
1026
1027          After this projection, the input_pass case indexes look like:
1028
1029            - 0: input_var.
1030            - 1: order_var.
1031            - 2 and up: cmd->n_group_vars group variables
1032            - 2 + cmd->n_group_vars and up: split variables
1033            - 2 + cmd->n_group_vars + n_split_vars: weight var
1034       */
1035       subcase_init_empty (&projection);
1036       subcase_add_var_always (&projection, input_var, SC_ASCEND);
1037       subcase_add_var_always (&projection, order_var, SC_ASCEND);
1038       subcase_add_vars_always (&projection,
1039                                cmd->group_vars, cmd->n_group_vars);
1040       subcase_add_vars_always (&projection, dict_get_split_vars (d),
1041                                dict_get_n_splits (d));
1042       if (weight_var != NULL)
1043         {
1044           subcase_add_var_always (&projection, weight_var, SC_ASCEND);
1045           weight_idx = 2 + cmd->n_group_vars + dict_get_n_splits (d);
1046         }
1047       else
1048         weight_idx = -1;
1049       input_pass = casereader_project (input_pass, &projection);
1050       subcase_uninit (&projection);
1051
1052       /* Prepare 'group_vars' as the set of grouping variables. */
1053       subcase_init_empty (&group_vars);
1054       for (j = 0; j < cmd->n_group_vars; j++)
1055         subcase_add_always (&group_vars,
1056                             j + 2, var_get_width (cmd->group_vars[j]),
1057                             SC_ASCEND);
1058
1059       /* Prepare 'rank_ordering' for sorting with the group variables as
1060          primary key and the input variable as secondary key. */
1061       subcase_clone (&rank_ordering, &group_vars);
1062       subcase_add (&rank_ordering, 0, 0, subcase_get_direction (&cmd->sc, i));
1063
1064       /* Group by split variables */
1065       subcase_init_empty (&split_vars);
1066       for (j = 0; j < dict_get_n_splits (d); j++)
1067         subcase_add_always (&split_vars, 2 + j + cmd->n_group_vars,
1068                             var_get_width (dict_get_split_vars (d)[j]),
1069                             SC_ASCEND);
1070       split_grouper = casegrouper_create_subcase (input_pass, &split_vars);
1071       subcase_uninit (&split_vars);
1072       while (casegrouper_get_next_group (split_grouper, &split_group))
1073         {
1074           struct casereader *ordered;
1075           struct casegrouper *by_grouper;
1076           struct casereader *by_group;
1077
1078           ordered = sort_execute (split_group, &rank_ordering);
1079           by_grouper = casegrouper_create_subcase (ordered, &group_vars);
1080           while (casegrouper_get_next_group (by_grouper, &by_group))
1081             rank_sorted_file (by_group, outputs[i], weight_idx, cmd);
1082           ok = casegrouper_destroy (by_grouper) && ok;
1083         }
1084       subcase_uninit (&group_vars);
1085       subcase_uninit (&rank_ordering);
1086
1087       ok = casegrouper_destroy (split_grouper) && ok;
1088     }
1089   ok = proc_commit (ds) && ok;
1090
1091   /* Re-fetch the dictionary and order variable, because if TEMPORARY was in
1092      effect then there's a new dictionary. */
1093   d = dataset_dict (ds);
1094   order_var = dict_lookup_var_assert (d, "$ORDER");
1095
1096   /* Merge the original data set with the ranks (which we already sorted on
1097      $ORDER). */
1098   trns = xmalloc (sizeof *trns);
1099   trns->order_case_idx = var_get_case_index (order_var);
1100   trns->input_vars = xnmalloc (cmd->n_vars, sizeof *trns->input_vars);
1101   trns->n_input_vars = cmd->n_vars;
1102   trns->n_funcs = cmd->n_rs;
1103   for (i = 0; i < trns->n_input_vars; i++)
1104     {
1105       struct rank_trns_input_var *iv = &trns->input_vars[i];
1106       int j;
1107
1108       iv->input = casewriter_make_reader (outputs[i]);
1109       iv->current = casereader_read (iv->input);
1110       iv->output_vars = xnmalloc (trns->n_funcs, sizeof *iv->output_vars);
1111       for (j = 0; j < trns->n_funcs; j++)
1112         {
1113           struct rank_spec *rs = &cmd->rs[j];
1114           struct variable *var;
1115
1116           var = dict_create_var_assert (d, rs->dest_names[i], 0);
1117           var_set_both_formats (var, &dest_format[rs->rfunc]);
1118           var_set_label (var, rs->dest_labels[i]);
1119
1120           iv->output_vars[j] = var;
1121         }
1122     }
1123   free (outputs);
1124
1125   add_transformation (ds, &rank_trns_class, trns);
1126
1127   /* Delete our sort key, which we don't need anymore. */
1128   dict_delete_var (d, order_var);
1129
1130   return ok;
1131 }