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