Update version number to 0.8.0.
[pspp] / src / language / stats / quick-cluster.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2011, 2012 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/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   int ngroups;                  /* Number of group. (Given by the user) */
62   int maxiter;                  /* Maximum iterations (Given by the user) */
63
64   const struct variable *wv;    /* Weighting variable. */
65
66   enum missing_type missing_type;
67   enum mv_class exclude;
68 };
69
70 /* Holds all of the information for the functions.  int n, holds the number of
71    observation and its default value is -1.  We set it in
72    kmeans_recalculate_centers in first invocation. */
73 struct Kmeans
74 {
75   gsl_matrix *centers;          /* Centers for groups. */
76   gsl_vector_long *num_elements_groups;
77
78   casenumber n;                 /* Number of observations (default -1). */
79
80   int lastiter;                 /* Iteration where it found the solution. */
81   int trials;                   /* If not convergence, how many times has
82                                    clustering done. */
83   gsl_matrix *initial_centers;  /* Initial random centers. */
84
85   gsl_permutation *group_order; /* Group order for reporting. */
86   struct caseproto *proto;
87   struct casereader *index_rdr; /* Group ids for each case. */
88 };
89
90 static struct Kmeans *kmeans_create (const struct qc *qc);
91
92 static void kmeans_randomize_centers (struct Kmeans *kmeans, const struct qc *qc);
93
94 static int kmeans_get_nearest_group (struct Kmeans *kmeans, struct ccase *c, const struct qc *);
95
96 static void kmeans_recalculate_centers (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *);
97
98 static int
99 kmeans_calculate_indexes_and_check_convergence (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *);
100
101 static void kmeans_order_groups (struct Kmeans *kmeans, const struct qc *);
102
103 static void kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader, const struct qc *);
104
105 static void quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *);
106
107 static void quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *);
108
109 static void quick_cluster_show_results (struct Kmeans *kmeans, const struct qc *);
110
111 int cmd_quick_cluster (struct lexer *lexer, struct dataset *ds);
112
113 static void kmeans_destroy (struct Kmeans *kmeans);
114
115 /* Creates and returns a struct of Kmeans with given casereader 'cs', parsed
116    variables 'variables', number of cases 'n', number of variables 'm', number
117    of clusters and amount of maximum iterations. */
118 static struct Kmeans *
119 kmeans_create (const struct qc *qc)
120 {
121   struct Kmeans *kmeans = xmalloc (sizeof (struct Kmeans));
122   kmeans->centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
123   kmeans->num_elements_groups = gsl_vector_long_alloc (qc->ngroups);
124   kmeans->n = 0;
125   kmeans->lastiter = 0;
126   kmeans->trials = 0;
127   kmeans->group_order = gsl_permutation_alloc (kmeans->centers->size1);
128   kmeans->initial_centers = NULL;
129
130   kmeans->proto = caseproto_create ();
131   kmeans->proto = caseproto_add_width (kmeans->proto, 0);
132   kmeans->index_rdr = NULL;
133   return (kmeans);
134 }
135
136 static void
137 kmeans_destroy (struct Kmeans *kmeans)
138 {
139   gsl_matrix_free (kmeans->centers);
140   gsl_matrix_free (kmeans->initial_centers);
141
142   gsl_vector_long_free (kmeans->num_elements_groups);
143
144   gsl_permutation_free (kmeans->group_order);
145
146   caseproto_unref (kmeans->proto);
147
148   casereader_destroy (kmeans->index_rdr);
149
150   free (kmeans);
151 }
152
153 /* Creates random centers using randomly selected cases from the data. */
154 static void
155 kmeans_randomize_centers (struct Kmeans *kmeans, const struct qc *qc)
156 {
157   int i, j;
158   for (i = 0; i < qc->ngroups; i++)
159     {
160       for (j = 0; j < qc->n_vars; j++)
161         {
162           if (i == j)
163             {
164               gsl_matrix_set (kmeans->centers, i, j, 1);
165             }
166           else
167             {
168               gsl_matrix_set (kmeans->centers, i, j, 0);
169             }
170         }
171     }
172   /* If it is the first iteration, the variable kmeans->initial_centers is NULL
173      and it is created once for reporting issues. In SPSS, initial centers are
174      shown in the reports but in PSPP it is not shown now. I am leaving it
175      here. */
176   if (!kmeans->initial_centers)
177     {
178       kmeans->initial_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
179       gsl_matrix_memcpy (kmeans->initial_centers, kmeans->centers);
180     }
181 }
182
183 static int
184 kmeans_get_nearest_group (struct Kmeans *kmeans, struct ccase *c, const struct qc *qc)
185 {
186   int result = -1;
187   int i, j;
188   double mindist = INFINITY;
189   for (i = 0; i < qc->ngroups; i++)
190     {
191       double dist = 0;
192       for (j = 0; j < qc->n_vars; j++)
193         {
194           const union value *val = case_data (c, qc->vars[j]);
195           if ( var_is_value_missing (qc->vars[j], val, qc->exclude))
196             continue;
197
198           dist += pow2 (gsl_matrix_get (kmeans->centers, i, j) - val->f);
199         }
200       if (dist < mindist)
201         {
202           mindist = dist;
203           result = i;
204         }
205     }
206   return (result);
207 }
208
209 /* Re-calculate the cluster centers. */
210 static void
211 kmeans_recalculate_centers (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
212 {
213   casenumber i = 0;
214   int v, j;
215   struct ccase *c;
216
217   struct casereader *cs = casereader_clone (reader);
218   struct casereader *cs_index = casereader_clone (kmeans->index_rdr);
219
220   gsl_matrix_set_all (kmeans->centers, 0.0);
221   for (; (c = casereader_read (cs)) != NULL; case_unref (c))
222     {
223       double weight = qc->wv ? case_data (c, qc->wv)->f : 1.0;
224       struct ccase *c_index = casereader_read (cs_index);
225       int index = case_data_idx (c_index, 0)->f;
226       for (v = 0; v < qc->n_vars; ++v)
227         {
228           const union value *val = case_data (c, qc->vars[v]);
229           double x = val->f * weight;
230           double curval;
231
232           if ( var_is_value_missing (qc->vars[v], val, qc->exclude))
233             continue;
234
235           curval = gsl_matrix_get (kmeans->centers, index, v);
236           gsl_matrix_set (kmeans->centers, index, v, curval + x);
237         }
238       i++;
239       case_unref (c_index);
240     }
241   casereader_destroy (cs);
242   casereader_destroy (cs_index);
243
244   /* Getting number of cases */
245   if (kmeans->n == 0)
246     kmeans->n = i;
247
248   /* We got sum of each center but we need averages.
249      We are dividing centers to numobs. This may be inefficient and
250      we should check it again. */
251   for (i = 0; i < qc->ngroups; i++)
252     {
253       casenumber numobs = kmeans->num_elements_groups->data[i];
254       for (j = 0; j < qc->n_vars; j++)
255         {
256           if (numobs > 0)
257             {
258               double *x = gsl_matrix_ptr (kmeans->centers, i, j);
259               *x /= numobs;
260             }
261           else
262             {
263               gsl_matrix_set (kmeans->centers, i, j, 0);
264             }
265         }
266     }
267 }
268
269 /* The variable index in struct Kmeans holds integer values that represents the
270    current groups of cases.  index[n]=a shows the nth case is belong to ath
271    cluster.  This function calculates these indexes and returns the number of
272    different cases of the new and old index variables.  If last two index
273    variables are equal, there is no any enhancement of clustering. */
274 static int
275 kmeans_calculate_indexes_and_check_convergence (struct Kmeans *kmeans, const struct casereader *reader, const struct qc *qc)
276 {
277   int totaldiff = 0;
278   struct ccase *c;
279   struct casereader *cs = casereader_clone (reader);
280
281   /* A casewriter into which we will write the indexes. */
282   struct casewriter *index_wtr = autopaging_writer_create (kmeans->proto);
283
284   gsl_vector_long_set_all (kmeans->num_elements_groups, 0);
285
286   for (; (c = casereader_read (cs)) != NULL; case_unref (c))
287     {
288       /* A case to hold the new index. */
289       struct ccase *index_case_new = case_create (kmeans->proto);
290       int bestindex = kmeans_get_nearest_group (kmeans, c, qc);
291       double weight = qc->wv ? case_data (c, qc->wv)->f : 1.0;
292       assert (bestindex < kmeans->num_elements_groups->size);
293       kmeans->num_elements_groups->data[bestindex] += weight;
294       if (kmeans->index_rdr)
295         {
296           /* A case from which the old index will be read. */
297           struct ccase *index_case_old = NULL;
298
299           /* Read the case from the index casereader. */
300           index_case_old = casereader_read (kmeans->index_rdr);
301
302           /* Set totaldiff, using the old_index. */
303           totaldiff += abs (case_data_idx (index_case_old, 0)->f - bestindex);
304
305           /* We have no use for the old case anymore, so unref it. */
306           case_unref (index_case_old);
307         }
308       else
309         {
310           /* If this is the first run, then assume index is zero. */
311           totaldiff += bestindex;
312         }
313
314       /* Set the value of the new inde.x */
315       case_data_rw_idx (index_case_new, 0)->f = bestindex;
316
317       /* and write the new index to the casewriter */
318       casewriter_write (index_wtr, index_case_new);
319     }
320   casereader_destroy (cs);
321   /* We have now read through the entire index_rdr, so it's of no use
322      anymore. */
323   casereader_destroy (kmeans->index_rdr);
324
325   /* Convert the writer into a reader, ready for the next iteration to read */
326   kmeans->index_rdr = casewriter_make_reader (index_wtr);
327
328   return (totaldiff);
329 }
330
331 static void
332 kmeans_order_groups (struct Kmeans *kmeans, const struct qc *qc)
333 {
334   gsl_vector *v = gsl_vector_alloc (qc->ngroups);
335   gsl_matrix_get_col (v, kmeans->centers, 0);
336   gsl_sort_vector_index (kmeans->group_order, v);
337   gsl_vector_free (v);
338 }
339
340 /* Main algorithm.
341    Does iterations, checks convergency. */
342 static void
343 kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader, const struct qc *qc)
344 {
345   int i;
346   bool redo;
347   int diffs;
348   bool show_warning1;
349
350   show_warning1 = true;
351 cluster:
352   redo = false;
353   kmeans_randomize_centers (kmeans, qc);
354   for (kmeans->lastiter = 0; kmeans->lastiter < qc->maxiter;
355        kmeans->lastiter++)
356     {
357       diffs = kmeans_calculate_indexes_and_check_convergence (kmeans, reader, qc);
358       kmeans_recalculate_centers (kmeans, reader, qc);
359       if (show_warning1 && qc->ngroups > kmeans->n)
360         {
361           msg (MW, _("Number of clusters may not be larger than the number "
362                      "of cases."));
363           show_warning1 = false;
364         }
365       if (diffs == 0)
366         break;
367     }
368
369   for (i = 0; i < qc->ngroups; i++)
370     {
371       if (kmeans->num_elements_groups->data[i] == 0)
372         {
373           kmeans->trials++;
374           if (kmeans->trials >= 3)
375             break;
376           redo = true;
377           break;
378         }
379     }
380   if (redo)
381     goto cluster;
382
383 }
384
385 /* Reports centers of clusters.
386    Initial parameter is optional for future use.
387    If initial is true, initial cluster centers are reported.  Otherwise,
388    resulted centers are reported. */
389 static void
390 quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *qc)
391 {
392   struct tab_table *t;
393   int nc, nr, currow;
394   int i, j;
395   nc = qc->ngroups + 1;
396   nr = qc->n_vars + 4;
397   t = tab_create (nc, nr);
398   tab_headers (t, 0, nc - 1, 0, 1);
399   currow = 0;
400   if (!initial)
401     {
402       tab_title (t, _("Final Cluster Centers"));
403     }
404   else
405     {
406       tab_title (t, _("Initial Cluster Centers"));
407     }
408   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
409   tab_joint_text (t, 1, 0, nc - 1, 0, TAB_CENTER, _("Cluster"));
410   tab_hline (t, TAL_1, 1, nc - 1, 2);
411   currow += 2;
412
413   for (i = 0; i < qc->ngroups; i++)
414     {
415       tab_text_format (t, (i + 1), currow, TAB_CENTER, "%d", (i + 1));
416     }
417   currow++;
418   tab_hline (t, TAL_1, 1, nc - 1, currow);
419   currow++;
420   for (i = 0; i < qc->n_vars; i++)
421     {
422       tab_text (t, 0, currow + i, TAB_LEFT,
423                 var_to_string (qc->vars[i]));
424     }
425
426   for (i = 0; i < qc->ngroups; i++)
427     {
428       for (j = 0; j < qc->n_vars; j++)
429         {
430           if (!initial)
431             {
432               tab_double (t, i + 1, j + 4, TAB_CENTER,
433                           gsl_matrix_get (kmeans->centers,
434                                           kmeans->group_order->data[i], j),
435                           var_get_print_format (qc->vars[j]));
436             }
437           else
438             {
439               tab_double (t, i + 1, j + 4, TAB_CENTER,
440                           gsl_matrix_get (kmeans->initial_centers,
441                                           kmeans->group_order->data[i], j),
442                           var_get_print_format (qc->vars[j]));
443             }
444         }
445     }
446   tab_submit (t);
447 }
448
449 /* Reports number of cases of each single cluster. */
450 static void
451 quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
452 {
453   struct tab_table *t;
454   int nc, nr;
455   int i, numelem;
456   long int total;
457   nc = 3;
458   nr = qc->ngroups + 1;
459   t = tab_create (nc, nr);
460   tab_headers (t, 0, nc - 1, 0, 0);
461   tab_title (t, _("Number of Cases in each Cluster"));
462   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
463   tab_text (t, 0, 0, TAB_LEFT, _("Cluster"));
464
465   total = 0;
466   for (i = 0; i < qc->ngroups; i++)
467     {
468       tab_text_format (t, 1, i, TAB_CENTER, "%d", (i + 1));
469       numelem =
470         kmeans->num_elements_groups->data[kmeans->group_order->data[i]];
471       tab_text_format (t, 2, i, TAB_CENTER, "%d", numelem);
472       total += numelem;
473     }
474
475   tab_text (t, 0, qc->ngroups, TAB_LEFT, _("Valid"));
476   tab_text_format (t, 2, qc->ngroups, TAB_LEFT, "%ld", total);
477   tab_submit (t);
478 }
479
480 /* Reports. */
481 static void
482 quick_cluster_show_results (struct Kmeans *kmeans, const struct qc *qc)
483 {
484   kmeans_order_groups (kmeans, qc);
485   /* Uncomment the line below for reporting initial centers. */
486   /* quick_cluster_show_centers (kmeans, true); */
487   quick_cluster_show_centers (kmeans, false, qc);
488   quick_cluster_show_number_cases (kmeans, qc);
489 }
490
491 int
492 cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
493 {
494   struct qc qc;
495   struct Kmeans *kmeans;
496   bool ok;
497   const struct dictionary *dict = dataset_dict (ds);
498   qc.ngroups = 2;
499   qc.maxiter = 2;
500   qc.missing_type = MISS_LISTWISE;
501   qc.exclude = MV_ANY;
502
503   if (!parse_variables_const (lexer, dict, &qc.vars, &qc.n_vars,
504                               PV_NO_DUPLICATE | PV_NUMERIC))
505     {
506       return (CMD_FAILURE);
507     }
508
509   while (lex_token (lexer) != T_ENDCMD)
510     {
511       lex_match (lexer, T_SLASH);
512
513       if (lex_match_id (lexer, "MISSING"))
514         {
515           lex_match (lexer, T_EQUALS);
516           while (lex_token (lexer) != T_ENDCMD
517                  && lex_token (lexer) != T_SLASH)
518             {
519               if (lex_match_id (lexer, "LISTWISE") || lex_match_id (lexer, "DEFAULT"))
520                 {
521                   qc.missing_type = MISS_LISTWISE;
522                 }
523               else if (lex_match_id (lexer, "PAIRWISE"))
524                 {
525                   qc.missing_type = MISS_PAIRWISE;
526                 }
527               else if (lex_match_id (lexer, "INCLUDE"))
528                 {
529                   qc.exclude = MV_SYSTEM;
530                 }
531               else if (lex_match_id (lexer, "EXCLUDE"))
532                 {
533                   qc.exclude = MV_ANY;
534                 }
535               else
536                 goto error;
537             }     
538         }
539       else if (lex_match_id (lexer, "CRITERIA"))
540         {
541           lex_match (lexer, T_EQUALS);
542           while (lex_token (lexer) != T_ENDCMD
543                  && lex_token (lexer) != T_SLASH)
544             {
545               if (lex_match_id (lexer, "CLUSTERS"))
546                 {
547                   if (lex_force_match (lexer, T_LPAREN))
548                     {
549                       lex_force_int (lexer);
550                       qc.ngroups = lex_integer (lexer);
551                       lex_get (lexer);
552                       lex_force_match (lexer, T_RPAREN);
553                     }
554                 }
555               else if (lex_match_id (lexer, "MXITER"))
556                 {
557                   if (lex_force_match (lexer, T_LPAREN))
558                     {
559                       lex_force_int (lexer);
560                       qc.maxiter = lex_integer (lexer);
561                       lex_get (lexer);
562                       lex_force_match (lexer, T_RPAREN);
563                     }
564                 }
565               else
566                 goto error;
567             }
568         }
569     }
570
571   qc.wv = dict_get_weight (dict);
572
573   {
574     struct casereader *group;
575     struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
576
577     while (casegrouper_get_next_group (grouper, &group))
578       {
579         if ( qc.missing_type == MISS_LISTWISE )
580           {
581             group  = casereader_create_filter_missing (group, qc.vars, qc.n_vars,
582                                                      qc.exclude,
583                                                      NULL,  NULL);
584           }
585
586         kmeans = kmeans_create (&qc);
587         kmeans_cluster (kmeans, group, &qc);
588         quick_cluster_show_results (kmeans, &qc);
589         kmeans_destroy (kmeans);
590         casereader_destroy (group);
591       }
592     ok = casegrouper_destroy (grouper);
593   }
594   ok = proc_commit (ds) && ok;
595
596   free (qc.vars);
597
598   return (ok);
599
600  error:
601   free (qc.vars);
602   return CMD_FAILURE;
603 }