pintos-os.org Git - pspp/blob - src/language/stats/quick-cluster.c

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