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/tab.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 tab_table *t;
 518   int nc, nr, currow;
 519   int i, j;
 520   nc = qc->ngroups + 1;
 521   nr = qc->n_vars + 4;
 522   t = tab_create (nc, nr);
 523   tab_headers (t, 0, nc - 1, 0, 1);
 524   currow = 0;
 525   if (!initial)
 526     {
 527       tab_title (t, _("Final Cluster Centers"));
 528     }
 529   else
 530     {
 531       tab_title (t, _("Initial Cluster Centers"));
 532     }
 533   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
 534   tab_joint_text (t, 1, 0, nc - 1, 0, TAB_CENTER, _("Cluster"));
 535   tab_hline (t, TAL_1, 1, nc - 1, 2);
 536   currow += 2;
 537
 538   for (i = 0; i < qc->ngroups; i++)
 539     {
 540       tab_text_format (t, (i + 1), currow, TAB_CENTER, "%d", (i + 1));
 541     }
 542   currow++;
 543   tab_hline (t, TAL_1, 1, nc - 1, currow);
 544   currow++;
 545   for (i = 0; i < qc->n_vars; i++)
 546     {
 547       tab_text (t, 0, currow + i, TAB_LEFT,
 548                 var_to_string (qc->vars[i]));
 549     }
 550
 551   for (i = 0; i < qc->ngroups; i++)
 552     {
 553       for (j = 0; j < qc->n_vars; j++)
 554         {
 555           if (!initial)
 556             {
 557               tab_double (t, i + 1, j + 4, TAB_CENTER,
 558                           gsl_matrix_get (kmeans->centers,
 559                                           kmeans->group_order->data[i], j),
 560                           var_get_print_format (qc->vars[j]), RC_OTHER);
 561             }
 562           else
 563             {
 564               tab_double (t, i + 1, j + 4, TAB_CENTER,
 565                           gsl_matrix_get (kmeans->initial_centers,
 566                                           kmeans->group_order->data[i], j),
 567                           var_get_print_format (qc->vars[j]), RC_OTHER);
 568             }
 569         }
 570     }
 571   tab_submit (t);
 572 }
 573
 574 /* Reports cluster membership for each case. */
 575 static void
 576 quick_cluster_show_membership (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
 577 {
 578   struct tab_table *t;
 579   int nc, nr, i;
 580
 581   struct ccase *c;
 582   struct casereader *cs = casereader_clone (reader);
 583   nc = 2;
 584   nr = kmeans->n + 1;
 585   t = tab_create (nc, nr);
 586   tab_headers (t, 0, nc - 1, 0, 0);
 587   tab_title (t, _("Cluster Membership"));
 588   tab_text (t, 0, 0, TAB_CENTER, _("Case Number"));
 589   tab_text (t, 1, 0, TAB_CENTER, _("Cluster"));
 590   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
 591   tab_hline (t, TAL_1, 0, nc - 1, 1);
 592
 593   gsl_permutation *ip = gsl_permutation_alloc (qc->ngroups);
 594   gsl_permutation_inverse (ip, kmeans->group_order);
 595
 596   for (i = 0; (c = casereader_read (cs)) != NULL; i++, case_unref (c))
 597     {
 598       int clust = -1;
 599       assert (i < kmeans->n);
 600       kmeans_get_nearest_group (kmeans, c, qc, &clust, NULL, NULL, NULL);
 601       clust = ip->data[clust];
 602       tab_text_format (t, 0, i+1, TAB_CENTER, "%d", (i + 1));
 603       tab_text_format (t, 1, i+1, TAB_CENTER, "%d", (clust + 1));
 604     }
 605   gsl_permutation_free (ip);
 606   assert (i == kmeans->n);
 607   tab_submit (t);
 608   casereader_destroy (cs);
 609 }
 610
 611
 612 /* Reports number of cases of each single cluster. */
 613 static void
 614 quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
 615 {
 616   struct tab_table *t;
 617   int nc, nr;
 618   int i, numelem;
 619   long int total;
 620   nc = 3;
 621   nr = qc->ngroups + 1;
 622   t = tab_create (nc, nr);
 623   tab_headers (t, 0, nc - 1, 0, 0);
 624   tab_title (t, _("Number of Cases in each Cluster"));
 625   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
 626   tab_text (t, 0, 0, TAB_LEFT, _("Cluster"));
 627
 628   total = 0;
 629   for (i = 0; i < qc->ngroups; i++)
 630     {
 631       tab_text_format (t, 1, i, TAB_CENTER, "%d", (i + 1));
 632       numelem =
 633         kmeans->num_elements_groups->data[kmeans->group_order->data[i]];
 634       tab_text_format (t, 2, i, TAB_CENTER, "%d", numelem);
 635       total += numelem;
 636     }
 637
 638   tab_text (t, 0, qc->ngroups, TAB_LEFT, _("Valid"));
 639   tab_text_format (t, 2, qc->ngroups, TAB_LEFT, "%ld", total);
 640   tab_submit (t);
 641 }
 642
 643 /* Reports. */
 644 static void
 645 quick_cluster_show_results (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
 646 {
 647   kmeans_order_groups (kmeans, qc); /* what does this do? */
 648
 649   if( qc->print_initial_clusters )
 650     quick_cluster_show_centers (kmeans, true, qc);
 651   quick_cluster_show_centers (kmeans, false, qc);
 652   quick_cluster_show_number_cases (kmeans, qc);
 653   if( qc->print_cluster_membership )
 654     quick_cluster_show_membership(kmeans, reader, qc);
 655 }
 656
 657 int
 658 cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
 659 {
 660   struct qc qc;
 661   struct Kmeans *kmeans;
 662   bool ok;
 663   const struct dictionary *dict = dataset_dict (ds);
 664   qc.ngroups = 2;
 665   qc.maxiter = 10;
 666   qc.epsilon = DBL_EPSILON;
 667   qc.missing_type = MISS_LISTWISE;
 668   qc.exclude = MV_ANY;
 669   qc.print_cluster_membership = false; /* default = do not output case cluster membership */
 670   qc.print_initial_clusters = false;   /* default = do not print initial clusters */
 671   qc.no_initial = false;               /* default = use well separated initial clusters */
 672   qc.no_update = false;               /* default = iterate until convergence or max iterations */
 673
 674   if (!parse_variables_const (lexer, dict, &qc.vars, &qc.n_vars,
 675                               PV_NO_DUPLICATE | PV_NUMERIC))
 676     {
 677       return (CMD_FAILURE);
 678     }
 679
 680   while (lex_token (lexer) != T_ENDCMD)
 681     {
 682       lex_match (lexer, T_SLASH);
 683
 684       if (lex_match_id (lexer, "MISSING"))
 685         {
 686           lex_match (lexer, T_EQUALS);
 687           while (lex_token (lexer) != T_ENDCMD
 688                  && lex_token (lexer) != T_SLASH)
 689             {
 690               if (lex_match_id (lexer, "LISTWISE") || lex_match_id (lexer, "DEFAULT"))
 691                 {
 692                   qc.missing_type = MISS_LISTWISE;
 693                 }
 694               else if (lex_match_id (lexer, "PAIRWISE"))
 695                 {
 696                   qc.missing_type = MISS_PAIRWISE;
 697                 }
 698               else if (lex_match_id (lexer, "INCLUDE"))
 699                 {
 700                   qc.exclude = MV_SYSTEM;
 701                 }
 702               else if (lex_match_id (lexer, "EXCLUDE"))
 703                 {
 704                   qc.exclude = MV_ANY;
 705                 }
 706               else
 707                 {
 708                   lex_error (lexer, NULL);
 709                   goto error;
 710                 }
 711             }
 712         }
 713       else if (lex_match_id (lexer, "PRINT"))
 714         {
 715           lex_match (lexer, T_EQUALS);
 716           while (lex_token (lexer) != T_ENDCMD
 717                  && lex_token (lexer) != T_SLASH)
 718             {
 719               if (lex_match_id (lexer, "CLUSTER"))
 720                 qc.print_cluster_membership = true;
 721               else if (lex_match_id (lexer, "INITIAL"))
 722                 qc.print_initial_clusters = true;
 723               else
 724                 {
 725                   lex_error (lexer, NULL);
 726                   goto error;
 727                 }
 728             }
 729         }
 730       else if (lex_match_id (lexer, "CRITERIA"))
 731         {
 732           lex_match (lexer, T_EQUALS);
 733           while (lex_token (lexer) != T_ENDCMD
 734                  && lex_token (lexer) != T_SLASH)
 735             {
 736               if (lex_match_id (lexer, "CLUSTERS"))
 737                 {
 738                   if (lex_force_match (lexer, T_LPAREN) &&
 739                       lex_force_int (lexer))
 740                     {
 741                       qc.ngroups = lex_integer (lexer);
 742                       if (qc.ngroups <= 0)
 743                         {
 744                           lex_error (lexer, _("The number of clusters must be positive"));
 745                           goto error;
 746                         }
 747                       lex_get (lexer);
 748                       if (!lex_force_match (lexer, T_RPAREN))
 749                         goto error;
 750                     }
 751                 }
 752               else if (lex_match_id (lexer, "CONVERGE"))
 753                 {
 754                   if (lex_force_match (lexer, T_LPAREN) &&
 755                       lex_force_num (lexer))
 756                     {
 757                       qc.epsilon = lex_number (lexer);
 758                       if (qc.epsilon <= 0)
 759                         {
 760                           lex_error (lexer, _("The convergence criterium must be positive"));
 761                           goto error;
 762                         }
 763                       lex_get (lexer);
 764                       if (!lex_force_match (lexer, T_RPAREN))
 765                         goto error;
 766                     }
 767                 }
 768               else if (lex_match_id (lexer, "MXITER"))
 769                 {
 770                   if (lex_force_match (lexer, T_LPAREN) &&
 771                       lex_force_int (lexer))
 772                     {
 773                       qc.maxiter = lex_integer (lexer);
 774                       if (qc.maxiter <= 0)
 775                         {
 776                           lex_error (lexer, _("The number of iterations must be positive"));
 777                           goto error;
 778                         }
 779                       lex_get (lexer);
 780                       if (!lex_force_match (lexer, T_RPAREN))
 781                         goto error;
 782                     }
 783                 }
 784               else if (lex_match_id (lexer, "NOINITIAL"))
 785                 {
 786                   qc.no_initial = true;
 787                 }
 788               else if (lex_match_id (lexer, "NOUPDATE"))
 789                 {
 790                   qc.no_update = true;
 791                 }
 792               else
 793                 {
 794                   lex_error (lexer, NULL);
 795                   goto error;
 796                 }
 797             }
 798         }
 799       else
 800         {
 801           lex_error (lexer, NULL);
 802           goto error;
 803         }
 804     }
 805
 806   qc.wv = dict_get_weight (dict);
 807
 808   {
 809     struct casereader *group;
 810     struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
 811
 812     while (casegrouper_get_next_group (grouper, &group))
 813       {
 814         if ( qc.missing_type == MISS_LISTWISE )
 815           {
 816             group  = casereader_create_filter_missing (group, qc.vars, qc.n_vars,
 817                                                        qc.exclude,
 818                                                        NULL,  NULL);
 819           }
 820
 821         kmeans = kmeans_create (&qc);
 822         kmeans_cluster (kmeans, group, &qc);
 823         quick_cluster_show_results (kmeans, group, &qc);
 824         kmeans_destroy (kmeans);
 825         casereader_destroy (group);
 826       }
 827     ok = casegrouper_destroy (grouper);
 828   }
 829   ok = proc_commit (ds) && ok;
 830
 831   free (qc.vars);
 832
 833   return (ok);
 834
 835  error:
 836   free (qc.vars);
 837   return CMD_FAILURE;
 838 }