GLM: Disable the interactions option until it actually gives the right results
[pspp] / src / language / stats / quick-cluster.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2011 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, heading_columns, currow;
394   int i, j;
395   nc = qc->ngroups + 1;
396   nr = qc->n_vars + 4;
397   heading_columns = 1;
398   t = tab_create (nc, nr);
399   tab_headers (t, 0, nc - 1, 0, 1);
400   currow = 0;
401   if (!initial)
402     {
403       tab_title (t, _("Final Cluster Centers"));
404     }
405   else
406     {
407       tab_title (t, _("Initial Cluster Centers"));
408     }
409   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
410   tab_joint_text (t, 1, 0, nc - 1, 0, TAB_CENTER, _("Cluster"));
411   tab_hline (t, TAL_1, 1, nc - 1, 2);
412   currow += 2;
413
414   for (i = 0; i < qc->ngroups; i++)
415     {
416       tab_text_format (t, (i + 1), currow, TAB_CENTER, "%d", (i + 1));
417     }
418   currow++;
419   tab_hline (t, TAL_1, 1, nc - 1, currow);
420   currow++;
421   for (i = 0; i < qc->n_vars; i++)
422     {
423       tab_text (t, 0, currow + i, TAB_LEFT,
424                 var_to_string (qc->vars[i]));
425     }
426
427   for (i = 0; i < qc->ngroups; i++)
428     {
429       for (j = 0; j < qc->n_vars; j++)
430         {
431           if (!initial)
432             {
433               tab_double (t, i + 1, j + 4, TAB_CENTER,
434                           gsl_matrix_get (kmeans->centers,
435                                           kmeans->group_order->data[i], j),
436                           var_get_print_format (qc->vars[j]));
437             }
438           else
439             {
440               tab_double (t, i + 1, j + 4, TAB_CENTER,
441                           gsl_matrix_get (kmeans->initial_centers,
442                                           kmeans->group_order->data[i], j),
443                           var_get_print_format (qc->vars[j]));
444             }
445         }
446     }
447   tab_submit (t);
448 }
449
450 /* Reports number of cases of each single cluster. */
451 static void
452 quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
453 {
454   struct tab_table *t;
455   int nc, nr;
456   int i, numelem;
457   long int total;
458   nc = 3;
459   nr = qc->ngroups + 1;
460   t = tab_create (nc, nr);
461   tab_headers (t, 0, nc - 1, 0, 0);
462   tab_title (t, _("Number of Cases in each Cluster"));
463   tab_box (t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, nc - 1, nr - 1);
464   tab_text (t, 0, 0, TAB_LEFT, _("Cluster"));
465
466   total = 0;
467   for (i = 0; i < qc->ngroups; i++)
468     {
469       tab_text_format (t, 1, i, TAB_CENTER, "%d", (i + 1));
470       numelem =
471         kmeans->num_elements_groups->data[kmeans->group_order->data[i]];
472       tab_text_format (t, 2, i, TAB_CENTER, "%d", numelem);
473       total += numelem;
474     }
475
476   tab_text (t, 0, qc->ngroups, TAB_LEFT, _("Valid"));
477   tab_text_format (t, 2, qc->ngroups, TAB_LEFT, "%ld", total);
478   tab_submit (t);
479 }
480
481 /* Reports. */
482 static void
483 quick_cluster_show_results (struct Kmeans *kmeans, const struct qc *qc)
484 {
485   kmeans_order_groups (kmeans, qc);
486   /* Uncomment the line below for reporting initial centers. */
487   /* quick_cluster_show_centers (kmeans, true); */
488   quick_cluster_show_centers (kmeans, false, qc);
489   quick_cluster_show_number_cases (kmeans, qc);
490 }
491
492 int
493 cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
494 {
495   struct qc qc;
496   struct Kmeans *kmeans;
497   bool ok;
498   const struct dictionary *dict = dataset_dict (ds);
499   qc.ngroups = 2;
500   qc.maxiter = 2;
501   qc.missing_type = MISS_LISTWISE;
502   qc.exclude = MV_ANY;
503
504   if (!parse_variables_const (lexer, dict, &qc.vars, &qc.n_vars,
505                               PV_NO_DUPLICATE | PV_NUMERIC))
506     {
507       return (CMD_FAILURE);
508     }
509
510   while (lex_token (lexer) != T_ENDCMD)
511     {
512       lex_match (lexer, T_SLASH);
513
514       if (lex_match_id (lexer, "MISSING"))
515         {
516           lex_match (lexer, T_EQUALS);
517           while (lex_token (lexer) != T_ENDCMD
518                  && lex_token (lexer) != T_SLASH)
519             {
520               if (lex_match_id (lexer, "LISTWISE") || lex_match_id (lexer, "DEFAULT"))
521                 {
522                   qc.missing_type = MISS_LISTWISE;
523                 }
524               else if (lex_match_id (lexer, "PAIRWISE"))
525                 {
526                   qc.missing_type = MISS_PAIRWISE;
527                 }
528               else if (lex_match_id (lexer, "INCLUDE"))
529                 {
530                   qc.exclude = MV_SYSTEM;
531                 }
532               else if (lex_match_id (lexer, "EXCLUDE"))
533                 {
534                   qc.exclude = MV_ANY;
535                 }
536               else
537                 goto error;
538             }     
539         }
540       else if (lex_match_id (lexer, "CRITERIA"))
541         {
542           lex_match (lexer, T_EQUALS);
543           while (lex_token (lexer) != T_ENDCMD
544                  && lex_token (lexer) != T_SLASH)
545             {
546               if (lex_match_id (lexer, "CLUSTERS"))
547                 {
548                   if (lex_force_match (lexer, T_LPAREN))
549                     {
550                       lex_force_int (lexer);
551                       qc.ngroups = lex_integer (lexer);
552                       lex_get (lexer);
553                       lex_force_match (lexer, T_RPAREN);
554                     }
555                 }
556               else if (lex_match_id (lexer, "MXITER"))
557                 {
558                   if (lex_force_match (lexer, T_LPAREN))
559                     {
560                       lex_force_int (lexer);
561                       qc.maxiter = lex_integer (lexer);
562                       lex_get (lexer);
563                       lex_force_match (lexer, T_RPAREN);
564                     }
565                 }
566               else
567                 goto error;
568             }
569         }
570     }
571
572   qc.wv = dict_get_weight (dict);
573
574   {
575     struct casereader *group;
576     struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
577
578     while (casegrouper_get_next_group (grouper, &group))
579       {
580         if ( qc.missing_type == MISS_LISTWISE )
581           {
582             group  = casereader_create_filter_missing (group, qc.vars, qc.n_vars,
583                                                      qc.exclude,
584                                                      NULL,  NULL);
585           }
586
587         kmeans = kmeans_create (&qc);
588         kmeans_cluster (kmeans, group, &qc);
589         quick_cluster_show_results (kmeans, &qc);
590         kmeans_destroy (kmeans);
591         casereader_destroy (group);
592       }
593     ok = casegrouper_destroy (grouper);
594   }
595   ok = proc_commit (ds) && ok;
596
597   free (qc.vars);
598
599   return (ok);
600
601  error:
602   free (qc.vars);
603   return CMD_FAILURE;
604 }