kmeans.sql_in

Go to the documentation of this file.

/* ----------------------------------------------------------------------- *//**
 *
 * @file kmeans.sql_in
 *
 * @brief Set of functions for k-means clustering.
 *
 * @sa For a brief introduction to k-means clustering, see the module
 *     description \ref grp_kmeans.
 *
 *//* ----------------------------------------------------------------------- */

m4_include(`SQLCommon.m4')

/**
@addtogroup grp_kmeans

@about

Clustering refers to the problem of partitioning a set of objects according to
some problem-dependent measure of <em>similarity</em>. In the k-means variant,
one is given \f$ n \f$ points \f$ x_1, \dots, x_n \in \mathbb R^d \f$, and the
goal is to position \f$ k \f$ centroids \f$ c_1, \dots, c_k \in \mathbb R^d \f$
so that the sum of \em distances between each point and its closest centroid is
minimized. Each centroid represents a cluster that consists of all points to
which this centroid is closest. Formally, we wish to minimize the following
objective function:
\f[
    (c_1, \dots, c_k) \mapsto \sum_{i=1}^n \min_{j=1}^k \operatorname{dist}(x_i, c_j)
\f]
In the most common case, \f$ \operatorname{dist} \f$ is the square of the
Euclidean distance.

This problem is computationally difficult (NP-hard), yet the
local-search heuristic proposed by Lloyd [4] performs reasonably well in
practice. In fact, it is so ubiquitous today that it is
often referred to as the <em>standard algorithm</em> or even just the
<em>k-means algorithm</em> [1]. It works as follows:

-# Seed the \f$ k \f$ centroids (see below)
-# Repeat until convergence:
 -# Assign each point to its closest centroid
 -# Move each centroid to a position that minimizes the sum of distances in this
    cluster
-# Convergence is achieved when no points change their assignments during step
   2a.

Since the objective function decreases in every step, this algorithm is
guaranteed to converge to a local optimum.

@implementation

Data points and predefined centroids (if used) are expected to be stored row-wise,
in a column of type <tt>\ref grp_svec "SVEC"</tt> (or any type convertible to
<tt>\ref grp_svec "SVEC"</tt>, like <tt>FLOAT[]</tt> or <tt>INTEGER[]</tt>).
Data points with non-finite values (NULL, NaN, infinity) in any component
will be skipped during analysis.

The following methods are available for the centroid seeding:
 - <strong>random selection</strong>:
   Select \f$ k \f$ centroids randomly among the input points.
 - <strong>kmeans++</strong> [2]:
   Start with a single centroid chosen randomly among the input points. Then
   iteratively choose new
   centroids from the input points until there is a total of \f$ k \f$
   centroids. The probability for picking a particular point is proportional to
   its minimum distance to any existing centroid.
   \n
   Intuitively, kmeans++ favors seedings where centroids are spread out over the
   whole range of the input points, while at the same time not being too
   susceptible to outliers [2].
 - <strong>user-specified set of initial centroids</strong>:
   See below for a description of the expected format of the set of initial
   centroids.

The following distance functions can be used (computation of barycenter/mean in parentheses):
 - <strong>\ref dist_norm1</strong>:  1-norm/Manhattan (element-wise median [Note that MADlib does not provide a median
   aggregate function for support and performance reasons.])
 - <strong>\ref dist_norm2</strong>: 2-norm/Euclidean (element-wise mean)
 - <strong>\ref squared_dist_norm2</strong>: squared Euclidean distance (element-wise mean)
 - <strong>\ref dist_angle</strong>: angle (element-wise mean of normalized points)
 - <strong>\ref dist_tanimoto</strong>: tanimoto (element-wise mean of normalized points [5])
 - <strong>user defined function</strong> with signature DOUBLE PRECISION[] x DOUBLE PRECISION[] -> DOUBLE PRECISION

The following aggregate functions for determining centroids can be used:
 - <strong>\ref avg</strong>: average
 - <strong>\ref normalized_avg</strong>: normalized average

The algorithm stops when one of the following conditions is met:
 - The fraction of updated points is smaller than the convergence threshold (default: 0.001).
 - The algorithm reaches the maximum number of allowed iterations (default: 20).

A popular method to assess the quality of the clustering is the
<em>silhouette coefficient</em>, a simplified version of which is provided as part of the k-means module.
Note that for large data sets, this computation is expensive.

@input
The <strong>source relation</strong> is expected to be of the following form (or
to be implicitly convertible into the following form):
<pre>{TABLE|VIEW} <em>rel_source</em> (
    ...
    <em>expr_points</em> FLOAT8[],
    ...
)</pre>
where:
 - <em>expr_points</em> is the name of a column with point coordinates.
   Types such as \c svec or <tt>INTEGER[]</tt> are implicitly converted to
   <tt>FLOAT8[]</tt>.

If kmeans is called with a set of initial centroids, the centroid relation is
expected to be of the following form:
<pre>{TABLE|VIEW} <em>rel_initial_centroids</em> (
    ...
    <em>expr_centroid</em> DOUBLE PRECISION[],
    ...
)</pre>
where:
 - <em>expr_centroid</em> is the name of a column with coordinates.

@usage
The k-means algorithm can be invoked in four possible ways:

- using <em>random</em> centroid seeding method for a
provided \f$ k \f$:
<pre>SELECT * FROM \ref kmeans_random(
  '<em>rel_source</em>', '<em>expr_point</em>', k,
  [ '<em>fn_dist</em>', '<em>agg_centroid</em>',
  <em>max_num_iterations</em>, <em>min_frac_reassigned</em> ]
);</pre>

- using <em>kmeans++</em> centroid seeding method for a
provided \f$ k \f$:
<pre>SELECT * FROM \ref kmeanspp(
  '<em>rel_source</em>', '<em>expr_point</em>', k,
  [ '<em>fn_dist</em>', '<em>agg_centroid</em>',
  <em>max_num_iterations</em>, <em>min_frac_reassigned</em> ]
);</pre>

- with a provided centroid set:
<pre>SELECT * FROM \ref kmeans(
  '<em>rel_source</em>', '<em>expr_point</em>',
  '<em>rel_initial_centroids</em>', '<em>expr_centroid</em>',
  [ '<em>fn_dist</em>', '<em>agg_centroid</em>',
  <em>max_num_iterations</em>, <em>min_frac_reassigned</em> ]
);</pre>
------------ OR ---------------
<pre>SELECT * FROM \ref kmeans(
  '<em>rel_source</em>', '<em>expr_point</em>',
  initial_centroids,
  [ '<em>fn_dist</em>', '<em>agg_centroid</em>',
  <em>max_num_iterations</em>, <em>min_frac_reassigned</em> ]
);</pre>
where:
 - <em>initial_centroids</em> is of type <tt>DOUBLE PRECISION[][]</tt>.

The output of the k-means module is a table that includes the final
centroid positions (DOUBLE PRECISION[][]), the objective function,
the fraction of reassigned points in the last iteration, and
the number of total iterations:
<pre>
            centroids             |   objective_fn   | frac_reassigned | num_iterations
----------------------------------+------------------+-----------------+----------------
             ...
</pre>

@examp

-#  Prepare some input data.
\code
sql> SELECT * FROM public.km_sample LIMIT 5;
                  points
-------------------------------------------
 {1,1}:{15.8822241332382,105.945462542586}
 {1,1}:{34.5065216883086,72.3126099305227}
 {1,1}:{22.5074400822632,95.3209559689276}
 {1,1}:{70.2589857042767,68.7395178806037}
 {1,1}:{30.9844257542863,25.3213323024102}
(5 rows)
\endcode
Note: the example <em>points</em> is type <tt>\ref grp_svec "SVEC"</tt>.
-#  Run k-means clustering using kmeans++ for centroid seeding:
\code
sql> SELECT * FROM madlib.kmeanspp('km_sample', 'points', 2, 'madlib.squared_dist_norm2', 'madlib.avg', 20, 0.001);
);
                                centroids                                |   objective_fn   | frac_reassigned | num_iterations
-------------------------------------------------------------------------+------------------+-----------------+----------------
 {{68.01668579784,48.9667382972952},{28.1452167573446,84.5992507653263}} | 586729.010675982 |           0.001 |              5
\endcode
-# Calculate the simplified silhouette coefficient:
\code
sql> SELECT * from madlib.simple_silhouette('km_test_svec','points',
     (select centroids from madlib.kmeanspp('km_test_svec','points',2,'madlib.squared_dist_norm2','madlib.avg',20,0.001)),
     'madlib.dist_norm2');
 simple_silhouette
-------------------
 0.611022970398174
\endcode

@literature

[1] Wikipedia, K-means Clustering,
    http://en.wikipedia.org/wiki/K-means_clustering

[2] David Arthur, Sergei Vassilvitskii: k-means++: the advantages of careful
    seeding, Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete
    Algorithms (SODA'07), pp. 1027-1035,
    http://www.stanford.edu/~darthur/kMeansPlusPlus.pdf

[3] E. R. Hruschka, L. N. C. Silva, R. J. G. B. Campello: Clustering
    Gene-Expression Data: A Hybrid Approach that Iterates Between k-Means and
    Evolutionary Search. In: Studies in Computational Intelligence - Hybrid
    Evolutionary Algorithms. pp. 313-335. Springer. 2007.

[4] Lloyd, Stuart: Least squares quantization in PCM. Technical Note, Bell
    Laboratories. Published much later in: IEEE Transactions on Information
    Theory 28(2), pp. 128-137. 1982.

[5] Leisch, Friedrich: A Toolbox for K-Centroids Cluster Analysis.  In: Computational
    Statistics and Data Analysis, 51(2). pp. 526-544. 2006.

@sa File kmeans.sql_in documenting the SQL functions.

@internal
@sa namespace kmeans (documenting the implementation in Python)
@endinternal
*/

/*
 * @brief k-Means return type
 *
 * A composite value:
 *  - <tt>centroids</tt> - Matrix containing the new \f$ l \leq k \f$
 *    repositioned centroids as columns. If this matrix has \f$ l < k \f$
 *    columns, one or more old centroids no longer were closest to any point.
 *  - <tt>old_centroid_its</tt> - The order of the centroids in
 *    <tt>centroid</tt> is not guaranteed to be consitent across iterations.
 *    In particular, if a centroid is no longer closest to any point it can be
 *    dropped and a new centroid is added afterwards. We therefore need to map
 *    positions in <tt>centroids</tt> to the respective positions in the
 *    previous iteration.
 *  - <tt>objective_fn</tt> - Value of the objective function, i.e.,
 *    \f$ \sum_{x \in P} \dist(x, C)^2 \f$ where
 *    \f$ P \f$ is the set of points, \f$ C \f$ is the set of centroids, and
 *    \f$ \dist(x, C) := \min_{c \in C} \operatorname{dist}(x, c) \f$.
 *  - <tt>frac_reassigned</tt> - Fraction of points that was assigned a
 *    different centroid in the current iteration.
 *  - <tt>num_iterations</tt> - Number of iterations performed (so far).
 */
CREATE TYPE MADLIB_SCHEMA.kmeans_result AS (
    centroids DOUBLE PRECISION[][],
    objective_fn DOUBLE PRECISION,
    frac_reassigned DOUBLE PRECISION,
    num_iterations INTEGER
);

/*
 * @brief k-Means inter-iteration state type
 *
 * A composite value like \ref{kmeans_result}. Additional fields:
 *  - <tt>old_centroid_its</tt> - The order of the centroids in
 *    <tt>centroid</tt> is not guaranteed to be consitent across iterations.
 *    In particular, if a centroid is no longer closest to any point it can be
 *    dropped and a new centroid is added afterwards. We therefore need to map
 *    positions in <tt>centroids</tt> to the respective positions in the
 *    previous iteration.
 *  - <tt>num_iterations</tt> - Number of iterations performed (so far).
 */
CREATE TYPE MADLIB_SCHEMA.kmeans_state AS (
    centroids DOUBLE PRECISION[][],
    old_centroid_ids INTEGER[],
    objective_fn DOUBLE PRECISION,
    frac_reassigned DOUBLE PRECISION
);

/**
 * @internal
 * @brief Execute a SQL command where $1, ..., $4 are substituted with the
  *     given arguments.
 */
CREATE FUNCTION MADLIB_SCHEMA.internal_execute_using_kmeans_args(
    sql VARCHAR, DOUBLE PRECISION[][], REGPROC, INTEGER, DOUBLE PRECISION
) RETURNS VOID
VOLATILE
CALLED ON NULL INPUT
LANGUAGE c
AS 'MODULE_PATHNAME', 'exec_sql_using';

CREATE FUNCTION MADLIB_SCHEMA.internal_compute_kmeans(
    rel_args VARCHAR,
    rel_state VARCHAR,
    rel_source VARCHAR,
    expr_point VARCHAR,
    agg_centroid VARCHAR)
RETURNS INTEGER
VOLATILE
LANGUAGE plpythonu
AS $$PythonFunction(kmeans, kmeans, compute_kmeans)$$;

/**
 * @brief Perform Lloyd's k-means local-search heuristic
 *
 * @param rel_source Name of the relation containing input points
 * @param expr_point Expression evaluating to point coordinates for each tuple
 * @param initial_centroids Matrix containing the initial centroids as columns
 * @param fn_dist Name of a function with signature
 *     <tt>DOUBLE PRECISION[] x DOUBLE PRECISION[] -> DOUBLE PRECISION</tt> that
 *     returns the distance between two points. The default is the
 *     \ref squared_dist_norm2(float8[],float8[]) "squared Euclidean distance".
 * @param agg_centroid Name of an aggregate function with signature
 *     <tt>DOUBLE PRECISION[] -> DOUBLE PRECISION[]</tt> that, for each group
 *     of points, returns a centroid. In order for Lloyd's local-search
 *     heuristic to provably converge and to return a local minimum, this
 *     centroid should minimize the sum of distances between each point in the
 *     group and the centroid. The default is the
 *     \ref avg(float8[]) "average (mean/barycenter in Euclidean space)",
 *     which satisfies this property if <tt>fn_dist = 'squared_dist_norm2'</tt>.
 * @param max_num_iterations Maximum number of iterations
 * @param min_frac_reassigned Fraction of reassigned points below which
 *     convergence is assumed and the algorithm terminates
 * @returns A composite value:
 *  - <tt>centroids</tt> - Matrix with \f$ k \f$ centroids as columns.
 *  - <tt>frac_reassigned</tt> - Fraction of points that were assigned a
 *    different centroid in the last iteration.
 *  - <tt>num_iterations</tt> - The number of iterations before the
 *    algorithm terminated
 */
CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    initial_centroids DOUBLE PRECISION[][],
    fn_dist VARCHAR /*+ DEFAULT 'squared_dist_norm2' */,
    agg_centroid VARCHAR /*+ DEFAULT 'avg' */,
    max_num_iterations INTEGER /*+ DEFAULT 20 */,
    min_frac_reassigned DOUBLE PRECISION /*+ DEFAULT 0.001 */
) RETURNS MADLIB_SCHEMA.kmeans_result AS $$
DECLARE
    theIteration INTEGER;
    theResult MADLIB_SCHEMA.kmeans_result;
    oldClientMinMessages VARCHAR;
    class_rel_source REGCLASS;
    proc_fn_dist REGPROCEDURE;
    proc_agg_centroid REGPROCEDURE;
    rel_filtered VARCHAR;
    num_points INTEGER;
    k INTEGER;
    centroids FLOAT8[];
BEGIN
    IF (array_upper(initial_centroids,1) IS NULL) THEN
    RAISE EXCEPTION 'No valid initial centroids given.';
    END IF;

    centroids := ARRAY(SELECT unnest(initial_centroids));
    IF (SELECT MADLIB_SCHEMA.svec_elsum(centroids)) >= 'Infinity'::float THEN
        RAISE EXCEPTION 'At least one initial centroid has non-finite values.';
    END IF;

    rel_filtered = MADLIB_SCHEMA.__filter_input_relation(rel_source, expr_point);
    class_rel_source := rel_filtered;
    proc_fn_dist := fn_dist
        || '(DOUBLE PRECISION[], DOUBLE PRECISION[])';
    IF (SELECT prorettype != 'DOUBLE PRECISION'::regtype OR proisagg = TRUE
        FROM pg_proc WHERE oid = proc_fn_dist) THEN
        RAISE EXCEPTION 'Distance function has wrong signature or is not a simple function.';
    END IF;
    proc_agg_centroid := agg_centroid || '(DOUBLE PRECISION[])';
    IF (SELECT prorettype != 'DOUBLE PRECISION[]'::regtype OR proisagg = FALSE
        FROM pg_proc WHERE oid = proc_agg_centroid) THEN
        RAISE EXCEPTION 'Mean aggregate has wrong signature or is not an aggregate.';
    END IF;
    IF (min_frac_reassigned < 0) OR (min_frac_reassigned > 1) THEN
        RAISE EXCEPTION 'Convergence threshold is not a valid value (must be a fraction between 0 and 1).';
    END IF;
    IF (max_num_iterations < 0) THEN
        RAISE EXCEPTION 'Number of iterations must be a non-negative integer.';
    END IF;

    -- Extra parameter check added so that ERROR output is more user-readable (doesn't include Python traceback)
    k := array_upper(initial_centroids,1);
    IF (k <= 0) THEN
        RAISE EXCEPTION 'Number of clusters k must be a positive integer.';
    END IF;
    IF (k > 32767) THEN
    RAISE EXCEPTION 'Number of clusters k must be <= 32767 (for results to be returned in a reasonable amount of time).';
    END IF;
    EXECUTE $sql$ SELECT count(*) FROM $sql$ || textin(regclassout(class_rel_source)) INTO num_points ;
    IF (num_points < k) THEN
    RAISE EXCEPTION 'Number of centroids is greater than number of points.';
    END IF;

    -- We first setup the argument table. Rationale: We want to avoid all data
    -- conversion between native types and Python code. Instead, we use Python
    -- as a pure driver layer.
    PERFORM MADLIB_SCHEMA.create_schema_pg_temp();
    oldClientMinMessages :=
        (SELECT setting FROM pg_settings WHERE name = 'client_min_messages');
    EXECUTE 'SET client_min_messages TO warning';

    -- Unfortunately, the EXECUTE USING syntax is only available starting
    -- PostgreSQL 8.4:
    -- http://www.postgresql.org/docs/8.4/static/plpgsql-statements.html#PLPGSQL-STATEMENTS-EXECUTING-DYN
    -- We therefore have to emulate.
    PERFORM MADLIB_SCHEMA.internal_execute_using_kmeans_args($sql$
        DROP TABLE IF EXISTS pg_temp._madlib_kmeans_args;
        CREATE TABLE pg_temp._madlib_kmeans_args AS
        SELECT
            $1 AS initial_centroids, array_upper($1, 1) AS k,
            $2 AS fn_dist, $3 AS max_num_iterations,
            $4 AS min_frac_reassigned;
        $sql$,
        initial_centroids, proc_fn_dist, max_num_iterations,
        min_frac_reassigned);
    EXECUTE 'SET client_min_messages TO ' || oldClientMinMessages;

    -- Perform acutal computation.
    -- Unfortunately, Greenplum and PostgreSQL <= 8.2 do not have conversion
    -- operators from regclass to varchar/text.
    theIteration := MADLIB_SCHEMA.internal_compute_kmeans('_madlib_kmeans_args',
            '_madlib_kmeans_state',
            textin(regclassout(class_rel_source)), expr_point,
            textin(regprocout(proc_agg_centroid)));

    -- Retrieve result from state table and return it
    EXECUTE
        $sql$
        SELECT (_state).centroids, (_state).objective_fn,
            (_state).frac_reassigned, NULL
        FROM _madlib_kmeans_state
        WHERE _iteration = $sql$ || theIteration || $sql$
        $sql$
        INTO theResult;
    -- The number of iterations are not updated in the C++ code. We do it here.
    IF NOT (theResult IS NULL) THEN
        theResult.num_iterations = theIteration;
    END IF;
    RETURN theResult;
END;
$$ LANGUAGE plpgsql VOLATILE;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    initial_centroids DOUBLE PRECISION[][],
    fn_dist VARCHAR,
    agg_centroid VARCHAR,
    max_num_iterations INTEGER
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans($1, $2, $3, $4, $5, $6, 0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    initial_centroids DOUBLE PRECISION[][],
    fn_dist VARCHAR,
    agg_centroid VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans($1, $2, $3, $4, $5, 20, 0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    initial_centroids DOUBLE PRECISION[][],
    fn_dist VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans($1, $2, $3, $4, 'MADLIB_SCHEMA.avg', 20,
        0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    initial_centroids DOUBLE PRECISION[][]
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans($1, $2, $3,
        'MADLIB_SCHEMA.squared_dist_norm2', 'MADLIB_SCHEMA.avg', 20, 0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.internal_execute_using_kmeanspp_seeding_args(
    sql VARCHAR, INTEGER, REGPROC, DOUBLE PRECISION[][]
) RETURNS VOID
VOLATILE
CALLED ON NULL INPUT
LANGUAGE c
AS 'MODULE_PATHNAME', 'exec_sql_using';

CREATE FUNCTION MADLIB_SCHEMA.internal_compute_kmeanspp_seeding(
    rel_args VARCHAR,
    rel_state VARCHAR,
    rel_source VARCHAR,
    expr_point VARCHAR)
RETURNS INTEGER
AS $$PythonFunction(kmeans, kmeans, compute_kmeanspp_seeding)$$
LANGUAGE plpythonu VOLATILE;

/**
 * @brief k-Means++ Seeding
 *
 * @param rel_source Name of the relation containing input points
 * @param expr_point Expression evaluating to point coordinates for each tuple
 * @param k Number of centroids
 * @param fn_dist Name of a function with signature
 *     <tt>DOUBLE PRECISION[] x DOUBLE PRECISION[] -> DOUBLE PRECISION</tt> that
 *     returns the distance between two points
 * @param initial_centroids A matrix containing up to \f$ k \f$ columns as
 *     columns. kmeanspp_seeding() proceeds exactly as if these centroids had
 *     already been generated in previous iterations. This parameter may be
 *     NULL in which all \f$ k \f$ centroids will be generated.
 * @returns A matrix containing \f$ k \f$ centroids as columns
 */
CREATE FUNCTION MADLIB_SCHEMA.kmeanspp_seeding(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR /*+ DEFAULT 'squared_dist_norm2' */,
    initial_centroids DOUBLE PRECISION[][] /*+ DEFAULT NULL */
) RETURNS DOUBLE PRECISION[][] AS $$
DECLARE
    theIteration INTEGER;
    theResult DOUBLE PRECISION[][];
    oldClientMinMessages VARCHAR;
    class_rel_source REGCLASS;
    proc_fn_dist REGPROCEDURE;
    num_points INTEGER;
    num_centroids INTEGER;
    rel_filtered VARCHAR;
BEGIN
    rel_filtered = MADLIB_SCHEMA.__filter_input_relation(rel_source, expr_point);
    class_rel_source := rel_filtered;

    IF (initial_centroids IS NOT NULL) THEN
    num_centroids := array_upper(initial_centroids,1);
    ELSE
    num_centroids := k;
    END IF;

    proc_fn_dist := fn_dist
        || '(DOUBLE PRECISION[], DOUBLE PRECISION[])';
    IF (SELECT prorettype != 'DOUBLE PRECISION'::regtype OR proisagg = TRUE
        FROM pg_proc WHERE oid = proc_fn_dist) THEN
        RAISE EXCEPTION 'Distance function has wrong signature or is not a simple function.';
    END IF;
    IF (k <= 0) THEN
        RAISE EXCEPTION 'Number of clusters k must be a positive integer.';
    END IF;
    IF (k > 32767) THEN
    RAISE EXCEPTION 'Number of clusters k must be <= 32767 (for results to be returned in a reasonable amount of time).';
    END IF;
    EXECUTE $sql$ SELECT count(*) FROM $sql$ || textin(regclassout(class_rel_source)) INTO num_points ;
    IF (num_points < k OR num_points < num_centroids) THEN
    RAISE EXCEPTION 'Number of centroids is greater than number of points.';
    END IF;
    IF (k < num_centroids) THEN
    RAISE WARNING 'Number of clusters k is less than number of supplied initial centroids. Number of final clusters will equal number of supplied initial centroids.';
    END IF;

    -- We first setup the argument table. Rationale: We want to avoid all data
    -- conversion between native types and Python code. Instead, we use Python
    -- as a pure driver layer.
    oldClientMinMessages :=
        (SELECT setting FROM pg_settings WHERE name = 'client_min_messages');
    EXECUTE 'SET client_min_messages TO warning';
    PERFORM MADLIB_SCHEMA.create_schema_pg_temp();
    -- Unfortunately, the EXECUTE USING syntax is only available starting
    -- PostgreSQL 8.4:
    -- http://www.postgresql.org/docs/8.4/static/plpgsql-statements.html#PLPGSQL-STATEMENTS-EXECUTING-DYN
    -- We therefore have to emulate.
    PERFORM MADLIB_SCHEMA.internal_execute_using_kmeanspp_seeding_args($sql$
        DROP TABLE IF EXISTS pg_temp._madlib_kmeanspp_args;
        CREATE TEMPORARY TABLE _madlib_kmeanspp_args AS
        SELECT $1 AS k, $2 AS fn_dist, $3 AS initial_centroids;
        $sql$,
        k, proc_fn_dist, initial_centroids);
    EXECUTE 'SET client_min_messages TO ' || oldClientMinMessages;

    -- Perform acutal computation.
    -- Unfortunately, Greenplum and PostgreSQL <= 8.2 do not have conversion
    -- operators from regclass to varchar/text.
    theIteration := (
        SELECT MADLIB_SCHEMA.internal_compute_kmeanspp_seeding(
            '_madlib_kmeanspp_args', '_madlib_kmeanspp_state',
            textin(regclassout(class_rel_source)), expr_point)
    );

    -- Retrieve result from state table and return it
    EXECUTE
        $sql$
        SELECT _state FROM _madlib_kmeanspp_state
        WHERE _iteration = $sql$ || theIteration || $sql$
        $sql$
        INTO theResult;
    RETURN theResult;
END;
$$ LANGUAGE plpgsql VOLATILE;

CREATE FUNCTION MADLIB_SCHEMA.kmeanspp_seeding(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR
) RETURNS DOUBLE PRECISION[][]
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3, $4, NULL)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeanspp_seeding(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER
) RETURNS DOUBLE PRECISION[][]
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3,
        'MADLIB_SCHEMA.squared_dist_norm2', NULL)
$$;

/**
 * @brief Run k-Means++.
 *
 * This is a shortcut for running k-means++. It is equivalent to
 * <pre>SELECT \ref kmeans(
    rel_source,
    expr_point,
    \ref kmeanspp_seeding(
        rel_source,
        expr_point,
        k,
        fn_dist
    ),
    fn_dist,
    agg_centroid,
    max_num_iterations,
    min_frac_reassigned
)</pre>
 */
CREATE FUNCTION MADLIB_SCHEMA.kmeanspp(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR /*+ DEFAULT 'squared_dist_norm2' */,
    agg_centroid VARCHAR /*+ DEFAULT 'avg' */,
    max_num_iterations INTEGER /*+ DEFAULT 20 */,
    min_frac_reassigned DOUBLE PRECISION /*+ DEFAULT 0.001 */
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3, $4),
        $4, $5, $6, $7);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeanspp(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR,
    agg_centroid VARCHAR,
    max_num_iterations INTEGER
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3, $4),
        $4, $5, $6, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeanspp(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR,
    agg_centroid VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3, $4),
        $4, $5, 20, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeanspp(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3, $4),
        $4, 'MADLIB_SCHEMA.avg', 20, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeanspp(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2,
        MADLIB_SCHEMA.kmeanspp_seeding($1, $2, $3,
            'MADLIB_SCHEMA.squared_dist_norm2'),
        'MADLIB_SCHEMA.squared_dist_norm2', 'MADLIB_SCHEMA.avg', 20, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.internal_execute_using_kmeans_random_seeding_args(
    sql VARCHAR, INTEGER, DOUBLE PRECISION[][]
) RETURNS VOID
VOLATILE
CALLED ON NULL INPUT
LANGUAGE c
AS 'MODULE_PATHNAME', 'exec_sql_using';

CREATE FUNCTION MADLIB_SCHEMA.internal_compute_kmeans_random_seeding(
    rel_args VARCHAR,
    rel_state VARCHAR,
    rel_source VARCHAR,
    expr_point VARCHAR)
RETURNS INTEGER
AS $$PythonFunction(kmeans, kmeans, compute_kmeans_random_seeding)$$
LANGUAGE plpythonu VOLATILE;

/**
 * @brief k-Means Random Seeding
 *
 * @param rel_source Name of the relation containing input points
 * @param expr_point Expression evaluating to point coordinates for each tuple
 * @param k Number of centroids
 * @param initial_centroids A matrix containing up to \f$ k \f$ columns as
 *     columns. kmeanspp_seeding() proceeds exactly as if these centroids had
 *     already been generated in previous iterations. This parameter may be
 *     NULL in which all \f$ k \f$ centroids will be generated.
 * @returns A matrix containing \f$ k \f$ centroids as columns
 */
CREATE FUNCTION MADLIB_SCHEMA.kmeans_random_seeding(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    initial_centroids DOUBLE PRECISION[][] /*+ DEFAULT NULL */
) RETURNS DOUBLE PRECISION[][] AS $$
DECLARE
    theIteration INTEGER;
    theResult DOUBLE PRECISION[][];
    oldClientMinMessages VARCHAR;
    class_rel_source REGCLASS;
    proc_fn_dist REGPROCEDURE;
    num_points INTEGER;
    num_centroids INTEGER;
    rel_filtered VARCHAR;
BEGIN
    rel_filtered = MADLIB_SCHEMA.__filter_input_relation(rel_source, expr_point);
    class_rel_source := rel_filtered;

    IF (initial_centroids IS NOT NULL) THEN
    num_centroids := array_upper(initial_centroids,1);
    ELSE
    num_centroids := k;
    END IF;

    IF (k <= 0) THEN
        RAISE EXCEPTION 'Number of clusters k must be a positive integer.';
    END IF;
    IF (k > 32767) THEN
    RAISE EXCEPTION 'Number of clusters k must be <= 32767 (for results to be returned in a reasonable amount of time).';
    END IF;
    EXECUTE $sql$ SELECT count(*) FROM $sql$ || textin(regclassout(class_rel_source)) INTO num_points;
    IF (num_points < k  OR num_points < num_centroids) THEN
    RAISE EXCEPTION 'Number of centroids is greater than number of points.';
    END IF;
    IF (k < num_centroids) THEN
    RAISE WARNING 'Number of clusters k is less than number of supplied initial centroids. Number of final clusters will equal number of supplied initial centroids.';
    END IF;

    -- We first setup the argument table. Rationale: We want to avoid all data
    -- conversion between native types and Python code. Instead, we use Python
    -- as a pure driver layer.
    oldClientMinMessages :=
        (SELECT setting FROM pg_settings WHERE name = 'client_min_messages');
    EXECUTE 'SET client_min_messages TO warning';
    PERFORM MADLIB_SCHEMA.create_schema_pg_temp();
    -- Unfortunately, the EXECUTE USING syntax is only available starting
    -- PostgreSQL 8.4:
    -- http://www.postgresql.org/docs/8.4/static/plpgsql-statements.html#PLPGSQL-STATEMENTS-EXECUTING-DYN
    -- We therefore have to emulate.
    PERFORM MADLIB_SCHEMA.internal_execute_using_kmeans_random_seeding_args($sql$
        DROP TABLE IF EXISTS pg_temp._madlib_kmeans_random_args;
        CREATE TEMPORARY TABLE _madlib_kmeans_random_args AS
        SELECT $1 AS k, $2 AS initial_centroids;
        $sql$,
        k, initial_centroids);
    EXECUTE 'SET client_min_messages TO ' || oldClientMinMessages;

    -- Perform acutal computation.
    -- Unfortunately, Greenplum and PostgreSQL <= 8.2 do not have conversion
    -- operators from regclass to varchar/text.
    theIteration := (
        SELECT MADLIB_SCHEMA.internal_compute_kmeans_random_seeding(
            '_madlib_kmeans_random_args', '_madlib_kmeans_random_state',
            textin(regclassout(class_rel_source)), expr_point)
    );

    -- Retrieve result from state table and return it
    EXECUTE
        $sql$
        SELECT _state FROM _madlib_kmeans_random_state
        WHERE _iteration = $sql$ || theIteration || $sql$
        $sql$
        INTO theResult;
    RETURN theResult;
END;
$$ LANGUAGE plpgsql VOLATILE;

CREATE FUNCTION MADLIB_SCHEMA.kmeans_random_seeding(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER
) RETURNS DOUBLE PRECISION[][]
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans_random_seeding($1, $2, $3, NULL)
$$;

/**
 * @brief Run k-Means with random seeding.
 *
 * This is a shortcut for running k-means with random seeding. It is equivalent
 * to
 * <pre>SELECT \ref kmeans(
    rel_source,
    expr_point,
    \ref kmeans_random_seeding(
        rel_source,
        expr_point,
        k
    ),
    fn_dist,
    agg_centroid,
    max_num_iterations,
    min_frac_reassigned
)</pre>
 */
CREATE FUNCTION MADLIB_SCHEMA.kmeans_random(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR /*+ DEFAULT 'squared_dist_norm2' */,
    agg_centroid VARCHAR /*+ DEFAULT 'avg' */,
    max_num_iterations INTEGER /*+ DEFAULT 20 */,
    min_frac_reassigned DOUBLE PRECISION /*+ DEFAULT 0.001 */
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeans_random_seeding($1, $2, $3),
        $4, $5, $6, $7);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans_random(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR,
    agg_centroid VARCHAR,
    max_num_iterations INTEGER
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeans_random_seeding($1, $2, $3),
        $4, $5, $6, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans_random(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR,
    agg_centroid VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2, MADLIB_SCHEMA.kmeans_random_seeding($1, $2, $3),
        $4, $5, 20, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans_random(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER,
    fn_dist VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2,
        MADLIB_SCHEMA.kmeans_random_seeding($1, $2, $3),
        $4, 'MADLIB_SCHEMA.avg', 20, 0.001);
    RETURN ret;
END
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans_random(
    rel_source VARCHAR,
    expr_point VARCHAR,
    k INTEGER
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    ret MADLIB_SCHEMA.kmeans_result;
BEGIN
    ret = MADLIB_SCHEMA.kmeans(
        $1, $2,
        MADLIB_SCHEMA.kmeans_random_seeding($1, $2, $3),
        'MADLIB_SCHEMA.squared_dist_norm2', 'MADLIB_SCHEMA.avg', 20, 0.001);
    RETURN ret;
END
$$;

/**
 * @internal
 * @brief Execute a SQL command where $1, ..., $6 are substituted with the
 *     given arguments.
 */
CREATE FUNCTION MADLIB_SCHEMA.internal_execute_using_kmeans_args(
    sql VARCHAR, rel_source VARCHAR, expr_point VARCHAR,
    fn_dist VARCHAR, agg_centroid VARCHAR, max_num_iterations INTEGER,
    min_frac_reassigned DOUBLE PRECISION
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
CALLED ON NULL INPUT
LANGUAGE c
AS 'MODULE_PATHNAME', 'exec_sql_using';

/**
 * @internal
 * @brief Filter out the invalid data points in the original input relation
 */
CREATE FUNCTION MADLIB_SCHEMA.__filter_input_relation(
    rel_source VARCHAR, expr_point VARCHAR)
RETURNS VARCHAR
AS $$
DECLARE
    oldClientMinMessages VARCHAR;
    rel_source_filtered VARCHAR;
BEGIN
    IF (SELECT position('.' in rel_source)) > 0 THEN
        rel_source_filtered := '_madlib_' || split_part(rel_source, '.', 2) || '_filtered';
    ELSE
    rel_source_filtered := '_madlib_' || rel_source || '_filtered';
    END IF;

    oldClientMinMessages :=
        (SELECT setting FROM pg_settings WHERE name = 'client_min_messages');
    EXECUTE 'SET client_min_messages TO warning';
    EXECUTE 'DROP VIEW IF EXISTS _madlib_'||rel_source_filtered||'_filtered';
    EXECUTE 'DROP VIEW IF EXISTS '||rel_source_filtered;
    EXECUTE 'CREATE TEMP VIEW '||rel_source_filtered||'
             AS SELECT * FROM '||rel_source||'
                    WHERE abs(
                              coalesce(
                                 MADLIB_SCHEMA.svec_elsum('||expr_point||'),
                                 ''Infinity''::FLOAT8
                              )
                             ) < ''Infinity''::FLOAT8';
    EXECUTE 'SET client_min_messages TO ' || oldClientMinMessages;
    RETURN rel_source_filtered;
    EXCEPTION
        WHEN undefined_function THEN
        RAISE EXCEPTION 'Point coordinates (%) are not a valid type
                        (SVEC, FLOAT[], or INTEGER[]).', expr_point;
END
$$
VOLATILE STRICT
LANGUAGE PLPGSQL;

/**
 * @brief Perform Lloyd's k-means local-search heuristic, but with initial
 *     centroids stored in a table
 *
 * This is a shortcut for running k-means with initial centroids stored in a
 * table (as opposed to an array of centroids). It is equivalent
 * to
 * <pre>SELECT \ref kmeans(
    rel_source,
    expr_point,
    (SELECT \ref matrix_agg($expr_centroid) FROM $rel_initial_centroids),
    fn_dist,
    agg_centroid,
    max_num_iterations,
    min_frac_reassigned
)</pre>
 * where <tt>$expr_centroid</tt> and <tt>$rel_initial_centroids</tt> denote
 * textual substituions.
 */
CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    rel_initial_centroids VARCHAR,
    expr_centroid VARCHAR,
    fn_dist VARCHAR /*+ DEFAULT 'squared_dist_norm2' */,
    agg_centroid VARCHAR /*+ DEFAULT 'avg' */,
    max_num_iterations INTEGER /*+ DEFAULT 20 */,
    min_frac_reassigned DOUBLE PRECISION /*+ DEFAULT 0.001 */
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    class_rel_initial_centroids REGCLASS;
    theResult MADLIB_SCHEMA.kmeans_result;
BEGIN
    class_rel_initial_centroids := rel_initial_centroids;
    SELECT * FROM MADLIB_SCHEMA.internal_execute_using_kmeans_args($sql$
        SELECT MADLIB_SCHEMA.kmeans(
            $1, $2,
            (
                SELECT MADLIB_SCHEMA.matrix_agg(($sql$ || expr_centroid || $sql$)::FLOAT8[])
                FROM $sql$ || textin(regclassout(class_rel_initial_centroids))
                    || $sql$
            ),
            $3, $4, $5, $6)
            $sql$,
        rel_source, expr_point,
        fn_dist, agg_centroid, max_num_iterations, min_frac_reassigned)
        INTO theResult;
    RETURN theResult;
END;
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    rel_initial_centroids VARCHAR,
    expr_centroid VARCHAR,
    fn_dist VARCHAR,
    agg_centroid VARCHAR,
    max_num_iterations INTEGER
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans(
        $1, $2,
        $3, $4, $5, $6, $7, 0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    rel_initial_centroids VARCHAR,
    expr_centroid VARCHAR,
    fn_dist VARCHAR,
    agg_centroid VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans(
        $1, $2,
        $3, $4, $5, $6, 20, 0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    rel_initial_centroids VARCHAR,
    expr_centroid VARCHAR,
    fn_dist VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans(
        $1, $2,
        $3, $4, $5, 'MADLIB_SCHEMA.avg', 20, 0.001)
$$;

CREATE FUNCTION MADLIB_SCHEMA.kmeans(
    rel_source VARCHAR,
    expr_point VARCHAR,
    rel_initial_centroids VARCHAR,
    expr_centroid VARCHAR
) RETURNS MADLIB_SCHEMA.kmeans_result
VOLATILE
STRICT
LANGUAGE sql AS $$
    SELECT MADLIB_SCHEMA.kmeans(
        $1, $2,
        $3, $4,
        'MADLIB_SCHEMA.squared_dist_norm2', 'MADLIB_SCHEMA.avg', 20, 0.001)
$$;


/**
 * @internal
 * @brief Execute a SQL command where $1, ..., $3 are substituted with the
 *     given arguments.
 */
CREATE FUNCTION MADLIB_SCHEMA.internal_execute_using_silhouette_args(
    sql VARCHAR, centroids DOUBLE PRECISION[][], fn_dist REGPROC
) RETURNS DOUBLE PRECISION
STABLE
CALLED ON NULL INPUT
LANGUAGE c
AS 'MODULE_PATHNAME', 'exec_sql_using';

/**
 * @brief Compute a simplified version of the silhouette coefficient
 *
 * @param rel_source Name of the relation containing input points
 * @param expr_point Expression evaluating to point coordinates \f$ x_i \f$ for
 *     each tuple
 * @param centroids Matrix \f$ M = (\vec{m_0} \dots \vec{m_{k-1}})
 *     \in \mathbb{R}^{d \times k} \f$ with \f$ k \f$ columns, where column
 *     \f$ i \f$ contains the position of centroid \f$ i \f$.
 * @param fn_dist Name of a function with signature
 *     <tt>DOUBLE PRECISION[] x DOUBLE PRECISION[] -> DOUBLE PRECISION</tt> that
 *     returns the distance between two points
 * @return For each point \f$ x_i \f$, let
 *     \f$ d_1( x_i ) \f$ and \f$ d_2( x_i ) \f$ be the distance to the closest
 *     and 2nd-closest centroid, respectively. If there is more than one
 *     closest centroids then \f$ d_1( x_i ) = d_2( x_i )\f$.
 *     The return value is the average, over all points \f$ x_i \f$, of
 *     \f[
 *         \frac{d_2( x_i ) - d_1(x_i)}{d_2(x_i)},
 *     \f]
 *     where 0/0 is interpreted as 0.
 *     Clearly, the simplified silhouette coefficient assumes values in
 *     \f$ [0,1] \f$.
 */
CREATE FUNCTION MADLIB_SCHEMA.simple_silhouette(
    rel_source VARCHAR,
    expr_point VARCHAR,
    centroids DOUBLE PRECISION[][],
    fn_dist VARCHAR /*+ DEFAULT 'dist_norm2' */
) RETURNS DOUBLE PRECISION
STABLE
STRICT
LANGUAGE plpgsql
AS $$
DECLARE
    class_rel_source REGCLASS;
    proc_fn_dist REGPROCEDURE;
    rel_filtered VARCHAR;
BEGIN
    IF (array_upper(centroids,1) IS NULL) THEN
    RAISE EXCEPTION 'No valid centroids given.';
    END IF;

    rel_filtered = MADLIB_SCHEMA.__filter_input_relation(rel_source, expr_point);
    class_rel_source := rel_filtered;
    proc_fn_dist := fn_dist
        || '(DOUBLE PRECISION[], DOUBLE PRECISION[])';
    IF (SELECT prorettype != 'DOUBLE PRECISION'::regtype OR proisagg = TRUE
        FROM pg_proc WHERE oid = proc_fn_dist) THEN
        RAISE EXCEPTION 'Distance function has wrong signature or is not a simple function.';
    END IF;

    RETURN MADLIB_SCHEMA.internal_execute_using_silhouette_args($sql$
        SELECT
            avg(CASE
                    WHEN distances[2] = 0 THEN 0
                    ELSE (distances[2] - distances[1]) / distances[2]
                END)
        FROM (
            SELECT
                (MADLIB_SCHEMA.closest_columns(
                    $1,
                    ($sql$ || expr_point || $sql$)::FLOAT8[],
                    2::INTEGER,
                    $2
                )).distances
            FROM
                $sql$ || textin(regclassout(class_rel_source)) || $sql$
        ) AS two_shortest_distances
        $sql$,
        centroids, proc_fn_dist);
END;
$$;

CREATE FUNCTION MADLIB_SCHEMA.simple_silhouette(
    rel_source VARCHAR,
    expr_point VARCHAR,
    centroids DOUBLE PRECISION[][]
) RETURNS DOUBLE PRECISION
STABLE
STRICT
LANGUAGE sql
AS $$
    SELECT MADLIB_SCHEMA.simple_silhouette($1, $2, $3,
        'MADLIB_SCHEMA.dist_norm2')
$$;