--- /dev/null
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2011, 2012, 2015, 2019 Free Software Foundation, Inc.
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include <config.h>
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_sort_vector.h>
+#include <gsl/gsl_statistics.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "data/case.h"
+#include "data/casegrouper.h"
+#include "data/casereader.h"
+#include "data/casewriter.h"
+#include "data/dataset.h"
+#include "data/dictionary.h"
+#include "data/format.h"
+#include "data/missing-values.h"
+#include "language/command.h"
+#include "language/lexer/lexer.h"
+#include "language/lexer/variable-parser.h"
+#include "libpspp/message.h"
+#include "libpspp/misc.h"
+#include "libpspp/assertion.h"
+#include "libpspp/str.h"
+#include "math/random.h"
+#include "output/pivot-table.h"
+#include "output/output-item.h"
+
+#include "gettext.h"
+#define _(msgid) gettext (msgid)
+#define N_(msgid) msgid
+
+enum missing_type
+ {
+ MISS_LISTWISE,
+ MISS_PAIRWISE,
+ };
+
+
+struct save_trans_data
+ {
+ /* A writer which contains the values (if any) to be appended to
+ each case in the active dataset */
+ struct casewriter *writer;
+
+ /* A reader created from the writer above. */
+ struct casereader *appending_reader;
+
+ /* The indices to be used to access values in the above,
+ reader/writer */
+ int membership_case_idx;
+ int distance_case_idx;
+
+ /* The variables created to hold the values appended to the dataset */
+ struct variable *membership;
+ struct variable *distance;
+ };
+
+
+struct qc
+ {
+ struct dataset *dataset;
+ struct dictionary *dict;
+
+ const struct variable **vars;
+ size_t n_vars;
+
+ double epsilon; /* The convergence criterion */
+
+ int ngroups; /* Number of group. (Given by the user) */
+ int maxiter; /* Maximum iterations (Given by the user) */
+ bool print_cluster_membership; /* true => print membership */
+ bool print_initial_clusters; /* true => print initial cluster */
+ bool initial; /* false => simplified initial cluster selection */
+ bool update; /* false => do not iterate */
+
+ const struct variable *wv; /* Weighting variable. */
+
+ enum missing_type missing_type;
+ enum mv_class exclude;
+
+ /* Which values are to be saved? */
+ bool save_membership;
+ bool save_distance;
+
+ /* The name of the new variable to contain the cluster of each case. */
+ char *var_membership;
+
+ /* The name of the new variable to contain the distance of each case
+ from its cluster centre. */
+ char *var_distance;
+
+ struct save_trans_data *save_trans_data;
+ };
+
+/* Holds all of the information for the functions. int n, holds the number of
+ observation and its default value is -1. We set it in
+ kmeans_recalculate_centers in first invocation. */
+struct Kmeans
+ {
+ gsl_matrix *centers; /* Centers for groups. */
+ gsl_matrix *updated_centers;
+ casenumber n;
+
+ gsl_vector_long *num_elements_groups;
+
+ gsl_matrix *initial_centers; /* Initial random centers. */
+ double convergence_criteria;
+ gsl_permutation *group_order; /* Group order for reporting. */
+ };
+
+static struct Kmeans *kmeans_create (const struct qc *);
+
+static void kmeans_get_nearest_group (const struct Kmeans *,
+ struct ccase *, const struct qc *,
+ int *, double *, int *, double *);
+
+static void kmeans_order_groups (struct Kmeans *, const struct qc *);
+
+static void kmeans_cluster (struct Kmeans *, struct casereader *,
+ const struct qc *);
+
+static void quick_cluster_show_centers (struct Kmeans *, bool initial,
+ const struct qc *);
+
+static void quick_cluster_show_membership (struct Kmeans *,
+ const struct casereader *,
+ struct qc *);
+
+static void quick_cluster_show_number_cases (struct Kmeans *,
+ const struct qc *);
+
+static void quick_cluster_show_results (struct Kmeans *,
+ const struct casereader *,
+ struct qc *);
+
+int cmd_quick_cluster (struct lexer *, struct dataset *);
+
+static void kmeans_destroy (struct Kmeans *);
+
+/* Creates and returns a struct of Kmeans with given casereader 'cs', parsed
+ variables 'variables', number of cases 'n', number of variables 'm', number
+ of clusters and amount of maximum iterations. */
+static struct Kmeans *
+kmeans_create (const struct qc *qc)
+{
+ struct Kmeans *kmeans = xmalloc (sizeof *kmeans);
+ *kmeans = (struct Kmeans) {
+ .centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars),
+ .updated_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars),
+ .num_elements_groups = gsl_vector_long_alloc (qc->ngroups),
+ .group_order = gsl_permutation_alloc (qc->ngroups),
+ };
+ return kmeans;
+}
+
+static void
+kmeans_destroy (struct Kmeans *kmeans)
+{
+ gsl_matrix_free (kmeans->centers);
+ gsl_matrix_free (kmeans->updated_centers);
+ gsl_matrix_free (kmeans->initial_centers);
+
+ gsl_vector_long_free (kmeans->num_elements_groups);
+
+ gsl_permutation_free (kmeans->group_order);
+
+ free (kmeans);
+}
+
+static double
+diff_matrix (const gsl_matrix *m1, const gsl_matrix *m2)
+{
+ double max_diff = -INFINITY;
+ for (size_t i = 0; i < m1->size1; ++i)
+ {
+ double diff = 0;
+ for (size_t j = 0; j < m1->size2; ++j)
+ diff += pow2 (gsl_matrix_get (m1,i,j) - gsl_matrix_get (m2,i,j));
+ if (diff > max_diff)
+ max_diff = diff;
+ }
+
+ return max_diff;
+}
+
+
+
+static double
+matrix_mindist (const gsl_matrix *m, int *mn, int *mm)
+{
+ double mindist = INFINITY;
+ for (size_t i = 0; i + 1 < m->size1; ++i)
+ for (size_t j = i + 1; j < m->size1; ++j)
+ {
+ double diff_sq = 0;
+ for (size_t k = 0; k < m->size2; ++k)
+ diff_sq += pow2 (gsl_matrix_get (m, j, k) - gsl_matrix_get (m, i, k));
+ if (diff_sq < mindist)
+ {
+ mindist = diff_sq;
+ if (mn)
+ *mn = i;
+ if (mm)
+ *mm = j;
+ }
+ }
+ return mindist;
+}
+
+/* Return the distance of C from the group whose index is WHICH */
+static double
+dist_from_case (const struct Kmeans *kmeans, const struct ccase *c,
+ const struct qc *qc, int which)
+{
+ double dist = 0;
+ for (size_t j = 0; j < qc->n_vars; j++)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ assert (!(var_is_value_missing (qc->vars[j], val) & qc->exclude));
+ dist += pow2 (gsl_matrix_get (kmeans->centers, which, j) - val->f);
+ }
+
+ return dist;
+}
+
+/* Return the minimum distance of the group WHICH and all other groups */
+static double
+min_dist_from (const struct Kmeans *kmeans, const struct qc *qc, int which)
+{
+ double mindist = INFINITY;
+ for (size_t i = 0; i < qc->ngroups; i++)
+ {
+ if (i == which)
+ continue;
+
+ double dist = 0;
+ for (size_t j = 0; j < qc->n_vars; j++)
+ dist += pow2 (gsl_matrix_get (kmeans->centers, i, j)
+ - gsl_matrix_get (kmeans->centers, which, j));
+
+ if (dist < mindist)
+ mindist = dist;
+ }
+
+ return mindist;
+}
+
+/* Calculate the initial cluster centers. */
+static void
+kmeans_initial_centers (struct Kmeans *kmeans,
+ const struct casereader *reader,
+ const struct qc *qc)
+{
+ int nc = 0;
+
+ struct casereader *cs = casereader_clone (reader);
+ struct ccase *c;
+ for (; (c = casereader_read (cs)) != NULL; case_unref (c))
+ {
+ bool missing = false;
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ if (var_is_value_missing (qc->vars[j], val) & qc->exclude)
+ {
+ missing = true;
+ break;
+ }
+
+ if (nc < qc->ngroups)
+ gsl_matrix_set (kmeans->centers, nc, j, val->f);
+ }
+ if (missing)
+ continue;
+
+ if (nc++ < qc->ngroups)
+ continue;
+
+ if (qc->initial)
+ {
+ int mn, mm;
+ double m = matrix_mindist (kmeans->centers, &mn, &mm);
+
+ int mq, mp;
+ double delta;
+ kmeans_get_nearest_group (kmeans, c, qc, &mq, &delta, &mp, NULL);
+ if (delta > m)
+ /* If the distance between C and the nearest group, is greater than the distance
+ between the two groups which are clostest to each
+ other, then one group must be replaced. */
+ {
+ /* Out of mn and mm, which is the clostest of the two groups to C ? */
+ int which = (dist_from_case (kmeans, c, qc, mn)
+ > dist_from_case (kmeans, c, qc, mm)) ? mm : mn;
+
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ gsl_matrix_set (kmeans->centers, which, j, val->f);
+ }
+ }
+ else if (dist_from_case (kmeans, c, qc, mp) > min_dist_from (kmeans, qc, mq))
+ /* If the distance between C and the second nearest group
+ (MP) is greater than the smallest distance between the
+ nearest group (MQ) and any other group, then replace
+ MQ with C. */
+ {
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ gsl_matrix_set (kmeans->centers, mq, j, val->f);
+ }
+ }
+ }
+ }
+
+ casereader_destroy (cs);
+
+ kmeans->convergence_criteria = qc->epsilon * matrix_mindist (kmeans->centers, NULL, NULL);
+
+ /* As it is the first iteration, the variable kmeans->initial_centers is NULL
+ and it is created once for reporting issues. */
+ kmeans->initial_centers = gsl_matrix_alloc (qc->ngroups, qc->n_vars);
+ gsl_matrix_memcpy (kmeans->initial_centers, kmeans->centers);
+}
+
+/* Return the index of the group which is nearest to the case C */
+static void
+kmeans_get_nearest_group (const struct Kmeans *kmeans, struct ccase *c,
+ const struct qc *qc, int *g_q, double *delta_q,
+ int *g_p, double *delta_p)
+{
+ int result0 = -1;
+ int result1 = -1;
+ double mindist0 = INFINITY;
+ double mindist1 = INFINITY;
+ for (size_t i = 0; i < qc->ngroups; i++)
+ {
+ double dist = 0;
+ for (size_t j = 0; j < qc->n_vars; j++)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ if (var_is_value_missing (qc->vars[j], val) & qc->exclude)
+ continue;
+
+ dist += pow2 (gsl_matrix_get (kmeans->centers, i, j) - val->f);
+ }
+
+ if (dist < mindist0)
+ {
+ mindist1 = mindist0;
+ result1 = result0;
+
+ mindist0 = dist;
+ result0 = i;
+ }
+ else if (dist < mindist1)
+ {
+ mindist1 = dist;
+ result1 = i;
+ }
+ }
+
+ if (delta_q)
+ *delta_q = mindist0;
+
+ if (g_q)
+ *g_q = result0;
+
+ if (delta_p)
+ *delta_p = mindist1;
+
+ if (g_p)
+ *g_p = result1;
+}
+
+static void
+kmeans_order_groups (struct Kmeans *kmeans, const struct qc *qc)
+{
+ gsl_vector *v = gsl_vector_alloc (qc->ngroups);
+ gsl_matrix_get_col (v, kmeans->centers, 0);
+ gsl_sort_vector_index (kmeans->group_order, v);
+ gsl_vector_free (v);
+}
+
+/* Main algorithm.
+ Does iterations, checks convergency. */
+static void
+kmeans_cluster (struct Kmeans *kmeans, struct casereader *reader,
+ const struct qc *qc)
+{
+ kmeans_initial_centers (kmeans, reader, qc);
+
+ gsl_matrix_memcpy (kmeans->updated_centers, kmeans->centers);
+ for (int xx = 0; xx < qc->maxiter; ++xx)
+ {
+ gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
+
+ kmeans->n = 0;
+ if (qc->update)
+ {
+ struct casereader *r = casereader_clone (reader);
+ struct ccase *c;
+ for (; (c = casereader_read (r)) != NULL; case_unref (c))
+ {
+ bool missing = false;
+ for (size_t j = 0; j < qc->n_vars; j++)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ if (var_is_value_missing (qc->vars[j], val) & qc->exclude)
+ missing = true;
+ }
+ if (missing)
+ continue;
+
+ double mindist = INFINITY;
+ int group = -1;
+ for (size_t g = 0; g < qc->ngroups; ++g)
+ {
+ double d = dist_from_case (kmeans, c, qc, g);
+
+ if (d < mindist)
+ {
+ mindist = d;
+ group = g;
+ }
+ }
+
+ long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
+ *n += qc->wv ? case_num (c, qc->wv) : 1.0;
+ kmeans->n++;
+
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ if (var_is_value_missing (qc->vars[j], val) & qc->exclude)
+ continue;
+ double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
+ *x += val->f * (qc->wv ? case_num (c, qc->wv) : 1.0);
+ }
+ }
+
+ casereader_destroy (r);
+ }
+
+ /* Divide the cluster sums by the number of items in each cluster */
+ for (size_t g = 0; g < qc->ngroups; ++g)
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
+ double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
+ *x /= n + 1; // Plus 1 for the initial centers
+ }
+ gsl_matrix_memcpy (kmeans->centers, kmeans->updated_centers);
+
+ kmeans->n = 0;
+ /* Step 3 */
+ gsl_vector_long_set_all (kmeans->num_elements_groups, 0.0);
+ gsl_matrix_set_all (kmeans->updated_centers, 0.0);
+ struct ccase *c;
+ struct casereader *cs = casereader_clone (reader);
+ for (; (c = casereader_read (cs)) != NULL; case_unref (c))
+ {
+ int group = -1;
+ kmeans_get_nearest_group (kmeans, c, qc, &group, NULL, NULL, NULL);
+
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ const union value *val = case_data (c, qc->vars[j]);
+ if (var_is_value_missing (qc->vars[j], val) & qc->exclude)
+ continue;
+
+ double *x = gsl_matrix_ptr (kmeans->updated_centers, group, j);
+ *x += val->f;
+ }
+
+ long *n = gsl_vector_long_ptr (kmeans->num_elements_groups, group);
+ *n += qc->wv ? case_num (c, qc->wv) : 1.0;
+ kmeans->n++;
+ }
+ casereader_destroy (cs);
+
+ /* Divide the cluster sums by the number of items in each cluster */
+ for (size_t g = 0; g < qc->ngroups; ++g)
+ for (size_t j = 0; j < qc->n_vars; ++j)
+ {
+ long n = gsl_vector_long_get (kmeans->num_elements_groups, g);
+ double *x = gsl_matrix_ptr (kmeans->updated_centers, g, j);
+ *x /= n;
+ }
+
+ double d = diff_matrix (kmeans->updated_centers, kmeans->centers);
+ if (d < kmeans->convergence_criteria)
+ break;
+
+ if (!qc->update)
+ break;
+ }
+}
+
+/* Reports centers of clusters.
+ Initial parameter is optional for future use.
+ If initial is true, initial cluster centers are reported. Otherwise,
+ resulted centers are reported. */
+static void
+quick_cluster_show_centers (struct Kmeans *kmeans, bool initial, const struct qc *qc)
+{
+ struct pivot_table *table
+ = pivot_table_create (initial
+ ? N_("Initial Cluster Centers")
+ : N_("Final Cluster Centers"));
+
+ struct pivot_dimension *clusters
+ = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Cluster"));
+
+ clusters->root->show_label = true;
+ for (size_t i = 0; i < qc->ngroups; i++)
+ pivot_category_create_leaf (clusters->root,
+ pivot_value_new_integer (i + 1));
+
+ struct pivot_dimension *variables
+ = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Variable"));
+
+ for (size_t i = 0; i < qc->n_vars; i++)
+ pivot_category_create_leaf (variables->root,
+ pivot_value_new_variable (qc->vars[i]));
+
+ const gsl_matrix *matrix = (initial
+ ? kmeans->initial_centers
+ : kmeans->centers);
+ for (size_t i = 0; i < qc->ngroups; i++)
+ for (size_t j = 0; j < qc->n_vars; j++)
+ {
+ double x = gsl_matrix_get (matrix, kmeans->group_order->data[i], j);
+ union value v = { .f = x };
+ pivot_table_put2 (table, i, j,
+ pivot_value_new_var_value (qc->vars[j], &v));
+ }
+
+ pivot_table_submit (table);
+}
+
+
+/* A transformation function which juxtaposes the dataset with the
+ (pre-prepared) dataset containing membership and/or distance
+ values. */
+static enum trns_result
+save_trans_func (void *aux, struct ccase **c, casenumber x UNUSED)
+{
+ const struct save_trans_data *std = aux;
+ struct ccase *ca = casereader_read (std->appending_reader);
+ if (ca == NULL)
+ return TRNS_CONTINUE;
+
+ *c = case_unshare (*c);
+
+ if (std->membership_case_idx >= 0)
+ *case_num_rw (*c, std->membership) = case_num_idx (ca, std->membership_case_idx);
+
+ if (std->distance_case_idx >= 0)
+ *case_num_rw (*c, std->distance) = case_num_idx (ca, std->distance_case_idx);
+
+ case_unref (ca);
+
+ return TRNS_CONTINUE;
+}
+
+/* Free the resources of the transformation. */
+static bool
+save_trans_destroy (void *aux)
+{
+ struct save_trans_data *std = aux;
+ casereader_destroy (std->appending_reader);
+ free (std);
+ return true;
+}
+
+/* Reports cluster membership for each case, and is requested saves the
+ membership and the distance of the case from the cluster centre. */
+static void
+quick_cluster_show_membership (struct Kmeans *kmeans,
+ const struct casereader *reader,
+ struct qc *qc)
+{
+ struct pivot_table *table = NULL;
+ struct pivot_dimension *cases = NULL;
+ if (qc->print_cluster_membership)
+ {
+ table = pivot_table_create (N_("Cluster Membership"));
+
+ pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Cluster"),
+ N_("Cluster"));
+
+ cases
+ = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Case Number"));
+
+ cases->root->show_label = true;
+ }
+
+ gsl_permutation *ip = gsl_permutation_alloc (qc->ngroups);
+ gsl_permutation_inverse (ip, kmeans->group_order);
+
+ struct caseproto *proto = caseproto_create ();
+ if (qc->save_membership || qc->save_distance)
+ {
+ /* Prepare data which may potentially be used in a
+ transformation appending new variables to the active
+ dataset. */
+ int idx = 0;
+ int membership_case_idx = -1;
+ if (qc->save_membership)
+ {
+ proto = caseproto_add_width (proto, 0);
+ membership_case_idx = idx++;
+ }
+
+ int distance_case_idx = -1;
+ if (qc->save_distance)
+ {
+ proto = caseproto_add_width (proto, 0);
+ distance_case_idx = idx++;
+ }
+
+ qc->save_trans_data = xmalloc (sizeof *qc->save_trans_data);
+ *qc->save_trans_data = (struct save_trans_data) {
+ .membership_case_idx = membership_case_idx,
+ .distance_case_idx = distance_case_idx,
+ .writer = autopaging_writer_create (proto),
+ };
+ }
+
+ struct casereader *cs = casereader_clone (reader);
+ struct ccase *c;
+ for (int i = 0; (c = casereader_read (cs)) != NULL; i++, case_unref (c))
+ {
+ assert (i < kmeans->n);
+ int clust;
+ kmeans_get_nearest_group (kmeans, c, qc, &clust, NULL, NULL, NULL);
+ int cluster = ip->data[clust];
+
+ if (qc->save_trans_data)
+ {
+ /* Calculate the membership and distance values. */
+ struct ccase *outc = case_create (proto);
+ if (qc->save_membership)
+ *case_num_rw_idx (outc, qc->save_trans_data->membership_case_idx) = cluster + 1;
+
+ if (qc->save_distance)
+ *case_num_rw_idx (outc, qc->save_trans_data->distance_case_idx)
+ = sqrt (dist_from_case (kmeans, c, qc, clust));
+
+ casewriter_write (qc->save_trans_data->writer, outc);
+ }
+
+ if (qc->print_cluster_membership)
+ {
+ /* Print the cluster membership to the table. */
+ int case_idx = pivot_category_create_leaf (cases->root,
+ pivot_value_new_integer (i + 1));
+ pivot_table_put2 (table, 0, case_idx,
+ pivot_value_new_integer (cluster + 1));
+ }
+ }
+
+ caseproto_unref (proto);
+ gsl_permutation_free (ip);
+
+ if (qc->print_cluster_membership)
+ pivot_table_submit (table);
+ casereader_destroy (cs);
+}
+
+
+/* Reports number of cases of each single cluster. */
+static void
+quick_cluster_show_number_cases (struct Kmeans *kmeans, const struct qc *qc)
+{
+ struct pivot_table *table
+ = pivot_table_create (N_("Number of Cases in each Cluster"));
+
+ pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
+ N_("Count"));
+
+ struct pivot_dimension *clusters
+ = pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Clusters"));
+
+ struct pivot_category *group
+ = pivot_category_create_group (clusters->root, N_("Cluster"));
+
+ long int total = 0;
+ for (int i = 0; i < qc->ngroups; i++)
+ {
+ int cluster_idx
+ = pivot_category_create_leaf (group, pivot_value_new_integer (i + 1));
+ int count = kmeans->num_elements_groups->data [kmeans->group_order->data[i]];
+ pivot_table_put2 (table, 0, cluster_idx, pivot_value_new_integer (count));
+ total += count;
+ }
+
+ int cluster_idx = pivot_category_create_leaf (clusters->root,
+ pivot_value_new_text (N_("Valid")));
+ pivot_table_put2 (table, 0, cluster_idx, pivot_value_new_integer (total));
+ pivot_table_submit (table);
+}
+
+/* Reports. */
+static void
+quick_cluster_show_results (struct Kmeans *kmeans, const struct casereader *reader,
+ struct qc *qc)
+{
+ kmeans_order_groups (kmeans, qc); /* what does this do? */
+
+ if (qc->print_initial_clusters)
+ quick_cluster_show_centers (kmeans, true, qc);
+ quick_cluster_show_centers (kmeans, false, qc);
+ quick_cluster_show_number_cases (kmeans, qc);
+
+ quick_cluster_show_membership (kmeans, reader, qc);
+}
+
+/* Parse the QUICK CLUSTER command and populate QC accordingly.
+ Returns false on error. */
+static bool
+quick_cluster_parse (struct lexer *lexer, struct qc *qc)
+{
+ if (!parse_variables_const (lexer, qc->dict, &qc->vars, &qc->n_vars,
+ PV_NO_DUPLICATE | PV_NUMERIC))
+ return false;
+
+ while (lex_token (lexer) != T_ENDCMD)
+ {
+ lex_match (lexer, T_SLASH);
+
+ if (lex_match_id (lexer, "MISSING"))
+ {
+ lex_match (lexer, T_EQUALS);
+ while (lex_token (lexer) != T_ENDCMD
+ && lex_token (lexer) != T_SLASH)
+ {
+ if (lex_match_id (lexer, "LISTWISE")
+ || lex_match_id (lexer, "DEFAULT"))
+ qc->missing_type = MISS_LISTWISE;
+ else if (lex_match_id (lexer, "PAIRWISE"))
+ qc->missing_type = MISS_PAIRWISE;
+ else if (lex_match_id (lexer, "INCLUDE"))
+ qc->exclude = MV_SYSTEM;
+ else if (lex_match_id (lexer, "EXCLUDE"))
+ qc->exclude = MV_ANY;
+ else
+ {
+ lex_error_expecting (lexer, "LISTWISE", "DEFAULT",
+ "PAIRWISE", "INCLUDE", "EXCLUDE");
+ return false;
+ }
+ }
+ }
+ else if (lex_match_id (lexer, "PRINT"))
+ {
+ lex_match (lexer, T_EQUALS);
+ while (lex_token (lexer) != T_ENDCMD
+ && lex_token (lexer) != T_SLASH)
+ {
+ if (lex_match_id (lexer, "CLUSTER"))
+ qc->print_cluster_membership = true;
+ else if (lex_match_id (lexer, "INITIAL"))
+ qc->print_initial_clusters = true;
+ else
+ {
+ lex_error_expecting (lexer, "CLUSTER", "INITIAL");
+ return false;
+ }
+ }
+ }
+ else if (lex_match_id (lexer, "SAVE"))
+ {
+ lex_match (lexer, T_EQUALS);
+ while (lex_token (lexer) != T_ENDCMD
+ && lex_token (lexer) != T_SLASH)
+ {
+ if (lex_match_id (lexer, "CLUSTER"))
+ {
+ qc->save_membership = true;
+ if (lex_match (lexer, T_LPAREN))
+ {
+ if (!lex_force_id (lexer))
+ return false;
+
+ free (qc->var_membership);
+ qc->var_membership = xstrdup (lex_tokcstr (lexer));
+ if (NULL != dict_lookup_var (qc->dict, qc->var_membership))
+ {
+ lex_error (lexer,
+ _("A variable called `%s' already exists."),
+ qc->var_membership);
+ free (qc->var_membership);
+ qc->var_membership = NULL;
+ return false;
+ }
+
+ lex_get (lexer);
+
+ if (!lex_force_match (lexer, T_RPAREN))
+ return false;
+ }
+ }
+ else if (lex_match_id (lexer, "DISTANCE"))
+ {
+ qc->save_distance = true;
+ if (lex_match (lexer, T_LPAREN))
+ {
+ if (!lex_force_id (lexer))
+ return false;
+
+ free (qc->var_distance);
+ qc->var_distance = xstrdup (lex_tokcstr (lexer));
+ if (NULL != dict_lookup_var (qc->dict, qc->var_distance))
+ {
+ lex_error (lexer,
+ _("A variable called `%s' already exists."),
+ qc->var_distance);
+ free (qc->var_distance);
+ qc->var_distance = NULL;
+ return false;
+ }
+
+ lex_get (lexer);
+
+ if (!lex_force_match (lexer, T_RPAREN))
+ return false;
+ }
+ }
+ else
+ {
+ lex_error_expecting (lexer, "CLUSTER", "DISTANCE");
+ return false;
+ }
+ }
+ }
+ else if (lex_match_id (lexer, "CRITERIA"))
+ {
+ lex_match (lexer, T_EQUALS);
+ while (lex_token (lexer) != T_ENDCMD
+ && lex_token (lexer) != T_SLASH)
+ {
+ if (lex_match_id (lexer, "CLUSTERS"))
+ {
+ if (!lex_force_match (lexer, T_LPAREN)
+ || !lex_force_int_range (lexer, "CLUSTERS", 1, INT_MAX))
+ return false;
+ qc->ngroups = lex_integer (lexer);
+ lex_get (lexer);
+ if (!lex_force_match (lexer, T_RPAREN))
+ return false;
+ }
+ else if (lex_match_id (lexer, "CONVERGE"))
+ {
+ if (!lex_force_match (lexer, T_LPAREN)
+ || !lex_force_num_range_open (lexer, "CONVERGE",
+ 0, DBL_MAX))
+ return false;
+ qc->epsilon = lex_number (lexer);
+ lex_get (lexer);
+ if (!lex_force_match (lexer, T_RPAREN))
+ return false;
+ }
+ else if (lex_match_id (lexer, "MXITER"))
+ {
+ if (!lex_force_match (lexer, T_LPAREN)
+ || !lex_force_int_range (lexer, "MXITER", 1, INT_MAX))
+ return false;
+ qc->maxiter = lex_integer (lexer);
+ lex_get (lexer);
+ if (!lex_force_match (lexer, T_RPAREN))
+ return false;
+ }
+ else if (lex_match_id (lexer, "NOINITIAL"))
+ qc->initial = false;
+ else if (lex_match_id (lexer, "NOUPDATE"))
+ qc->update = false;
+ else
+ {
+ lex_error_expecting (lexer, "CLUSTERS", "CONVERGE", "MXITER",
+ "NOINITIAL", "NOUPDATE");
+ return false;
+ }
+ }
+ }
+ else
+ {
+ lex_error_expecting (lexer, "MISSING", "PRINT", "SAVE", "CRITERIA");
+ return false;
+ }
+ }
+ return true;
+}
+
+int
+cmd_quick_cluster (struct lexer *lexer, struct dataset *ds)
+{
+ struct qc qc = {
+ .dataset = ds,
+ .dict = dataset_dict (ds),
+ .ngroups = 2,
+ .maxiter = 10,
+ .epsilon = DBL_EPSILON,
+ .missing_type = MISS_LISTWISE,
+ .exclude = MV_ANY,
+ .initial = true,
+ .update = true,
+ };
+
+ if (!quick_cluster_parse (lexer, &qc))
+ goto error;
+
+ qc.wv = dict_get_weight (qc.dict);
+
+ struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), qc.dict);
+ struct casereader *group;
+ while (casegrouper_get_next_group (grouper, &group))
+ {
+ if (qc.missing_type == MISS_LISTWISE)
+ group = casereader_create_filter_missing (group, qc.vars, qc.n_vars,
+ qc.exclude, NULL, NULL);
+
+ struct Kmeans *kmeans = kmeans_create (&qc);
+ kmeans_cluster (kmeans, group, &qc);
+ quick_cluster_show_results (kmeans, group, &qc);
+ kmeans_destroy (kmeans);
+ casereader_destroy (group);
+ }
+ bool ok = casegrouper_destroy (grouper);
+ ok = proc_commit (ds) && ok;
+
+ /* If requested, set a transformation to append the cluster and
+ distance values to the current dataset. */
+ if (qc.save_trans_data)
+ {
+ struct save_trans_data *std = qc.save_trans_data;
+
+ std->appending_reader = casewriter_make_reader (std->writer);
+
+ if (qc.save_membership)
+ {
+ /* Invent a variable name if necessary. */
+ int idx = 0;
+ struct string name;
+ ds_init_empty (&name);
+ while (qc.var_membership == NULL)
+ {
+ ds_clear (&name);
+ ds_put_format (&name, "QCL_%d", idx++);
+
+ if (!dict_lookup_var (qc.dict, ds_cstr (&name)))
+ {
+ qc.var_membership = strdup (ds_cstr (&name));
+ break;
+ }
+ }
+ ds_destroy (&name);
+
+ std->membership = dict_create_var_assert (qc.dict, qc.var_membership, 0);
+ }
+
+ if (qc.save_distance)
+ {
+ /* Invent a variable name if necessary. */
+ int idx = 0;
+ struct string name;
+ ds_init_empty (&name);
+ while (qc.var_distance == NULL)
+ {
+ ds_clear (&name);
+ ds_put_format (&name, "QCL_%d", idx++);
+
+ if (!dict_lookup_var (qc.dict, ds_cstr (&name)))
+ {
+ qc.var_distance = strdup (ds_cstr (&name));
+ break;
+ }
+ }
+ ds_destroy (&name);
+
+ std->distance = dict_create_var_assert (qc.dict, qc.var_distance, 0);
+ }
+
+ static const struct trns_class trns_class = {
+ .name = "QUICK CLUSTER",
+ .execute = save_trans_func,
+ .destroy = save_trans_destroy,
+ };
+ add_transformation (qc.dataset, &trns_class, std);
+ }
+
+ free (qc.var_distance);
+ free (qc.var_membership);
+ free (qc.vars);
+ return ok;
+
+ error:
+ free (qc.var_distance);
+ free (qc.var_membership);
+ free (qc.vars);
+ return CMD_FAILURE;
+}
projects / pspp / blobdiff