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