X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fquick-cluster.c;h=0cd08d02a9268224436bac93339ffc876c9669a7;hb=88cf14d5e42bb262b245098114bb6eb8f83fc87d;hp=014406f098aa53adf20ef12136e30d800b751cdb;hpb=fe8dc2171009e90d2335f159d05f7e6660e24780;p=pspp

diff --git a/src/language/stats/quick-cluster.c b/src/language/stats/quick-cluster.c
index 014406f098..0cd08d02a9 100644
--- a/src/language/stats/quick-cluster.c
+++ b/src/language/stats/quick-cluster.c
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2011 Free Software Foundation, Inc.
+   Copyright (C) 2011, 2012, 2015, 2019 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,60 +36,121 @@
 #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"
-#include "output/text-item.h"
+#include "output/pivot-table.h"
+#include "output/output-item.h"
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
 #define N_(msgid) msgid
 
+enum missing_type
+  {
+    MISS_LISTWISE,
+    MISS_PAIRWISE,
+  };
+
+
+struct save_trans_data
+{
+  /* A writer which contains the values (if any) to be appended to
+     each case in the active dataset   */
+  struct casewriter *writer;
+
+  /* A reader created from the writer above. */
+  struct casereader *appending_reader;
+
+  /* The indices to be used to access values in the above,
+     reader/writer  */
+  int CASE_IDX_MEMBERSHIP;
+  int CASE_IDX_DISTANCE;
+
+  /* The variables created to hold the values appended to the dataset  */
+  struct variable *membership;
+  struct variable *distance;
+};
+
+
+#define SAVE_MEMBERSHIP 0x1
+#define SAVE_DISTANCE   0x2
+
+struct qc
+{
+  struct dataset *dataset;
+  struct dictionary *dict;
+
+  const struct variable **vars;
+  size_t n_vars;
+
+  double epsilon;               /* The convergence criterion */
+
+  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;
+
+  /* Which values are to be saved?  */
+  int save_values;
+
+  /* The name of the new variable to contain the cluster of each case.  */
+  char *var_membership;
+
+  /* The name of the new variable to contain the distance of each case
+     from its cluster centre.  */
+  char *var_distance;
+
+  struct save_trans_data *save_trans_data;
+};
+
 /* Holds all of the information for the functions.  int n, holds the number of
    observation and its default value is -1.  We set it in
    kmeans_recalculate_centers in first invocation. */
 struct Kmeans
 {
   gsl_matrix *centers;		/* Centers for groups. */
+  gsl_matrix *updated_centers;
+  casenumber n;
+
   gsl_vector_long *num_elements_groups;
-  int ngroups;			/* Number of group. (Given by the user) */
-  casenumber n;			/* Number of observations (default -1). */
-  int m;			/* Number of variables. (Given by the user) */
-  int maxiter;			/* Maximum iterations (Given by the user) */
-  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. */
-  const struct variable **variables;
+  double convergence_criteria;
   gsl_permutation *group_order;	/* Group order for reporting. */
-  struct casereader *original_casereader;
-  struct caseproto *proto;
-  struct casereader *index_rdr;	/* Group ids for each case. */
-  const struct variable *wv;	/* Weighting variable. */
 };
 
-static struct Kmeans *kmeans_create (struct casereader *cs,
-				     const struct variable **variables,
-				     int m, int ngroups, int maxiter);
+static struct Kmeans *kmeans_create (const struct qc *qc);
 
-static void kmeans_randomize_centers (struct Kmeans *kmeans);
+static void kmeans_get_nearest_group (const struct Kmeans *kmeans,
+				      struct ccase *c, const struct qc *,
+				      int *, double *, int *, double *);
 
-static int kmeans_get_nearest_group (struct Kmeans *kmeans, struct ccase *c);
+static void kmeans_order_groups (struct Kmeans *kmeans, const struct qc *);
 
-static void kmeans_recalculate_centers (struct Kmeans *kmeans);
+static void kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader,
+			    const struct qc *);
 
-static int
-kmeans_calculate_indexes_and_check_convergence (struct Kmeans *kmeans);
+static void quick_cluster_show_centers (struct Kmeans *kmeans, bool initial,
+					const struct qc *);
 
-static void kmeans_order_groups (struct Kmeans *kmeans);
+static void quick_cluster_show_membership (struct Kmeans *kmeans,
+					   const struct casereader *reader,
+					   struct qc *);
 
-static void kmeans_cluster (struct Kmeans *kmeans);
+static void quick_cluster_show_number_cases (struct Kmeans *kmeans,
+					     const struct qc *);
 
-static void quick_cluster_show_centers (struct Kmeans *kmeans, bool initial);
-
-static void quick_cluster_show_number_cases (struct Kmeans *kmeans);
-
-static void quick_cluster_show_results (struct Kmeans *kmeans);
+static void quick_cluster_show_results (struct Kmeans *kmeans,
+					const struct casereader *reader,
+					struct qc *);
 
 int cmd_quick_cluster (struct lexer *lexer, struct dataset *ds);
 
@@ -100,26 +160,15 @@ static void kmeans_destroy (struct Kmeans *kmeans);
    variables 'variables', number of cases 'n', number of variables 'm', number
    of clusters and amount of maximum iterations. */
 static struct Kmeans *
-kmeans_create (struct casereader *cs, const struct variable **variables,
-	       int m, int ngroups, int maxiter)
+kmeans_create (const struct qc *qc)
 {
   struct Kmeans *kmeans = xmalloc (sizeof (struct Kmeans));
-  kmeans->centers = gsl_matrix_alloc (ngroups, m);
-  kmeans->num_elements_groups = gsl_vector_long_alloc (ngroups);
-  kmeans->ngroups = ngroups;
-  kmeans->n = 0;
-  kmeans->m = m;
-  kmeans->maxiter = maxiter;
-  kmeans->lastiter = 0;
-  kmeans->trials = 0;
-  kmeans->variables = variables;
+  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->group_order = gsl_permutation_alloc (kmeans->centers->size1);
-  kmeans->original_casereader = cs;
   kmeans->initial_centers = NULL;
 
-  kmeans->proto = caseproto_create ();
-  kmeans->proto = caseproto_add_width (kmeans->proto, 0);
-  kmeans->index_rdr = NULL;
   return (kmeans);
 }
 
@@ -127,267 +176,392 @@ 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);
+}
 
-  /*
-     These reader and writer were already destroyed.
-     free (kmeans->original_casereader);
-     free (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)
+
+
+static double
+matrix_mindist (const gsl_matrix *m, int *mn, int *mm)
 {
   int i, j;
-  for (i = 0; i < kmeans->ngroups; i++)
+  double mindist = INFINITY;
+  for (i = 0; i < m->size1 - 1; ++i)
     {
-      for (j = 0; j < kmeans->m; 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 (kmeans->ngroups, kmeans->m);
-      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)
+/* 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;
-  for (i = 0; i < kmeans->ngroups; i++)
+  int j, i;
+
+  double mindist = INFINITY;
+  for (i = 0; i < qc->ngroups; i++)
     {
-      dist = 0;
-      for (j = 0; j < kmeans->m; j++)
+      if (i == which)
+	continue;
+
+      double dist = 0;
+      for (j = 0; j < qc->n_vars; j++)
 	{
-	  x = case_data (c, kmeans->variables[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)
+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;
-  double weight;
-
-  i = 0;
-  cs = casereader_clone (kmeans->original_casereader);
-  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))
     {
-      c_index = casereader_read (cs_index);
-      index = case_data_idx (c_index, 0)->f;
-      for (v = 0; v < kmeans->m; ++v)
+      bool missing = false;
+      for (j = 0; j < qc->n_vars; ++j)
 	{
-	  if (kmeans->wv)
-	    {
-	      weight = case_data (c, kmeans->wv)->f;
-	    }
-	  else
+	  const union value *val = case_data (c, qc->vars[j]);
+	  if (var_is_value_missing (qc->vars[j], val) & qc->exclude)
 	    {
-	      weight = 1.0;
+	      missing = true;
+	      break;
 	    }
-	  x = case_data (c, kmeans->variables[v])->f * weight;
-	  curval = gsl_matrix_get (kmeans->centers, index, v);
-	  gsl_matrix_set (kmeans->centers, index, v, curval + x);
+
+	  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 < kmeans->ngroups; i++)
-    {
-      casenumber numobs = kmeans->num_elements_groups->data[i];
-      for (j = 0; j < kmeans->m; 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)
-{
-  int totaldiff = 0;
-  double weight;
-  struct ccase *c;
-  struct casereader *cs = casereader_clone (kmeans->original_casereader);
+  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);
-      if (kmeans->wv)
-	{
-	  weight = (casenumber) case_data (c, kmeans->wv)->f;
-	}
-      else
+      double dist = 0;
+      for (j = 0; j < qc->n_vars; j++)
 	{
-	  weight = 1.0;
-	}
-      kmeans->num_elements_groups->data[bestindex] += weight;
-      if (kmeans->index_rdr)
-	{
-	  /* 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;
+
+  if (g_q)
+    *g_q = result0;
 
-      /* 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);
 
-  /* Convert the writer into a reader, ready for the next iteration to read */
-  kmeans->index_rdr = casewriter_make_reader (index_wtr);
+  if (delta_p)
+    *delta_p = mindist1;
 
-  return (totaldiff);
+  if (g_p)
+    *g_p = result1;
 }
 
+
+
 static void
-kmeans_order_groups (struct Kmeans *kmeans)
+kmeans_order_groups (struct Kmeans *kmeans, const struct qc *qc)
 {
-  gsl_vector *v = gsl_vector_alloc (kmeans->ngroups);
+  gsl_vector *v = gsl_vector_alloc (qc->ngroups);
   gsl_matrix_get_col (v, kmeans->centers, 0);
   gsl_sort_vector_index (kmeans->group_order, v);
+  gsl_vector_free (v);
 }
 
 /* Main algorithm.
    Does iterations, checks convergency. */
 static void
-kmeans_cluster (struct Kmeans *kmeans)
+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);
-  for (kmeans->lastiter = 0; kmeans->lastiter < kmeans->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);
-      kmeans_recalculate_centers (kmeans);
-      if (show_warning1 && kmeans->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_num (c, qc->wv) : 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_num (c, qc->wv) : 1.0);
+		}
+	    }
+
+	  casereader_destroy (r);
 	}
-      if (diffs == 0)
-	break;
-    }
 
-  for (i = 0; i < kmeans->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_num (c, qc->wv) : 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.
@@ -395,130 +569,349 @@ cluster:
    If initial is true, initial cluster centers are reported.  Otherwise,
    resulted centers are reported. */
 static void
-quick_cluster_show_centers (struct Kmeans *kmeans, bool initial)
+quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *qc)
 {
-  struct tab_table *t;
-  int nc, nr, heading_columns, currow;
-  int i, j;
-  nc = kmeans->ngroups + 1;
-  nr = kmeans->m + 4;
-  heading_columns = 1;
-  t = tab_create (nc, nr);
-  tab_headers (t, 0, nc - 1, 0, 1);
-  currow = 0;
-  if (!initial)
-    {
-      tab_title (t, _("Final Cluster Centers"));
-    }
-  else
-    {
-      tab_title (t, _("Initial Cluster Centers"));
-    }
-  tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
-  tab_joint_text (t, 1, 0, nc - 1, 0, TAB_CENTER, _("Cluster"));
-  tab_hline (t, TAL_1, 1, nc - 1, 2);
-  currow += 2;
+  struct pivot_table *table
+    = pivot_table_create (initial
+			  ? N_("Initial Cluster Centers")
+			  : N_("Final Cluster Centers"));
+
+  struct pivot_dimension *clusters
+    = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Cluster"));
+
+  clusters->root->show_label = true;
+  for (size_t i = 0; i < qc->ngroups; i++)
+    pivot_category_create_leaf (clusters->root,
+                                pivot_value_new_integer (i + 1));
+
+  struct pivot_dimension *variables
+    = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Variable"));
+
+  for (size_t i = 0; i < qc->n_vars; i++)
+    pivot_category_create_leaf (variables->root,
+                                pivot_value_new_variable (qc->vars[i]));
+
+  const gsl_matrix *matrix = (initial
+                              ? kmeans->initial_centers
+                              : kmeans->centers);
+  for (size_t i = 0; i < qc->ngroups; i++)
+    for (size_t j = 0; j < qc->n_vars; j++)
+      {
+        double x = gsl_matrix_get (matrix, kmeans->group_order->data[i], j);
+        union value v = { .f = x };
+        pivot_table_put2 (table, i, j,
+                          pivot_value_new_var_value (qc->vars[j], &v));
+      }
+
+  pivot_table_submit (table);
+}
+
+
+/* A transformation function which juxtaposes the dataset with the
+   (pre-prepared) dataset containing membership and/or distance
+   values.  */
+static enum trns_result
+save_trans_func (void *aux, struct ccase **c, casenumber x UNUSED)
+{
+  const struct save_trans_data *std = aux;
+  struct ccase *ca  = casereader_read (std->appending_reader);
+  if (ca == NULL)
+    return TRNS_CONTINUE;
+
+  *c = case_unshare (*c);
 
-  for (i = 0; i < kmeans->ngroups; i++)
+  if (std->CASE_IDX_MEMBERSHIP >= 0)
+    *case_num_rw (*c, std->membership) = case_num_idx (ca, std->CASE_IDX_MEMBERSHIP);
+
+  if (std->CASE_IDX_DISTANCE >= 0)
+    *case_num_rw (*c, std->distance) = case_num_idx (ca, std->CASE_IDX_DISTANCE);
+
+  case_unref (ca);
+
+  return TRNS_CONTINUE;
+}
+
+
+/* Free the resources of the transformation.  */
+static bool
+save_trans_destroy (void *aux)
+{
+  struct save_trans_data *std = aux;
+  casereader_destroy (std->appending_reader);
+  free (std);
+  return true;
+}
+
+
+/* Reports cluster membership for each case, and is requested
+saves the membership and the distance of the case from the cluster
+centre.  */
+static void
+quick_cluster_show_membership (struct Kmeans *kmeans,
+			       const struct casereader *reader,
+			       struct qc *qc)
+{
+  struct pivot_table *table = NULL;
+  struct pivot_dimension *cases = NULL;
+  if (qc->print_cluster_membership)
     {
-      tab_text_format (t, (i + 1), currow, TAB_CENTER, "%d", (i + 1));
+      table = pivot_table_create (N_("Cluster Membership"));
+
+      pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Cluster"),
+			      N_("Cluster"));
+
+      cases
+	= pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Case Number"));
+
+      cases->root->show_label = true;
     }
-  currow++;
-  tab_hline (t, TAL_1, 1, nc - 1, currow);
-  currow++;
-  for (i = 0; i < kmeans->m; i++)
+
+  gsl_permutation *ip = gsl_permutation_alloc (qc->ngroups);
+  gsl_permutation_inverse (ip, kmeans->group_order);
+
+  struct caseproto *proto = caseproto_create ();
+  if (qc->save_values)
     {
-      tab_text (t, 0, currow + i, TAB_LEFT,
-		var_to_string (kmeans->variables[i]));
+      /* Prepare data which may potentially be used in a
+	 transformation appending new variables to the active
+	 dataset.  */
+      qc->save_trans_data = xzalloc (sizeof *qc->save_trans_data);
+      qc->save_trans_data->CASE_IDX_MEMBERSHIP = -1;
+      qc->save_trans_data->CASE_IDX_DISTANCE = -1;
+      qc->save_trans_data->writer = autopaging_writer_create (proto);
+
+      int idx = 0;
+      if (qc->save_values & SAVE_MEMBERSHIP)
+	{
+	  proto = caseproto_add_width (proto, 0);
+	  qc->save_trans_data->CASE_IDX_MEMBERSHIP = idx++;
+	}
+
+      if (qc->save_values & SAVE_DISTANCE)
+	{
+	  proto = caseproto_add_width (proto, 0);
+	  qc->save_trans_data->CASE_IDX_DISTANCE = idx++;
+	}
     }
 
-  for (i = 0; i < kmeans->ngroups; i++)
+  struct casereader *cs = casereader_clone (reader);
+  struct ccase *c;
+  for (int i = 0; (c = casereader_read (cs)) != NULL; i++, case_unref (c))
     {
-      for (j = 0; j < kmeans->m; j++)
+      assert (i < kmeans->n);
+      int clust;
+      kmeans_get_nearest_group (kmeans, c, qc, &clust, NULL, NULL, NULL);
+      int cluster = ip->data[clust];
+
+      if (qc->save_trans_data)
+      {
+	/* Calculate the membership and distance values.  */
+	struct ccase *outc = case_create (proto);
+	if (qc->save_values & SAVE_MEMBERSHIP)
+	  *case_num_rw_idx (outc, qc->save_trans_data->CASE_IDX_MEMBERSHIP) = cluster + 1;
+
+	if (qc->save_values & SAVE_DISTANCE)
+	  *case_num_rw_idx (outc, qc->save_trans_data->CASE_IDX_DISTANCE)
+	    = sqrt (dist_from_case (kmeans, c, qc, clust));
+
+	casewriter_write (qc->save_trans_data->writer, outc);
+      }
+
+      if (qc->print_cluster_membership)
 	{
-	  if (!initial)
-	    {
-	      tab_double (t, i + 1, j + 4, TAB_CENTER,
-			  gsl_matrix_get (kmeans->centers,
-					  kmeans->group_order->data[i], j),
-			  var_get_print_format (kmeans->variables[j]));
-	    }
-	  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 (kmeans->variables[j]));
-	    }
+	  /* Print the cluster membership to the table.  */
+	  int case_idx = pivot_category_create_leaf (cases->root,
+						 pivot_value_new_integer (i + 1));
+	  pivot_table_put2 (table, 0, case_idx,
+			    pivot_value_new_integer (cluster + 1));
 	}
     }
-  tab_submit (t);
+
+  caseproto_unref (proto);
+  gsl_permutation_free (ip);
+
+  if (qc->print_cluster_membership)
+    pivot_table_submit (table);
+  casereader_destroy (cs);
 }
 
+
 /* Reports number of cases of each single cluster. */
 static void
-quick_cluster_show_number_cases (struct Kmeans *kmeans)
+quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
 {
-  struct tab_table *t;
-  int nc, nr;
-  int i, numelem;
-  long int total;
-  nc = 3;
-  nr = kmeans->ngroups + 1;
-  t = tab_create (nc, nr);
-  tab_headers (t, 0, nc - 1, 0, 0);
-  tab_title (t, _("Number of Cases in each Cluster"));
-  tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
-  tab_text (t, 0, 0, TAB_LEFT, _("Cluster"));
-
-  total = 0;
-  for (i = 0; i < kmeans->ngroups; i++)
+  struct pivot_table *table
+    = pivot_table_create (N_("Number of Cases in each Cluster"));
+
+  pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
+                          N_("Count"));
+
+  struct pivot_dimension *clusters
+    = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Clusters"));
+
+  struct pivot_category *group
+    = pivot_category_create_group (clusters->root, N_("Cluster"));
+
+  long int total = 0;
+  for (int i = 0; i < qc->ngroups; i++)
     {
-      tab_text_format (t, 1, i, TAB_CENTER, "%d", (i + 1));
-      numelem =
-	kmeans->num_elements_groups->data[kmeans->group_order->data[i]];
-      tab_text_format (t, 2, i, TAB_CENTER, "%d", numelem);
-      total += numelem;
+      int cluster_idx
+	= pivot_category_create_leaf (group, pivot_value_new_integer (i + 1));
+      int count = kmeans->num_elements_groups->data [kmeans->group_order->data[i]];
+      pivot_table_put2 (table, 0, cluster_idx, pivot_value_new_integer (count));
+      total += count;
     }
 
-  tab_text (t, 0, kmeans->ngroups, TAB_LEFT, _("Valid"));
-  tab_text_format (t, 2, kmeans->ngroups, TAB_LEFT, "%ld", total);
-  tab_submit (t);
+  int cluster_idx = pivot_category_create_leaf (clusters->root,
+						pivot_value_new_text (N_("Valid")));
+  pivot_table_put2 (table, 0, cluster_idx, pivot_value_new_integer (total));
+  pivot_table_submit (table);
 }
 
 /* Reports. */
 static void
-quick_cluster_show_results (struct Kmeans *kmeans)
+quick_cluster_show_results (struct Kmeans *kmeans, const struct casereader *reader,
+			    struct qc *qc)
 {
-  kmeans_order_groups (kmeans);
-  /* Uncomment the line below for reporting initial centers. */
-  /* quick_cluster_show_centers (kmeans, true); */
-  quick_cluster_show_centers (kmeans, false);
-  quick_cluster_show_number_cases (kmeans);
+  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);
+
+  quick_cluster_show_membership (kmeans, reader, qc);
 }
 
-int
-cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
+/* Parse the QUICK CLUSTER command and populate QC accordingly.
+   Returns false on error.  */
+static bool
+quick_cluster_parse (struct lexer *lexer, struct qc *qc)
 {
-  struct Kmeans *kmeans;
-  bool ok;
-  const struct dictionary *dict = dataset_dict (ds);
-  const struct variable **variables;
-  struct casereader *cs;
-  int groups = 2;
-  int maxiter = 2;
-  size_t p;
-
-  if (!parse_variables_const (lexer, dict, &variables, &p,
+  if (!parse_variables_const (lexer, qc->dict, &qc->vars, &qc->n_vars,
 			      PV_NO_DUPLICATE | PV_NUMERIC))
     {
-      msg (ME, _("Variables cannot be parsed"));
       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);
+		  return false;
+		}
+	    }
+	}
+      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);
+		  return false;
+		}
+	    }
+	}
+      else if (lex_match_id (lexer, "SAVE"))
+	{
+	  lex_match (lexer, T_EQUALS);
+	  while (lex_token (lexer) != T_ENDCMD
+		 && lex_token (lexer) != T_SLASH)
+	    {
+	      if (lex_match_id (lexer, "CLUSTER"))
+		{
+		  qc->save_values |= SAVE_MEMBERSHIP;
+		  if (lex_match (lexer, T_LPAREN))
+		    {
+		      if (!lex_force_id (lexer))
+			return false;
+
+		      free (qc->var_membership);
+		      qc->var_membership = xstrdup (lex_tokcstr (lexer));
+		      if (NULL != dict_lookup_var (qc->dict, qc->var_membership))
+			{
+			  lex_error (lexer,
+				     _("A variable called `%s' already exists."),
+				     qc->var_membership);
+			  free (qc->var_membership);
+			  qc->var_membership = NULL;
+			  return false;
+			}
+
+		      lex_get (lexer);
+
+		      if (!lex_force_match (lexer, T_RPAREN))
+			return false;
+		    }
+		}
+	      else if (lex_match_id (lexer, "DISTANCE"))
+		{
+		  qc->save_values |= SAVE_DISTANCE;
+		  if (lex_match (lexer, T_LPAREN))
+		    {
+		      if (!lex_force_id (lexer))
+			return false;
+
+		      free (qc->var_distance);
+		      qc->var_distance = xstrdup (lex_tokcstr (lexer));
+		      if (NULL != dict_lookup_var (qc->dict, qc->var_distance))
+			{
+			  lex_error (lexer,
+				     _("A variable called `%s' already exists."),
+				     qc->var_distance);
+			  free (qc->var_distance);
+			  qc->var_distance = NULL;
+			  return false;
+			}
+
+		      lex_get (lexer);
+
+		      if (!lex_force_match (lexer, T_RPAREN))
+			return false;
+		    }
+		}
+	      else
+		{
+		  lex_error (lexer, _("Expecting %s or %s."),
+			     "CLUSTER", "DISTANCE");
+		  return false;
+		}
+	    }
+	}
+      else if (lex_match_id (lexer, "CRITERIA"))
 	{
 	  lex_match (lexer, T_EQUALS);
 	  while (lex_token (lexer) != T_ENDCMD
@@ -526,40 +919,174 @@ 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_range (lexer, "CLUSTERS", 1, INT_MAX))
+		    {
+		      qc->ngroups = lex_integer (lexer);
+		      lex_get (lexer);
+		      if (!lex_force_match (lexer, T_RPAREN))
+			return false;
+		    }
+		}
+	      else if (lex_match_id (lexer, "CONVERGE"))
+		{
+		  if (lex_force_match (lexer, T_LPAREN) &&
+		      lex_force_num_range_open (lexer, "CONVERGE", 0, DBL_MAX))
 		    {
-		      lex_force_int (lexer);
-		      groups = lex_integer (lexer);
+		      qc->epsilon = lex_number (lexer);
 		      lex_get (lexer);
-		      lex_force_match (lexer, T_RPAREN);
+		      if (!lex_force_match (lexer, T_RPAREN))
+			return false;
 		    }
 		}
 	      else if (lex_match_id (lexer, "MXITER"))
 		{
-		  if (lex_force_match (lexer, T_LPAREN))
+		  if (lex_force_match (lexer, T_LPAREN) &&
+		      lex_force_int_range (lexer, "MXITER", 1, INT_MAX))
 		    {
-		      lex_force_int (lexer);
-		      maxiter = lex_integer (lexer);
+		      qc->maxiter = lex_integer (lexer);
 		      lex_get (lexer);
-		      lex_force_match (lexer, T_RPAREN);
+		      if (!lex_force_match (lexer, T_RPAREN))
+			return false;
 		    }
 		}
+	      else if (lex_match_id (lexer, "NOINITIAL"))
+		{
+		  qc->no_initial = true;
+		}
+	      else if (lex_match_id (lexer, "NOUPDATE"))
+		{
+		  qc->no_update = true;
+		}
 	      else
-                return CMD_FAILURE;
+		{
+		  lex_error (lexer, NULL);
+		  return false;
+		}
 	    }
 	}
+      else
+        {
+          lex_error (lexer, NULL);
+          return false;
+        }
     }
+  return true;
+}
+
+int
+cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
+{
+  struct qc qc;
+  struct Kmeans *kmeans;
+  bool ok;
+  memset (&qc, 0, sizeof qc);
+  qc.dataset = ds;
+  qc.dict =  dataset_dict (ds);
+  qc.ngroups = 2;
+  qc.maxiter = 10;
+  qc.epsilon = DBL_EPSILON;
+  qc.missing_type = MISS_LISTWISE;
+  qc.exclude = MV_ANY;
+
+
+  if (!quick_cluster_parse (lexer, &qc))
+    goto error;
+
+  qc.wv = dict_get_weight (qc.dict);
+
+  {
+    struct casereader *group;
+    struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), qc.dict);
+
+    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, group, &qc);
+	kmeans_destroy (kmeans);
+	casereader_destroy (group);
+      }
+    ok = casegrouper_destroy (grouper);
+  }
+  ok = proc_commit (ds) && ok;
+
+
+  /* If requested, set a transformation to append the cluster and
+     distance values to the current dataset.  */
+  if (qc.save_trans_data)
+    {
+      struct save_trans_data *std = qc.save_trans_data;
+      std->appending_reader = casewriter_make_reader (std->writer);
+      std->writer = NULL;
 
-  cs = proc_open (ds);
+      if (qc.save_values & SAVE_MEMBERSHIP)
+	{
+	  /* Invent a variable name if necessary.  */
+	  int idx = 0;
+	  struct string name;
+	  ds_init_empty (&name);
+	  while (qc.var_membership == NULL)
+	    {
+	      ds_clear (&name);
+	      ds_put_format (&name, "QCL_%d", idx++);
 
-  kmeans = kmeans_create (cs, variables, p, groups, maxiter);
+	      if (!dict_lookup_var (qc.dict, ds_cstr (&name)))
+		{
+		  qc.var_membership = strdup (ds_cstr (&name));
+		  break;
+		}
+	    }
+	  ds_destroy (&name);
 
-  kmeans->wv = dict_get_weight (dict);
-  kmeans_cluster (kmeans);
-  quick_cluster_show_results (kmeans);
-  ok = proc_commit (ds);
+	  std->membership = dict_create_var_assert (qc.dict, qc.var_membership, 0);
+	}
 
-  kmeans_destroy (kmeans);
+      if (qc.save_values & SAVE_DISTANCE)
+	{
+	  /* Invent a variable name if necessary.  */
+	  int idx = 0;
+	  struct string name;
+	  ds_init_empty (&name);
+	  while (qc.var_distance == NULL)
+	    {
+	      ds_clear (&name);
+	      ds_put_format (&name, "QCL_%d", idx++);
 
+	      if (!dict_lookup_var (qc.dict, ds_cstr (&name)))
+		{
+		  qc.var_distance = strdup (ds_cstr (&name));
+		  break;
+		}
+	    }
+	  ds_destroy (&name);
+
+	  std->distance = dict_create_var_assert (qc.dict, qc.var_distance, 0);
+	}
+
+      static const struct trns_class trns_class = {
+        .name = "QUICK CLUSTER",
+        .execute = save_trans_func,
+        .destroy = save_trans_destroy,
+      };
+      add_transformation (qc.dataset, &trns_class, std);
+    }
+
+  free (qc.var_distance);
+  free (qc.var_membership);
+  free (qc.vars);
   return (ok);
+
+ error:
+  free (qc.var_distance);
+  free (qc.var_membership);
+  free (qc.vars);
+  return CMD_FAILURE;
 }