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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2011, 2012, 2015 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 <gsl/gsl_matrix.h>
  20 #include <gsl/gsl_permutation.h>
  21 #include <gsl/gsl_sort_vector.h>
  22 #include <gsl/gsl_statistics.h>
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25
  26 #include "data/case.h"
  27 #include "data/casegrouper.h"
  28 #include "data/casereader.h"
  29 #include "data/casewriter.h"
  30 #include "data/dataset.h"
  31 #include "data/dictionary.h"
  32 #include "data/format.h"
  33 #include "data/missing-values.h"
  34 #include "language/command.h"
  35 #include "language/lexer/lexer.h"
  36 #include "language/lexer/variable-parser.h"
  37 #include "libpspp/message.h"
  38 #include "libpspp/misc.h"
  39 #include "libpspp/assertion.h"
  40 #include "libpspp/str.h"
  41 #include "math/random.h"
  42 #include "output/pivot-table.h"
  43 #include "output/text-item.h"
  44
  45 #include "gettext.h"
  46 #define _(msgid) gettext (msgid)
  47 #define N_(msgid) msgid
  48
  49 enum missing_type
  50   {
  51     MISS_LISTWISE,
  52     MISS_PAIRWISE,
  53   };
  54
  55
  56 struct qc
  57 {
  58   const struct variable **vars;
  59   size_t n_vars;
  60
  61   double epsilon;               /* The convergence criterion */
  62
  63   int ngroups;                  /* Number of group. (Given by the user) */
  64   int maxiter;                  /* Maximum iterations (Given by the user) */
  65   bool print_cluster_membership; /* true => print membership */
  66   bool print_initial_clusters;   /* true => print initial cluster */
  67   bool no_initial;              /* true => simplified initial cluster selection */
  68   bool no_update;               /* true => do not iterate  */
  69
  70   const struct variable *wv;    /* Weighting variable. */
  71
  72   enum missing_type missing_type;
  73   enum mv_class exclude;
  74 };
  75
  76 /* Holds all of the information for the functions.  int n, holds the number of
  77    observation and its default value is -1.  We set it in
  78    kmeans_recalculate_centers in first invocation. */
  79 struct Kmeans
  80 {
  81   gsl_matrix *centers;          /* Centers for groups. */
  82   gsl_matrix *updated_centers;
  83   casenumber n;
  84
  85   gsl_vector_long *num_elements_groups;
  86
  87   gsl_matrix *initial_centers;  /* Initial random centers. */
  88   double convergence_criteria;
  89   gsl_permutation *group_order; /* Group order for reporting. */
  90 };
  91
  92 static struct Kmeans *kmeans_create (const struct qc *qc);
  93
  94 static void kmeans_get_nearest_group (const struct Kmeans *kmeans,
  95                                       struct ccase *c, const struct qc *,
  96                                       int *, double *, int *, double *);
  97
  98 static void kmeans_order_groups (struct Kmeans *kmeans, const struct qc *);
  99
 100 static void kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader,
 101                             const struct qc *);
 102
 103 static void quick_cluster_show_centers (struct Kmeans *kmeans, bool initial,
 104                                         const struct qc *);
 105
 106 static void quick_cluster_show_membership (struct Kmeans *kmeans,
 107                                            const struct casereader *reader,
 108                                            const struct qc *);
 109
 110 static void quick_cluster_show_number_cases (struct Kmeans *kmeans,
 111                                              const struct qc *);
 112
 113 static void quick_cluster_show_results (struct Kmeans *kmeans,
 114                                         const struct casereader *reader,
 115                                         const struct qc *);
 116
 117 int cmd_quick_cluster (struct lexer *lexer, struct dataset *ds);
 118
 119 static void kmeans_destroy (struct Kmeans *kmeans);
 120
 121 /* Creates and returns a struct of Kmeans with given casereader 'cs', parsed
 122    variables 'variables', number of cases 'n', number of variables 'm', number
 123    of clusters and amount of maximum iterations. */
 124 static struct Kmeans *
 125 kmeans_create (const struct qc *qc)
 126 {
 127   struct Kmeans *kmeans = xmalloc (sizeof (struct Kmeans));
 128   kmeans->centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
 129   kmeans->updated_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
 130   kmeans->num_elements_groups = gsl_vector_long_alloc (qc->ngroups);
 131   kmeans->group_order = gsl_permutation_alloc (kmeans->centers->size1);
 132   kmeans->initial_centers = NULL;
 133
 134   return (kmeans);
 135 }
 136
 137 static void
 138 kmeans_destroy (struct Kmeans *kmeans)
 139 {
 140   gsl_matrix_free (kmeans->centers);
 141   gsl_matrix_free (kmeans->updated_centers);
 142   gsl_matrix_free (kmeans->initial_centers);
 143
 144   gsl_vector_long_free (kmeans->num_elements_groups);
 145
 146   gsl_permutation_free (kmeans->group_order);
 147
 148   free (kmeans);
 149 }
 150
 151 static double
 152 diff_matrix (const gsl_matrix *m1, const gsl_matrix *m2)
 153 {
 154   int i,j;
 155   double max_diff = -INFINITY;
 156   for (i = 0; i < m1->size1; ++i)
 157     {
 158       double diff = 0;
 159       for (j = 0; j < m1->size2; ++j)
 160         {
 161           diff += pow2 (gsl_matrix_get (m1,i,j) - gsl_matrix_get (m2,i,j) );
 162         }
 163       if (diff > max_diff)
 164         max_diff = diff;
 165     }
 166
 167   return max_diff;
 168 }
 169
 170
 171
 172 static double
 173 matrix_mindist (const gsl_matrix *m, int *mn, int *mm)
 174 {
 175   int i, j;
 176   double mindist = INFINITY;
 177   for (i = 0; i < m->size1 - 1; ++i)
 178     {
 179       for (j = i + 1; j < m->size1; ++j)
 180         {
 181           int k;
 182           double diff_sq = 0;
 183           for (k = 0; k < m->size2; ++k)
 184             {
 185               diff_sq += pow2 (gsl_matrix_get (m, j, k) - gsl_matrix_get (m, i, k));
 186             }
 187           if (diff_sq < mindist)
 188             {
 189               mindist = diff_sq;
 190               if (mn)
 191                 *mn = i;
 192               if (mm)
 193                 *mm = j;
 194             }
 195         }
 196     }
 197
 198   return mindist;
 199 }
 200
 201
 202 /* Return the distance of C from the group whose index is WHICH */
 203 static double
 204 dist_from_case (const struct Kmeans *kmeans, const struct ccase *c,
 205                 const struct qc *qc, int which)
 206 {
 207   int j;
 208   double dist = 0;
 209   for (j = 0; j < qc->n_vars; j++)
 210     {
 211       const union value *val = case_data (c, qc->vars[j]);
 212       if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
 213         NOT_REACHED ();
 214
 215       dist += pow2 (gsl_matrix_get (kmeans->centers, which, j) - val->f);
 216     }
 217
 218   return dist;
 219 }
 220
 221 /* Return the minimum distance of the group WHICH and all other groups */
 222 static double
 223 min_dist_from (const struct Kmeans *kmeans, const struct qc *qc, int which)
 224 {
 225   int j, i;
 226
 227   double mindist = INFINITY;
 228   for (i = 0; i < qc->ngroups; i++)
 229     {
 230       if (i == which)
 231         continue;
 232
 233       double dist = 0;
 234       for (j = 0; j < qc->n_vars; j++)
 235         {
 236           dist += pow2 (gsl_matrix_get (kmeans->centers, i, j)
 237                         - gsl_matrix_get (kmeans->centers, which, j));
 238         }
 239
 240       if (dist < mindist)
 241         {
 242           mindist = dist;
 243         }
 244     }
 245
 246   return mindist;
 247 }
 248
 249
 250
 251 /* Calculate the initial cluster centers. */
 252 static void
 253 kmeans_initial_centers (struct Kmeans *kmeans,
 254                         const struct casereader *reader,
 255                         const struct qc *qc)
 256 {
 257   struct ccase *c;
 258   int nc = 0, j;
 259
 260   struct casereader *cs = casereader_clone (reader);
 261   for (; (c = casereader_read (cs)) != NULL; case_unref (c))
 262     {
 263       bool missing = false;
 264       for (j = 0; j < qc->n_vars; ++j)
 265         {
 266           const union value *val = case_data (c, qc->vars[j]);
 267           if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
 268             {
 269               missing = true;
 270               break;
 271             }
 272
 273           if (nc < qc->ngroups)
 274             gsl_matrix_set (kmeans->centers, nc, j, val->f);
 275         }
 276
 277       if (missing)
 278         continue;
 279
 280       if (nc++ < qc->ngroups)
 281         continue;
 282
 283       if (!qc->no_initial)
 284         {
 285           int mq, mp;
 286           double delta;
 287
 288           int mn, mm;
 289           double m = matrix_mindist (kmeans->centers, &mn, &mm);
 290
 291           kmeans_get_nearest_group (kmeans, c, qc, &mq, &delta, &mp, NULL);
 292           if (delta > m)
 293             /* If the distance between C and the nearest group, is greater than the distance
 294                between the two  groups which are clostest to each
 295                other, then one group must be replaced.  */
 296             {
 297               /* Out of mn and mm, which is the clostest of the two groups to C ? */
 298               int which = (dist_from_case (kmeans, c, qc, mn)
 299                            > dist_from_case (kmeans, c, qc, mm)) ? mm : mn;
 300
 301               for (j = 0; j < qc->n_vars; ++j)
 302                 {
 303                   const union value *val = case_data (c, qc->vars[j]);
 304                   gsl_matrix_set (kmeans->centers, which, j, val->f);
 305                 }
 306             }
 307           else if (dist_from_case (kmeans, c, qc, mp) > min_dist_from (kmeans, qc, mq))
 308             /* If the distance between C and the second nearest group
 309                (MP) is greater than the smallest distance between the
 310                nearest group (MQ) and any other group, then replace
 311                MQ with C.  */
 312             {
 313               for (j = 0; j < qc->n_vars; ++j)
 314                 {
 315                   const union value *val = case_data (c, qc->vars[j]);
 316                   gsl_matrix_set (kmeans->centers, mq, j, val->f);
 317                 }
 318             }
 319         }
 320     }
 321
 322   casereader_destroy (cs);
 323
 324   kmeans->convergence_criteria = qc->epsilon * matrix_mindist (kmeans->centers, NULL, NULL);
 325
 326   /* As it is the first iteration, the variable kmeans->initial_centers is NULL
 327      and it is created once for reporting issues. */
 328   kmeans->initial_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
 329   gsl_matrix_memcpy (kmeans->initial_centers, kmeans->centers);
 330 }
 331
 332
 333 /* Return the index of the group which is nearest to the case C */
 334 static void
 335 kmeans_get_nearest_group (const struct Kmeans *kmeans, struct ccase *c,
 336                           const struct qc *qc, int *g_q, double *delta_q,
 337                           int *g_p, double *delta_p)
 338 {
 339   int result0 = -1;
 340   int result1 = -1;
 341   int i, j;
 342   double mindist0 = INFINITY;
 343   double mindist1 = INFINITY;
 344   for (i = 0; i < qc->ngroups; i++)
 345     {
 346       double dist = 0;
 347       for (j = 0; j < qc->n_vars; j++)
 348         {
 349           const union value *val = case_data (c, qc->vars[j]);
 350           if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
 351             continue;
 352
 353           dist += pow2 (gsl_matrix_get (kmeans->centers, i, j) - val->f);
 354         }
 355
 356       if (dist < mindist0)
 357         {
 358           mindist1 = mindist0;
 359           result1 = result0;
 360
 361           mindist0 = dist;
 362           result0 = i;
 363         }
 364       else if (dist < mindist1)
 365         {
 366           mindist1 = dist;
 367           result1 = i;
 368         }
 369     }
 370
 371   if (delta_q)
 372     *delta_q = mindist0;
 373
 374   if (g_q)
 375     *g_q = result0;
 376
 377
 378   if (delta_p)
 379     *delta_p = mindist1;
 380
 381   if (g_p)
 382     *g_p = result1;
 383 }
 384
 385
 386
 387 static void
 388 kmeans_order_groups (struct Kmeans *kmeans, const struct qc *qc)
 389 {
 390   gsl_vector *v = gsl_vector_alloc (qc->ngroups);
 391   gsl_matrix_get_col (v, kmeans->centers, 0);
 392   gsl_sort_vector_index (kmeans->group_order, v);
 393   gsl_vector_free (v);
 394 }
 395
 396 /* Main algorithm.
 397    Does iterations, checks convergency. */
 398 static void
 399 kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader,
 400                 const struct qc *qc)
 401 {
 402   int j;
 403
 404   kmeans_initial_centers (kmeans, reader, qc);
 405
 406   gsl_matrix_memcpy (kmeans->updated_centers, kmeans->centers);
 407
 408
 409   for (int xx = 0 ; xx < qc->maxiter ; ++xx)
 410     {
 411       gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
 412
 413       kmeans->n = 0;
 414       if (!qc->no_update)
 415         {
 416           struct casereader *r = casereader_clone (reader);
 417           struct ccase *c;
 418           for (; (c = casereader_read (r)) != NULL; case_unref (c))
 419             {
 420               int group = -1;
 421               int g;
 422               bool missing = false;
 423
 424               for (j = 0; j < qc->n_vars; j++)
 425                 {
 426                   const union value *val = case_data (c, qc->vars[j]);
 427                   if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
 428                     missing = true;
 429                 }
 430
 431               if (missing)
 432                 continue;
 433
 434               double mindist = INFINITY;
 435               for (g = 0; g < qc->ngroups; ++g)
 436                 {
 437                   double d = dist_from_case (kmeans, c, qc, g);
 438
 439                   if (d < mindist)
 440                     {
 441                       mindist = d;
 442                       group = g;
 443                     }
 444                 }
 445
 446               long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
 447               *n += qc->wv ? case_data (c, qc->wv)->f : 1.0;
 448               kmeans->n++;
 449
 450               for (j = 0; j < qc->n_vars; ++j)
 451                 {
 452                   const union value *val = case_data (c, qc->vars[j]);
 453                   if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
 454                     continue;
 455                   double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
 456                   *x += val->f * (qc->wv ? case_data (c, qc->wv)->f : 1.0);
 457                 }
 458             }
 459
 460           casereader_destroy (r);
 461         }
 462
 463       int g;
 464
 465       /* Divide the cluster sums by the number of items in each cluster */
 466       for (g = 0; g < qc->ngroups; ++g)
 467         {
 468           for (j = 0; j < qc->n_vars; ++j)
 469             {
 470               long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
 471               double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
 472               *x /= n + 1;  // Plus 1 for the initial centers
 473             }
 474         }
 475
 476
 477       gsl_matrix_memcpy (kmeans->centers, kmeans->updated_centers);
 478
 479       {
 480         kmeans->n = 0;
 481         /* Step 3 */
 482         gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
 483         gsl_matrix_set_all (kmeans->updated_centers, 0.0);
 484         struct ccase *c;
 485         struct casereader *cs = casereader_clone (reader);
 486         for (; (c = casereader_read (cs)) != NULL; case_unref (c))
 487           {
 488             int group = -1;
 489             kmeans_get_nearest_group (kmeans, c, qc, &group, NULL, NULL, NULL);
 490
 491             for (j = 0; j < qc->n_vars; ++j)
 492               {
 493                 const union value *val = case_data (c, qc->vars[j]);
 494                 if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
 495                   continue;
 496
 497                 double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
 498                 *x += val->f;
 499               }
 500
 501             long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
 502             *n += qc->wv ? case_data (c, qc->wv)->f : 1.0;
 503             kmeans->n++;
 504           }
 505         casereader_destroy (cs);
 506
 507
 508         /* Divide the cluster sums by the number of items in each cluster */
 509         for (g = 0; g < qc->ngroups; ++g)
 510           {
 511             for (j = 0; j < qc->n_vars; ++j)
 512               {
 513                 long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
 514                 double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
 515                 *x /= n ;
 516               }
 517           }
 518
 519         double d = diff_matrix (kmeans->updated_centers, kmeans->centers);
 520         if (d < kmeans->convergence_criteria)
 521           break;
 522       }
 523
 524       if (qc->no_update)
 525         break;
 526     }
 527 }
 528
 529 /* Reports centers of clusters.
 530    Initial parameter is optional for future use.
 531    If initial is true, initial cluster centers are reported.  Otherwise,
 532    resulted centers are reported. */
 533 static void
 534 quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *qc)
 535 {
 536   struct pivot_table *table
 537     = pivot_table_create (initial
 538                           ? N_("Initial Cluster Centers")
 539                           : N_("Final Cluster Centers"));
 540
 541   struct pivot_dimension *clusters
 542     = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Cluster"));
 543
 544   clusters->root->show_label = true;
 545   for (size_t i = 0; i < qc->ngroups; i++)
 546     pivot_category_create_leaf (clusters->root,
 547                                 pivot_value_new_integer (i + 1));
 548
 549   struct pivot_dimension *variables
 550     = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Variable"));
 551
 552   for (size_t i = 0; i < qc->n_vars; i++)
 553     pivot_category_create_leaf (variables->root,
 554                                 pivot_value_new_variable (qc->vars[i]));
 555
 556   const gsl_matrix *matrix = (initial
 557                               ? kmeans->initial_centers
 558                               : kmeans->centers);
 559   for (size_t i = 0; i < qc->ngroups; i++)
 560     for (size_t j = 0; j < qc->n_vars; j++)
 561       {
 562         double x = gsl_matrix_get (matrix, kmeans->group_order->data[i], j);
 563         union value v = { .f = x };
 564         pivot_table_put2 (table, i, j,
 565                           pivot_value_new_var_value (qc->vars[j], &v));
 566       }
 567
 568   pivot_table_submit (table);
 569 }
 570
 571 /* Reports cluster membership for each case. */
 572 static void
 573 quick_cluster_show_membership (struct Kmeans *kmeans,
 574                                const struct casereader *reader,
 575                                const struct qc *qc)
 576 {
 577   struct pivot_table *table = pivot_table_create (N_("Cluster Membership"));
 578
 579   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Cluster"),
 580                           N_("Cluster"));
 581
 582   struct pivot_dimension *cases
 583     = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Case Number"));
 584
 585   cases->root->show_label = true;
 586
 587   gsl_permutation *ip = gsl_permutation_alloc (qc->ngroups);
 588   gsl_permutation_inverse (ip, kmeans->group_order);
 589
 590   struct casereader *cs = casereader_clone (reader);
 591   struct ccase *c;
 592   for (int i = 0; (c = casereader_read (cs)) != NULL; i++, case_unref (c))
 593     {
 594       assert (i < kmeans->n);
 595       int clust;
 596       kmeans_get_nearest_group (kmeans, c, qc, &clust, NULL, NULL, NULL);
 597       int cluster = ip->data[clust];
 598
 599       int case_idx = pivot_category_create_leaf (cases->root,
 600                                                  pivot_value_new_integer (i + 1));
 601       pivot_table_put2 (table, 0, case_idx,
 602                         pivot_value_new_integer (cluster + 1));
 603     }
 604
 605   gsl_permutation_free (ip);
 606   pivot_table_submit (table);
 607   casereader_destroy (cs);
 608 }
 609
 610
 611 /* Reports number of cases of each single cluster. */
 612 static void
 613 quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
 614 {
 615   struct pivot_table *table
 616     = pivot_table_create (N_("Number of Cases in each Cluster"));
 617
 618   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 619                           N_("Count"));
 620
 621   struct pivot_dimension *clusters
 622     = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Clusters"));
 623
 624   struct pivot_category *group
 625     = pivot_category_create_group (clusters->root, N_("Cluster"));
 626
 627   long int total = 0;
 628   for (int i = 0; i < qc->ngroups; i++)
 629     {
 630       int cluster_idx
 631         = pivot_category_create_leaf (group, pivot_value_new_integer (i + 1));
 632       int count = kmeans->num_elements_groups->data [kmeans->group_order->data[i]];
 633       pivot_table_put2 (table, 0, cluster_idx, pivot_value_new_integer (count));
 634       total += count;
 635     }
 636
 637   int cluster_idx = pivot_category_create_leaf (clusters->root,
 638                                                 pivot_value_new_text (N_("Valid")));
 639   pivot_table_put2 (table, 0, cluster_idx, pivot_value_new_integer (total));
 640   pivot_table_submit (table);
 641 }
 642
 643 /* Reports. */
 644 static void
 645 quick_cluster_show_results (struct Kmeans *kmeans, const struct casereader *reader,
 646                             const struct qc *qc)
 647 {
 648   kmeans_order_groups (kmeans, qc); /* what does this do? */
 649
 650   if (qc->print_initial_clusters)
 651     quick_cluster_show_centers (kmeans, true, qc);
 652   quick_cluster_show_centers (kmeans, false, qc);
 653   quick_cluster_show_number_cases (kmeans, qc);
 654   if (qc->print_cluster_membership)
 655     quick_cluster_show_membership (kmeans, reader, qc);
 656 }
 657
 658 int
 659 cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
 660 {
 661   struct qc qc;
 662   struct Kmeans *kmeans;
 663   bool ok;
 664   const struct dictionary *dict = dataset_dict (ds);
 665   qc.ngroups = 2;
 666   qc.maxiter = 10;
 667   qc.epsilon = DBL_EPSILON;
 668   qc.missing_type = MISS_LISTWISE;
 669   qc.exclude = MV_ANY;
 670   qc.print_cluster_membership = false; /* default = do not output case cluster membership */
 671   qc.print_initial_clusters = false;   /* default = do not print initial clusters */
 672   qc.no_initial = false;               /* default = use well separated initial clusters */
 673   qc.no_update = false;               /* default = iterate until convergence or max iterations */
 674
 675   if (!parse_variables_const (lexer, dict, &qc.vars, &qc.n_vars,
 676                               PV_NO_DUPLICATE | PV_NUMERIC))
 677     {
 678       return (CMD_FAILURE);
 679     }
 680
 681   while (lex_token (lexer) != T_ENDCMD)
 682     {
 683       lex_match (lexer, T_SLASH);
 684
 685       if (lex_match_id (lexer, "MISSING"))
 686         {
 687           lex_match (lexer, T_EQUALS);
 688           while (lex_token (lexer) != T_ENDCMD
 689                  && lex_token (lexer) != T_SLASH)
 690             {
 691               if (lex_match_id (lexer, "LISTWISE")
 692                   || lex_match_id (lexer, "DEFAULT"))
 693                 {
 694                   qc.missing_type = MISS_LISTWISE;
 695                 }
 696               else if (lex_match_id (lexer, "PAIRWISE"))
 697                 {
 698                   qc.missing_type = MISS_PAIRWISE;
 699                 }
 700               else if (lex_match_id (lexer, "INCLUDE"))
 701                 {
 702                   qc.exclude = MV_SYSTEM;
 703                 }
 704               else if (lex_match_id (lexer, "EXCLUDE"))
 705                 {
 706                   qc.exclude = MV_ANY;
 707                 }
 708               else
 709                 {
 710                   lex_error (lexer, NULL);
 711                   goto error;
 712                 }
 713             }
 714         }
 715       else if (lex_match_id (lexer, "PRINT"))
 716         {
 717           lex_match (lexer, T_EQUALS);
 718           while (lex_token (lexer) != T_ENDCMD
 719                  && lex_token (lexer) != T_SLASH)
 720             {
 721               if (lex_match_id (lexer, "CLUSTER"))
 722                 qc.print_cluster_membership = true;
 723               else if (lex_match_id (lexer, "INITIAL"))
 724                 qc.print_initial_clusters = true;
 725               else
 726                 {
 727                   lex_error (lexer, NULL);
 728                   goto error;
 729                 }
 730             }
 731         }
 732       else if (lex_match_id (lexer, "CRITERIA"))
 733         {
 734           lex_match (lexer, T_EQUALS);
 735           while (lex_token (lexer) != T_ENDCMD
 736                  && lex_token (lexer) != T_SLASH)
 737             {
 738               if (lex_match_id (lexer, "CLUSTERS"))
 739                 {
 740                   if (lex_force_match (lexer, T_LPAREN) &&
 741                       lex_force_int (lexer))
 742                     {
 743                       qc.ngroups = lex_integer (lexer);
 744                       if (qc.ngroups <= 0)
 745                         {
 746                           lex_error (lexer, _("The number of clusters must be positive"));
 747                           goto error;
 748                         }
 749                       lex_get (lexer);
 750                       if (!lex_force_match (lexer, T_RPAREN))
 751                         goto error;
 752                     }
 753                 }
 754               else if (lex_match_id (lexer, "CONVERGE"))
 755                 {
 756                   if (lex_force_match (lexer, T_LPAREN) &&
 757                       lex_force_num (lexer))
 758                     {
 759                       qc.epsilon = lex_number (lexer);
 760                       if (qc.epsilon <= 0)
 761                         {
 762                           lex_error (lexer, _("The convergence criterion must be positive"));
 763                           goto error;
 764                         }
 765                       lex_get (lexer);
 766                       if (!lex_force_match (lexer, T_RPAREN))
 767                         goto error;
 768                     }
 769                 }
 770               else if (lex_match_id (lexer, "MXITER"))
 771                 {
 772                   if (lex_force_match (lexer, T_LPAREN) &&
 773                       lex_force_int (lexer))
 774                     {
 775                       qc.maxiter = lex_integer (lexer);
 776                       if (qc.maxiter <= 0)
 777                         {
 778                           lex_error (lexer, _("The number of iterations must be positive"));
 779                           goto error;
 780                         }
 781                       lex_get (lexer);
 782                       if (!lex_force_match (lexer, T_RPAREN))
 783                         goto error;
 784                     }
 785                 }
 786               else if (lex_match_id (lexer, "NOINITIAL"))
 787                 {
 788                   qc.no_initial = true;
 789                 }
 790               else if (lex_match_id (lexer, "NOUPDATE"))
 791                 {
 792                   qc.no_update = true;
 793                 }
 794               else
 795                 {
 796                   lex_error (lexer, NULL);
 797                   goto error;
 798                 }
 799             }
 800         }
 801       else
 802         {
 803           lex_error (lexer, NULL);
 804           goto error;
 805         }
 806     }
 807
 808   qc.wv = dict_get_weight (dict);
 809
 810   {
 811     struct casereader *group;
 812     struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
 813
 814     while (casegrouper_get_next_group (grouper, &group))
 815       {
 816         if ( qc.missing_type == MISS_LISTWISE )
 817           {
 818             group = casereader_create_filter_missing (group, qc.vars, qc.n_vars,
 819                                                        qc.exclude,
 820                                                        NULL,  NULL);
 821           }
 822
 823         kmeans = kmeans_create (&qc);
 824         kmeans_cluster (kmeans, group, &qc);
 825         quick_cluster_show_results (kmeans, group, &qc);
 826         kmeans_destroy (kmeans);
 827         casereader_destroy (group);
 828       }
 829     ok = casegrouper_destroy (grouper);
 830   }
 831   ok = proc_commit (ds) && ok;
 832
 833   free (qc.vars);
 834
 835   return (ok);
 836
 837  error:
 838   free (qc.vars);
 839   return CMD_FAILURE;
 840 }