# coding=utf-8
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import plpy

from utilities.control import MinWarning
from utilities.utilities import _assert
from utilities.utilities import unique_string
from utilities.utilities import add_postfix
from utilities.utilities import INTEGER, NUMERIC, ONLY_ARRAY
from utilities.utilities import is_valid_psql_type
from utilities.utilities import is_platform_pg
from utilities.utilities import num_features
from utilities.utilities import get_seg_number
from utilities.validate_args import input_tbl_valid, output_tbl_valid
from utilities.validate_args import is_var_valid
from utilities.validate_args import cols_in_tbl_valid
from utilities.validate_args import get_expr_type
from utilities.validate_args import get_algorithm_name
from graph.wcc import wcc

from math import log
from math import ceil
from math import sqrt
from time import time
import numpy as np
from collections import deque

import utilities.debug as DEBUG
DEBUG.plpy_info_enabled = False
DEBUG.plpy_execute_enabled = False
DEBUG.timings_enabled = False

try:
    from rtree import index
except ImportError:
    RTREE_ENABLED=0
else:
    RTREE_ENABLED=1

METHOD_BRUTE_FORCE = 'brute_force'
METHOD_OPTIMIZED = 'optimized'
DEFAULT_MIN_SAMPLES = 5
DEFAULT_METRIC = 'squared_dist_norm2'

def dbscan(schema_madlib, source_table, output_table, id_column, expr_point,
           eps, min_samples, metric, algorithm, max_segmentation_depth, **kwargs):

    with MinWarning("warning"):
        # algorithm=None is handled in get_algorithm_name()
        min_samples = DEFAULT_MIN_SAMPLES if not min_samples else min_samples
        metric = DEFAULT_METRIC if not metric else metric
        num_segs = get_seg_number()

        algorithm = get_algorithm_name(algorithm, METHOD_OPTIMIZED,
            [METHOD_BRUTE_FORCE, METHOD_OPTIMIZED], 'DBSCAN')

        if max_segmentation_depth is None:
            # Default to num_segs
            max_depth = num_segs
        else:
            max_depth = max_segmentation_depth
            if algorithm != METHOD_OPTIMIZED:
                plpy.warning("Ignoring max_segmentation_depth={} param, "
                          "N/A to algorithm={}".format(max_depth, algorithm))

        _validate_dbscan(schema_madlib, source_table, output_table, id_column,
                         expr_point, eps, min_samples, metric, algorithm, max_depth)

        dist_by_src = ''  if is_platform_pg() else 'DISTRIBUTED BY (__src__)'
        dist_by_id = ''  if is_platform_pg() else 'DISTRIBUTED BY (id)'
        dist_by_reach_id = ''  if is_platform_pg() else 'DISTRIBUTED BY (__reachable_id__)'
        dist_by_seg_id = '' if is_platform_pg() else 'DISTRIBUTED BY (seg_id)'
        distributed_by = '' if is_platform_pg() else 'DISTRIBUTED BY (__dist_id__)'

        core_points_table = unique_string(desp='core_points_table')
        core_edge_table = unique_string(desp='core_edge_table')
        distance_table = unique_string(desp='distance_table')

        source_view = unique_string(desp='source_view')
        plpy.execute("DROP VIEW IF EXISTS {0}".format(source_view))
        sql = """
            CREATE VIEW {source_view} AS
            SELECT ({id_column})::BIGINT AS id, ({expr_point})::DOUBLE PRECISION[] AS point
            FROM {source_table}
            """.format(**locals())
        plpy.execute(sql)

        if algorithm == METHOD_OPTIMIZED:
            dbscan_opt = DBSCAN_optimized(schema_madlib, source_view, output_table,
                                          eps, min_samples, metric, max_depth)
            dbscan_opt.run()
        else:
            # Calculate pairwise distances
            sql = """
                CREATE TEMP TABLE {distance_table} AS
                SELECT __src__, __dest__ FROM (
                    SELECT  __t1__.id AS __src__,
                            __t2__.id AS __dest__,
                            {schema_madlib}.{metric}(
                                __t1__.point, __t2__.point) AS __dist__
                    FROM {source_view} AS __t1__, {source_view} AS __t2__)q1
                WHERE __dist__ < {eps}
                {dist_by_src}
                """.format(**locals())
            plpy.execute(sql)

            # Find core points
            sql = """
                CREATE TEMP TABLE {core_points_table} AS
                SELECT * FROM (
                    SELECT __src__ AS id, count(*) AS __count__
                    FROM {distance_table} GROUP BY __src__) q1
                WHERE __count__ >= {min_samples}
                {dist_by_id}
                """.format(**locals())
            plpy.execute(sql)

            # Find the connections between core points to form the clusters
            sql = """
                CREATE TEMP TABLE {core_edge_table} AS
                SELECT __src__, __dest__
                FROM {distance_table} AS __t1__, (SELECT array_agg(id) AS arr
                                                  FROM {core_points_table}) __t2__
                WHERE __t1__.__src__ = ANY(arr) AND __t1__.__dest__ = ANY(arr)
                {dist_by_src}
            """.format(**locals())
            plpy.execute(sql)

            sql = """
                SELECT count(*) FROM {core_points_table}
                """.format(**locals())
            core_count = plpy.execute(sql)[0]['count']
            _assert(core_count != 0, "DBSCAN: Cannot find any core points/clusters.")

            # Start snowflake creation
            if is_platform_pg():
                sql = """
                    SELECT count(*) FROM {core_edge_table}
                    """.format(**locals())
                count = plpy.execute(sql)[0]['count']

                counts_list = [int(count)]
                seg_id = 0
                group_by_clause = ''

            else:
                sql = """
                    SELECT count(*), gp_segment_id FROM {core_edge_table} GROUP BY gp_segment_id
                    """.format(**locals())
                count_result = plpy.execute(sql)
                seg_num = get_seg_number()
                counts_list = [0]*seg_num
                for i in count_result:
                    counts_list[int(i['gp_segment_id'])] = int(i['count'])
                seg_id = 'gp_segment_id'
                group_by_clause = 'GROUP BY gp_segment_id'

            snowflake_table = unique_string(desp='snowflake_table')
            sf_edge_table = unique_string(desp='sf_edge_table')

            plpy.execute("DROP TABLE IF EXISTS {0}, {1}".format(
                snowflake_table, sf_edge_table))
            sql = """
                CREATE TEMP TABLE {snowflake_table} AS
                SELECT {seg_id}::INTEGER AS seg_id,
                       {schema_madlib}.build_snowflake_table( __src__,
                                                __dest__,
                                                ARRAY{counts_list},
                                                {seg_id}
                                               ) AS __sf__
                FROM {core_edge_table} {group_by_clause}
                {dist_by_seg_id}
            """.format(**locals())
            plpy.execute(sql)

            sql = """
                CREATE TEMP TABLE {sf_edge_table} AS
                SELECT seg_id, (edge).src AS __src__, (edge).dest AS __dest__
                FROM (
                    SELECT seg_id,
                           {schema_madlib}.__dbscan_unpack_edge(__sf__) AS edge
                    FROM {snowflake_table}
                    ) q1
                {dist_by_seg_id}
                """.format(**locals())
            plpy.execute(sql)

            # Run wcc to get the min id for each cluster
            wcc(schema_madlib, core_points_table, 'id', sf_edge_table,
                'src=__src__, dest=__dest__', output_table, None)

            plpy.execute("""
                ALTER TABLE {0}
                ADD COLUMN is_core_point BOOLEAN,
                ADD COLUMN point DOUBLE PRECISION[];
                """.format(output_table)
            )
            plpy.execute("""
                ALTER TABLE {0}
                RENAME COLUMN component_id TO cluster_id;
                """.format(output_table)
            )
            plpy.execute("""
                UPDATE {0}
                SET is_core_point = TRUE
                """.format(output_table)
            )

            # Find reachable points
            reachable_points_table = unique_string(desp='reachable_points_table')
            plpy.execute("DROP TABLE IF EXISTS {0}".format(reachable_points_table))
            sql = """
                CREATE TEMP TABLE {reachable_points_table} AS
                    SELECT array_agg(__src__) AS __src_list__,
                           __dest__ AS __reachable_id__
                    FROM {distance_table} AS __t1__,
                         (SELECT array_agg(id) AS __arr__
                          FROM {core_points_table}) __t2__
                    WHERE __src__ = ANY(__arr__) AND __dest__ != ALL(__arr__)
                    GROUP BY __dest__
                    {dist_by_reach_id}
                """.format(**locals())
            plpy.execute(sql)

            sql = """
                INSERT INTO {output_table}
                SELECT  __reachable_id__,
                        cluster_id,
                        FALSE AS is_core_point,
                        NULL AS point
                FROM {reachable_points_table} AS __t1__ JOIN
                     {output_table} AS __t2__
                     ON (__src_list__[1] = id)
                """.format(**locals())
            plpy.execute(sql)

            # Add features of points to the output table to use them for prediction
            sql = """
                UPDATE {output_table} AS __t1__
                SET point = __t2__.point
                FROM {source_view} AS __t2__
                WHERE __t1__.id = __t2__.id
            """.format(**locals())
            plpy.execute(sql)

            plpy.execute("DROP TABLE IF EXISTS {0}, {1}, {2}, {3}, {4}, {5}".format(
                             distance_table, core_points_table, core_edge_table,
                             reachable_points_table, snowflake_table, sf_edge_table))

        plpy.execute("DROP VIEW IF EXISTS {0}".format(source_view))

        # -- Update the cluster ids to be consecutive --
        sql = """
            UPDATE {output_table} AS __t1__
            SET cluster_id = new_id-1
            FROM (
                SELECT cluster_id, row_number() OVER(ORDER BY cluster_id) AS new_id
                FROM {output_table}
                GROUP BY cluster_id) __t2__
            WHERE __t1__.cluster_id = __t2__.cluster_id
        """.format(**locals())
        plpy.execute(sql)

        # -- Generate output summary table ---

        output_summary_table = add_postfix(output_table, '_summary')

        plpy.execute("DROP TABLE IF EXISTS {0}".format(output_summary_table))
        sql = """
            CREATE TABLE {output_summary_table} AS
            SELECT  $${id_column}$$::VARCHAR AS id_column,
                    $${expr_point}$$::VARCHAR AS expr_point,
                    {eps}::DOUBLE PRECISION AS eps,
                    {min_samples}::BIGINT AS min_points,
                    '{metric}'::VARCHAR AS metric
            """.format(**locals())
        plpy.execute(sql)

class DBSCAN(object):
    def __init__(self, schema_madlib, source_table, output_table,
                 eps, min_samples, metric):
        """
            Args:
                @param schema_madlib        Name of the Madlib Schema
                @param source_table         Training data table
                @param output_table         Output table of points labelled by cluster_id
                @param eps                  The eps value defined by the user
                @param min_samples          min_samples defined by user
                @param metric               metric selected by user

        """
        self.schema_madlib = schema_madlib
        self.source_table = source_table
        self.output_table = output_table
        self.n_dims = num_features(source_table, 'point')
        self.distributed_by = '' if is_platform_pg() else 'DISTRIBUTED BY (__dist_id__)'
        self.min_samples = min_samples
        self.metric = metric
        # If squared_dist_norm2 is used, we assume eps is set for the squared
        # distance.  That means the border width must be sqrt(eps)
        self.eps = sqrt(eps) if metric == DEFAULT_METRIC else eps

class DBSCAN_optimized(DBSCAN):
    def __init__(self, schema_madlib, source_table, output_table,
                 eps, min_samples, metric, max_depth=None):
        """
            Args:
                @param schema_madlib        Name of the Madlib Schema
                @param source_table         Training data table
                @param output_table         Output table of points labelled by cluster_id
                @param eps                  The eps value defined by the user
                @param min_samples          min_samples defined by user
                @param metric               metric selected by user
        """

        DBSCAN.__init__(self, schema_madlib, source_table, output_table,
                        eps, min_samples, metric)
        self.max_depth = max_depth
        self.unique_str = unique_string('DBSCAN_opt')
        self.num_segs = get_seg_number()
        res = plpy.execute("SELECT COUNT(*) FROM {source_table}".format(**locals()))
        self.num_points = res[0]['count']

    def run(self):
        unique_str = self.unique_str
        distributed_by = self.distributed_by
        schema_madlib = self.schema_madlib

        self.dist_map = 'dist_map' + unique_str

        segmented_source = 'segmented_source' + unique_str
        self.segmented_source = segmented_source
        self.build_optimized_tree()

        rowcount_table = 'rowcounts' + unique_str
        leaf_clusters_table = 'leaf_clusters' + unique_str

        count_rows_sql = """
            DROP TABLE IF EXISTS {rowcount_table};
            CREATE TEMP TABLE {rowcount_table} AS
            WITH i AS (
                SELECT
                    count(*) AS num_rows,
                    dist_id,
                    __dist_id__
                FROM internal_{segmented_source}
                GROUP BY dist_id, __dist_id__
            ), e AS (
                SELECT
                    count(*) AS num_rows,
                    dist_id,
                    __dist_id__
                FROM external_{segmented_source}
                GROUP BY dist_id, __dist_id__
            )
            SELECT
                i.num_rows AS num_internal_rows,
                COALESCE(e.num_rows,0) AS num_external_rows,
                i.dist_id,
                __dist_id__
            FROM i LEFT JOIN e
                USING (dist_id, __dist_id__)
            {distributed_by}
            """.format(**locals())
        plpy.execute(count_rows_sql)

        # Run dbscan in parallel on each leaf
        dbscan_leaves_sql = """
            CREATE TEMP TABLE {leaf_clusters_table} AS
            WITH input AS (
                    SELECT *
                    FROM internal_{segmented_source}
                UNION ALL
                    SELECT *
                    FROM external_{segmented_source}
            ), segmented_source AS (
                SELECT
                    ROW(
                        id,
                        point,
                        FALSE,
                        0,
                        leaf_id,
                        dist_id
                    )::{schema_madlib}.__dbscan_record AS rec,
                    __dist_id__
                FROM input
            )
            SELECT (output).*, __dist_id__ FROM (
                SELECT
                    {schema_madlib}.__dbscan_leaf(
                                                  rec,
                                                  {self.eps},
                                                  {self.min_samples},
                                                  '{self.metric}',
                                                  num_internal_rows,
                                                  num_external_rows
                                                 ) AS output,
                    __dist_id__
                FROM segmented_source JOIN {rowcount_table}
                    USING (__dist_id__)
            ) a {distributed_by}
        """.format(**locals())
        DEBUG.plpy_execute(dbscan_leaves_sql, report_segment_tracebacks=True)

        plpy.execute("""
            DROP TABLE IF EXISTS
            internal_{segmented_source},
            external_{segmented_source}
        """.format(**locals()))

        cluster_map = 'cluster_map' + unique_str
        noise_point = label.NOISE_POINT

        # Replace local cluster_id's in each leaf with
        #   global cluster_id's which are unique across all leaves
        gen_cluster_map_sql = """
            CREATE TEMP TABLE {cluster_map} AS
                SELECT
                    dist_id,
                    __dist_id__,
                    cluster_id AS local_id,
                    ROW_NUMBER() OVER(
                        ORDER BY dist_id, cluster_id
                    ) AS global_id
                FROM {leaf_clusters_table}
                    WHERE dist_id = leaf_id AND cluster_id != {noise_point}
                GROUP BY __dist_id__, dist_id, cluster_id
                {distributed_by}
        """.format(**locals())
        plpy.execute(gen_cluster_map_sql)

        globalize_cluster_ids_sql = """
            UPDATE {leaf_clusters_table} lc
                SET cluster_id = global_id
            FROM {cluster_map} cm WHERE
                cm.dist_id = lc.dist_id AND
                cluster_id = local_id;
        """.format(**locals())
        DEBUG.plpy.execute(globalize_cluster_ids_sql)

        intercluster_edges = 'intercluster_edges' + unique_str
        dist_by_src = '' if is_platform_pg() else 'DISTRIBUTED BY (src)'

        # Build intercluster table of edges connecting clusters
        #  on different leaves
        intercluster_edge_sql = """
            CREATE TEMP TABLE {intercluster_edges} AS
                SELECT
                    internal.cluster_id AS src,
                    external.cluster_id AS dest
                FROM {leaf_clusters_table} internal
                JOIN {leaf_clusters_table} external
                    USING (id)
                WHERE internal.is_core_point
                    AND internal.leaf_id  = internal.dist_id
                    AND external.leaf_id != external.dist_id
                GROUP BY src,dest
                {dist_by_src}
        """.format(**locals())
        DEBUG.plpy.execute(intercluster_edge_sql)

        res = plpy.execute("SELECT count(*) FROM {intercluster_edges}".format(**locals()))

        if int(res[0]['count']) > 0:
            wcc_table = 'wcc' + unique_str
            plpy.execute("""
            DROP TABLE IF EXISTS
                {wcc_table}, {wcc_table}_summary,
                {self.output_table}
            """.format(**locals()))

            # Find connected components of intercluster_edge_table
            wcc(schema_madlib, cluster_map, 'global_id', intercluster_edges,
            'src=src,dest=dest', wcc_table, None)

           # Rename each cluster_id to its component id in the intercluster graph
            merge_clusters_sql = """
            UPDATE {leaf_clusters_table} lc
                SET cluster_id = wcc.component_id
            FROM {wcc_table} wcc
                WHERE lc.cluster_id = wcc.global_id
            """.format(**locals())
            plpy.execute(merge_clusters_sql)

            plpy.execute("""
            DROP TABLE IF EXISTS
                {wcc_table}, {wcc_table}_summary
            """.format(**locals()))

        add_crossborder_points_sql = """
            CREATE TABLE {self.output_table} AS
            SELECT
                id, point, cluster_id, is_core_point
            FROM (
                SELECT
                    id,
                    point,
                    is_core_point,
                    leaf_id = dist_id AS is_internal,
                    dist_id,
                    CASE WHEN cluster_id = {noise_point}
                        THEN MAX(cluster_id) OVER(PARTITION BY id)
                        ELSE cluster_id
                    END
                FROM {leaf_clusters_table}
            ) x WHERE is_internal AND cluster_id != {noise_point}
        """.format(**locals())
        DEBUG.plpy.execute(add_crossborder_points_sql)

        plpy.execute("""
            DROP TABLE IF EXISTS {rowcount_table},
            {leaf_clusters_table}, {cluster_map},
            {intercluster_edges}
        """.format(**locals()))

    def build_optimized_tree(self):
        eps = self.eps
        num_segs = self.num_segs
        num_points = self.num_points
        target = num_points / num_segs
        source_table = self.source_table
        leaf_bounds_table = 'leaf_bounds' + self.unique_str
        max_depth = self.max_depth
        n_dims = self.n_dims
        dist_map = self.dist_map
        distributed_by = self.distributed_by
        internal_points = 'internal_' + self.segmented_source
        external_points = 'external_' + self.segmented_source

        next_cut = 'next_cut' + self.unique_str
        prev_node = 'prev_node' + self.unique_str

        first_cut_sql = """
            DROP TABLE IF EXISTS {prev_node};
            CREATE TEMP TABLE {prev_node} AS
            WITH world_bounds AS (
                SELECT
                    i,
                    min(point[i]) AS lower,
                    max(point[i])+{eps} AS upper,
                    max(point[i])+{eps} - min(point[i]) AS size,
                    floor((max(point[i])+{eps} - min(point[i])) / {eps})::INT AS eps_bins
                FROM generate_series(1,{n_dims}) AS i, {source_table}
                GROUP BY i
            ),
            density_map AS (
                SELECT i, eps_bin,
                    eps_bins,
                    COALESCE(density, 0) AS density
                FROM (
                    SELECT i, eps_bins,
                        generate_series(0, eps_bins - 1) AS eps_bin
                    FROM world_bounds
                ) g LEFT JOIN (
                    SELECT i,
                        floor((point[i] - lower)/{eps})::INT AS eps_bin,
                        eps_bins,
                        COUNT(*) AS density
                    FROM
                        world_bounds,
                        {source_table}
                    GROUP BY i, eps_bin, eps_bins
                ) a
                USING (i, eps_bin, eps_bins)
            ),
            loss_table AS (
                SELECT i, eps_bin,
                    left_internal,
                    {num_points} as points_in_node,
                    left_external,
                    right_external,
                    {self.schema_madlib}.__dbscan_segmentation_loss(
                        left_internal,
                        {num_points} - left_internal,
                        left_external,
                        right_external,
                        {num_segs}::BIGINT,
                        V,
                        {eps}::DOUBLE PRECISION,
                        {n_dims}::BIGINT,
                        eps_bin::DOUBLE PRECISION,
                        eps_bins::DOUBLE PRECISION
                    ) AS losses
                FROM (
                    SELECT i, eps_bin, eps_bins,
                        {num_segs}::BIGINT AS num_segs,
                        COALESCE(density, 0) AS right_external,
                        COALESCE(LEAD(density) OVER(PARTITION BY i ORDER BY eps_bin), 0) AS left_external,
                        SUM(density) OVER(PARTITION BY i ORDER BY eps_bin)::BIGINT AS left_internal
                    FROM (
                        SELECT i, generate_series(0, eps_bins - 1) AS eps_bin FROM world_bounds
                    ) g LEFT JOIN density_map USING(i, eps_bin)
                ) params,
                (
                    SELECT {self.schema_madlib}.prod(size) AS V
                    FROM world_bounds
                ) AS volume
            ),
            optimize AS (
                SELECT i, lower, cutoff, upper,
                       left_avail, right_avail,
                       left_internal, right_internal
                FROM (
                    SELECT
                        i,
			CASE WHEN (losses).right_avail = 0
			   THEN
                               points_in_node
                           ELSE
                               left_internal
                        END AS left_internal,
			CASE WHEN (losses).right_avail = 0
			   THEN
                               0::BIGINT
                           ELSE
                               points_in_node - left_internal
                        END AS right_internal,
                        (losses).left_avail AS left_avail,
                        (losses).right_avail AS right_avail,
                        GREATEST((losses).left_loss,(losses).right_loss) AS total_loss,
                        min(
                            GREATEST((losses).left_loss,(losses).right_loss)
                        ) OVER() AS total_min,
                        wb.lower,
                        CASE WHEN (losses).right_avail = 0
                        THEN
                            wb.upper
                        ELSE
                            wb.lower + {eps}*(eps_bin+1)
                        END AS cutoff,
                        wb.upper
                    FROM loss_table JOIN world_bounds wb USING(i)
                ) a WHERE total_loss = total_min LIMIT 1
            )
            SELECT s.id,
                ARRAY[i] AS coords,
                CASE WHEN s.point[i] < cutoff
                    THEN left_avail
                    ELSE right_avail
                END AS avail_segs,
                CASE WHEN point[i] < cutoff
                    THEN left_internal
                    ELSE right_internal
                END AS points_in_node,
                CASE WHEN s.point[i] < cutoff
                    THEN ARRAY[lower]
                    ELSE ARRAY[cutoff]
                END AS lower_bounds,
                CASE WHEN s.point[i] < cutoff
                    THEN ARRAY[cutoff]
                    ELSE ARRAY[upper]
                END AS upper_bounds,
                (s.point[i] >= cutoff)::INT AS node_id
            FROM optimize, {source_table} s
        """.format(**locals())
        plpy.execute(first_cut_sql)

        segment_source_sql = """
            CREATE TEMP TABLE {next_cut} AS
            WITH node_bounds AS (
                SELECT
                    i,
                    node_id,
                    avail_segs,
                    points_in_node,
                    min(point[i]) AS lower,
                    max(point[i])+{eps} AS upper,
                    max(point[i])+{eps} - min(point[i]) AS size,
                    floor((max(point[i])+{eps} - min(point[i])) / {eps})::INT AS eps_bins
                FROM generate_series(1, {n_dims}) AS i,
                    {prev_node} JOIN {source_table} USING (id)
                GROUP BY 1,2,3,4
            ), density_map AS (
                SELECT i, node_id, eps_bin,
                    COALESCE( density, 0) AS density
                FROM (
                    SELECT node_id, i,
                        generate_series(0, eps_bins - 1) AS eps_bin
                    FROM node_bounds
                ) g LEFT JOIN (
                    SELECT node_id, i,
                        floor((point[i] - lower)/{eps})::INT AS eps_bin,
                        COUNT(*) AS density
                    FROM
                        node_bounds JOIN {prev_node} USING (node_id)
                        JOIN {source_table} USING (id)
                    GROUP BY node_id, i, eps_bin
                ) a
                USING(i, node_id, eps_bin)
            ), params AS (
                SELECT i, node_id, eps_bin, eps_bins,
                    COALESCE(density, 0) AS right_external,
                    COALESCE(LEAD(density) OVER(PARTITION BY node_id,i ORDER BY eps_bin), 0) AS left_external,
                    SUM(density) OVER(PARTITION BY node_id,i ORDER BY eps_bin)::BIGINT AS left_internal
                FROM (
                    SELECT i, node_id, eps_bins, generate_series(0, eps_bins - 1) AS eps_bin FROM node_bounds
                ) g LEFT JOIN density_map USING(i, node_id, eps_bin)
            ), volume AS (
                SELECT
                    node_id,
                    avail_segs,
                    points_in_node,
                    {self.schema_madlib}.prod(size) AS V
                FROM node_bounds
                WHERE i=1
                GROUP BY 1,2,3
            ), loss_table AS (
                SELECT i, node_id, eps_bin,
                    left_internal,
                    points_in_node - left_internal AS right_internal,
                    left_external,
                    right_external,
                    {self.schema_madlib}.__dbscan_segmentation_loss(
                        left_internal,
                        points_in_node,
                        left_external,
                        right_external,
                        avail_segs,
                        V,
                        {eps}::DOUBLE PRECISION,
                        {n_dims}::BIGINT,
                        eps_bin::DOUBLE PRECISION,
                        eps_bins::DOUBLE PRECISION
                    ) AS losses
                FROM params JOIN volume USING (node_id)
            ), optimize AS (
                SELECT * FROM (
                    SELECT
                        node_id,
                        i,
                        (losses).left_avail,
                        (losses).right_avail,
			CASE WHEN (losses).right_avail = 0
			   THEN
                               points_in_node
                           ELSE
                               left_internal
                        END AS left_internal,
			CASE WHEN (losses).right_avail = 0
			   THEN
                               0::BIGINT
                           ELSE
                               points_in_node - left_internal
                        END AS right_internal,
                        node_bounds.lower,
                        CASE WHEN (losses).right_avail = 0
                            THEN
                                node_bounds.upper
                            ELSE
                                node_bounds.lower + {eps}*(eps_bin+1)
                            END AS cutoff,
                        node_bounds.upper,
                        ROW_NUMBER() OVER(
                            PARTITION BY node_id
                            ORDER BY GREATEST((losses).left_loss,(losses).right_loss)
                        ) AS loss_rank
                    FROM loss_table JOIN node_bounds USING(i, node_id)
                ) a WHERE loss_rank = 1
            )
            SELECT id, coords || i AS coords,
                CASE WHEN point[i] < opt.cutoff
                    THEN left_avail
                    ELSE right_avail
                END AS avail_segs,
                CASE WHEN point[i] < opt.cutoff
                    THEN left_internal
                    ELSE right_internal
                END AS points_in_node,
                lower_bounds ||
                    CASE WHEN point[i] < opt.cutoff
                    THEN
                        lower
                    ELSE
                        opt.cutoff
                    END AS lower_bounds,
                    opt.cutoff,
                upper_bounds ||
                    CASE WHEN point[i] < opt.cutoff
                    THEN
                        opt.cutoff
                    ELSE
                        upper
                    END AS upper_bounds,
                (node_id << 1) | (point[i] >= opt.cutoff)::INT AS node_id
            FROM {source_table}
            JOIN {prev_node}
                USING (id)
            JOIN optimize opt
                USING (node_id)
        """.format(**locals())

        for i in range(max_depth):
            DEBUG.plpy.info(segment_source_sql)
            plpy.execute(segment_source_sql)
            plpy.execute("""
                DROP TABLE IF EXISTS {prev_node};
                ALTER TABLE {next_cut} RENAME TO {prev_node}
            """.format(**locals()))

        if is_platform_pg():
            # For postgres, there should only be 1 node,
            #  so we assign it dist_id = __dist_id__ = 0
            plpy.execute("""
                CREATE TABLE {dist_map} AS
                SELECT 0::BIGINT AS dist_id, 0::BIGINT as __dist_id__
            """.format(**locals()))
        else:
            # Create implicit many-to-1 map from __dist_id__ to seg_id
            temp_dist_table = 'dist_table' + self.unique_str
            create_dist_table_sql = """
                CREATE TEMP TABLE {temp_dist_table} AS
                    SELECT generate_series(1,10000) AS __dist_id__
                {distributed_by}
            """.format(**locals())
            plpy.execute(create_dist_table_sql)

            # Create a 1-1 map from dist_id to __dist_id__,
            # where dist_id's are re-ordered so they are evenly
            # distributed among segments.  This should work even
            # if there are gaps in the dist_id assignments (some
            # dist_id's missing)
            leaves_sql = """
                     SELECT
                         node_id AS dist_id,
                         ROW_NUMBER() OVER() AS leaf_num
                     FROM {prev_node}
                     GROUP BY node_id
            """.format(**locals())
            leaves_cnt = plpy.execute("""
                WITH leaves as ({leaves_sql})
                SELECT COUNT(*) cnt from leaves
                """.format(**locals()))[0]['cnt']
            create_dist_map_sql = """
                CREATE TEMP TABLE {dist_map} AS
                WITH leaves as (
                    {leaves_sql}
                ), segmap AS (
                     SELECT
                         gp_segment_id AS seg_id,
                         ROW_NUMBER() OVER(
                             ORDER BY ctid, gp_segment_id
                         ) AS leaf_num,
                         __dist_id__
                     FROM {temp_dist_table}
                     LIMIT {leaves_cnt}
                )
                SELECT dist_id, __dist_id__
                FROM leaves JOIN segmap
                     USING (leaf_num)
                {distributed_by}
            """.format(**locals())
            plpy.execute(create_dist_map_sql)

            plpy.execute("""
                DROP TABLE IF EXISTS {temp_dist_table}
            """.format(**locals()))

        gen_internal_points_sql = """
            CREATE TEMP TABLE {internal_points} AS
            SELECT id, point, node_id AS leaf_id, node_id AS dist_id, __dist_id__
            FROM {source_table}
            JOIN {prev_node}
                USING (id)
            JOIN {dist_map}
                ON node_id = dist_id
            {distributed_by}
        """.format(**locals())
        DEBUG.plpy.execute(gen_internal_points_sql)

        plpy.execute("""
            DROP TABLE IF EXISTS {prev_node}
        """.format(**locals()))

        generate_leaf_bounds_sql = """
            CREATE TEMP TABLE {leaf_bounds_table} AS
            SELECT
               array_agg(min ORDER BY i) AS leaf_lower,
               array_agg(max ORDER BY i) AS leaf_upper,
               array_agg(min - {eps} ORDER BY i) AS dist_lower,
               array_agg(max + {eps} ORDER BY i) AS dist_upper,
               leaf_id,
               __dist_id__
            FROM (
                SELECT
                    leaf_id,
                    i,
                    min(point[i]),
                    max(point[i]),
                count(*)
                FROM {internal_points}, generate_series(1, {n_dims}) AS i
                GROUP BY leaf_id, i
            ) x
            JOIN {dist_map}
                ON leaf_id = dist_id
            GROUP BY leaf_id, __dist_id__
            {distributed_by}
        """.format(**locals())
        DEBUG.plpy.execute(generate_leaf_bounds_sql)

        gen_external_points_sql = """
            CREATE TEMP TABLE {external_points} AS
            SELECT i.id, i.point, i.leaf_id, b.leaf_id AS dist_id, b.__dist_id__
            FROM {leaf_bounds_table} b JOIN {internal_points} i
                ON b.leaf_id != i.leaf_id
            WHERE 0 <= all({self.schema_madlib}.array_sub(point, dist_lower)) AND
                  0 < all({self.schema_madlib}.array_sub(dist_upper, point))
            {distributed_by}
        """.format(**locals())
        plpy.execute(gen_external_points_sql)

        plpy.execute("""
            DROP TABLE IF EXISTS
                {leaf_bounds_table}, {dist_map}
        """.format(**locals()))

def dbscan_predict(schema_madlib, dbscan_table, source_table, id_column,
    expr_point, output_table, **kwargs):

    with MinWarning("warning"):

        _validate_dbscan_predict(schema_madlib, dbscan_table, source_table, id_column,
                                 expr_point, output_table)

        dbscan_summary_table = add_postfix(dbscan_table, '_summary')
        summary = plpy.execute("SELECT * FROM {0}".format(dbscan_summary_table))[0]

        eps = summary['eps']
        metric = summary['metric']
        sql = """
            CREATE TABLE {output_table} AS
            WITH src AS (
                SELECT ({id_column}) AS id, ({expr_point}) AS point FROM {source_table}
            )
            SELECT __q1__.id, cluster_id, distance
            FROM (
                SELECT __t2__.id, cluster_id,
                       min({schema_madlib}.{metric}(__t1__.point,
                                                __t2__.point)) as distance
                FROM {dbscan_table} AS __t1__, src AS __t2__
                WHERE is_core_point = TRUE
                GROUP BY __t2__.id, cluster_id
            ) __q1__ WHERE distance <= {eps}
            """.format(**locals())
        result = plpy.execute(sql)

class label:
    UNLABELLED = 0
    NOISE_POINT = -1

pstats = None
class dbscan_perfstats:
    def __init__(self):
        self.range_query_calls = 0
        self.range_query_time = 0.0
        self.range_query_candidates = 0
        self.range_query_neighbors = 0

    def show(self, msg):
        DEBUG.plpy.info("{}_stats: range_query total calls = {}"
             .format(msg, self.range_query_calls))

        range_query_calls = self.range_query_calls if self.range_query_calls else 0.00001

        DEBUG.plpy.info("{}_stats: range_query total time = {}s"
             .format(msg, self.range_query_time))
        DEBUG.plpy.info("{}_stats: range_query total candidates = {}"
            .format(msg, self.range_query_candidates))
        DEBUG.plpy.info("{}_stats: range_query total neighbors returned = {}"
            .format(msg, self.range_query_neighbors))
        DEBUG.plpy.info("{}_stats: range_query avg time = {}s"
            .format(msg, self.range_query_time / range_query_calls))
        DEBUG.plpy.info("{}_stats: range_query avg candidates returned = {}"
            .format(msg, self.range_query_candidates * 1.0 / range_query_calls))
        DEBUG.plpy.info("{}_stats: range_query avg neighbors returned = {}"
            .format(msg, self.range_query_neighbors * 1.0 / range_query_calls))
    	self.__init__()  # Clear counters and start over

def within_range_tf(point, candidate_ids, candidate_points,
                    metric_cutoff, norm=2):
    sd = tf.squared_difference(tfcandidates, tfpoint)
    distances = tf.math.reduce_sum(sd, 1)
    output = tf.where(distances <= 2)
    res = output.eval(session=sess)
    res.shape = res.shape[0]
    return res

def within_range_np_linalg(point, candidate_ids, candidate_points,
                           eps, norm=2):
    distances = np.linalg.norm(candidate_points - point, ord=norm, axis=1)
    return candidate_ids[distances <= eps]

def within_range_np(point, candidate_ids, candidate_points,
                    metric_cutoff, norm=2):
    if norm==1:
        distances = np.sum(candidate_points, axis=1)
    else:
        distances = np.sum(np.square(candidate_points - point), axis=1)
    return candidate_ids[distances <= metric_cutoff]

class dbscan_record(object):
    def __init__(self, obj, eps=None):
            if eps is not None:
                self._init_from_dict(obj, eps)
            else:
                self._clone(obj)

    @staticmethod
    def from_dict(d, eps):
        if d['leaf_id'] == d['dist_id']:
            return dbscan_internal_record(d, eps)
        else:
            return dbscan_external_record(d, eps)

    def _init_from_dict(self, d, eps):
        """
            dbscan_record dict constructor
            The only reason for passing eps here is so we
            can pre-compute the point's neighborhood for later
        """
        self.id = d['id']
        self.point = d['point']
        self.is_core_point = False
        self.cluster_id = label.UNLABELLED
        self.leaf_id = d['leaf_id']
        self.dist_id = d['dist_id']

        # private members (for convenience, not for output)
        p = np.array(self.point)
        self._np_point = p
        self._bounds = np.concatenate([p, p])
        self._neighborhood = np.concatenate([p - eps, p + eps])

    def _clone(self, rec):
        """
            dbscan_record copy constructor - we call this
            just before outputing any row.  For external
            points, we have to clone it each time so that
            we can output multiple rows.  For internal points,
            it's just a quick way of cleaning it up (removing
            the extra private vars).
        """
        # Copy only public members
        self.id = rec.id
        self.point = rec.point
        self.is_core_point = rec.is_core_point
        self.cluster_id = rec.cluster_id
        self.leaf_id = rec.leaf_id
        self.dist_id = rec.dist_id

    def __repr__(self):
        fmt = 'id={},point[{}],cluster_id={},core={},leaf_id={},dist_id={}'
        rep = fmt.format(self.id, len(self.point), self.cluster_id,
                         self.is_core_point, self.leaf_id, self.dist_id)
        return 'dbscan_record({})'.format(rep)

    def __eq__(self, other):
        if not isinstance(other, dbscan_record):
            return False
        return self.__dict__ == other.__dict__

class dbscan_internal_record(dbscan_record):
    def __init__(self, obj, eps=None):
        dbscan_record.__init__(self, obj, eps)

    def _init_from_dict(self, d, eps):
        """
            private members specific to internal records
        """
        dbscan_record._init_from_dict(self, d, eps)
        self._external_neighbors = []
        self._internal_neighbor_count = 0
        self._possible_border_points = []

    def not_core_point(self, min_samples):
        """
            Computes self.is_core_point from internal and
            external neighbor counts, and returns its negation
        """
        ext_neighbors = len(self._external_neighbors)
        int_neighbors = self._internal_neighbor_count
        num_neighbors = ext_neighbors + int_neighbors
        self.is_core_point = (num_neighbors + 1 >= min_samples)
        return not self.is_core_point

class dbscan_external_record(dbscan_record):
    def __init__(self, obj, eps=None):
        dbscan_record.__init__(self, obj, eps)

    def _init_from_dict(self, d, eps):
        """
            private members specific to external records
        """
        dbscan_record._init_from_dict(self, d, eps)
        self._nearby_clusters = set()

class LeafStorage:
    """
    LeafStorage

    This is a per-leaf storage class for dbscan_leaf().
    There are 3 containers within LeafStorage:

    rtree - a spatial index of the id's of the leaf's internal points
    records - a hash table for looking up dbscan_internal_records by id
    external_records - an ordered list of dbscan_external_records

    Methods:

    add_point() - adds a point to LeafStorage
    range_query() - returns a list of internal point id's within eps of a point
    update_neighbors() - updates internal neighbor counts based on a point and
                         list of neighboring internal points
    label_point() - assigns a point to a specific cluster, or marks it as noise

    yield_external_records() - clones and yields a copy of each dbscan_external_record
                               based on its list of nearby_clusters
    """
    def __init__(self, dimension, dtype, eps, metric, min_samples):
        self.n_dims = dimension
        self.dtype = dtype
        self.eps = eps
        self.metric = metric
        self.min_samples = min_samples
        props = index.Property(dimension=dimension)
        self.props = props
        self.rtree = index.Index(properties=props)
        self.records = {}
        self.external_records = []

    def add_point(self, rec):
        if rec.leaf_id == rec.dist_id:
            self.rtree.add(rec.id, rec._bounds)
            self.records[rec.id] = rec
        else:
            self.external_records.append(rec)

    def range_query(self, db_rec):
        """
        Returns a list of records of all neighboring
        points within eps distance from point.  Returns None if the
        number of results will be < min_results.  Also returns None
        if number of results > max_results > 0  The special value
        max_results = 0 means unlimited results (return all).
        """
        global pstats
        pstats.range_query_calls += 1
        start_time = time()

        eps = self.eps
        metric = self.metric
        records = self.records
        rtree = self.rtree

        if not metric in ('squared_dist_norm2', 'dist_norm2', 'dist_norm1'):
            plpy.error("Sorry, optimized method does not support {} metric yet".format(metric))

        norm = 1 if metric == 'dist_norm1' else 2
        metric_cutoff = eps if metric == 'dist_norm1' else eps*eps

        if isinstance(db_rec, dbscan_internal_record):
            # As long as we update the neighbor counts of
            #  of a core point's neighbors before returning,
            #  we can safely remove it from the rtree for
            #  future searches.  If this is a non-core point,
            #  we still delete it, but add its id to the list of
            #  _possible_border_points in the db_rec of each
            #  of its unlabelled neighbors
            self.rtree.delete(db_rec.id, db_rec._bounds)

        DEBUG.start_timing('rtree_count')
        n = rtree.count(db_rec._neighborhood)
        DEBUG.print_timing('rtree_count')

        DEBUG.start_timing('rtree_intersection')
        ids = rtree.intersection(db_rec._neighborhood, objects=False)
        DEBUG.print_timing('rtree_intersection')

        candidate_ids = np.empty(n, dtype='int64')
        candidate_points = np.empty([n, self.n_dims], dtype=self.dtype)
        i = 0
        for id in ids:
            candidate_ids[i] = id
            candidate_points[i] = records[id]._np_point
            i += 1

        # Avoid catastrophic issues which are extremeley difficult to debug
        assert i == n

        pstats.range_query_candidates += n

        DEBUG.start_timing('within_range_np_linalg')
        neighbors = within_range_np_linalg(
            db_rec._np_point,
            candidate_ids,
            candidate_points,
            eps,
            norm
        )
        DEBUG.print_timing('within_range_np_linalg')

        pstats.range_query_neighbors += len(neighbors)
        end_time = time()
        pstats.range_query_time += (end_time - start_time)

        if isinstance(db_rec, dbscan_internal_record):
            # Necessary because we are deleting it from the rtree
            #  (won't show up in future searches)
            self.update_neighbors(db_rec, neighbors)
            neighbors = np.concatenate([neighbors, np.array(db_rec._possible_border_points, dtype='int64')])
        elif isinstance(db_rec, dbscan_external_record):
            for id in neighbors:
                records[id]._external_neighbors.append(db_rec)
        else:
            raise Exception("range_query() received unexpected type {}".format(type(db_rec)))

        return neighbors

    def update_neighbors(self, db_rec, neighbors):
        db_rec._internal_neighbor_count += len(neighbors)

        not_core_point = db_rec.not_core_point(self.min_samples)

        for n_id in neighbors:
            n_rec = self.records[n_id]
            n_rec._internal_neighbor_count += 1
            if not_core_point and n_rec.cluster_id == label.UNLABELLED:
                n_rec._possible_border_points.append(db_rec.id)


    def label_point(self, rec, cluster):
        """
            Labels a point as a member of a particular cluster, or if
            cluster = None as a NOISE_POINT.
        """
        rec.cluster_id = cluster

        if isinstance(rec, dbscan_internal_record):
            if cluster != label.NOISE_POINT:
                # Mark nearby external points as near this cluster
                for ext_neighbor in rec._external_neighbors:
                    ext_neighbor._nearby_clusters.add(cluster)

    def yield_external_records(self):
        external_records = self.external_records

        for rec in external_records:
            for cluster_id in rec._nearby_clusters:
                rec_clone = dbscan_external_record(rec)
                self.label_point(rec_clone, cluster_id)
                yield rec_clone

DEBUG_WITH_TRACEBACKS = False

# Just for debugging... so we receive the traceback
def dbscan_leaf(db_rec, eps, min_samples, metric, num_internal_points,
                num_external_points, SD, **kwargs):
    res = _dbscan_leaf(db_rec, eps, min_samples, metric, num_internal_points,
        num_external_points, SD)

    if DEBUG_WITH_TRACEBACKS:
        ret = []
        for i, db_rec in enumerate(res):
            ret.append(db_rec)  # Returning a generator breaks the traceback forwarding
        return ret
    else:
        return res

def _dbscan_leaf(db_rec, eps, min_samples, metric, num_internal_points,
                num_external_points, SD):
    """
        Performs classic dbscan algorithm in parallel on each leaf
    """
    global pstats

    num_points = num_internal_points + num_external_points

    db_rec = dbscan_record.from_dict(db_rec, eps)
    id = db_rec.id
    dist_id = db_rec.dist_id

    if dist_id in SD:
        leaf = SD[dist_id]
    else:
        pstats = dbscan_perfstats()
        DEBUG.start_timing("dbscan_leaf")
        DEBUG.start_timing("dbscan_leaf_load_rtree")
        leaf = LeafStorage( len(db_rec.point),
                            db_rec._np_point.dtype,
                            eps,
                            metric,
                            min_samples
                          )
        SD[dist_id] = leaf

    records = leaf.records
    external_records = leaf.external_records

    leaf.add_point(db_rec)

    if len(records) == num_internal_points and \
       len(external_records) == num_external_points:
        DEBUG.print_timing("dbscan_leaf_load_rtree")
        DEBUG.start_timing("dbscan_leaf_external")

        # Find all internal neighbors of external points
        for ext_rec in external_records:
            DEBUG.start_timing("external_rangequery")
            # We don't need to save the result, we only care about
            #  the side effect of adding this point to the list of
            #  neighbor._ext_neighbors for each internal neighbor found
            _ = leaf.range_query(ext_rec)
            DEBUG.print_timing("external_rangequery")

        DEBUG.print_timing("dbscan_leaf_external")
        pstats.show("external")
        pstats = dbscan_perfstats()

        # Classify internal points into clusters (main dbscan algorithm),
        #  and yield resulting internal records
        DEBUG.start_timing("dbscan_leaf_internal")
        cluster = 0
        for rec in records.values():
            if rec.cluster_id != label.UNLABELLED:
                continue

            DEBUG.start_timing('internal_rangequery')
            neighbors = leaf.range_query(rec)
            DEBUG.print_timing('internal_rangequery')

            if not rec.is_core_point: # rec.is_core_point field will be set inside range_query(rec)
                leaf.label_point(rec, label.NOISE_POINT)
                continue

            cluster += 1

            leaf.label_point(rec, cluster)
            yield dbscan_record(rec)

            queue = deque([neighbors])
            DEBUG.start_timing('internal_per_cluster')
            while len(queue) > 0:
                for next_id in queue.popleft():
                    rec2 = records[next_id]
                    if rec2.cluster_id != label.NOISE_POINT and \
                       rec2.cluster_id != label.UNLABELLED:
                        # We delete all points from the unlabelled
                        # rtree as soon as we assign them to a cluster.
                        # But we will still get here if the same point
                        # gets returned by more than one range query
                        # and added to the queue before the first reference
                        # to it gets popped and processed.  The point will
                        # be labelled as a part of this cluster the first time
                        # it gets popped.  We can ignore any duplicate references
                        # to it that get popped after that.
                        continue

                    DEBUG.start_timing('internal_rangequery')
                    neighbors = leaf.range_query(rec2)
                    DEBUG.print_timing('internal_rangequery')

                    if not rec2.is_core_point:  # rec2.is_core_point field will be set inside range_query(rec2)
                        # Label as border point
                        leaf.label_point(rec2, cluster)
                        yield dbscan_record(rec2)
                        continue

                    queue.append(neighbors)
                    # Label as core point
                    leaf.label_point(rec2, cluster)
                    yield dbscan_record(rec2)

            DEBUG.print_timing('internal_per_cluster')

        DEBUG.print_timing("dbscan_leaf_internal")
        pstats.show("internal")

        # Return rows for external records (one for each
        #  cluster an external point is close to)
        DEBUG.start_timing("yield_external_records")
        for ext_rec in leaf.yield_external_records():
           yield ext_rec
        DEBUG.print_timing("yield_external_records")

        # Clean up for next time, if there is one
        del SD[dist_id]

        DEBUG.print_timing("dbscan_leaf")

# The snowflake table is used to reduce the size of the edge table.
# Note that the sole purpose of the edge table is finding the connected
# components. Which means removing some of the edges is fine as long as the
# component is intact.

# We call it snowflake because the end result will look like a point in the
# middle with a bunch of edges coming out of it to the other connected points.

# Example:
# Edges: [1,2] [1,3] [1,4] [2,3] [2,4] [3,1] [3,4] [5,6] [6,7] [7,8] [7,5]
# Result: 1 and 5 are snowflake cores
# Edges: [1,2] [1,3] [1,4] [5,6] [5,7] [5,8]

# This is a proto-wcc operation (creates connected components) but we still
# need to run the actual wcc to combine snowflakes from different segments.

def sf_transition(state, src, dest, n_rows, gp_segment_id, **kwargs):

    SD = kwargs['SD']
    if not state:
        data = []
        SD['counter'] = 0
    else:
        data = SD['data']

    counter = SD['counter']

    data.append([src,dest])
    ret = [[-1,-1],[-1,-1]]

    my_n_rows = n_rows[gp_segment_id]

    if len(data) == my_n_rows:

        cl_ids = {}
        clid_counter = 1
        cl_counts = {}
        for i in data:

            key1 = i[0]
            key2 = i[1]

            cl_id_k1 = cl_ids[key1] if key1 in cl_ids else None
            cl_id_k2 = cl_ids[key2] if key2 in cl_ids else None

            if not cl_id_k1 and not cl_id_k2:
                cl_ids[key1] = clid_counter
                cl_ids[key2] = clid_counter
                cl_counts[clid_counter] = 2
                clid_counter += 1
            elif cl_id_k1 and not cl_id_k2:
                cl_ids[key2] = cl_id_k1
                cl_counts[cl_ids[key1]] += 1
            elif not cl_id_k1 and cl_id_k2:
                cl_ids[key1] = cl_id_k2
                cl_counts[cl_ids[key2]] += 1
            else:
                if cl_id_k1 != cl_id_k2:
                    if cl_counts[cl_id_k1] > cl_counts[cl_id_k2]:
                        for chkey,chvalue in cl_ids.items():
                            if chvalue == cl_id_k2:
                                cl_ids[chkey] = cl_id_k1
                        cl_counts[cl_id_k1] += cl_counts[cl_id_k2]
                        cl_counts[cl_id_k2] = 0
                    else:
                        for chkey,chvalue in cl_ids.items():
                            if chvalue == cl_id_k1:
                                cl_ids[chkey] = cl_id_k2
                        cl_counts[cl_id_k2] += cl_counts[cl_id_k1]
                        cl_counts[cl_id_k1] = 0

        if cl_ids:

            running_cl_id = -1
            running_sf_center = -1
            ret = []
            for vertex_id, vertex_cl in sorted(cl_ids.items(), key=lambda item: item[1]):

                # Check if we are still in the same snowflake
                if vertex_cl != running_cl_id:

                    running_sf_center = vertex_id
                    running_cl_id = vertex_cl
                else:
                    ret.append([running_sf_center,vertex_id])

    SD['data'] = data

    return ret

def sf_merge(state1, state2, **kwargs):

    if state1:
        return state1
    else:
        return state2

def sf_final(state, **kwargs):

    return state

def _validate_dbscan(schema_madlib, source_table, output_table, id_column,
    expr_point, eps, min_samples, metric, algorithm, max_depth):

    input_tbl_valid(source_table, 'dbscan')
    output_tbl_valid(output_table, 'dbscan')
    output_summary_table = add_postfix(output_table, '_summary')
    output_tbl_valid(output_summary_table, 'dbscan')

    _assert(is_var_valid(source_table, id_column),
            "dbscan error: {0} is an invalid column or "
            "expression for id param".format(id_column))

    id_col_type = get_expr_type(id_column, source_table)
    _assert(is_valid_psql_type(id_col_type, INTEGER),
            "dbscan Error: unique id column or expression '{0}' in input table must"
            " evaluate to an integer type.".format(id_column))

    _assert(is_var_valid(source_table, expr_point),
            "dbscan error: {0} is an invalid column name or "
            "expression for expr_point param".format(expr_point))

    point_col_type = get_expr_type(expr_point, source_table)
    _assert(is_valid_psql_type(point_col_type, NUMERIC | ONLY_ARRAY),
            "dbscan Error: Feature column or expression '{0}' in input table is not"
            " a numeric array.".format(expr_point))

    _assert(eps > 0, "dbscan Error: eps has to be a positive number")

    _assert(min_samples > 0, "dbscan Error: min_samples has to be a positive number")

    fn_dist_list = ['dist_norm1', 'dist_norm2', 'squared_dist_norm2',
                    'dist_angle', 'dist_tanimoto']
    _assert(metric in fn_dist_list,
        "dbscan Error: metric has to be one of the madlib defined distance functions")

    _assert(algorithm == METHOD_BRUTE_FORCE or RTREE_ENABLED == 1,
        "dbscan Error: Cannot use rtree without the necessary python module: rtree")
    _assert(max_depth is None or max_depth >= 0,
        "dbscan Error: max_segmentation_depth must be a non-negative integer or NULL")

def _validate_dbscan_predict(schema_madlib, dbscan_table, source_table,
    id_column, expr_point, output_table):

    input_tbl_valid(source_table, 'dbscan')
    input_tbl_valid(dbscan_table, 'dbscan')
    dbscan_summary_table = add_postfix(dbscan_table, '_summary')
    input_tbl_valid(dbscan_summary_table, 'dbscan')
    output_tbl_valid(output_table, 'dbscan')

    _assert(is_var_valid(source_table, id_column),
            "dbscan error: {0} is an invalid column name or "
            "expression for id param".format(id_column))

    id_col_type = get_expr_type(id_column, source_table)
    _assert(is_valid_psql_type(id_col_type, INTEGER),
            "dbscan Error: unique id column or expression '{0}' in input table must"
            " evaluate to an integer type.".format(id_column))

    _assert(is_var_valid(source_table, expr_point),
            "dbscan error: {0} is an invalid column name or "
            "expression for expr_point param".format(expr_point))

    point_col_type = get_expr_type(expr_point, source_table)
    _assert(is_valid_psql_type(point_col_type, NUMERIC | ONLY_ARRAY),
            "dbscan Error: Feature column or expression '{0}' in train table is not"
            " a numeric array.".format(expr_point))

def dbscan_help(schema_madlib, message=None, **kwargs):
    """
    Help function for dbscan

    Args:
        @param schema_madlib
        @param message: string, Help message string
        @param kwargs

    Returns:
        String. Help/usage information
    """
    if message is not None and \
            message.lower() in ("usage", "help", "?"):
        help_string = """
-----------------------------------------------------------------------
                            USAGE
-----------------------------------------------------------------------
SELECT {schema_madlib}.dbscan(
    source_table,       -- Name of the training data table
    output_table,       -- Name of the output table
    id_column,          -- A column or expression of columns specifying a unique id for rows in the source_table
    expr_point,         -- Column name or expression for data points
    eps,                -- The minimum radius of a cluster
    min_samples,        -- The minimum size of a cluster
    metric,             -- The name of the function to use to calculate the
                        -- distance
    algorithm,          -- The algorithm to use for dbscan: brute_force or optimized (default=optimized)
    max_segmentation_depth  -- This has no effect on method=brute.  For method=rtree, this
                               limits the number of splits to be considered while dividing
                               the dataset into smaller spacial regions for optimal parallel performance.
                               The default ( NULL or ommitted ) value means it will be unlimited,
                               but note that it will never break the data up into more regions than
                               segments in a greenplum cluster, and in some cases may choose fewer.
                               Setting this to non-NULL will decrease up-front overhead, and can
                               be useful if you expect to find only a small number of densely populated
                               clusters, or if the dataset itself is small enough where initial planning
                               overhead matters more than scaling performance.  This parameter will
                               not affect the output, only the performance.
    );

-----------------------------------------------------------------------
                            OUTPUT
-----------------------------------------------------------------------
The output of the dbscan function is a table with the following columns:

id                  The id of each point which has been assigned to a cluster
cluster_id          The cluster_id's assigned to each point
is_core_point       Boolean column that indicates if the point is core or not
point               The coordinates of each classified point
"""
    else:
        help_string = """
----------------------------------------------------------------------------
                                SUMMARY
----------------------------------------------------------------------------
DBSCAN is a density-based clustering algorithm. Given a set of points in
some space, it groups together points that are closely packed together
(points with many nearby neighbors), marking as outliers points that lie
alone in low-density regions (whose nearest neighbors are too far away).
--
For an overview on usage, run:
SELECT {schema_madlib}.dbscan('usage');
SELECT {schema_madlib}.dbscan_predict('usage');
"""
    return help_string.format(schema_madlib=schema_madlib)
# ------------------------------------------------------------------------------

def dbscan_predict_help(schema_madlib, message=None, **kwargs):
    """
    Help function for dbscan

    Args:
        @param schema_madlib
        @param message: string, Help message string
        @param kwargs

    Returns:
        String. Help/usage information
    """
    if message is not None and \
            message.lower() in ("usage", "help", "?"):
        help_string = """
-----------------------------------------------------------------------
                            USAGE
-----------------------------------------------------------------------
SELECT {schema_madlib}.dbscan_predict(
    dbscan_table,       -- Name of the tdbscan output table
    new_point           -- Double precision array representing the point
                        -- for prediction
    );

-----------------------------------------------------------------------
                            OUTPUT
-----------------------------------------------------------------------
The output of the dbscan_predict is an integer indicating the cluster_id
of given point
"""
    else:
        help_string = """
----------------------------------------------------------------------------
                                SUMMARY
----------------------------------------------------------------------------
DBSCAN is a density-based clustering algorithm. Given a set of points in
some space, it groups together points that are closely packed together
(points with many nearby neighbors), marking as outliers points that lie
alone in low-density regions (whose nearest neighbors are too far away).
--
For an overview on usage, run:
SELECT {schema_madlib}.dbscan('usage');
SELECT {schema_madlib}.dbscan_predict('usage');
"""
    return help_string.format(schema_madlib=schema_madlib)
# ------------------------------------------------------------------------------