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