Fix some typos (found by codespell)
[pspp] / src / language / stats / quick-cluster.c
index 9adcc64243a70c50539909304f758d09488a4020..b171951557f95c1a464aacc9e1dfcfabc8576dd6 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2011 Free Software Foundation, Inc.
+   Copyright (C) 2011, 2012, 2015 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -20,7 +20,6 @@
 #include <gsl/gsl_permutation.h>
 #include <gsl/gsl_sort_vector.h>
 #include <gsl/gsl_statistics.h>
-#include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 
@@ -37,6 +36,7 @@
 #include "language/lexer/variable-parser.h"
 #include "libpspp/message.h"
 #include "libpspp/misc.h"
+#include "libpspp/assertion.h"
 #include "libpspp/str.h"
 #include "math/random.h"
 #include "output/tab.h"
 #define _(msgid) gettext (msgid)
 #define N_(msgid) msgid
 
+enum missing_type
+  {
+    MISS_LISTWISE,
+    MISS_PAIRWISE,
+  };
+
+
 struct qc
 {
   const struct variable **vars;
   size_t n_vars;
 
+  double epsilon;               /* The convergence criterium */
+
   int ngroups;                 /* Number of group. (Given by the user) */
   int maxiter;                 /* Maximum iterations (Given by the user) */
+  bool print_cluster_membership; /* true => print membership */
+  bool print_initial_clusters;   /* true => print initial cluster */
+  bool no_initial;              /* true => simplified initial cluster selection */
+  bool no_update;               /* true => do not iterate  */
 
   const struct variable *wv;   /* Weighting variable. */
+
+  enum missing_type missing_type;
+  enum mv_class exclude;
 };
 
 /* Holds all of the information for the functions.  int n, holds the number of
@@ -63,30 +79,19 @@ struct qc
 struct Kmeans
 {
   gsl_matrix *centers;         /* Centers for groups. */
-  gsl_vector_long *num_elements_groups;
+  gsl_matrix *updated_centers;
+  casenumber n;
 
-  casenumber n;                        /* Number of observations (default -1). */
+  gsl_vector_long *num_elements_groups;
 
-  int lastiter;                        /* Iteration where it found the solution. */
-  int trials;                  /* If not convergence, how many times has
-                                   clustering done. */
   gsl_matrix *initial_centers; /* Initial random centers. */
-
+  double convergence_criteria;
   gsl_permutation *group_order;        /* Group order for reporting. */
-  struct caseproto *proto;
-  struct casereader *index_rdr;        /* Group ids for each case. */
 };
 
 static struct Kmeans *kmeans_create (const struct qc *qc);
 
-static void kmeans_randomize_centers (struct Kmeans *kmeans, const struct qc *qc);
-
-static int kmeans_get_nearest_group (struct Kmeans *kmeans, struct ccase *c, const struct qc *);
-
-static void kmeans_recalculate_centers (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *);
-
-static int
-kmeans_calculate_indexes_and_check_convergence (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *);
+static void kmeans_get_nearest_group (const struct Kmeans *kmeans, struct ccase *c, const struct qc *, int *, double *, int *, double *);
 
 static void kmeans_order_groups (struct Kmeans *kmeans, const struct qc *);
 
@@ -94,9 +99,11 @@ static void kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader, co
 
 static void quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *);
 
+static void quick_cluster_show_membership (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *);
+
 static void quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *);
 
-static void quick_cluster_show_results (struct Kmeans *kmeans, const struct qc *);
+static void quick_cluster_show_results (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *);
 
 int cmd_quick_cluster (struct lexer *lexer, struct dataset *ds);
 
@@ -110,16 +117,11 @@ kmeans_create (const struct qc *qc)
 {
   struct Kmeans *kmeans = xmalloc (sizeof (struct Kmeans));
   kmeans->centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
+  kmeans->updated_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
   kmeans->num_elements_groups = gsl_vector_long_alloc (qc->ngroups);
-  kmeans->n = 0;
-  kmeans->lastiter = 0;
-  kmeans->trials = 0;
   kmeans->group_order = gsl_permutation_alloc (kmeans->centers->size1);
   kmeans->initial_centers = NULL;
 
-  kmeans->proto = caseproto_create ();
-  kmeans->proto = caseproto_add_width (kmeans->proto, 0);
-  kmeans->index_rdr = NULL;
   return (kmeans);
 }
 
@@ -127,196 +129,243 @@ static void
 kmeans_destroy (struct Kmeans *kmeans)
 {
   gsl_matrix_free (kmeans->centers);
+  gsl_matrix_free (kmeans->updated_centers);
   gsl_matrix_free (kmeans->initial_centers);
 
   gsl_vector_long_free (kmeans->num_elements_groups);
 
   gsl_permutation_free (kmeans->group_order);
 
-  caseproto_unref (kmeans->proto);
+  free (kmeans);
+}
 
-  casereader_destroy (kmeans->index_rdr);
+static double
+diff_matrix (const gsl_matrix *m1, const gsl_matrix *m2)
+{
+  int i,j;
+  double max_diff = -INFINITY;
+  for (i = 0; i < m1->size1; ++i)
+    {
+      double diff = 0;
+      for (j = 0; j < m1->size2; ++j)
+       {
+         diff += pow2 (gsl_matrix_get (m1,i,j) - gsl_matrix_get (m2,i,j) );
+       }
+      if (diff > max_diff)
+       max_diff = diff;
+    }
 
-  free (kmeans);
+  return max_diff;
 }
 
-/* Creates random centers using randomly selected cases from the data. */
-static void
-kmeans_randomize_centers (struct Kmeans *kmeans, const struct qc *qc)
+
+
+static double
+matrix_mindist (const gsl_matrix *m, int *mn, int *mm)
 {
   int i, j;
-  for (i = 0; i < qc->ngroups; i++)
+  double mindist = INFINITY;
+  for (i = 0; i < m->size1 - 1; ++i)
     {
-      for (j = 0; j < qc->n_vars; j++)
+      for (j = i + 1; j < m->size1; ++j)
        {
-         if (i == j)
+         int k;
+         double diff_sq = 0;
+         for (k = 0; k < m->size2; ++k)
            {
-             gsl_matrix_set (kmeans->centers, i, j, 1);
+             diff_sq += pow2 (gsl_matrix_get (m, j, k) - gsl_matrix_get (m, i, k));
            }
-         else
+         if (diff_sq < mindist)
            {
-             gsl_matrix_set (kmeans->centers, i, j, 0);
+             mindist = diff_sq;
+             if (mn)
+               *mn = i;
+             if (mm)
+               *mm = j;
            }
        }
     }
-  /* If it is the first iteration, the variable kmeans->initial_centers is NULL
-     and it is created once for reporting issues. In SPSS, initial centers are
-     shown in the reports but in PSPP it is not shown now. I am leaving it
-     here. */
-  if (!kmeans->initial_centers)
+
+  return mindist;
+}
+
+
+/* Return the distance of C from the group whose index is WHICH */
+static double
+dist_from_case (const struct Kmeans *kmeans, const struct ccase *c, const struct qc *qc, int which)
+{
+  int j;
+  double dist = 0;
+  for (j = 0; j < qc->n_vars; j++)
     {
-      kmeans->initial_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
-      gsl_matrix_memcpy (kmeans->initial_centers, kmeans->centers);
+      const union value *val = case_data (c, qc->vars[j]);
+      if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+       NOT_REACHED ();
+
+      dist += pow2 (gsl_matrix_get (kmeans->centers, which, j) - val->f);
     }
+
+  return dist;
 }
 
-static int
-kmeans_get_nearest_group (struct Kmeans *kmeans, struct ccase *c, const struct qc *qc)
+/* Return the minimum distance of the group WHICH and all other groups */
+static double
+min_dist_from (const struct Kmeans *kmeans, const struct qc *qc, int which)
 {
-  int result = -1;
-  double x;
-  int i, j;
-  double dist;
-  double mindist;
-  mindist = INFINITY;
+  int j, i;
+
+  double mindist = INFINITY;
   for (i = 0; i < qc->ngroups; i++)
     {
-      dist = 0;
+      if (i == which)
+       continue;
+
+      double dist = 0;
       for (j = 0; j < qc->n_vars; j++)
        {
-         x = case_data (c, qc->vars[j])->f;
-         dist += pow2 (gsl_matrix_get (kmeans->centers, i, j) - x);
+         dist += pow2 (gsl_matrix_get (kmeans->centers, i, j) - gsl_matrix_get (kmeans->centers, which, j));
        }
+
       if (dist < mindist)
        {
          mindist = dist;
-         result = i;
        }
     }
-  return (result);
+
+  return mindist;
 }
 
-/* Re-calculate the cluster centers. */
+
+
+/* Calculate the initial cluster centers. */
 static void
-kmeans_recalculate_centers (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
+kmeans_initial_centers (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
 {
-  casenumber i;
-  int v, j;
-  double x, curval;
   struct ccase *c;
-  struct ccase *c_index;
-  struct casereader *cs;
-  struct casereader *cs_index;
-  int index;
-
-  i = 0;
-  cs = casereader_clone (reader);
-  cs_index = casereader_clone (kmeans->index_rdr);
+  int nc = 0, j;
 
-  gsl_matrix_set_all (kmeans->centers, 0.0);
+  struct casereader *cs = casereader_clone (reader);
   for (; (c = casereader_read (cs)) != NULL; case_unref (c))
     {
-      double weight = qc->wv ? case_data (c, qc->wv)->f : 1.0;
-      c_index = casereader_read (cs_index);
-      index = case_data_idx (c_index, 0)->f;
-      for (v = 0; v < qc->n_vars; ++v)
+      bool missing = false;
+      for (j = 0; j < qc->n_vars; ++j)
        {
-         x = case_data (c, qc->vars[v])->f * weight;
-         curval = gsl_matrix_get (kmeans->centers, index, v);
-         gsl_matrix_set (kmeans->centers, index, v, curval + x);
+         const union value *val = case_data (c, qc->vars[j]);
+         if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+           {
+             missing = true;
+             break;
+           }
+
+         if (nc < qc->ngroups)
+           gsl_matrix_set (kmeans->centers, nc, j, val->f);
        }
-      i++;
-      case_unref (c_index);
-    }
-  casereader_destroy (cs);
-  casereader_destroy (cs_index);
 
-  /* Getting number of cases */
-  if (kmeans->n == 0)
-    kmeans->n = i;
+      if (missing)
+       continue;
 
-  /* We got sum of each center but we need averages.
-     We are dividing centers to numobs. This may be inefficient and
-     we should check it again. */
-  for (i = 0; i < qc->ngroups; i++)
-    {
-      casenumber numobs = kmeans->num_elements_groups->data[i];
-      for (j = 0; j < qc->n_vars; j++)
+      if (nc++ < qc->ngroups)
+       continue;
+
+      if (!qc->no_initial)
        {
-         if (numobs > 0)
+         int mq, mp;
+         double delta;
+
+         int mn, mm;
+         double m = matrix_mindist (kmeans->centers, &mn, &mm);
+
+         kmeans_get_nearest_group (kmeans, c, qc, &mq, &delta, &mp, NULL);
+         if (delta > m)
+           /* If the distance between C and the nearest group, is greater than the distance
+              between the two  groups which are clostest to each other, then one group must be replaced */
            {
-             double *x = gsl_matrix_ptr (kmeans->centers, i, j);
-             *x /= numobs;
+             /* Out of mn and mm, which is the clostest of the two groups to C ? */
+             int which = (dist_from_case (kmeans, c, qc, mn) > dist_from_case (kmeans, c, qc, mm)) ? mm : mn;
+
+             for (j = 0; j < qc->n_vars; ++j)
+               {
+                 const union value *val = case_data (c, qc->vars[j]);
+                 gsl_matrix_set (kmeans->centers, which, j, val->f);
+               }
            }
-         else
+         else if (dist_from_case (kmeans, c, qc, mp) > min_dist_from (kmeans, qc, mq))
+           /* If the distance between C and the second nearest group (MP) is greater than the
+              smallest distance between the nearest group (MQ) and any other group, then replace
+              MQ with C */
            {
-             gsl_matrix_set (kmeans->centers, i, j, 0);
+             for (j = 0; j < qc->n_vars; ++j)
+               {
+                 const union value *val = case_data (c, qc->vars[j]);
+                 gsl_matrix_set (kmeans->centers, mq, j, val->f);
+               }
            }
        }
     }
-}
 
-/* The variable index in struct Kmeans holds integer values that represents the
-   current groups of cases.  index[n]=a shows the nth case is belong to ath
-   cluster.  This function calculates these indexes and returns the number of
-   different cases of the new and old index variables.  If last two index
-   variables are equal, there is no any enhancement of clustering. */
-static int
-kmeans_calculate_indexes_and_check_convergence (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
-{
-  int totaldiff = 0;
-  struct ccase *c;
-  struct casereader *cs = casereader_clone (reader);
+  casereader_destroy (cs);
 
-  /* A casewriter into which we will write the indexes. */
-  struct casewriter *index_wtr = autopaging_writer_create (kmeans->proto);
+  kmeans->convergence_criteria = qc->epsilon * matrix_mindist (kmeans->centers, NULL, NULL);
 
-  gsl_vector_long_set_all (kmeans->num_elements_groups, 0);
+  /* As it is the first iteration, the variable kmeans->initial_centers is NULL
+     and it is created once for reporting issues. */
+  kmeans->initial_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
+  gsl_matrix_memcpy (kmeans->initial_centers, kmeans->centers);
+}
 
-  for (; (c = casereader_read (cs)) != NULL; case_unref (c))
+
+/* Return the index of the group which is nearest to the case C */
+static void
+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)
+{
+  int result0 = -1;
+  int result1 = -1;
+  int i, j;
+  double mindist0 = INFINITY;
+  double mindist1 = INFINITY;
+  for (i = 0; i < qc->ngroups; i++)
     {
-      /* A case to hold the new index. */
-      struct ccase *index_case_new = case_create (kmeans->proto);
-      int bestindex = kmeans_get_nearest_group (kmeans, c, qc);
-      double weight = qc->wv ? case_data (c, qc->wv)->f : 1.0;
-      kmeans->num_elements_groups->data[bestindex] += weight;
-      if (kmeans->index_rdr)
+      double dist = 0;
+      for (j = 0; j < qc->n_vars; j++)
        {
-         /* A case from which the old index will be read. */
-         struct ccase *index_case_old = NULL;
+         const union value *val = case_data (c, qc->vars[j]);
+         if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+           continue;
 
-         /* Read the case from the index casereader. */
-         index_case_old = casereader_read (kmeans->index_rdr);
+         dist += pow2 (gsl_matrix_get (kmeans->centers, i, j) - val->f);
+       }
 
-         /* Set totaldiff, using the old_index. */
-         totaldiff += abs (case_data_idx (index_case_old, 0)->f - bestindex);
+      if (dist < mindist0)
+       {
+         mindist1 = mindist0;
+         result1 = result0;
 
-         /* We have no use for the old case anymore, so unref it. */
-         case_unref (index_case_old);
+         mindist0 = dist;
+         result0 = i;
        }
-      else
+      else if (dist < mindist1)
        {
-         /* If this is the first run, then assume index is zero. */
-         totaldiff += bestindex;
+         mindist1 = dist;
+         result1 = i;
        }
+    }
 
-      /* Set the value of the new inde.x */
-      case_data_rw_idx (index_case_new, 0)->f = bestindex;
+  if (delta_q)
+    *delta_q = mindist0;
 
-      /* and write the new index to the casewriter */
-      casewriter_write (index_wtr, index_case_new);
-    }
-  casereader_destroy (cs);
-  /* We have now read through the entire index_rdr, so it's of no use
-     anymore. */
-  casereader_destroy (kmeans->index_rdr);
+  if (g_q)
+    *g_q = result0;
 
-  /* Convert the writer into a reader, ready for the next iteration to read */
-  kmeans->index_rdr = casewriter_make_reader (index_wtr);
 
-  return (totaldiff);
+  if (delta_p)
+    *delta_p = mindist1;
+
+  if (g_p)
+    *g_p = result1;
 }
 
+
+
 static void
 kmeans_order_groups (struct Kmeans *kmeans, const struct qc *qc)
 {
@@ -331,44 +380,131 @@ kmeans_order_groups (struct Kmeans *kmeans, const struct qc *qc)
 static void
 kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader, const struct qc *qc)
 {
-  int i;
-  bool redo;
-  int diffs;
-  bool show_warning1;
-
-  show_warning1 = true;
-cluster:
-  redo = false;
-  kmeans_randomize_centers (kmeans, qc);
-  for (kmeans->lastiter = 0; kmeans->lastiter < qc->maxiter;
-       kmeans->lastiter++)
+  int j;
+
+  kmeans_initial_centers (kmeans, reader, qc);
+
+  gsl_matrix_memcpy (kmeans->updated_centers, kmeans->centers);
+
+
+  for (int xx = 0 ; xx < qc->maxiter ; ++xx)
     {
-      diffs = kmeans_calculate_indexes_and_check_convergence (kmeans, reader, qc);
-      kmeans_recalculate_centers (kmeans, reader, qc);
-      if (show_warning1 && qc->ngroups > kmeans->n)
+      gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
+
+      kmeans->n = 0;
+      if (!qc->no_update)
        {
-         msg (MW, _("Number of clusters may not be larger than the number "
-                     "of cases."));
-         show_warning1 = false;
+         struct casereader *r = casereader_clone (reader);
+         struct ccase *c;
+         for (; (c = casereader_read (r)) != NULL; case_unref (c))
+           {
+             int group = -1;
+             int g;
+             bool missing = false;
+
+             for (j = 0; j < qc->n_vars; j++)
+               {
+                 const union value *val = case_data (c, qc->vars[j]);
+                 if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+                   missing = true;
+               }
+
+             if (missing)
+               continue;
+
+             double mindist = INFINITY;
+             for (g = 0; g < qc->ngroups; ++g)
+               {
+                 double d = dist_from_case (kmeans, c, qc, g);
+
+                 if (d < mindist)
+                   {
+                     mindist = d;
+                     group = g;
+                   }
+               }
+
+             long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
+             *n += qc->wv ? case_data (c, qc->wv)->f : 1.0;
+             kmeans->n++;
+
+             for (j = 0; j < qc->n_vars; ++j)
+               {
+                 const union value *val = case_data (c, qc->vars[j]);
+                 if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+                   continue;
+                 double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
+                 *x += val->f * (qc->wv ? case_data (c, qc->wv)->f : 1.0);
+               }
+           }
+
+         casereader_destroy (r);
        }
-      if (diffs == 0)
-       break;
-    }
 
-  for (i = 0; i < qc->ngroups; i++)
-    {
-      if (kmeans->num_elements_groups->data[i] == 0)
+      int g;
+
+      /* Divide the cluster sums by the number of items in each cluster */
+      for (g = 0; g < qc->ngroups; ++g)
        {
-         kmeans->trials++;
-         if (kmeans->trials >= 3)
-           break;
-         redo = true;
-         break;
+         for (j = 0; j < qc->n_vars; ++j)
+           {
+             long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
+             double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
+             *x /= n + 1;  // Plus 1 for the initial centers
+           }
        }
-    }
-  if (redo)
-    goto cluster;
 
+
+      gsl_matrix_memcpy (kmeans->centers, kmeans->updated_centers);
+
+      {
+       kmeans->n = 0;
+       /* Step 3 */
+       gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
+       gsl_matrix_set_all (kmeans->updated_centers, 0.0);
+       struct ccase *c;
+       struct casereader *cs = casereader_clone (reader);
+       for (; (c = casereader_read (cs)) != NULL; case_unref (c))
+         {
+           int group = -1;
+           kmeans_get_nearest_group (kmeans, c, qc, &group, NULL, NULL, NULL);
+
+           for (j = 0; j < qc->n_vars; ++j)
+             {
+               const union value *val = case_data (c, qc->vars[j]);
+               if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
+                 continue;
+
+               double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
+               *x += val->f;
+             }
+
+           long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
+           *n += qc->wv ? case_data (c, qc->wv)->f : 1.0;
+           kmeans->n++;
+         }
+       casereader_destroy (cs);
+
+
+       /* Divide the cluster sums by the number of items in each cluster */
+       for (g = 0; g < qc->ngroups; ++g)
+         {
+           for (j = 0; j < qc->n_vars; ++j)
+             {
+               long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
+               double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
+               *x /= n ;
+             }
+         }
+
+       double d = diff_matrix (kmeans->updated_centers, kmeans->centers);
+       if (d < kmeans->convergence_criteria)
+         break;
+      }
+
+      if (qc->no_update)
+       break;
+    }
 }
 
 /* Reports centers of clusters.
@@ -379,11 +515,10 @@ static void
 quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *qc)
 {
   struct tab_table *t;
-  int nc, nr, heading_columns, currow;
+  int nc, nr, currow;
   int i, j;
   nc = qc->ngroups + 1;
   nr = qc->n_vars + 4;
-  heading_columns = 1;
   t = tab_create (nc, nr);
   tab_headers (t, 0, nc - 1, 0, 1);
   currow = 0;
@@ -422,20 +557,58 @@ quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc
              tab_double (t, i + 1, j + 4, TAB_CENTER,
                          gsl_matrix_get (kmeans->centers,
                                          kmeans->group_order->data[i], j),
-                         var_get_print_format (qc->vars[j]));
+                         var_get_print_format (qc->vars[j]), RC_OTHER);
            }
          else
            {
              tab_double (t, i + 1, j + 4, TAB_CENTER,
                          gsl_matrix_get (kmeans->initial_centers,
                                          kmeans->group_order->data[i], j),
-                         var_get_print_format (qc->vars[j]));
+                         var_get_print_format (qc->vars[j]), RC_OTHER);
            }
        }
     }
   tab_submit (t);
 }
 
+/* Reports cluster membership for each case. */
+static void
+quick_cluster_show_membership (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
+{
+  struct tab_table *t;
+  int nc, nr, i;
+
+  struct ccase *c;
+  struct casereader *cs = casereader_clone (reader);
+  nc = 2;
+  nr = kmeans->n + 1;
+  t = tab_create (nc, nr);
+  tab_headers (t, 0, nc - 1, 0, 0);
+  tab_title (t, _("Cluster Membership"));
+  tab_text (t, 0, 0, TAB_CENTER, _("Case Number"));
+  tab_text (t, 1, 0, TAB_CENTER, _("Cluster"));
+  tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
+  tab_hline (t, TAL_1, 0, nc - 1, 1);
+
+  gsl_permutation *ip = gsl_permutation_alloc (qc->ngroups);
+  gsl_permutation_inverse (ip, kmeans->group_order);
+
+  for (i = 0; (c = casereader_read (cs)) != NULL; i++, case_unref (c))
+    {
+      int clust = -1;
+      assert (i < kmeans->n);
+      kmeans_get_nearest_group (kmeans, c, qc, &clust, NULL, NULL, NULL);
+      clust = ip->data[clust];
+      tab_text_format (t, 0, i+1, TAB_CENTER, "%d", (i + 1));
+      tab_text_format (t, 1, i+1, TAB_CENTER, "%d", (clust + 1));
+    }
+  gsl_permutation_free (ip);
+  assert (i == kmeans->n);
+  tab_submit (t);
+  casereader_destroy (cs);
+}
+
+
 /* Reports number of cases of each single cluster. */
 static void
 quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
@@ -469,13 +642,16 @@ quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
 
 /* Reports. */
 static void
-quick_cluster_show_results (struct Kmeans *kmeans, const struct qc *qc)
+quick_cluster_show_results (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
 {
-  kmeans_order_groups (kmeans, qc);
-  /* Uncomment the line below for reporting initial centers. */
-  /* quick_cluster_show_centers (kmeans, true); */
+  kmeans_order_groups (kmeans, qc); /* what does this do? */
+
+  if( qc->print_initial_clusters )
+    quick_cluster_show_centers (kmeans, true, qc);
   quick_cluster_show_centers (kmeans, false, qc);
   quick_cluster_show_number_cases (kmeans, qc);
+  if( qc->print_cluster_membership )
+    quick_cluster_show_membership(kmeans, reader, qc);
 }
 
 int
@@ -486,7 +662,14 @@ cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
   bool ok;
   const struct dictionary *dict = dataset_dict (ds);
   qc.ngroups = 2;
-  qc.maxiter = 2;
+  qc.maxiter = 10;
+  qc.epsilon = DBL_EPSILON;
+  qc.missing_type = MISS_LISTWISE;
+  qc.exclude = MV_ANY;
+  qc.print_cluster_membership = false; /* default = do not output case cluster membership */
+  qc.print_initial_clusters = false;   /* default = do not print initial clusters */
+  qc.no_initial = false;               /* default = use well separated initial clusters */
+  qc.no_update = false;               /* default = iterate until convergence or max iterations */
 
   if (!parse_variables_const (lexer, dict, &qc.vars, &qc.n_vars,
                              PV_NO_DUPLICATE | PV_NUMERIC))
@@ -494,9 +677,57 @@ cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
       return (CMD_FAILURE);
     }
 
-  if (lex_match (lexer, T_SLASH))
+  while (lex_token (lexer) != T_ENDCMD)
     {
-      if (lex_match_id (lexer, "CRITERIA"))
+      lex_match (lexer, T_SLASH);
+
+      if (lex_match_id (lexer, "MISSING"))
+       {
+         lex_match (lexer, T_EQUALS);
+         while (lex_token (lexer) != T_ENDCMD
+                && lex_token (lexer) != T_SLASH)
+           {
+             if (lex_match_id (lexer, "LISTWISE") || lex_match_id (lexer, "DEFAULT"))
+               {
+                 qc.missing_type = MISS_LISTWISE;
+               }
+             else if (lex_match_id (lexer, "PAIRWISE"))
+               {
+                 qc.missing_type = MISS_PAIRWISE;
+               }
+             else if (lex_match_id (lexer, "INCLUDE"))
+               {
+                 qc.exclude = MV_SYSTEM;
+               }
+             else if (lex_match_id (lexer, "EXCLUDE"))
+               {
+                 qc.exclude = MV_ANY;
+               }
+             else
+               {
+                 lex_error (lexer, NULL);
+                 goto error;
+               }
+           }
+       }
+      else if (lex_match_id (lexer, "PRINT"))
+       {
+         lex_match (lexer, T_EQUALS);
+         while (lex_token (lexer) != T_ENDCMD
+                && lex_token (lexer) != T_SLASH)
+           {
+             if (lex_match_id (lexer, "CLUSTER"))
+                qc.print_cluster_membership = true;
+             else if (lex_match_id (lexer, "INITIAL"))
+               qc.print_initial_clusters = true;
+             else
+               {
+                 lex_error (lexer, NULL);
+                 goto error;
+               }
+           }
+       }
+      else if (lex_match_id (lexer, "CRITERIA"))
        {
          lex_match (lexer, T_EQUALS);
          while (lex_token (lexer) != T_ENDCMD
@@ -504,28 +735,72 @@ cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
            {
              if (lex_match_id (lexer, "CLUSTERS"))
                {
-                 if (lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN) &&
+                     lex_force_int (lexer))
                    {
-                     lex_force_int (lexer);
                      qc.ngroups = lex_integer (lexer);
+                     if (qc.ngroups <= 0)
+                       {
+                         lex_error (lexer, _("The number of clusters must be positive"));
+                         goto error;
+                       }
+                     lex_get (lexer);
+                     if (!lex_force_match (lexer, T_RPAREN))
+                       goto error;
+                   }
+               }
+             else if (lex_match_id (lexer, "CONVERGE"))
+               {
+                 if (lex_force_match (lexer, T_LPAREN) &&
+                     lex_force_num (lexer))
+                   {
+                     qc.epsilon = lex_number (lexer);
+                     if (qc.epsilon <= 0)
+                       {
+                         lex_error (lexer, _("The convergence criterium must be positive"));
+                         goto error;
+                       }
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (!lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
              else if (lex_match_id (lexer, "MXITER"))
                {
-                 if (lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN) &&
+                     lex_force_int (lexer))
                    {
-                     lex_force_int (lexer);
                      qc.maxiter = lex_integer (lexer);
+                     if (qc.maxiter <= 0)
+                       {
+                         lex_error (lexer, _("The number of iterations must be positive"));
+                         goto error;
+                       }
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (!lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
+             else if (lex_match_id (lexer, "NOINITIAL"))
+               {
+                 qc.no_initial = true;
+               }
+             else if (lex_match_id (lexer, "NOUPDATE"))
+               {
+                 qc.no_update = true;
+               }
              else
-                goto error;
+               {
+                 lex_error (lexer, NULL);
+                 goto error;
+               }
            }
        }
+      else
+        {
+          lex_error (lexer, NULL);
+          goto error;
+        }
     }
 
   qc.wv = dict_get_weight (dict);
@@ -536,9 +811,16 @@ cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
 
     while (casegrouper_get_next_group (grouper, &group))
       {
+       if ( qc.missing_type == MISS_LISTWISE )
+         {
+           group  = casereader_create_filter_missing (group, qc.vars, qc.n_vars,
+                                                      qc.exclude,
+                                                      NULL,  NULL);
+         }
+
        kmeans = kmeans_create (&qc);
        kmeans_cluster (kmeans, group, &qc);
-       quick_cluster_show_results (kmeans, &qc);
+       quick_cluster_show_results (kmeans, group, &qc);
        kmeans_destroy (kmeans);
        casereader_destroy (group);
       }