release(camera-ai): publish 0.4.12 catalog portability fix
This commit is contained in:
KINDNICK
@@ -1,41 +1,41 @@
|
||||
# Fast inner loop for DBSCAN.
|
||||
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
from libcpp.vector cimport vector
|
||||
|
||||
from ..utils._typedefs cimport uint8_t, intp_t
|
||||
|
||||
|
||||
def dbscan_inner(const uint8_t[::1] is_core,
|
||||
object[:] neighborhoods,
|
||||
intp_t[::1] labels):
|
||||
cdef intp_t i, label_num = 0, v
|
||||
cdef intp_t[:] neighb
|
||||
cdef vector[intp_t] stack
|
||||
|
||||
for i in range(labels.shape[0]):
|
||||
if labels[i] != -1 or not is_core[i]:
|
||||
continue
|
||||
|
||||
# Depth-first search starting from i, ending at the non-core points.
|
||||
# This is very similar to the classic algorithm for computing connected
|
||||
# components, the difference being that we label non-core points as
|
||||
# part of a cluster (component), but don't expand their neighborhoods.
|
||||
while True:
|
||||
if labels[i] == -1:
|
||||
labels[i] = label_num
|
||||
if is_core[i]:
|
||||
neighb = neighborhoods[i]
|
||||
for i in range(neighb.shape[0]):
|
||||
v = neighb[i]
|
||||
if labels[v] == -1:
|
||||
stack.push_back(v)
|
||||
|
||||
if stack.size() == 0:
|
||||
break
|
||||
i = stack.back()
|
||||
stack.pop_back()
|
||||
|
||||
label_num += 1
|
||||
# Fast inner loop for DBSCAN.
|
||||
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
from libcpp.vector cimport vector
|
||||
|
||||
from ..utils._typedefs cimport uint8_t, intp_t
|
||||
|
||||
|
||||
def dbscan_inner(const uint8_t[::1] is_core,
|
||||
object[:] neighborhoods,
|
||||
intp_t[::1] labels):
|
||||
cdef intp_t i, label_num = 0, v
|
||||
cdef intp_t[:] neighb
|
||||
cdef vector[intp_t] stack
|
||||
|
||||
for i in range(labels.shape[0]):
|
||||
if labels[i] != -1 or not is_core[i]:
|
||||
continue
|
||||
|
||||
# Depth-first search starting from i, ending at the non-core points.
|
||||
# This is very similar to the classic algorithm for computing connected
|
||||
# components, the difference being that we label non-core points as
|
||||
# part of a cluster (component), but don't expand their neighborhoods.
|
||||
while True:
|
||||
if labels[i] == -1:
|
||||
labels[i] = label_num
|
||||
if is_core[i]:
|
||||
neighb = neighborhoods[i]
|
||||
for i in range(neighb.shape[0]):
|
||||
v = neighb[i]
|
||||
if labels[v] == -1:
|
||||
stack.push_back(v)
|
||||
|
||||
if stack.size() == 0:
|
||||
break
|
||||
i = stack.back()
|
||||
stack.pop_back()
|
||||
|
||||
label_num += 1
|
||||
|
||||
@@ -1,274 +1,274 @@
|
||||
# Minimum spanning tree single linkage implementation for hdbscan
|
||||
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
# Redistribution and use in source and binary forms, with or without
|
||||
# modification, are permitted provided that the following conditions are met:
|
||||
|
||||
# 1. Redistributions of source code must retain the above copyright notice,
|
||||
# this list of conditions and the following disclaimer.
|
||||
|
||||
# 2. Redistributions in binary form must reproduce the above copyright notice,
|
||||
# this list of conditions and the following disclaimer in the documentation
|
||||
# and/or other materials provided with the distribution.
|
||||
|
||||
# 3. Neither the name of the copyright holder nor the names of its contributors
|
||||
# may be used to endorse or promote products derived from this software without
|
||||
# specific prior written permission.
|
||||
|
||||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
# POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
cimport numpy as cnp
|
||||
from libc.float cimport DBL_MAX
|
||||
|
||||
import numpy as np
|
||||
from ...metrics._dist_metrics cimport DistanceMetric64
|
||||
from ...cluster._hierarchical_fast cimport UnionFind
|
||||
from ...cluster._hdbscan._tree cimport HIERARCHY_t
|
||||
from ...cluster._hdbscan._tree import HIERARCHY_dtype
|
||||
from ...utils._typedefs cimport intp_t, float64_t, int64_t, uint8_t
|
||||
|
||||
cnp.import_array()
|
||||
|
||||
cdef extern from "numpy/arrayobject.h":
|
||||
intp_t * PyArray_SHAPE(cnp.PyArrayObject *)
|
||||
|
||||
# Numpy structured dtype representing a single ordered edge in Prim's algorithm
|
||||
MST_edge_dtype = np.dtype([
|
||||
("current_node", np.int64),
|
||||
("next_node", np.int64),
|
||||
("distance", np.float64),
|
||||
])
|
||||
|
||||
# Packed shouldn't make a difference since they're all 8-byte quantities,
|
||||
# but it's included just to be safe.
|
||||
ctypedef packed struct MST_edge_t:
|
||||
int64_t current_node
|
||||
int64_t next_node
|
||||
float64_t distance
|
||||
|
||||
cpdef cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst_from_mutual_reachability(
|
||||
cnp.ndarray[float64_t, ndim=2] mutual_reachability
|
||||
):
|
||||
"""Compute the Minimum Spanning Tree (MST) representation of the mutual-
|
||||
reachability graph using Prim's algorithm.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mutual_reachability : ndarray of shape (n_samples, n_samples)
|
||||
Array of mutual-reachabilities between samples.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mst : ndarray of shape (n_samples - 1,), dtype=MST_edge_dtype
|
||||
The MST representation of the mutual-reachability graph. The MST is
|
||||
represented as a collection of edges.
|
||||
"""
|
||||
cdef:
|
||||
# Note: we utilize ndarray's over memory-views to make use of numpy
|
||||
# binary indexing and sub-selection below.
|
||||
cnp.ndarray[int64_t, ndim=1, mode='c'] current_labels
|
||||
cnp.ndarray[float64_t, ndim=1, mode='c'] min_reachability, left, right
|
||||
cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst
|
||||
|
||||
cnp.ndarray[uint8_t, mode='c'] label_filter
|
||||
|
||||
int64_t n_samples = PyArray_SHAPE(<cnp.PyArrayObject*> mutual_reachability)[0]
|
||||
int64_t current_node, new_node_index, new_node, i
|
||||
|
||||
mst = np.empty(n_samples - 1, dtype=MST_edge_dtype)
|
||||
current_labels = np.arange(n_samples, dtype=np.int64)
|
||||
current_node = 0
|
||||
min_reachability = np.full(n_samples, fill_value=np.inf, dtype=np.float64)
|
||||
for i in range(0, n_samples - 1):
|
||||
label_filter = current_labels != current_node
|
||||
current_labels = current_labels[label_filter]
|
||||
left = min_reachability[label_filter]
|
||||
right = mutual_reachability[current_node][current_labels]
|
||||
min_reachability = np.minimum(left, right)
|
||||
|
||||
new_node_index = np.argmin(min_reachability)
|
||||
new_node = current_labels[new_node_index]
|
||||
mst[i].current_node = current_node
|
||||
mst[i].next_node = new_node
|
||||
mst[i].distance = min_reachability[new_node_index]
|
||||
current_node = new_node
|
||||
|
||||
return mst
|
||||
|
||||
|
||||
cpdef cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst_from_data_matrix(
|
||||
const float64_t[:, ::1] raw_data,
|
||||
const float64_t[::1] core_distances,
|
||||
DistanceMetric64 dist_metric,
|
||||
float64_t alpha=1.0
|
||||
):
|
||||
"""Compute the Minimum Spanning Tree (MST) representation of the mutual-
|
||||
reachability graph generated from the provided `raw_data` and
|
||||
`core_distances` using Prim's algorithm.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
raw_data : ndarray of shape (n_samples, n_features)
|
||||
Input array of data samples.
|
||||
|
||||
core_distances : ndarray of shape (n_samples,)
|
||||
An array containing the core-distance calculated for each corresponding
|
||||
sample.
|
||||
|
||||
dist_metric : DistanceMetric
|
||||
The distance metric to use when calculating pairwise distances for
|
||||
determining mutual-reachability.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mst : ndarray of shape (n_samples - 1,), dtype=MST_edge_dtype
|
||||
The MST representation of the mutual-reachability graph. The MST is
|
||||
represented as a collection of edges.
|
||||
"""
|
||||
|
||||
cdef:
|
||||
uint8_t[::1] in_tree
|
||||
float64_t[::1] min_reachability
|
||||
int64_t[::1] current_sources
|
||||
cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst
|
||||
|
||||
int64_t current_node, source_node, new_node, next_node_source
|
||||
int64_t i, j, n_samples, num_features
|
||||
|
||||
float64_t current_node_core_dist, new_reachability, mutual_reachability_distance
|
||||
float64_t next_node_min_reach, pair_distance, next_node_core_dist
|
||||
|
||||
n_samples = raw_data.shape[0]
|
||||
num_features = raw_data.shape[1]
|
||||
|
||||
mst = np.empty(n_samples - 1, dtype=MST_edge_dtype)
|
||||
|
||||
in_tree = np.zeros(n_samples, dtype=np.uint8)
|
||||
min_reachability = np.full(n_samples, fill_value=np.inf, dtype=np.float64)
|
||||
current_sources = np.ones(n_samples, dtype=np.int64)
|
||||
|
||||
current_node = 0
|
||||
|
||||
# The following loop dynamically updates minimum reachability node-by-node,
|
||||
# avoiding unnecessary computation where possible.
|
||||
for i in range(0, n_samples - 1):
|
||||
|
||||
in_tree[current_node] = 1
|
||||
|
||||
current_node_core_dist = core_distances[current_node]
|
||||
|
||||
new_reachability = DBL_MAX
|
||||
source_node = 0
|
||||
new_node = 0
|
||||
|
||||
for j in range(n_samples):
|
||||
if in_tree[j]:
|
||||
continue
|
||||
|
||||
next_node_min_reach = min_reachability[j]
|
||||
next_node_source = current_sources[j]
|
||||
|
||||
pair_distance = dist_metric.dist(
|
||||
&raw_data[current_node, 0],
|
||||
&raw_data[j, 0],
|
||||
num_features
|
||||
)
|
||||
|
||||
pair_distance /= alpha
|
||||
|
||||
next_node_core_dist = core_distances[j]
|
||||
mutual_reachability_distance = max(
|
||||
current_node_core_dist,
|
||||
next_node_core_dist,
|
||||
pair_distance
|
||||
)
|
||||
|
||||
# If MRD(i, j) is smaller than node j's min_reachability, we update
|
||||
# node j's min_reachability for future reference.
|
||||
if mutual_reachability_distance < next_node_min_reach:
|
||||
min_reachability[j] = mutual_reachability_distance
|
||||
current_sources[j] = current_node
|
||||
|
||||
# If MRD(i, j) is also smaller than node i's current
|
||||
# min_reachability, we update and set their edge as the current
|
||||
# MST edge candidate.
|
||||
if mutual_reachability_distance < new_reachability:
|
||||
new_reachability = mutual_reachability_distance
|
||||
source_node = current_node
|
||||
new_node = j
|
||||
|
||||
# If the node j is closer to another node already in the tree, we
|
||||
# make their edge the current MST candidate edge.
|
||||
elif next_node_min_reach < new_reachability:
|
||||
new_reachability = next_node_min_reach
|
||||
source_node = next_node_source
|
||||
new_node = j
|
||||
|
||||
mst[i].current_node = source_node
|
||||
mst[i].next_node = new_node
|
||||
mst[i].distance = new_reachability
|
||||
current_node = new_node
|
||||
|
||||
return mst
|
||||
|
||||
cpdef cnp.ndarray[HIERARCHY_t, ndim=1, mode="c"] make_single_linkage(const MST_edge_t[::1] mst):
|
||||
"""Construct a single-linkage tree from an MST.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mst : ndarray of shape (n_samples - 1,), dtype=MST_edge_dtype
|
||||
The MST representation of the mutual-reachability graph. The MST is
|
||||
represented as a collection of edges.
|
||||
|
||||
Returns
|
||||
-------
|
||||
single_linkage : ndarray of shape (n_samples - 1,), dtype=HIERARCHY_dtype
|
||||
The single-linkage tree tree (dendrogram) built from the MST. Each
|
||||
of the array represents the following:
|
||||
|
||||
- left node/cluster
|
||||
- right node/cluster
|
||||
- distance
|
||||
- new cluster size
|
||||
"""
|
||||
cdef:
|
||||
cnp.ndarray[HIERARCHY_t, ndim=1, mode="c"] single_linkage
|
||||
|
||||
# Note mst.shape[0] is one fewer than the number of samples
|
||||
int64_t n_samples = mst.shape[0] + 1
|
||||
intp_t current_node_cluster, next_node_cluster
|
||||
int64_t current_node, next_node, i
|
||||
float64_t distance
|
||||
UnionFind U = UnionFind(n_samples)
|
||||
|
||||
single_linkage = np.zeros(n_samples - 1, dtype=HIERARCHY_dtype)
|
||||
|
||||
for i in range(n_samples - 1):
|
||||
|
||||
current_node = mst[i].current_node
|
||||
next_node = mst[i].next_node
|
||||
distance = mst[i].distance
|
||||
|
||||
current_node_cluster = U.fast_find(current_node)
|
||||
next_node_cluster = U.fast_find(next_node)
|
||||
|
||||
single_linkage[i].left_node = current_node_cluster
|
||||
single_linkage[i].right_node = next_node_cluster
|
||||
single_linkage[i].value = distance
|
||||
single_linkage[i].cluster_size = U.size[current_node_cluster] + U.size[next_node_cluster]
|
||||
|
||||
U.union(current_node_cluster, next_node_cluster)
|
||||
|
||||
return single_linkage
|
||||
# Minimum spanning tree single linkage implementation for hdbscan
|
||||
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
# Redistribution and use in source and binary forms, with or without
|
||||
# modification, are permitted provided that the following conditions are met:
|
||||
|
||||
# 1. Redistributions of source code must retain the above copyright notice,
|
||||
# this list of conditions and the following disclaimer.
|
||||
|
||||
# 2. Redistributions in binary form must reproduce the above copyright notice,
|
||||
# this list of conditions and the following disclaimer in the documentation
|
||||
# and/or other materials provided with the distribution.
|
||||
|
||||
# 3. Neither the name of the copyright holder nor the names of its contributors
|
||||
# may be used to endorse or promote products derived from this software without
|
||||
# specific prior written permission.
|
||||
|
||||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
# POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
cimport numpy as cnp
|
||||
from libc.float cimport DBL_MAX
|
||||
|
||||
import numpy as np
|
||||
from ...metrics._dist_metrics cimport DistanceMetric64
|
||||
from ...cluster._hierarchical_fast cimport UnionFind
|
||||
from ...cluster._hdbscan._tree cimport HIERARCHY_t
|
||||
from ...cluster._hdbscan._tree import HIERARCHY_dtype
|
||||
from ...utils._typedefs cimport intp_t, float64_t, int64_t, uint8_t
|
||||
|
||||
cnp.import_array()
|
||||
|
||||
cdef extern from "numpy/arrayobject.h":
|
||||
intp_t * PyArray_SHAPE(cnp.PyArrayObject *)
|
||||
|
||||
# Numpy structured dtype representing a single ordered edge in Prim's algorithm
|
||||
MST_edge_dtype = np.dtype([
|
||||
("current_node", np.int64),
|
||||
("next_node", np.int64),
|
||||
("distance", np.float64),
|
||||
])
|
||||
|
||||
# Packed shouldn't make a difference since they're all 8-byte quantities,
|
||||
# but it's included just to be safe.
|
||||
ctypedef packed struct MST_edge_t:
|
||||
int64_t current_node
|
||||
int64_t next_node
|
||||
float64_t distance
|
||||
|
||||
cpdef cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst_from_mutual_reachability(
|
||||
cnp.ndarray[float64_t, ndim=2] mutual_reachability
|
||||
):
|
||||
"""Compute the Minimum Spanning Tree (MST) representation of the mutual-
|
||||
reachability graph using Prim's algorithm.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mutual_reachability : ndarray of shape (n_samples, n_samples)
|
||||
Array of mutual-reachabilities between samples.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mst : ndarray of shape (n_samples - 1,), dtype=MST_edge_dtype
|
||||
The MST representation of the mutual-reachability graph. The MST is
|
||||
represented as a collection of edges.
|
||||
"""
|
||||
cdef:
|
||||
# Note: we utilize ndarray's over memory-views to make use of numpy
|
||||
# binary indexing and sub-selection below.
|
||||
cnp.ndarray[int64_t, ndim=1, mode='c'] current_labels
|
||||
cnp.ndarray[float64_t, ndim=1, mode='c'] min_reachability, left, right
|
||||
cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst
|
||||
|
||||
cnp.ndarray[uint8_t, mode='c'] label_filter
|
||||
|
||||
int64_t n_samples = PyArray_SHAPE(<cnp.PyArrayObject*> mutual_reachability)[0]
|
||||
int64_t current_node, new_node_index, new_node, i
|
||||
|
||||
mst = np.empty(n_samples - 1, dtype=MST_edge_dtype)
|
||||
current_labels = np.arange(n_samples, dtype=np.int64)
|
||||
current_node = 0
|
||||
min_reachability = np.full(n_samples, fill_value=np.inf, dtype=np.float64)
|
||||
for i in range(0, n_samples - 1):
|
||||
label_filter = current_labels != current_node
|
||||
current_labels = current_labels[label_filter]
|
||||
left = min_reachability[label_filter]
|
||||
right = mutual_reachability[current_node][current_labels]
|
||||
min_reachability = np.minimum(left, right)
|
||||
|
||||
new_node_index = np.argmin(min_reachability)
|
||||
new_node = current_labels[new_node_index]
|
||||
mst[i].current_node = current_node
|
||||
mst[i].next_node = new_node
|
||||
mst[i].distance = min_reachability[new_node_index]
|
||||
current_node = new_node
|
||||
|
||||
return mst
|
||||
|
||||
|
||||
cpdef cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst_from_data_matrix(
|
||||
const float64_t[:, ::1] raw_data,
|
||||
const float64_t[::1] core_distances,
|
||||
DistanceMetric64 dist_metric,
|
||||
float64_t alpha=1.0
|
||||
):
|
||||
"""Compute the Minimum Spanning Tree (MST) representation of the mutual-
|
||||
reachability graph generated from the provided `raw_data` and
|
||||
`core_distances` using Prim's algorithm.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
raw_data : ndarray of shape (n_samples, n_features)
|
||||
Input array of data samples.
|
||||
|
||||
core_distances : ndarray of shape (n_samples,)
|
||||
An array containing the core-distance calculated for each corresponding
|
||||
sample.
|
||||
|
||||
dist_metric : DistanceMetric
|
||||
The distance metric to use when calculating pairwise distances for
|
||||
determining mutual-reachability.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mst : ndarray of shape (n_samples - 1,), dtype=MST_edge_dtype
|
||||
The MST representation of the mutual-reachability graph. The MST is
|
||||
represented as a collection of edges.
|
||||
"""
|
||||
|
||||
cdef:
|
||||
uint8_t[::1] in_tree
|
||||
float64_t[::1] min_reachability
|
||||
int64_t[::1] current_sources
|
||||
cnp.ndarray[MST_edge_t, ndim=1, mode='c'] mst
|
||||
|
||||
int64_t current_node, source_node, new_node, next_node_source
|
||||
int64_t i, j, n_samples, num_features
|
||||
|
||||
float64_t current_node_core_dist, new_reachability, mutual_reachability_distance
|
||||
float64_t next_node_min_reach, pair_distance, next_node_core_dist
|
||||
|
||||
n_samples = raw_data.shape[0]
|
||||
num_features = raw_data.shape[1]
|
||||
|
||||
mst = np.empty(n_samples - 1, dtype=MST_edge_dtype)
|
||||
|
||||
in_tree = np.zeros(n_samples, dtype=np.uint8)
|
||||
min_reachability = np.full(n_samples, fill_value=np.inf, dtype=np.float64)
|
||||
current_sources = np.ones(n_samples, dtype=np.int64)
|
||||
|
||||
current_node = 0
|
||||
|
||||
# The following loop dynamically updates minimum reachability node-by-node,
|
||||
# avoiding unnecessary computation where possible.
|
||||
for i in range(0, n_samples - 1):
|
||||
|
||||
in_tree[current_node] = 1
|
||||
|
||||
current_node_core_dist = core_distances[current_node]
|
||||
|
||||
new_reachability = DBL_MAX
|
||||
source_node = 0
|
||||
new_node = 0
|
||||
|
||||
for j in range(n_samples):
|
||||
if in_tree[j]:
|
||||
continue
|
||||
|
||||
next_node_min_reach = min_reachability[j]
|
||||
next_node_source = current_sources[j]
|
||||
|
||||
pair_distance = dist_metric.dist(
|
||||
&raw_data[current_node, 0],
|
||||
&raw_data[j, 0],
|
||||
num_features
|
||||
)
|
||||
|
||||
pair_distance /= alpha
|
||||
|
||||
next_node_core_dist = core_distances[j]
|
||||
mutual_reachability_distance = max(
|
||||
current_node_core_dist,
|
||||
next_node_core_dist,
|
||||
pair_distance
|
||||
)
|
||||
|
||||
# If MRD(i, j) is smaller than node j's min_reachability, we update
|
||||
# node j's min_reachability for future reference.
|
||||
if mutual_reachability_distance < next_node_min_reach:
|
||||
min_reachability[j] = mutual_reachability_distance
|
||||
current_sources[j] = current_node
|
||||
|
||||
# If MRD(i, j) is also smaller than node i's current
|
||||
# min_reachability, we update and set their edge as the current
|
||||
# MST edge candidate.
|
||||
if mutual_reachability_distance < new_reachability:
|
||||
new_reachability = mutual_reachability_distance
|
||||
source_node = current_node
|
||||
new_node = j
|
||||
|
||||
# If the node j is closer to another node already in the tree, we
|
||||
# make their edge the current MST candidate edge.
|
||||
elif next_node_min_reach < new_reachability:
|
||||
new_reachability = next_node_min_reach
|
||||
source_node = next_node_source
|
||||
new_node = j
|
||||
|
||||
mst[i].current_node = source_node
|
||||
mst[i].next_node = new_node
|
||||
mst[i].distance = new_reachability
|
||||
current_node = new_node
|
||||
|
||||
return mst
|
||||
|
||||
cpdef cnp.ndarray[HIERARCHY_t, ndim=1, mode="c"] make_single_linkage(const MST_edge_t[::1] mst):
|
||||
"""Construct a single-linkage tree from an MST.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mst : ndarray of shape (n_samples - 1,), dtype=MST_edge_dtype
|
||||
The MST representation of the mutual-reachability graph. The MST is
|
||||
represented as a collection of edges.
|
||||
|
||||
Returns
|
||||
-------
|
||||
single_linkage : ndarray of shape (n_samples - 1,), dtype=HIERARCHY_dtype
|
||||
The single-linkage tree tree (dendrogram) built from the MST. Each
|
||||
of the array represents the following:
|
||||
|
||||
- left node/cluster
|
||||
- right node/cluster
|
||||
- distance
|
||||
- new cluster size
|
||||
"""
|
||||
cdef:
|
||||
cnp.ndarray[HIERARCHY_t, ndim=1, mode="c"] single_linkage
|
||||
|
||||
# Note mst.shape[0] is one fewer than the number of samples
|
||||
int64_t n_samples = mst.shape[0] + 1
|
||||
intp_t current_node_cluster, next_node_cluster
|
||||
int64_t current_node, next_node, i
|
||||
float64_t distance
|
||||
UnionFind U = UnionFind(n_samples)
|
||||
|
||||
single_linkage = np.zeros(n_samples - 1, dtype=HIERARCHY_dtype)
|
||||
|
||||
for i in range(n_samples - 1):
|
||||
|
||||
current_node = mst[i].current_node
|
||||
next_node = mst[i].next_node
|
||||
distance = mst[i].distance
|
||||
|
||||
current_node_cluster = U.fast_find(current_node)
|
||||
next_node_cluster = U.fast_find(next_node)
|
||||
|
||||
single_linkage[i].left_node = current_node_cluster
|
||||
single_linkage[i].right_node = next_node_cluster
|
||||
single_linkage[i].value = distance
|
||||
single_linkage[i].cluster_size = U.size[current_node_cluster] + U.size[next_node_cluster]
|
||||
|
||||
U.union(current_node_cluster, next_node_cluster)
|
||||
|
||||
return single_linkage
|
||||
|
||||
+210
-210
@@ -1,210 +1,210 @@
|
||||
# mutual reachability distance computations
|
||||
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
# Redistribution and use in source and binary forms, with or without
|
||||
# modification, are permitted provided that the following conditions are met:
|
||||
|
||||
# 1. Redistributions of source code must retain the above copyright notice,
|
||||
# this list of conditions and the following disclaimer.
|
||||
|
||||
# 2. Redistributions in binary form must reproduce the above copyright notice,
|
||||
# this list of conditions and the following disclaimer in the documentation
|
||||
# and/or other materials provided with the distribution.
|
||||
|
||||
# 3. Neither the name of the copyright holder nor the names of its contributors
|
||||
# may be used to endorse or promote products derived from this software without
|
||||
# specific prior written permission.
|
||||
|
||||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
# POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
cimport numpy as cnp
|
||||
|
||||
import numpy as np
|
||||
from scipy.sparse import issparse
|
||||
from cython cimport floating, integral
|
||||
from libc.math cimport isfinite, INFINITY
|
||||
from ...utils._typedefs cimport intp_t
|
||||
cnp.import_array()
|
||||
|
||||
|
||||
def mutual_reachability_graph(
|
||||
distance_matrix, min_samples=5, max_distance=0.0
|
||||
):
|
||||
"""Compute the weighted adjacency matrix of the mutual reachability graph.
|
||||
|
||||
The mutual reachability distance used to build the graph is defined as::
|
||||
|
||||
max(d_core(x_p), d_core(x_q), d(x_p, x_q))
|
||||
|
||||
and the core distance `d_core` is defined as the distance between a point
|
||||
`x_p` and its k-th nearest neighbor.
|
||||
|
||||
Note that all computations are done in-place.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
distance_matrix : {ndarray, sparse matrix} of shape (n_samples, n_samples)
|
||||
Array of distances between samples. If sparse, the array must be in
|
||||
`CSR` format.
|
||||
|
||||
min_samples : int, default=5
|
||||
The parameter `k` used to calculate the distance between a point
|
||||
`x_p` and its k-th nearest neighbor.
|
||||
|
||||
max_distance : float, default=0.0
|
||||
The distance which `np.inf` is replaced with. When the true mutual-
|
||||
reachability distance is measured to be infinite, it is instead
|
||||
truncated to `max_dist`. Only used when `distance_matrix` is a sparse
|
||||
matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mututal_reachability_graph: {ndarray, sparse matrix} of shape \
|
||||
(n_samples, n_samples)
|
||||
Weighted adjacency matrix of the mutual reachability graph.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Campello, R. J., Moulavi, D., & Sander, J. (2013, April).
|
||||
Density-based clustering based on hierarchical density estimates.
|
||||
In Pacific-Asia Conference on Knowledge Discovery and Data Mining
|
||||
(pp. 160-172). Springer Berlin Heidelberg.
|
||||
"""
|
||||
further_neighbor_idx = min_samples - 1
|
||||
if issparse(distance_matrix):
|
||||
if distance_matrix.format != "csr":
|
||||
raise ValueError(
|
||||
"Only sparse CSR matrices are supported for `distance_matrix`."
|
||||
)
|
||||
_sparse_mutual_reachability_graph(
|
||||
distance_matrix.data,
|
||||
distance_matrix.indices,
|
||||
distance_matrix.indptr,
|
||||
distance_matrix.shape[0],
|
||||
further_neighbor_idx=further_neighbor_idx,
|
||||
max_distance=max_distance,
|
||||
)
|
||||
else:
|
||||
_dense_mutual_reachability_graph(
|
||||
distance_matrix, further_neighbor_idx=further_neighbor_idx
|
||||
)
|
||||
return distance_matrix
|
||||
|
||||
|
||||
def _dense_mutual_reachability_graph(
|
||||
floating[:, :] distance_matrix,
|
||||
intp_t further_neighbor_idx,
|
||||
):
|
||||
"""Dense implementation of mutual reachability graph.
|
||||
|
||||
The computation is done in-place, i.e. the distance matrix is modified
|
||||
directly.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
distance_matrix : ndarray of shape (n_samples, n_samples)
|
||||
Array of distances between samples.
|
||||
|
||||
further_neighbor_idx : int
|
||||
The index of the furthest neighbor to use to define the core distances.
|
||||
"""
|
||||
cdef:
|
||||
intp_t i, j, n_samples = distance_matrix.shape[0]
|
||||
floating mutual_reachability_distance
|
||||
floating[::1] core_distances
|
||||
|
||||
# We assume that the distance matrix is symmetric. We choose to sort every
|
||||
# row to have the same implementation than the sparse case that requires
|
||||
# CSR matrix.
|
||||
core_distances = np.ascontiguousarray(
|
||||
np.partition(
|
||||
distance_matrix, further_neighbor_idx, axis=1
|
||||
)[:, further_neighbor_idx]
|
||||
)
|
||||
|
||||
with nogil:
|
||||
# TODO: Update w/ prange with thread count based on
|
||||
# _openmp_effective_n_threads
|
||||
for i in range(n_samples):
|
||||
for j in range(n_samples):
|
||||
mutual_reachability_distance = max(
|
||||
core_distances[i],
|
||||
core_distances[j],
|
||||
distance_matrix[i, j],
|
||||
)
|
||||
distance_matrix[i, j] = mutual_reachability_distance
|
||||
|
||||
|
||||
def _sparse_mutual_reachability_graph(
|
||||
cnp.ndarray[floating, ndim=1, mode="c"] data,
|
||||
cnp.ndarray[integral, ndim=1, mode="c"] indices,
|
||||
cnp.ndarray[integral, ndim=1, mode="c"] indptr,
|
||||
intp_t n_samples,
|
||||
intp_t further_neighbor_idx,
|
||||
floating max_distance,
|
||||
):
|
||||
"""Sparse implementation of mutual reachability graph.
|
||||
|
||||
The computation is done in-place, i.e. the distance matrix is modified
|
||||
directly. This implementation only accepts `CSR` format sparse matrices.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
distance_matrix : sparse matrix of shape (n_samples, n_samples)
|
||||
Sparse matrix of distances between samples. The sparse format should
|
||||
be `CSR`.
|
||||
|
||||
further_neighbor_idx : int
|
||||
The index of the furthest neighbor to use to define the core distances.
|
||||
|
||||
max_distance : float
|
||||
The distance which `np.inf` is replaced with. When the true mutual-
|
||||
reachability distance is measured to be infinite, it is instead
|
||||
truncated to `max_dist`. Only used when `distance_matrix` is a sparse
|
||||
matrix.
|
||||
"""
|
||||
cdef:
|
||||
integral i, col_ind, row_ind
|
||||
floating mutual_reachability_distance
|
||||
floating[:] core_distances
|
||||
floating[:] row_data
|
||||
|
||||
if floating is float:
|
||||
dtype = np.float32
|
||||
else:
|
||||
dtype = np.float64
|
||||
|
||||
core_distances = np.empty(n_samples, dtype=dtype)
|
||||
|
||||
for i in range(n_samples):
|
||||
row_data = data[indptr[i]:indptr[i + 1]]
|
||||
if further_neighbor_idx < row_data.size:
|
||||
core_distances[i] = np.partition(
|
||||
row_data, further_neighbor_idx
|
||||
)[further_neighbor_idx]
|
||||
else:
|
||||
core_distances[i] = INFINITY
|
||||
|
||||
with nogil:
|
||||
for row_ind in range(n_samples):
|
||||
for i in range(indptr[row_ind], indptr[row_ind + 1]):
|
||||
col_ind = indices[i]
|
||||
mutual_reachability_distance = max(
|
||||
core_distances[row_ind], core_distances[col_ind], data[i]
|
||||
)
|
||||
if isfinite(mutual_reachability_distance):
|
||||
data[i] = mutual_reachability_distance
|
||||
elif max_distance > 0:
|
||||
data[i] = max_distance
|
||||
# mutual reachability distance computations
|
||||
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
# Redistribution and use in source and binary forms, with or without
|
||||
# modification, are permitted provided that the following conditions are met:
|
||||
|
||||
# 1. Redistributions of source code must retain the above copyright notice,
|
||||
# this list of conditions and the following disclaimer.
|
||||
|
||||
# 2. Redistributions in binary form must reproduce the above copyright notice,
|
||||
# this list of conditions and the following disclaimer in the documentation
|
||||
# and/or other materials provided with the distribution.
|
||||
|
||||
# 3. Neither the name of the copyright holder nor the names of its contributors
|
||||
# may be used to endorse or promote products derived from this software without
|
||||
# specific prior written permission.
|
||||
|
||||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
# POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
cimport numpy as cnp
|
||||
|
||||
import numpy as np
|
||||
from scipy.sparse import issparse
|
||||
from cython cimport floating, integral
|
||||
from libc.math cimport isfinite, INFINITY
|
||||
from ...utils._typedefs cimport intp_t
|
||||
cnp.import_array()
|
||||
|
||||
|
||||
def mutual_reachability_graph(
|
||||
distance_matrix, min_samples=5, max_distance=0.0
|
||||
):
|
||||
"""Compute the weighted adjacency matrix of the mutual reachability graph.
|
||||
|
||||
The mutual reachability distance used to build the graph is defined as::
|
||||
|
||||
max(d_core(x_p), d_core(x_q), d(x_p, x_q))
|
||||
|
||||
and the core distance `d_core` is defined as the distance between a point
|
||||
`x_p` and its k-th nearest neighbor.
|
||||
|
||||
Note that all computations are done in-place.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
distance_matrix : {ndarray, sparse matrix} of shape (n_samples, n_samples)
|
||||
Array of distances between samples. If sparse, the array must be in
|
||||
`CSR` format.
|
||||
|
||||
min_samples : int, default=5
|
||||
The parameter `k` used to calculate the distance between a point
|
||||
`x_p` and its k-th nearest neighbor.
|
||||
|
||||
max_distance : float, default=0.0
|
||||
The distance which `np.inf` is replaced with. When the true mutual-
|
||||
reachability distance is measured to be infinite, it is instead
|
||||
truncated to `max_dist`. Only used when `distance_matrix` is a sparse
|
||||
matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mututal_reachability_graph: {ndarray, sparse matrix} of shape \
|
||||
(n_samples, n_samples)
|
||||
Weighted adjacency matrix of the mutual reachability graph.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Campello, R. J., Moulavi, D., & Sander, J. (2013, April).
|
||||
Density-based clustering based on hierarchical density estimates.
|
||||
In Pacific-Asia Conference on Knowledge Discovery and Data Mining
|
||||
(pp. 160-172). Springer Berlin Heidelberg.
|
||||
"""
|
||||
further_neighbor_idx = min_samples - 1
|
||||
if issparse(distance_matrix):
|
||||
if distance_matrix.format != "csr":
|
||||
raise ValueError(
|
||||
"Only sparse CSR matrices are supported for `distance_matrix`."
|
||||
)
|
||||
_sparse_mutual_reachability_graph(
|
||||
distance_matrix.data,
|
||||
distance_matrix.indices,
|
||||
distance_matrix.indptr,
|
||||
distance_matrix.shape[0],
|
||||
further_neighbor_idx=further_neighbor_idx,
|
||||
max_distance=max_distance,
|
||||
)
|
||||
else:
|
||||
_dense_mutual_reachability_graph(
|
||||
distance_matrix, further_neighbor_idx=further_neighbor_idx
|
||||
)
|
||||
return distance_matrix
|
||||
|
||||
|
||||
def _dense_mutual_reachability_graph(
|
||||
floating[:, :] distance_matrix,
|
||||
intp_t further_neighbor_idx,
|
||||
):
|
||||
"""Dense implementation of mutual reachability graph.
|
||||
|
||||
The computation is done in-place, i.e. the distance matrix is modified
|
||||
directly.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
distance_matrix : ndarray of shape (n_samples, n_samples)
|
||||
Array of distances between samples.
|
||||
|
||||
further_neighbor_idx : int
|
||||
The index of the furthest neighbor to use to define the core distances.
|
||||
"""
|
||||
cdef:
|
||||
intp_t i, j, n_samples = distance_matrix.shape[0]
|
||||
floating mutual_reachability_distance
|
||||
floating[::1] core_distances
|
||||
|
||||
# We assume that the distance matrix is symmetric. We choose to sort every
|
||||
# row to have the same implementation than the sparse case that requires
|
||||
# CSR matrix.
|
||||
core_distances = np.ascontiguousarray(
|
||||
np.partition(
|
||||
distance_matrix, further_neighbor_idx, axis=1
|
||||
)[:, further_neighbor_idx]
|
||||
)
|
||||
|
||||
with nogil:
|
||||
# TODO: Update w/ prange with thread count based on
|
||||
# _openmp_effective_n_threads
|
||||
for i in range(n_samples):
|
||||
for j in range(n_samples):
|
||||
mutual_reachability_distance = max(
|
||||
core_distances[i],
|
||||
core_distances[j],
|
||||
distance_matrix[i, j],
|
||||
)
|
||||
distance_matrix[i, j] = mutual_reachability_distance
|
||||
|
||||
|
||||
def _sparse_mutual_reachability_graph(
|
||||
cnp.ndarray[floating, ndim=1, mode="c"] data,
|
||||
cnp.ndarray[integral, ndim=1, mode="c"] indices,
|
||||
cnp.ndarray[integral, ndim=1, mode="c"] indptr,
|
||||
intp_t n_samples,
|
||||
intp_t further_neighbor_idx,
|
||||
floating max_distance,
|
||||
):
|
||||
"""Sparse implementation of mutual reachability graph.
|
||||
|
||||
The computation is done in-place, i.e. the distance matrix is modified
|
||||
directly. This implementation only accepts `CSR` format sparse matrices.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
distance_matrix : sparse matrix of shape (n_samples, n_samples)
|
||||
Sparse matrix of distances between samples. The sparse format should
|
||||
be `CSR`.
|
||||
|
||||
further_neighbor_idx : int
|
||||
The index of the furthest neighbor to use to define the core distances.
|
||||
|
||||
max_distance : float
|
||||
The distance which `np.inf` is replaced with. When the true mutual-
|
||||
reachability distance is measured to be infinite, it is instead
|
||||
truncated to `max_dist`. Only used when `distance_matrix` is a sparse
|
||||
matrix.
|
||||
"""
|
||||
cdef:
|
||||
integral i, col_ind, row_ind
|
||||
floating mutual_reachability_distance
|
||||
floating[:] core_distances
|
||||
floating[:] row_data
|
||||
|
||||
if floating is float:
|
||||
dtype = np.float32
|
||||
else:
|
||||
dtype = np.float64
|
||||
|
||||
core_distances = np.empty(n_samples, dtype=dtype)
|
||||
|
||||
for i in range(n_samples):
|
||||
row_data = data[indptr[i]:indptr[i + 1]]
|
||||
if further_neighbor_idx < row_data.size:
|
||||
core_distances[i] = np.partition(
|
||||
row_data, further_neighbor_idx
|
||||
)[further_neighbor_idx]
|
||||
else:
|
||||
core_distances[i] = INFINITY
|
||||
|
||||
with nogil:
|
||||
for row_ind in range(n_samples):
|
||||
for i in range(indptr[row_ind], indptr[row_ind + 1]):
|
||||
col_ind = indices[i]
|
||||
mutual_reachability_distance = max(
|
||||
core_distances[row_ind], core_distances[col_ind], data[i]
|
||||
)
|
||||
if isfinite(mutual_reachability_distance):
|
||||
data[i] = mutual_reachability_distance
|
||||
elif max_distance > 0:
|
||||
data[i] = max_distance
|
||||
|
||||
@@ -1,49 +1,49 @@
|
||||
# Copyright (c) 2015, Leland McInnes
|
||||
# All rights reserved.
|
||||
|
||||
# Redistribution and use in source and binary forms, with or without
|
||||
# modification, are permitted provided that the following conditions are met:
|
||||
|
||||
# 1. Redistributions of source code must retain the above copyright notice,
|
||||
# this list of conditions and the following disclaimer.
|
||||
|
||||
# 2. Redistributions in binary form must reproduce the above copyright notice,
|
||||
# this list of conditions and the following disclaimer in the documentation
|
||||
# and/or other materials provided with the distribution.
|
||||
|
||||
# 3. Neither the name of the copyright holder nor the names of its contributors
|
||||
# may be used to endorse or promote products derived from this software without
|
||||
# specific prior written permission.
|
||||
|
||||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
# POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
from ...utils._typedefs cimport intp_t, float64_t, uint8_t
|
||||
cimport numpy as cnp
|
||||
|
||||
# This corresponds to the scipy.cluster.hierarchy format
|
||||
ctypedef packed struct HIERARCHY_t:
|
||||
intp_t left_node
|
||||
intp_t right_node
|
||||
float64_t value
|
||||
intp_t cluster_size
|
||||
|
||||
# Effectively an edgelist encoding a parent/child pair, along with a value and
|
||||
# the corresponding cluster_size in each row providing a tree structure.
|
||||
ctypedef packed struct CONDENSED_t:
|
||||
intp_t parent
|
||||
intp_t child
|
||||
float64_t value
|
||||
intp_t cluster_size
|
||||
|
||||
cdef extern from "numpy/arrayobject.h":
|
||||
intp_t * PyArray_SHAPE(cnp.PyArrayObject *)
|
||||
# Copyright (c) 2015, Leland McInnes
|
||||
# All rights reserved.
|
||||
|
||||
# Redistribution and use in source and binary forms, with or without
|
||||
# modification, are permitted provided that the following conditions are met:
|
||||
|
||||
# 1. Redistributions of source code must retain the above copyright notice,
|
||||
# this list of conditions and the following disclaimer.
|
||||
|
||||
# 2. Redistributions in binary form must reproduce the above copyright notice,
|
||||
# this list of conditions and the following disclaimer in the documentation
|
||||
# and/or other materials provided with the distribution.
|
||||
|
||||
# 3. Neither the name of the copyright holder nor the names of its contributors
|
||||
# may be used to endorse or promote products derived from this software without
|
||||
# specific prior written permission.
|
||||
|
||||
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
||||
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
|
||||
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
||||
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
||||
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
||||
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
||||
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
||||
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
||||
# POSSIBILITY OF SUCH DAMAGE.
|
||||
|
||||
from ...utils._typedefs cimport intp_t, float64_t, uint8_t
|
||||
cimport numpy as cnp
|
||||
|
||||
# This corresponds to the scipy.cluster.hierarchy format
|
||||
ctypedef packed struct HIERARCHY_t:
|
||||
intp_t left_node
|
||||
intp_t right_node
|
||||
float64_t value
|
||||
intp_t cluster_size
|
||||
|
||||
# Effectively an edgelist encoding a parent/child pair, along with a value and
|
||||
# the corresponding cluster_size in each row providing a tree structure.
|
||||
ctypedef packed struct CONDENSED_t:
|
||||
intp_t parent
|
||||
intp_t child
|
||||
float64_t value
|
||||
intp_t cluster_size
|
||||
|
||||
cdef extern from "numpy/arrayobject.h":
|
||||
intp_t * PyArray_SHAPE(cnp.PyArrayObject *)
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,15 +1,15 @@
|
||||
cluster_hdbscan_extension_metadata = {
|
||||
'_linkage': {'sources': [cython_gen.process('_linkage.pyx'), metrics_cython_tree]},
|
||||
'_reachability': {'sources': [cython_gen.process('_reachability.pyx')]},
|
||||
'_tree': {'sources': [cython_gen.process('_tree.pyx')]}
|
||||
}
|
||||
|
||||
foreach ext_name, ext_dict : cluster_hdbscan_extension_metadata
|
||||
py.extension_module(
|
||||
ext_name,
|
||||
ext_dict.get('sources'),
|
||||
dependencies: [np_dep],
|
||||
subdir: 'sklearn/cluster/_hdbscan',
|
||||
install: true
|
||||
)
|
||||
endforeach
|
||||
cluster_hdbscan_extension_metadata = {
|
||||
'_linkage': {'sources': [cython_gen.process('_linkage.pyx'), metrics_cython_tree]},
|
||||
'_reachability': {'sources': [cython_gen.process('_reachability.pyx')]},
|
||||
'_tree': {'sources': [cython_gen.process('_tree.pyx')]}
|
||||
}
|
||||
|
||||
foreach ext_name, ext_dict : cluster_hdbscan_extension_metadata
|
||||
py.extension_module(
|
||||
ext_name,
|
||||
ext_dict.get('sources'),
|
||||
dependencies: [np_dep],
|
||||
subdir: 'sklearn/cluster/_hdbscan',
|
||||
install: true
|
||||
)
|
||||
endforeach
|
||||
|
||||
@@ -1,9 +1,9 @@
|
||||
from ..utils._typedefs cimport intp_t
|
||||
|
||||
cdef class UnionFind:
|
||||
cdef intp_t next_label
|
||||
cdef intp_t[:] parent
|
||||
cdef intp_t[:] size
|
||||
|
||||
cdef void union(self, intp_t m, intp_t n) noexcept
|
||||
cdef intp_t fast_find(self, intp_t n) noexcept
|
||||
from ..utils._typedefs cimport intp_t
|
||||
|
||||
cdef class UnionFind:
|
||||
cdef intp_t next_label
|
||||
cdef intp_t[:] parent
|
||||
cdef intp_t[:] size
|
||||
|
||||
cdef void union(self, intp_t m, intp_t n) noexcept
|
||||
cdef intp_t fast_find(self, intp_t n) noexcept
|
||||
|
||||
+507
-507
File diff suppressed because it is too large
Load Diff
@@ -1,48 +1,48 @@
|
||||
from cython cimport floating
|
||||
|
||||
|
||||
cdef floating _euclidean_dense_dense(
|
||||
const floating*,
|
||||
const floating*,
|
||||
int,
|
||||
bint
|
||||
) noexcept nogil
|
||||
|
||||
cdef floating _euclidean_sparse_dense(
|
||||
const floating[::1],
|
||||
const int[::1],
|
||||
const floating[::1],
|
||||
floating,
|
||||
bint
|
||||
) noexcept nogil
|
||||
|
||||
cpdef void _relocate_empty_clusters_dense(
|
||||
const floating[:, ::1],
|
||||
const floating[::1],
|
||||
const floating[:, ::1],
|
||||
floating[:, ::1],
|
||||
floating[::1],
|
||||
const int[::1]
|
||||
)
|
||||
|
||||
cpdef void _relocate_empty_clusters_sparse(
|
||||
const floating[::1],
|
||||
const int[::1],
|
||||
const int[::1],
|
||||
const floating[::1],
|
||||
const floating[:, ::1],
|
||||
floating[:, ::1],
|
||||
floating[::1],
|
||||
const int[::1]
|
||||
)
|
||||
|
||||
cdef void _average_centers(
|
||||
floating[:, ::1],
|
||||
const floating[::1]
|
||||
)
|
||||
|
||||
cdef void _center_shift(
|
||||
const floating[:, ::1],
|
||||
const floating[:, ::1],
|
||||
floating[::1]
|
||||
)
|
||||
from cython cimport floating
|
||||
|
||||
|
||||
cdef floating _euclidean_dense_dense(
|
||||
const floating*,
|
||||
const floating*,
|
||||
int,
|
||||
bint
|
||||
) noexcept nogil
|
||||
|
||||
cdef floating _euclidean_sparse_dense(
|
||||
const floating[::1],
|
||||
const int[::1],
|
||||
const floating[::1],
|
||||
floating,
|
||||
bint
|
||||
) noexcept nogil
|
||||
|
||||
cpdef void _relocate_empty_clusters_dense(
|
||||
const floating[:, ::1],
|
||||
const floating[::1],
|
||||
const floating[:, ::1],
|
||||
floating[:, ::1],
|
||||
floating[::1],
|
||||
const int[::1]
|
||||
)
|
||||
|
||||
cpdef void _relocate_empty_clusters_sparse(
|
||||
const floating[::1],
|
||||
const int[::1],
|
||||
const int[::1],
|
||||
const floating[::1],
|
||||
const floating[:, ::1],
|
||||
floating[:, ::1],
|
||||
floating[::1],
|
||||
const int[::1]
|
||||
)
|
||||
|
||||
cdef void _average_centers(
|
||||
floating[:, ::1],
|
||||
const floating[::1]
|
||||
)
|
||||
|
||||
cdef void _center_shift(
|
||||
const floating[:, ::1],
|
||||
const floating[:, ::1],
|
||||
floating[::1]
|
||||
)
|
||||
|
||||
@@ -1,328 +1,328 @@
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
import numpy as np
|
||||
from cython cimport floating
|
||||
from cython.parallel cimport prange
|
||||
from libc.math cimport sqrt
|
||||
|
||||
from ..utils.extmath import row_norms
|
||||
|
||||
|
||||
# Number of samples per data chunk defined as a global constant.
|
||||
CHUNK_SIZE = 256
|
||||
|
||||
|
||||
cdef floating _euclidean_dense_dense(
|
||||
const floating* a, # IN
|
||||
const floating* b, # IN
|
||||
int n_features,
|
||||
bint squared
|
||||
) noexcept nogil:
|
||||
"""Euclidean distance between a dense and b dense"""
|
||||
cdef:
|
||||
int i
|
||||
int n = n_features // 4
|
||||
int rem = n_features % 4
|
||||
floating result = 0
|
||||
|
||||
# We manually unroll the loop for better cache optimization.
|
||||
for i in range(n):
|
||||
result += (
|
||||
(a[0] - b[0]) * (a[0] - b[0]) +
|
||||
(a[1] - b[1]) * (a[1] - b[1]) +
|
||||
(a[2] - b[2]) * (a[2] - b[2]) +
|
||||
(a[3] - b[3]) * (a[3] - b[3])
|
||||
)
|
||||
a += 4
|
||||
b += 4
|
||||
|
||||
for i in range(rem):
|
||||
result += (a[i] - b[i]) * (a[i] - b[i])
|
||||
|
||||
return result if squared else sqrt(result)
|
||||
|
||||
|
||||
def _euclidean_dense_dense_wrapper(
|
||||
const floating[::1] a,
|
||||
const floating[::1] b,
|
||||
bint squared
|
||||
):
|
||||
"""Wrapper of _euclidean_dense_dense for testing purpose"""
|
||||
return _euclidean_dense_dense(&a[0], &b[0], a.shape[0], squared)
|
||||
|
||||
|
||||
cdef floating _euclidean_sparse_dense(
|
||||
const floating[::1] a_data, # IN
|
||||
const int[::1] a_indices, # IN
|
||||
const floating[::1] b, # IN
|
||||
floating b_squared_norm,
|
||||
bint squared
|
||||
) noexcept nogil:
|
||||
"""Euclidean distance between a sparse and b dense"""
|
||||
cdef:
|
||||
int nnz = a_indices.shape[0]
|
||||
int i
|
||||
floating tmp, bi
|
||||
floating result = 0.0
|
||||
|
||||
for i in range(nnz):
|
||||
bi = b[a_indices[i]]
|
||||
tmp = a_data[i] - bi
|
||||
result += tmp * tmp - bi * bi
|
||||
|
||||
result += b_squared_norm
|
||||
|
||||
if result < 0:
|
||||
result = 0.0
|
||||
|
||||
return result if squared else sqrt(result)
|
||||
|
||||
|
||||
def _euclidean_sparse_dense_wrapper(
|
||||
const floating[::1] a_data,
|
||||
const int[::1] a_indices,
|
||||
const floating[::1] b,
|
||||
floating b_squared_norm,
|
||||
bint squared
|
||||
):
|
||||
"""Wrapper of _euclidean_sparse_dense for testing purpose"""
|
||||
return _euclidean_sparse_dense(
|
||||
a_data, a_indices, b, b_squared_norm, squared)
|
||||
|
||||
|
||||
cpdef floating _inertia_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers, # IN
|
||||
const int[::1] labels, # IN
|
||||
int n_threads,
|
||||
int single_label=-1,
|
||||
):
|
||||
"""Compute inertia for dense input data
|
||||
|
||||
Sum of squared distance between each sample and its assigned center.
|
||||
|
||||
If single_label is >= 0, the inertia is computed only for that label.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_features = X.shape[1]
|
||||
int i, j
|
||||
|
||||
floating sq_dist = 0.0
|
||||
floating inertia = 0.0
|
||||
|
||||
for i in prange(n_samples, nogil=True, num_threads=n_threads,
|
||||
schedule='static'):
|
||||
j = labels[i]
|
||||
if single_label < 0 or single_label == j:
|
||||
sq_dist = _euclidean_dense_dense(&X[i, 0], ¢ers[j, 0],
|
||||
n_features, True)
|
||||
inertia += sq_dist * sample_weight[i]
|
||||
|
||||
return inertia
|
||||
|
||||
|
||||
cpdef floating _inertia_sparse(
|
||||
X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers, # IN
|
||||
const int[::1] labels, # IN
|
||||
int n_threads,
|
||||
int single_label=-1,
|
||||
):
|
||||
"""Compute inertia for sparse input data
|
||||
|
||||
Sum of squared distance between each sample and its assigned center.
|
||||
|
||||
If single_label is >= 0, the inertia is computed only for that label.
|
||||
"""
|
||||
cdef:
|
||||
floating[::1] X_data = X.data
|
||||
int[::1] X_indices = X.indices
|
||||
int[::1] X_indptr = X.indptr
|
||||
|
||||
int n_samples = X.shape[0]
|
||||
int i, j
|
||||
|
||||
floating sq_dist = 0.0
|
||||
floating inertia = 0.0
|
||||
|
||||
floating[::1] centers_squared_norms = row_norms(centers, squared=True)
|
||||
|
||||
for i in prange(n_samples, nogil=True, num_threads=n_threads,
|
||||
schedule='static'):
|
||||
j = labels[i]
|
||||
if single_label < 0 or single_label == j:
|
||||
sq_dist = _euclidean_sparse_dense(
|
||||
X_data[X_indptr[i]: X_indptr[i + 1]],
|
||||
X_indices[X_indptr[i]: X_indptr[i + 1]],
|
||||
centers[j], centers_squared_norms[j], True)
|
||||
inertia += sq_dist * sample_weight[i]
|
||||
|
||||
return inertia
|
||||
|
||||
|
||||
cpdef void _relocate_empty_clusters_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # INOUT
|
||||
floating[::1] weight_in_clusters, # INOUT
|
||||
const int[::1] labels # IN
|
||||
):
|
||||
"""Relocate centers which have no sample assigned to them."""
|
||||
cdef:
|
||||
int[::1] empty_clusters = np.where(np.equal(weight_in_clusters, 0))[0].astype(np.int32)
|
||||
int n_empty = empty_clusters.shape[0]
|
||||
|
||||
if n_empty == 0:
|
||||
return
|
||||
|
||||
cdef:
|
||||
int n_features = X.shape[1]
|
||||
|
||||
floating[::1] distances = ((np.asarray(X) - np.asarray(centers_old)[labels])**2).sum(axis=1)
|
||||
int[::1] far_from_centers = np.argpartition(distances, -n_empty)[:-n_empty-1:-1].astype(np.int32)
|
||||
|
||||
int new_cluster_id, old_cluster_id, far_idx, idx, k
|
||||
floating weight
|
||||
|
||||
if np.max(distances) == 0:
|
||||
# Happens when there are more clusters than non-duplicate samples. Relocating
|
||||
# is pointless in this case.
|
||||
return
|
||||
|
||||
for idx in range(n_empty):
|
||||
|
||||
new_cluster_id = empty_clusters[idx]
|
||||
|
||||
far_idx = far_from_centers[idx]
|
||||
weight = sample_weight[far_idx]
|
||||
|
||||
old_cluster_id = labels[far_idx]
|
||||
|
||||
for k in range(n_features):
|
||||
centers_new[old_cluster_id, k] -= X[far_idx, k] * weight
|
||||
centers_new[new_cluster_id, k] = X[far_idx, k] * weight
|
||||
|
||||
weight_in_clusters[new_cluster_id] = weight
|
||||
weight_in_clusters[old_cluster_id] -= weight
|
||||
|
||||
|
||||
cpdef void _relocate_empty_clusters_sparse(
|
||||
const floating[::1] X_data, # IN
|
||||
const int[::1] X_indices, # IN
|
||||
const int[::1] X_indptr, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # INOUT
|
||||
floating[::1] weight_in_clusters, # INOUT
|
||||
const int[::1] labels # IN
|
||||
):
|
||||
"""Relocate centers which have no sample assigned to them."""
|
||||
cdef:
|
||||
int[::1] empty_clusters = np.where(np.equal(weight_in_clusters, 0))[0].astype(np.int32)
|
||||
int n_empty = empty_clusters.shape[0]
|
||||
|
||||
if n_empty == 0:
|
||||
return
|
||||
|
||||
cdef:
|
||||
int n_samples = X_indptr.shape[0] - 1
|
||||
int i, j, k
|
||||
|
||||
floating[::1] distances = np.zeros(n_samples, dtype=X_data.base.dtype)
|
||||
floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)
|
||||
|
||||
for i in range(n_samples):
|
||||
j = labels[i]
|
||||
distances[i] = _euclidean_sparse_dense(
|
||||
X_data[X_indptr[i]: X_indptr[i + 1]],
|
||||
X_indices[X_indptr[i]: X_indptr[i + 1]],
|
||||
centers_old[j], centers_squared_norms[j], True)
|
||||
|
||||
if np.max(distances) == 0:
|
||||
# Happens when there are more clusters than non-duplicate samples. Relocating
|
||||
# is pointless in this case.
|
||||
return
|
||||
|
||||
cdef:
|
||||
int[::1] far_from_centers = np.argpartition(distances, -n_empty)[:-n_empty-1:-1].astype(np.int32)
|
||||
|
||||
int new_cluster_id, old_cluster_id, far_idx, idx
|
||||
floating weight
|
||||
|
||||
for idx in range(n_empty):
|
||||
|
||||
new_cluster_id = empty_clusters[idx]
|
||||
|
||||
far_idx = far_from_centers[idx]
|
||||
weight = sample_weight[far_idx]
|
||||
|
||||
old_cluster_id = labels[far_idx]
|
||||
|
||||
for k in range(X_indptr[far_idx], X_indptr[far_idx + 1]):
|
||||
centers_new[old_cluster_id, X_indices[k]] -= X_data[k] * weight
|
||||
centers_new[new_cluster_id, X_indices[k]] = X_data[k] * weight
|
||||
|
||||
weight_in_clusters[new_cluster_id] = weight
|
||||
weight_in_clusters[old_cluster_id] -= weight
|
||||
|
||||
|
||||
cdef void _average_centers(
|
||||
floating[:, ::1] centers, # INOUT
|
||||
const floating[::1] weight_in_clusters # IN
|
||||
):
|
||||
"""Average new centers wrt weights."""
|
||||
cdef:
|
||||
int n_clusters = centers.shape[0]
|
||||
int n_features = centers.shape[1]
|
||||
int j, k
|
||||
floating alpha
|
||||
int argmax_weight = np.argmax(weight_in_clusters)
|
||||
|
||||
for j in range(n_clusters):
|
||||
if weight_in_clusters[j] > 0:
|
||||
alpha = 1.0 / weight_in_clusters[j]
|
||||
for k in range(n_features):
|
||||
centers[j, k] *= alpha
|
||||
else:
|
||||
# For convenience, we avoid setting empty clusters at the origin but place
|
||||
# them at the location of the biggest cluster.
|
||||
for k in range(n_features):
|
||||
centers[j, k] = centers[argmax_weight, k]
|
||||
|
||||
|
||||
cdef void _center_shift(
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
const floating[:, ::1] centers_new, # IN
|
||||
floating[::1] center_shift # OUT
|
||||
):
|
||||
"""Compute shift between old and new centers."""
|
||||
cdef:
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
int j
|
||||
|
||||
for j in range(n_clusters):
|
||||
center_shift[j] = _euclidean_dense_dense(
|
||||
¢ers_new[j, 0], ¢ers_old[j, 0], n_features, False)
|
||||
|
||||
|
||||
def _is_same_clustering(
|
||||
const int[::1] labels1,
|
||||
const int[::1] labels2,
|
||||
n_clusters
|
||||
):
|
||||
"""Check if two arrays of labels are the same up to a permutation of the labels"""
|
||||
cdef int[::1] mapping = np.full(fill_value=-1, shape=(n_clusters,), dtype=np.int32)
|
||||
cdef int i
|
||||
|
||||
for i in range(labels1.shape[0]):
|
||||
if mapping[labels1[i]] == -1:
|
||||
mapping[labels1[i]] = labels2[i]
|
||||
elif mapping[labels1[i]] != labels2[i]:
|
||||
return False
|
||||
return True
|
||||
# Authors: The scikit-learn developers
|
||||
# SPDX-License-Identifier: BSD-3-Clause
|
||||
|
||||
import numpy as np
|
||||
from cython cimport floating
|
||||
from cython.parallel cimport prange
|
||||
from libc.math cimport sqrt
|
||||
|
||||
from ..utils.extmath import row_norms
|
||||
|
||||
|
||||
# Number of samples per data chunk defined as a global constant.
|
||||
CHUNK_SIZE = 256
|
||||
|
||||
|
||||
cdef floating _euclidean_dense_dense(
|
||||
const floating* a, # IN
|
||||
const floating* b, # IN
|
||||
int n_features,
|
||||
bint squared
|
||||
) noexcept nogil:
|
||||
"""Euclidean distance between a dense and b dense"""
|
||||
cdef:
|
||||
int i
|
||||
int n = n_features // 4
|
||||
int rem = n_features % 4
|
||||
floating result = 0
|
||||
|
||||
# We manually unroll the loop for better cache optimization.
|
||||
for i in range(n):
|
||||
result += (
|
||||
(a[0] - b[0]) * (a[0] - b[0]) +
|
||||
(a[1] - b[1]) * (a[1] - b[1]) +
|
||||
(a[2] - b[2]) * (a[2] - b[2]) +
|
||||
(a[3] - b[3]) * (a[3] - b[3])
|
||||
)
|
||||
a += 4
|
||||
b += 4
|
||||
|
||||
for i in range(rem):
|
||||
result += (a[i] - b[i]) * (a[i] - b[i])
|
||||
|
||||
return result if squared else sqrt(result)
|
||||
|
||||
|
||||
def _euclidean_dense_dense_wrapper(
|
||||
const floating[::1] a,
|
||||
const floating[::1] b,
|
||||
bint squared
|
||||
):
|
||||
"""Wrapper of _euclidean_dense_dense for testing purpose"""
|
||||
return _euclidean_dense_dense(&a[0], &b[0], a.shape[0], squared)
|
||||
|
||||
|
||||
cdef floating _euclidean_sparse_dense(
|
||||
const floating[::1] a_data, # IN
|
||||
const int[::1] a_indices, # IN
|
||||
const floating[::1] b, # IN
|
||||
floating b_squared_norm,
|
||||
bint squared
|
||||
) noexcept nogil:
|
||||
"""Euclidean distance between a sparse and b dense"""
|
||||
cdef:
|
||||
int nnz = a_indices.shape[0]
|
||||
int i
|
||||
floating tmp, bi
|
||||
floating result = 0.0
|
||||
|
||||
for i in range(nnz):
|
||||
bi = b[a_indices[i]]
|
||||
tmp = a_data[i] - bi
|
||||
result += tmp * tmp - bi * bi
|
||||
|
||||
result += b_squared_norm
|
||||
|
||||
if result < 0:
|
||||
result = 0.0
|
||||
|
||||
return result if squared else sqrt(result)
|
||||
|
||||
|
||||
def _euclidean_sparse_dense_wrapper(
|
||||
const floating[::1] a_data,
|
||||
const int[::1] a_indices,
|
||||
const floating[::1] b,
|
||||
floating b_squared_norm,
|
||||
bint squared
|
||||
):
|
||||
"""Wrapper of _euclidean_sparse_dense for testing purpose"""
|
||||
return _euclidean_sparse_dense(
|
||||
a_data, a_indices, b, b_squared_norm, squared)
|
||||
|
||||
|
||||
cpdef floating _inertia_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers, # IN
|
||||
const int[::1] labels, # IN
|
||||
int n_threads,
|
||||
int single_label=-1,
|
||||
):
|
||||
"""Compute inertia for dense input data
|
||||
|
||||
Sum of squared distance between each sample and its assigned center.
|
||||
|
||||
If single_label is >= 0, the inertia is computed only for that label.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_features = X.shape[1]
|
||||
int i, j
|
||||
|
||||
floating sq_dist = 0.0
|
||||
floating inertia = 0.0
|
||||
|
||||
for i in prange(n_samples, nogil=True, num_threads=n_threads,
|
||||
schedule='static'):
|
||||
j = labels[i]
|
||||
if single_label < 0 or single_label == j:
|
||||
sq_dist = _euclidean_dense_dense(&X[i, 0], ¢ers[j, 0],
|
||||
n_features, True)
|
||||
inertia += sq_dist * sample_weight[i]
|
||||
|
||||
return inertia
|
||||
|
||||
|
||||
cpdef floating _inertia_sparse(
|
||||
X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers, # IN
|
||||
const int[::1] labels, # IN
|
||||
int n_threads,
|
||||
int single_label=-1,
|
||||
):
|
||||
"""Compute inertia for sparse input data
|
||||
|
||||
Sum of squared distance between each sample and its assigned center.
|
||||
|
||||
If single_label is >= 0, the inertia is computed only for that label.
|
||||
"""
|
||||
cdef:
|
||||
floating[::1] X_data = X.data
|
||||
int[::1] X_indices = X.indices
|
||||
int[::1] X_indptr = X.indptr
|
||||
|
||||
int n_samples = X.shape[0]
|
||||
int i, j
|
||||
|
||||
floating sq_dist = 0.0
|
||||
floating inertia = 0.0
|
||||
|
||||
floating[::1] centers_squared_norms = row_norms(centers, squared=True)
|
||||
|
||||
for i in prange(n_samples, nogil=True, num_threads=n_threads,
|
||||
schedule='static'):
|
||||
j = labels[i]
|
||||
if single_label < 0 or single_label == j:
|
||||
sq_dist = _euclidean_sparse_dense(
|
||||
X_data[X_indptr[i]: X_indptr[i + 1]],
|
||||
X_indices[X_indptr[i]: X_indptr[i + 1]],
|
||||
centers[j], centers_squared_norms[j], True)
|
||||
inertia += sq_dist * sample_weight[i]
|
||||
|
||||
return inertia
|
||||
|
||||
|
||||
cpdef void _relocate_empty_clusters_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # INOUT
|
||||
floating[::1] weight_in_clusters, # INOUT
|
||||
const int[::1] labels # IN
|
||||
):
|
||||
"""Relocate centers which have no sample assigned to them."""
|
||||
cdef:
|
||||
int[::1] empty_clusters = np.where(np.equal(weight_in_clusters, 0))[0].astype(np.int32)
|
||||
int n_empty = empty_clusters.shape[0]
|
||||
|
||||
if n_empty == 0:
|
||||
return
|
||||
|
||||
cdef:
|
||||
int n_features = X.shape[1]
|
||||
|
||||
floating[::1] distances = ((np.asarray(X) - np.asarray(centers_old)[labels])**2).sum(axis=1)
|
||||
int[::1] far_from_centers = np.argpartition(distances, -n_empty)[:-n_empty-1:-1].astype(np.int32)
|
||||
|
||||
int new_cluster_id, old_cluster_id, far_idx, idx, k
|
||||
floating weight
|
||||
|
||||
if np.max(distances) == 0:
|
||||
# Happens when there are more clusters than non-duplicate samples. Relocating
|
||||
# is pointless in this case.
|
||||
return
|
||||
|
||||
for idx in range(n_empty):
|
||||
|
||||
new_cluster_id = empty_clusters[idx]
|
||||
|
||||
far_idx = far_from_centers[idx]
|
||||
weight = sample_weight[far_idx]
|
||||
|
||||
old_cluster_id = labels[far_idx]
|
||||
|
||||
for k in range(n_features):
|
||||
centers_new[old_cluster_id, k] -= X[far_idx, k] * weight
|
||||
centers_new[new_cluster_id, k] = X[far_idx, k] * weight
|
||||
|
||||
weight_in_clusters[new_cluster_id] = weight
|
||||
weight_in_clusters[old_cluster_id] -= weight
|
||||
|
||||
|
||||
cpdef void _relocate_empty_clusters_sparse(
|
||||
const floating[::1] X_data, # IN
|
||||
const int[::1] X_indices, # IN
|
||||
const int[::1] X_indptr, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # INOUT
|
||||
floating[::1] weight_in_clusters, # INOUT
|
||||
const int[::1] labels # IN
|
||||
):
|
||||
"""Relocate centers which have no sample assigned to them."""
|
||||
cdef:
|
||||
int[::1] empty_clusters = np.where(np.equal(weight_in_clusters, 0))[0].astype(np.int32)
|
||||
int n_empty = empty_clusters.shape[0]
|
||||
|
||||
if n_empty == 0:
|
||||
return
|
||||
|
||||
cdef:
|
||||
int n_samples = X_indptr.shape[0] - 1
|
||||
int i, j, k
|
||||
|
||||
floating[::1] distances = np.zeros(n_samples, dtype=X_data.base.dtype)
|
||||
floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)
|
||||
|
||||
for i in range(n_samples):
|
||||
j = labels[i]
|
||||
distances[i] = _euclidean_sparse_dense(
|
||||
X_data[X_indptr[i]: X_indptr[i + 1]],
|
||||
X_indices[X_indptr[i]: X_indptr[i + 1]],
|
||||
centers_old[j], centers_squared_norms[j], True)
|
||||
|
||||
if np.max(distances) == 0:
|
||||
# Happens when there are more clusters than non-duplicate samples. Relocating
|
||||
# is pointless in this case.
|
||||
return
|
||||
|
||||
cdef:
|
||||
int[::1] far_from_centers = np.argpartition(distances, -n_empty)[:-n_empty-1:-1].astype(np.int32)
|
||||
|
||||
int new_cluster_id, old_cluster_id, far_idx, idx
|
||||
floating weight
|
||||
|
||||
for idx in range(n_empty):
|
||||
|
||||
new_cluster_id = empty_clusters[idx]
|
||||
|
||||
far_idx = far_from_centers[idx]
|
||||
weight = sample_weight[far_idx]
|
||||
|
||||
old_cluster_id = labels[far_idx]
|
||||
|
||||
for k in range(X_indptr[far_idx], X_indptr[far_idx + 1]):
|
||||
centers_new[old_cluster_id, X_indices[k]] -= X_data[k] * weight
|
||||
centers_new[new_cluster_id, X_indices[k]] = X_data[k] * weight
|
||||
|
||||
weight_in_clusters[new_cluster_id] = weight
|
||||
weight_in_clusters[old_cluster_id] -= weight
|
||||
|
||||
|
||||
cdef void _average_centers(
|
||||
floating[:, ::1] centers, # INOUT
|
||||
const floating[::1] weight_in_clusters # IN
|
||||
):
|
||||
"""Average new centers wrt weights."""
|
||||
cdef:
|
||||
int n_clusters = centers.shape[0]
|
||||
int n_features = centers.shape[1]
|
||||
int j, k
|
||||
floating alpha
|
||||
int argmax_weight = np.argmax(weight_in_clusters)
|
||||
|
||||
for j in range(n_clusters):
|
||||
if weight_in_clusters[j] > 0:
|
||||
alpha = 1.0 / weight_in_clusters[j]
|
||||
for k in range(n_features):
|
||||
centers[j, k] *= alpha
|
||||
else:
|
||||
# For convenience, we avoid setting empty clusters at the origin but place
|
||||
# them at the location of the biggest cluster.
|
||||
for k in range(n_features):
|
||||
centers[j, k] = centers[argmax_weight, k]
|
||||
|
||||
|
||||
cdef void _center_shift(
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
const floating[:, ::1] centers_new, # IN
|
||||
floating[::1] center_shift # OUT
|
||||
):
|
||||
"""Compute shift between old and new centers."""
|
||||
cdef:
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
int j
|
||||
|
||||
for j in range(n_clusters):
|
||||
center_shift[j] = _euclidean_dense_dense(
|
||||
¢ers_new[j, 0], ¢ers_old[j, 0], n_features, False)
|
||||
|
||||
|
||||
def _is_same_clustering(
|
||||
const int[::1] labels1,
|
||||
const int[::1] labels2,
|
||||
n_clusters
|
||||
):
|
||||
"""Check if two arrays of labels are the same up to a permutation of the labels"""
|
||||
cdef int[::1] mapping = np.full(fill_value=-1, shape=(n_clusters,), dtype=np.int32)
|
||||
cdef int i
|
||||
|
||||
for i in range(labels1.shape[0]):
|
||||
if mapping[labels1[i]] == -1:
|
||||
mapping[labels1[i]] = labels2[i]
|
||||
elif mapping[labels1[i]] != labels2[i]:
|
||||
return False
|
||||
return True
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,420 +1,420 @@
|
||||
# Licence: BSD 3 clause
|
||||
|
||||
from cython cimport floating
|
||||
from cython.parallel import prange, parallel
|
||||
from libc.stdlib cimport malloc, calloc, free
|
||||
from libc.string cimport memset
|
||||
from libc.float cimport DBL_MAX, FLT_MAX
|
||||
|
||||
from ..utils._openmp_helpers cimport omp_lock_t
|
||||
from ..utils._openmp_helpers cimport omp_init_lock
|
||||
from ..utils._openmp_helpers cimport omp_destroy_lock
|
||||
from ..utils._openmp_helpers cimport omp_set_lock
|
||||
from ..utils._openmp_helpers cimport omp_unset_lock
|
||||
from ..utils.extmath import row_norms
|
||||
from ..utils._cython_blas cimport _gemm
|
||||
from ..utils._cython_blas cimport RowMajor, Trans, NoTrans
|
||||
from ._k_means_common import CHUNK_SIZE
|
||||
from ._k_means_common cimport _relocate_empty_clusters_dense
|
||||
from ._k_means_common cimport _relocate_empty_clusters_sparse
|
||||
from ._k_means_common cimport _average_centers, _center_shift
|
||||
|
||||
|
||||
def lloyd_iter_chunked_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_in_clusters, # OUT
|
||||
int[::1] labels, # OUT
|
||||
floating[::1] center_shift, # OUT
|
||||
int n_threads,
|
||||
bint update_centers=True):
|
||||
"""Single iteration of K-means lloyd algorithm with dense input.
|
||||
|
||||
Update labels and centers (inplace), for one iteration, distributed
|
||||
over data chunks.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : ndarray of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration. `centers_new` can be `None` if
|
||||
`update_centers` is False.
|
||||
|
||||
weight_in_clusters : ndarray of shape (n_clusters,), dtype=floating
|
||||
Placeholder for the sums of the weights of every observation assigned
|
||||
to each center. `weight_in_clusters` can be `None` if `update_centers`
|
||||
is False.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
center_shift : ndarray of shape (n_clusters,), dtype=floating
|
||||
Distance between old and new centers.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
|
||||
update_centers : bool
|
||||
- If True, the labels and the new centers will be computed, i.e. runs
|
||||
the E-step and the M-step of the algorithm.
|
||||
- If False, only the labels will be computed, i.e runs the E-step of
|
||||
the algorithm. This is useful especially when calling predict on a
|
||||
fitted model.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_features = X.shape[1]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
|
||||
if n_samples == 0:
|
||||
# An empty array was passed, do nothing and return early (before
|
||||
# attempting to compute n_chunks). This can typically happen when
|
||||
# calling the prediction function of a bisecting k-means model with a
|
||||
# large fraction of outliers.
|
||||
return
|
||||
|
||||
cdef:
|
||||
# hard-coded number of samples per chunk. Appeared to be close to
|
||||
# optimal in all situations.
|
||||
int n_samples_chunk = CHUNK_SIZE if n_samples > CHUNK_SIZE else n_samples
|
||||
int n_chunks = n_samples // n_samples_chunk
|
||||
int n_samples_rem = n_samples % n_samples_chunk
|
||||
int chunk_idx
|
||||
int start, end
|
||||
|
||||
int j, k
|
||||
|
||||
floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)
|
||||
|
||||
floating *centers_new_chunk
|
||||
floating *weight_in_clusters_chunk
|
||||
floating *pairwise_distances_chunk
|
||||
|
||||
omp_lock_t lock
|
||||
|
||||
# count remainder chunk in total number of chunks
|
||||
n_chunks += n_samples != n_chunks * n_samples_chunk
|
||||
|
||||
# number of threads should not be bigger than number of chunks
|
||||
n_threads = min(n_threads, n_chunks)
|
||||
|
||||
if update_centers:
|
||||
memset(¢ers_new[0, 0], 0, n_clusters * n_features * sizeof(floating))
|
||||
memset(&weight_in_clusters[0], 0, n_clusters * sizeof(floating))
|
||||
omp_init_lock(&lock)
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
# thread local buffers
|
||||
centers_new_chunk = <floating*> calloc(n_clusters * n_features, sizeof(floating))
|
||||
weight_in_clusters_chunk = <floating*> calloc(n_clusters, sizeof(floating))
|
||||
pairwise_distances_chunk = <floating*> malloc(n_samples_chunk * n_clusters * sizeof(floating))
|
||||
|
||||
for chunk_idx in prange(n_chunks, schedule='static'):
|
||||
start = chunk_idx * n_samples_chunk
|
||||
if chunk_idx == n_chunks - 1 and n_samples_rem > 0:
|
||||
end = start + n_samples_rem
|
||||
else:
|
||||
end = start + n_samples_chunk
|
||||
|
||||
_update_chunk_dense(
|
||||
X[start: end],
|
||||
sample_weight[start: end],
|
||||
centers_old,
|
||||
centers_squared_norms,
|
||||
labels[start: end],
|
||||
centers_new_chunk,
|
||||
weight_in_clusters_chunk,
|
||||
pairwise_distances_chunk,
|
||||
update_centers)
|
||||
|
||||
# reduction from local buffers.
|
||||
if update_centers:
|
||||
# The lock is necessary to avoid race conditions when aggregating
|
||||
# info from different thread-local buffers.
|
||||
omp_set_lock(&lock)
|
||||
for j in range(n_clusters):
|
||||
weight_in_clusters[j] += weight_in_clusters_chunk[j]
|
||||
for k in range(n_features):
|
||||
centers_new[j, k] += centers_new_chunk[j * n_features + k]
|
||||
|
||||
omp_unset_lock(&lock)
|
||||
|
||||
free(centers_new_chunk)
|
||||
free(weight_in_clusters_chunk)
|
||||
free(pairwise_distances_chunk)
|
||||
|
||||
if update_centers:
|
||||
omp_destroy_lock(&lock)
|
||||
_relocate_empty_clusters_dense(
|
||||
X, sample_weight, centers_old, centers_new, weight_in_clusters, labels
|
||||
)
|
||||
|
||||
_average_centers(centers_new, weight_in_clusters)
|
||||
_center_shift(centers_old, centers_new, center_shift)
|
||||
|
||||
|
||||
cdef void _update_chunk_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
const floating[::1] centers_squared_norms, # IN
|
||||
int[::1] labels, # OUT
|
||||
floating *centers_new, # OUT
|
||||
floating *weight_in_clusters, # OUT
|
||||
floating *pairwise_distances, # OUT
|
||||
bint update_centers) noexcept nogil:
|
||||
"""K-means combined EM step for one dense data chunk.
|
||||
|
||||
Compute the partial contribution of a single data chunk to the labels and
|
||||
centers.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = labels.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
|
||||
floating sq_dist, min_sq_dist
|
||||
int i, j, k, label
|
||||
|
||||
# Instead of computing the full pairwise squared distances matrix,
|
||||
# ||X - C||² = ||X||² - 2 X.C^T + ||C||², we only need to store
|
||||
# the - 2 X.C^T + ||C||² term since the argmin for a given sample only
|
||||
# depends on the centers.
|
||||
# pairwise_distances = ||C||²
|
||||
for i in range(n_samples):
|
||||
for j in range(n_clusters):
|
||||
pairwise_distances[i * n_clusters + j] = centers_squared_norms[j]
|
||||
|
||||
# pairwise_distances += -2 * X.dot(C.T)
|
||||
_gemm(RowMajor, NoTrans, Trans, n_samples, n_clusters, n_features,
|
||||
-2.0, &X[0, 0], n_features, ¢ers_old[0, 0], n_features,
|
||||
1.0, pairwise_distances, n_clusters)
|
||||
|
||||
for i in range(n_samples):
|
||||
min_sq_dist = pairwise_distances[i * n_clusters]
|
||||
label = 0
|
||||
for j in range(1, n_clusters):
|
||||
sq_dist = pairwise_distances[i * n_clusters + j]
|
||||
if sq_dist < min_sq_dist:
|
||||
min_sq_dist = sq_dist
|
||||
label = j
|
||||
labels[i] = label
|
||||
|
||||
if update_centers:
|
||||
weight_in_clusters[label] += sample_weight[i]
|
||||
for k in range(n_features):
|
||||
centers_new[label * n_features + k] += X[i, k] * sample_weight[i]
|
||||
|
||||
|
||||
def lloyd_iter_chunked_sparse(
|
||||
X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_in_clusters, # OUT
|
||||
int[::1] labels, # OUT
|
||||
floating[::1] center_shift, # OUT
|
||||
int n_threads,
|
||||
bint update_centers=True):
|
||||
"""Single iteration of K-means lloyd algorithm with sparse input.
|
||||
|
||||
Update labels and centers (inplace), for one iteration, distributed
|
||||
over data chunks.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : sparse matrix of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster. Must be in CSR format.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration. `centers_new` can be `None` if
|
||||
`update_centers` is False.
|
||||
|
||||
weight_in_clusters : ndarray of shape (n_clusters,), dtype=floating
|
||||
Placeholder for the sums of the weights of every observation assigned
|
||||
to each center. `weight_in_clusters` can be `None` if `update_centers`
|
||||
is False.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
center_shift : ndarray of shape (n_clusters,), dtype=floating
|
||||
Distance between old and new centers.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
|
||||
update_centers : bool
|
||||
- If True, the labels and the new centers will be computed, i.e. runs
|
||||
the E-step and the M-step of the algorithm.
|
||||
- If False, only the labels will be computed, i.e runs the E-step of
|
||||
the algorithm. This is useful especially when calling predict on a
|
||||
fitted model.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_features = X.shape[1]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
|
||||
if n_samples == 0:
|
||||
# An empty array was passed, do nothing and return early (before
|
||||
# attempting to compute n_chunks). This can typically happen when
|
||||
# calling the prediction function of a bisecting k-means model with a
|
||||
# large fraction of outliers.
|
||||
return
|
||||
|
||||
cdef:
|
||||
# Choose same as for dense. Does not have the same impact since with
|
||||
# sparse data the pairwise distances matrix is not precomputed.
|
||||
# However, splitting in chunks is necessary to get parallelism.
|
||||
int n_samples_chunk = CHUNK_SIZE if n_samples > CHUNK_SIZE else n_samples
|
||||
int n_chunks = n_samples // n_samples_chunk
|
||||
int n_samples_rem = n_samples % n_samples_chunk
|
||||
int chunk_idx
|
||||
int start = 0, end = 0
|
||||
|
||||
int j, k
|
||||
|
||||
floating[::1] X_data = X.data
|
||||
int[::1] X_indices = X.indices
|
||||
int[::1] X_indptr = X.indptr
|
||||
|
||||
floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)
|
||||
|
||||
floating *centers_new_chunk
|
||||
floating *weight_in_clusters_chunk
|
||||
|
||||
omp_lock_t lock
|
||||
|
||||
# count remainder chunk in total number of chunks
|
||||
n_chunks += n_samples != n_chunks * n_samples_chunk
|
||||
|
||||
# number of threads should not be bigger than number of chunks
|
||||
n_threads = min(n_threads, n_chunks)
|
||||
|
||||
if update_centers:
|
||||
memset(¢ers_new[0, 0], 0, n_clusters * n_features * sizeof(floating))
|
||||
memset(&weight_in_clusters[0], 0, n_clusters * sizeof(floating))
|
||||
omp_init_lock(&lock)
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
# thread local buffers
|
||||
centers_new_chunk = <floating*> calloc(n_clusters * n_features, sizeof(floating))
|
||||
weight_in_clusters_chunk = <floating*> calloc(n_clusters, sizeof(floating))
|
||||
|
||||
for chunk_idx in prange(n_chunks, schedule='static'):
|
||||
start = chunk_idx * n_samples_chunk
|
||||
if chunk_idx == n_chunks - 1 and n_samples_rem > 0:
|
||||
end = start + n_samples_rem
|
||||
else:
|
||||
end = start + n_samples_chunk
|
||||
|
||||
_update_chunk_sparse(
|
||||
X_data[X_indptr[start]: X_indptr[end]],
|
||||
X_indices[X_indptr[start]: X_indptr[end]],
|
||||
X_indptr[start: end+1],
|
||||
sample_weight[start: end],
|
||||
centers_old,
|
||||
centers_squared_norms,
|
||||
labels[start: end],
|
||||
centers_new_chunk,
|
||||
weight_in_clusters_chunk,
|
||||
update_centers)
|
||||
|
||||
# reduction from local buffers.
|
||||
if update_centers:
|
||||
# The lock is necessary to avoid race conditions when aggregating
|
||||
# info from different thread-local buffers.
|
||||
omp_set_lock(&lock)
|
||||
for j in range(n_clusters):
|
||||
weight_in_clusters[j] += weight_in_clusters_chunk[j]
|
||||
for k in range(n_features):
|
||||
centers_new[j, k] += centers_new_chunk[j * n_features + k]
|
||||
omp_unset_lock(&lock)
|
||||
|
||||
free(centers_new_chunk)
|
||||
free(weight_in_clusters_chunk)
|
||||
|
||||
if update_centers:
|
||||
omp_destroy_lock(&lock)
|
||||
_relocate_empty_clusters_sparse(
|
||||
X_data, X_indices, X_indptr, sample_weight,
|
||||
centers_old, centers_new, weight_in_clusters, labels)
|
||||
|
||||
_average_centers(centers_new, weight_in_clusters)
|
||||
_center_shift(centers_old, centers_new, center_shift)
|
||||
|
||||
|
||||
cdef void _update_chunk_sparse(
|
||||
const floating[::1] X_data, # IN
|
||||
const int[::1] X_indices, # IN
|
||||
const int[::1] X_indptr, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
const floating[::1] centers_squared_norms, # IN
|
||||
int[::1] labels, # OUT
|
||||
floating *centers_new, # OUT
|
||||
floating *weight_in_clusters, # OUT
|
||||
bint update_centers) noexcept nogil:
|
||||
"""K-means combined EM step for one sparse data chunk.
|
||||
|
||||
Compute the partial contribution of a single data chunk to the labels and
|
||||
centers.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = labels.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
|
||||
floating sq_dist, min_sq_dist
|
||||
int i, j, k, label
|
||||
floating max_floating = FLT_MAX if floating is float else DBL_MAX
|
||||
int s = X_indptr[0]
|
||||
|
||||
# XXX Precompute the pairwise distances matrix is not worth for sparse
|
||||
# currently. Should be tested when BLAS (sparse x dense) matrix
|
||||
# multiplication is available.
|
||||
for i in range(n_samples):
|
||||
min_sq_dist = max_floating
|
||||
label = 0
|
||||
|
||||
for j in range(n_clusters):
|
||||
sq_dist = 0.0
|
||||
for k in range(X_indptr[i] - s, X_indptr[i + 1] - s):
|
||||
sq_dist += centers_old[j, X_indices[k]] * X_data[k]
|
||||
|
||||
# Instead of computing the full squared distance with each cluster,
|
||||
# ||X - C||² = ||X||² - 2 X.C^T + ||C||², we only need to compute
|
||||
# the - 2 X.C^T + ||C||² term since the argmin for a given sample
|
||||
# only depends on the centers C.
|
||||
sq_dist = centers_squared_norms[j] -2 * sq_dist
|
||||
if sq_dist < min_sq_dist:
|
||||
min_sq_dist = sq_dist
|
||||
label = j
|
||||
|
||||
labels[i] = label
|
||||
|
||||
if update_centers:
|
||||
weight_in_clusters[label] += sample_weight[i]
|
||||
for k in range(X_indptr[i] - s, X_indptr[i + 1] - s):
|
||||
centers_new[label * n_features + X_indices[k]] += X_data[k] * sample_weight[i]
|
||||
# Licence: BSD 3 clause
|
||||
|
||||
from cython cimport floating
|
||||
from cython.parallel import prange, parallel
|
||||
from libc.stdlib cimport malloc, calloc, free
|
||||
from libc.string cimport memset
|
||||
from libc.float cimport DBL_MAX, FLT_MAX
|
||||
|
||||
from ..utils._openmp_helpers cimport omp_lock_t
|
||||
from ..utils._openmp_helpers cimport omp_init_lock
|
||||
from ..utils._openmp_helpers cimport omp_destroy_lock
|
||||
from ..utils._openmp_helpers cimport omp_set_lock
|
||||
from ..utils._openmp_helpers cimport omp_unset_lock
|
||||
from ..utils.extmath import row_norms
|
||||
from ..utils._cython_blas cimport _gemm
|
||||
from ..utils._cython_blas cimport RowMajor, Trans, NoTrans
|
||||
from ._k_means_common import CHUNK_SIZE
|
||||
from ._k_means_common cimport _relocate_empty_clusters_dense
|
||||
from ._k_means_common cimport _relocate_empty_clusters_sparse
|
||||
from ._k_means_common cimport _average_centers, _center_shift
|
||||
|
||||
|
||||
def lloyd_iter_chunked_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_in_clusters, # OUT
|
||||
int[::1] labels, # OUT
|
||||
floating[::1] center_shift, # OUT
|
||||
int n_threads,
|
||||
bint update_centers=True):
|
||||
"""Single iteration of K-means lloyd algorithm with dense input.
|
||||
|
||||
Update labels and centers (inplace), for one iteration, distributed
|
||||
over data chunks.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : ndarray of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration. `centers_new` can be `None` if
|
||||
`update_centers` is False.
|
||||
|
||||
weight_in_clusters : ndarray of shape (n_clusters,), dtype=floating
|
||||
Placeholder for the sums of the weights of every observation assigned
|
||||
to each center. `weight_in_clusters` can be `None` if `update_centers`
|
||||
is False.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
center_shift : ndarray of shape (n_clusters,), dtype=floating
|
||||
Distance between old and new centers.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
|
||||
update_centers : bool
|
||||
- If True, the labels and the new centers will be computed, i.e. runs
|
||||
the E-step and the M-step of the algorithm.
|
||||
- If False, only the labels will be computed, i.e runs the E-step of
|
||||
the algorithm. This is useful especially when calling predict on a
|
||||
fitted model.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_features = X.shape[1]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
|
||||
if n_samples == 0:
|
||||
# An empty array was passed, do nothing and return early (before
|
||||
# attempting to compute n_chunks). This can typically happen when
|
||||
# calling the prediction function of a bisecting k-means model with a
|
||||
# large fraction of outliers.
|
||||
return
|
||||
|
||||
cdef:
|
||||
# hard-coded number of samples per chunk. Appeared to be close to
|
||||
# optimal in all situations.
|
||||
int n_samples_chunk = CHUNK_SIZE if n_samples > CHUNK_SIZE else n_samples
|
||||
int n_chunks = n_samples // n_samples_chunk
|
||||
int n_samples_rem = n_samples % n_samples_chunk
|
||||
int chunk_idx
|
||||
int start, end
|
||||
|
||||
int j, k
|
||||
|
||||
floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)
|
||||
|
||||
floating *centers_new_chunk
|
||||
floating *weight_in_clusters_chunk
|
||||
floating *pairwise_distances_chunk
|
||||
|
||||
omp_lock_t lock
|
||||
|
||||
# count remainder chunk in total number of chunks
|
||||
n_chunks += n_samples != n_chunks * n_samples_chunk
|
||||
|
||||
# number of threads should not be bigger than number of chunks
|
||||
n_threads = min(n_threads, n_chunks)
|
||||
|
||||
if update_centers:
|
||||
memset(¢ers_new[0, 0], 0, n_clusters * n_features * sizeof(floating))
|
||||
memset(&weight_in_clusters[0], 0, n_clusters * sizeof(floating))
|
||||
omp_init_lock(&lock)
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
# thread local buffers
|
||||
centers_new_chunk = <floating*> calloc(n_clusters * n_features, sizeof(floating))
|
||||
weight_in_clusters_chunk = <floating*> calloc(n_clusters, sizeof(floating))
|
||||
pairwise_distances_chunk = <floating*> malloc(n_samples_chunk * n_clusters * sizeof(floating))
|
||||
|
||||
for chunk_idx in prange(n_chunks, schedule='static'):
|
||||
start = chunk_idx * n_samples_chunk
|
||||
if chunk_idx == n_chunks - 1 and n_samples_rem > 0:
|
||||
end = start + n_samples_rem
|
||||
else:
|
||||
end = start + n_samples_chunk
|
||||
|
||||
_update_chunk_dense(
|
||||
X[start: end],
|
||||
sample_weight[start: end],
|
||||
centers_old,
|
||||
centers_squared_norms,
|
||||
labels[start: end],
|
||||
centers_new_chunk,
|
||||
weight_in_clusters_chunk,
|
||||
pairwise_distances_chunk,
|
||||
update_centers)
|
||||
|
||||
# reduction from local buffers.
|
||||
if update_centers:
|
||||
# The lock is necessary to avoid race conditions when aggregating
|
||||
# info from different thread-local buffers.
|
||||
omp_set_lock(&lock)
|
||||
for j in range(n_clusters):
|
||||
weight_in_clusters[j] += weight_in_clusters_chunk[j]
|
||||
for k in range(n_features):
|
||||
centers_new[j, k] += centers_new_chunk[j * n_features + k]
|
||||
|
||||
omp_unset_lock(&lock)
|
||||
|
||||
free(centers_new_chunk)
|
||||
free(weight_in_clusters_chunk)
|
||||
free(pairwise_distances_chunk)
|
||||
|
||||
if update_centers:
|
||||
omp_destroy_lock(&lock)
|
||||
_relocate_empty_clusters_dense(
|
||||
X, sample_weight, centers_old, centers_new, weight_in_clusters, labels
|
||||
)
|
||||
|
||||
_average_centers(centers_new, weight_in_clusters)
|
||||
_center_shift(centers_old, centers_new, center_shift)
|
||||
|
||||
|
||||
cdef void _update_chunk_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
const floating[::1] centers_squared_norms, # IN
|
||||
int[::1] labels, # OUT
|
||||
floating *centers_new, # OUT
|
||||
floating *weight_in_clusters, # OUT
|
||||
floating *pairwise_distances, # OUT
|
||||
bint update_centers) noexcept nogil:
|
||||
"""K-means combined EM step for one dense data chunk.
|
||||
|
||||
Compute the partial contribution of a single data chunk to the labels and
|
||||
centers.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = labels.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
|
||||
floating sq_dist, min_sq_dist
|
||||
int i, j, k, label
|
||||
|
||||
# Instead of computing the full pairwise squared distances matrix,
|
||||
# ||X - C||² = ||X||² - 2 X.C^T + ||C||², we only need to store
|
||||
# the - 2 X.C^T + ||C||² term since the argmin for a given sample only
|
||||
# depends on the centers.
|
||||
# pairwise_distances = ||C||²
|
||||
for i in range(n_samples):
|
||||
for j in range(n_clusters):
|
||||
pairwise_distances[i * n_clusters + j] = centers_squared_norms[j]
|
||||
|
||||
# pairwise_distances += -2 * X.dot(C.T)
|
||||
_gemm(RowMajor, NoTrans, Trans, n_samples, n_clusters, n_features,
|
||||
-2.0, &X[0, 0], n_features, ¢ers_old[0, 0], n_features,
|
||||
1.0, pairwise_distances, n_clusters)
|
||||
|
||||
for i in range(n_samples):
|
||||
min_sq_dist = pairwise_distances[i * n_clusters]
|
||||
label = 0
|
||||
for j in range(1, n_clusters):
|
||||
sq_dist = pairwise_distances[i * n_clusters + j]
|
||||
if sq_dist < min_sq_dist:
|
||||
min_sq_dist = sq_dist
|
||||
label = j
|
||||
labels[i] = label
|
||||
|
||||
if update_centers:
|
||||
weight_in_clusters[label] += sample_weight[i]
|
||||
for k in range(n_features):
|
||||
centers_new[label * n_features + k] += X[i, k] * sample_weight[i]
|
||||
|
||||
|
||||
def lloyd_iter_chunked_sparse(
|
||||
X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_in_clusters, # OUT
|
||||
int[::1] labels, # OUT
|
||||
floating[::1] center_shift, # OUT
|
||||
int n_threads,
|
||||
bint update_centers=True):
|
||||
"""Single iteration of K-means lloyd algorithm with sparse input.
|
||||
|
||||
Update labels and centers (inplace), for one iteration, distributed
|
||||
over data chunks.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : sparse matrix of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster. Must be in CSR format.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration. `centers_new` can be `None` if
|
||||
`update_centers` is False.
|
||||
|
||||
weight_in_clusters : ndarray of shape (n_clusters,), dtype=floating
|
||||
Placeholder for the sums of the weights of every observation assigned
|
||||
to each center. `weight_in_clusters` can be `None` if `update_centers`
|
||||
is False.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
center_shift : ndarray of shape (n_clusters,), dtype=floating
|
||||
Distance between old and new centers.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
|
||||
update_centers : bool
|
||||
- If True, the labels and the new centers will be computed, i.e. runs
|
||||
the E-step and the M-step of the algorithm.
|
||||
- If False, only the labels will be computed, i.e runs the E-step of
|
||||
the algorithm. This is useful especially when calling predict on a
|
||||
fitted model.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_features = X.shape[1]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
|
||||
if n_samples == 0:
|
||||
# An empty array was passed, do nothing and return early (before
|
||||
# attempting to compute n_chunks). This can typically happen when
|
||||
# calling the prediction function of a bisecting k-means model with a
|
||||
# large fraction of outliers.
|
||||
return
|
||||
|
||||
cdef:
|
||||
# Choose same as for dense. Does not have the same impact since with
|
||||
# sparse data the pairwise distances matrix is not precomputed.
|
||||
# However, splitting in chunks is necessary to get parallelism.
|
||||
int n_samples_chunk = CHUNK_SIZE if n_samples > CHUNK_SIZE else n_samples
|
||||
int n_chunks = n_samples // n_samples_chunk
|
||||
int n_samples_rem = n_samples % n_samples_chunk
|
||||
int chunk_idx
|
||||
int start = 0, end = 0
|
||||
|
||||
int j, k
|
||||
|
||||
floating[::1] X_data = X.data
|
||||
int[::1] X_indices = X.indices
|
||||
int[::1] X_indptr = X.indptr
|
||||
|
||||
floating[::1] centers_squared_norms = row_norms(centers_old, squared=True)
|
||||
|
||||
floating *centers_new_chunk
|
||||
floating *weight_in_clusters_chunk
|
||||
|
||||
omp_lock_t lock
|
||||
|
||||
# count remainder chunk in total number of chunks
|
||||
n_chunks += n_samples != n_chunks * n_samples_chunk
|
||||
|
||||
# number of threads should not be bigger than number of chunks
|
||||
n_threads = min(n_threads, n_chunks)
|
||||
|
||||
if update_centers:
|
||||
memset(¢ers_new[0, 0], 0, n_clusters * n_features * sizeof(floating))
|
||||
memset(&weight_in_clusters[0], 0, n_clusters * sizeof(floating))
|
||||
omp_init_lock(&lock)
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
# thread local buffers
|
||||
centers_new_chunk = <floating*> calloc(n_clusters * n_features, sizeof(floating))
|
||||
weight_in_clusters_chunk = <floating*> calloc(n_clusters, sizeof(floating))
|
||||
|
||||
for chunk_idx in prange(n_chunks, schedule='static'):
|
||||
start = chunk_idx * n_samples_chunk
|
||||
if chunk_idx == n_chunks - 1 and n_samples_rem > 0:
|
||||
end = start + n_samples_rem
|
||||
else:
|
||||
end = start + n_samples_chunk
|
||||
|
||||
_update_chunk_sparse(
|
||||
X_data[X_indptr[start]: X_indptr[end]],
|
||||
X_indices[X_indptr[start]: X_indptr[end]],
|
||||
X_indptr[start: end+1],
|
||||
sample_weight[start: end],
|
||||
centers_old,
|
||||
centers_squared_norms,
|
||||
labels[start: end],
|
||||
centers_new_chunk,
|
||||
weight_in_clusters_chunk,
|
||||
update_centers)
|
||||
|
||||
# reduction from local buffers.
|
||||
if update_centers:
|
||||
# The lock is necessary to avoid race conditions when aggregating
|
||||
# info from different thread-local buffers.
|
||||
omp_set_lock(&lock)
|
||||
for j in range(n_clusters):
|
||||
weight_in_clusters[j] += weight_in_clusters_chunk[j]
|
||||
for k in range(n_features):
|
||||
centers_new[j, k] += centers_new_chunk[j * n_features + k]
|
||||
omp_unset_lock(&lock)
|
||||
|
||||
free(centers_new_chunk)
|
||||
free(weight_in_clusters_chunk)
|
||||
|
||||
if update_centers:
|
||||
omp_destroy_lock(&lock)
|
||||
_relocate_empty_clusters_sparse(
|
||||
X_data, X_indices, X_indptr, sample_weight,
|
||||
centers_old, centers_new, weight_in_clusters, labels)
|
||||
|
||||
_average_centers(centers_new, weight_in_clusters)
|
||||
_center_shift(centers_old, centers_new, center_shift)
|
||||
|
||||
|
||||
cdef void _update_chunk_sparse(
|
||||
const floating[::1] X_data, # IN
|
||||
const int[::1] X_indices, # IN
|
||||
const int[::1] X_indptr, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
const floating[::1] centers_squared_norms, # IN
|
||||
int[::1] labels, # OUT
|
||||
floating *centers_new, # OUT
|
||||
floating *weight_in_clusters, # OUT
|
||||
bint update_centers) noexcept nogil:
|
||||
"""K-means combined EM step for one sparse data chunk.
|
||||
|
||||
Compute the partial contribution of a single data chunk to the labels and
|
||||
centers.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = labels.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
|
||||
floating sq_dist, min_sq_dist
|
||||
int i, j, k, label
|
||||
floating max_floating = FLT_MAX if floating is float else DBL_MAX
|
||||
int s = X_indptr[0]
|
||||
|
||||
# XXX Precompute the pairwise distances matrix is not worth for sparse
|
||||
# currently. Should be tested when BLAS (sparse x dense) matrix
|
||||
# multiplication is available.
|
||||
for i in range(n_samples):
|
||||
min_sq_dist = max_floating
|
||||
label = 0
|
||||
|
||||
for j in range(n_clusters):
|
||||
sq_dist = 0.0
|
||||
for k in range(X_indptr[i] - s, X_indptr[i + 1] - s):
|
||||
sq_dist += centers_old[j, X_indices[k]] * X_data[k]
|
||||
|
||||
# Instead of computing the full squared distance with each cluster,
|
||||
# ||X - C||² = ||X||² - 2 X.C^T + ||C||², we only need to compute
|
||||
# the - 2 X.C^T + ||C||² term since the argmin for a given sample
|
||||
# only depends on the centers C.
|
||||
sq_dist = centers_squared_norms[j] -2 * sq_dist
|
||||
if sq_dist < min_sq_dist:
|
||||
min_sq_dist = sq_dist
|
||||
label = j
|
||||
|
||||
labels[i] = label
|
||||
|
||||
if update_centers:
|
||||
weight_in_clusters[label] += sample_weight[i]
|
||||
for k in range(X_indptr[i] - s, X_indptr[i + 1] - s):
|
||||
centers_new[label * n_features + X_indices[k]] += X_data[k] * sample_weight[i]
|
||||
|
||||
+218
-218
@@ -1,218 +1,218 @@
|
||||
from cython cimport floating
|
||||
from cython.parallel cimport parallel, prange
|
||||
from libc.stdlib cimport malloc, free
|
||||
|
||||
|
||||
def _minibatch_update_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int n_threads):
|
||||
"""Update of the centers for dense MiniBatchKMeans.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : ndarray of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration.
|
||||
|
||||
weight_sums : ndarray of shape (n_clusters,), dtype=floating
|
||||
Current sums of the accumulated weights for each center.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int cluster_idx
|
||||
|
||||
int *indices
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
indices = <int*> malloc(n_samples * sizeof(int))
|
||||
|
||||
for cluster_idx in prange(n_clusters, schedule="static"):
|
||||
update_center_dense(cluster_idx, X, sample_weight,
|
||||
centers_old, centers_new, weight_sums, labels,
|
||||
indices)
|
||||
|
||||
free(indices)
|
||||
|
||||
|
||||
cdef void update_center_dense(
|
||||
int cluster_idx,
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int *indices) noexcept nogil: # TMP
|
||||
"""Update of a single center for dense MinibatchKMeans"""
|
||||
cdef:
|
||||
int n_samples = sample_weight.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
floating alpha
|
||||
int n_indices
|
||||
int k, sample_idx, feature_idx
|
||||
|
||||
floating wsum = 0
|
||||
|
||||
# indices = np.where(labels == cluster_idx)[0]
|
||||
k = 0
|
||||
for sample_idx in range(n_samples):
|
||||
if labels[sample_idx] == cluster_idx:
|
||||
indices[k] = sample_idx
|
||||
wsum += sample_weight[sample_idx]
|
||||
k += 1
|
||||
n_indices = k
|
||||
|
||||
if wsum > 0:
|
||||
# Undo the previous count-based scaling for this cluster center
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx] * weight_sums[cluster_idx]
|
||||
|
||||
# Update cluster with new point members
|
||||
for k in range(n_indices):
|
||||
sample_idx = indices[k]
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] += X[sample_idx, feature_idx] * sample_weight[sample_idx]
|
||||
|
||||
# Update the count statistics for this center
|
||||
weight_sums[cluster_idx] += wsum
|
||||
|
||||
# Rescale to compute mean of all points (old and new)
|
||||
alpha = 1 / weight_sums[cluster_idx]
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] *= alpha
|
||||
else:
|
||||
# No sample was assigned to this cluster in this batch of data
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx]
|
||||
|
||||
|
||||
def _minibatch_update_sparse(
|
||||
X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int n_threads):
|
||||
"""Update of the centers for sparse MiniBatchKMeans.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : sparse matrix of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster. Must be in CSR format.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration.
|
||||
|
||||
weight_sums : ndarray of shape (n_clusters,), dtype=floating
|
||||
Current sums of the accumulated weights for each center.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
"""
|
||||
cdef:
|
||||
floating[::1] X_data = X.data
|
||||
int[::1] X_indices = X.indices
|
||||
int[::1] X_indptr = X.indptr
|
||||
int n_samples = X.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int cluster_idx
|
||||
|
||||
int *indices
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
indices = <int*> malloc(n_samples * sizeof(int))
|
||||
|
||||
for cluster_idx in prange(n_clusters, schedule="static"):
|
||||
update_center_sparse(cluster_idx, X_data, X_indices, X_indptr,
|
||||
sample_weight, centers_old, centers_new,
|
||||
weight_sums, labels, indices)
|
||||
|
||||
free(indices)
|
||||
|
||||
|
||||
cdef void update_center_sparse(
|
||||
int cluster_idx,
|
||||
const floating[::1] X_data, # IN
|
||||
const int[::1] X_indices, # IN
|
||||
const int[::1] X_indptr, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int *indices) noexcept nogil: # TMP
|
||||
"""Update of a single center for sparse MinibatchKMeans"""
|
||||
cdef:
|
||||
int n_samples = sample_weight.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
floating alpha
|
||||
int n_indices
|
||||
int k, sample_idx, feature_idx
|
||||
|
||||
floating wsum = 0
|
||||
|
||||
# indices = np.where(labels == cluster_idx)[0]
|
||||
k = 0
|
||||
for sample_idx in range(n_samples):
|
||||
if labels[sample_idx] == cluster_idx:
|
||||
indices[k] = sample_idx
|
||||
wsum += sample_weight[sample_idx]
|
||||
k += 1
|
||||
n_indices = k
|
||||
|
||||
if wsum > 0:
|
||||
# Undo the previous count-based scaling for this cluster center:
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx] * weight_sums[cluster_idx]
|
||||
|
||||
# Update cluster with new point members
|
||||
for k in range(n_indices):
|
||||
sample_idx = indices[k]
|
||||
for feature_idx in range(X_indptr[sample_idx], X_indptr[sample_idx + 1]):
|
||||
centers_new[cluster_idx, X_indices[feature_idx]] += X_data[feature_idx] * sample_weight[sample_idx]
|
||||
|
||||
# Update the count statistics for this center
|
||||
weight_sums[cluster_idx] += wsum
|
||||
|
||||
# Rescale to compute mean of all points (old and new)
|
||||
alpha = 1 / weight_sums[cluster_idx]
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] *= alpha
|
||||
else:
|
||||
# No sample was assigned to this cluster in this batch of data
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx]
|
||||
from cython cimport floating
|
||||
from cython.parallel cimport parallel, prange
|
||||
from libc.stdlib cimport malloc, free
|
||||
|
||||
|
||||
def _minibatch_update_dense(
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int n_threads):
|
||||
"""Update of the centers for dense MiniBatchKMeans.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : ndarray of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration.
|
||||
|
||||
weight_sums : ndarray of shape (n_clusters,), dtype=floating
|
||||
Current sums of the accumulated weights for each center.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
"""
|
||||
cdef:
|
||||
int n_samples = X.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int cluster_idx
|
||||
|
||||
int *indices
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
indices = <int*> malloc(n_samples * sizeof(int))
|
||||
|
||||
for cluster_idx in prange(n_clusters, schedule="static"):
|
||||
update_center_dense(cluster_idx, X, sample_weight,
|
||||
centers_old, centers_new, weight_sums, labels,
|
||||
indices)
|
||||
|
||||
free(indices)
|
||||
|
||||
|
||||
cdef void update_center_dense(
|
||||
int cluster_idx,
|
||||
const floating[:, ::1] X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int *indices) noexcept nogil: # TMP
|
||||
"""Update of a single center for dense MinibatchKMeans"""
|
||||
cdef:
|
||||
int n_samples = sample_weight.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
floating alpha
|
||||
int n_indices
|
||||
int k, sample_idx, feature_idx
|
||||
|
||||
floating wsum = 0
|
||||
|
||||
# indices = np.where(labels == cluster_idx)[0]
|
||||
k = 0
|
||||
for sample_idx in range(n_samples):
|
||||
if labels[sample_idx] == cluster_idx:
|
||||
indices[k] = sample_idx
|
||||
wsum += sample_weight[sample_idx]
|
||||
k += 1
|
||||
n_indices = k
|
||||
|
||||
if wsum > 0:
|
||||
# Undo the previous count-based scaling for this cluster center
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx] * weight_sums[cluster_idx]
|
||||
|
||||
# Update cluster with new point members
|
||||
for k in range(n_indices):
|
||||
sample_idx = indices[k]
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] += X[sample_idx, feature_idx] * sample_weight[sample_idx]
|
||||
|
||||
# Update the count statistics for this center
|
||||
weight_sums[cluster_idx] += wsum
|
||||
|
||||
# Rescale to compute mean of all points (old and new)
|
||||
alpha = 1 / weight_sums[cluster_idx]
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] *= alpha
|
||||
else:
|
||||
# No sample was assigned to this cluster in this batch of data
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx]
|
||||
|
||||
|
||||
def _minibatch_update_sparse(
|
||||
X, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int n_threads):
|
||||
"""Update of the centers for sparse MiniBatchKMeans.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
X : sparse matrix of shape (n_samples, n_features), dtype=floating
|
||||
The observations to cluster. Must be in CSR format.
|
||||
|
||||
sample_weight : ndarray of shape (n_samples,), dtype=floating
|
||||
The weights for each observation in X.
|
||||
|
||||
centers_old : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers before previous iteration, placeholder for the centers after
|
||||
previous iteration.
|
||||
|
||||
centers_new : ndarray of shape (n_clusters, n_features), dtype=floating
|
||||
Centers after previous iteration, placeholder for the new centers
|
||||
computed during this iteration.
|
||||
|
||||
weight_sums : ndarray of shape (n_clusters,), dtype=floating
|
||||
Current sums of the accumulated weights for each center.
|
||||
|
||||
labels : ndarray of shape (n_samples,), dtype=int
|
||||
labels assignment.
|
||||
|
||||
n_threads : int
|
||||
The number of threads to be used by openmp.
|
||||
"""
|
||||
cdef:
|
||||
floating[::1] X_data = X.data
|
||||
int[::1] X_indices = X.indices
|
||||
int[::1] X_indptr = X.indptr
|
||||
int n_samples = X.shape[0]
|
||||
int n_clusters = centers_old.shape[0]
|
||||
int cluster_idx
|
||||
|
||||
int *indices
|
||||
|
||||
with nogil, parallel(num_threads=n_threads):
|
||||
indices = <int*> malloc(n_samples * sizeof(int))
|
||||
|
||||
for cluster_idx in prange(n_clusters, schedule="static"):
|
||||
update_center_sparse(cluster_idx, X_data, X_indices, X_indptr,
|
||||
sample_weight, centers_old, centers_new,
|
||||
weight_sums, labels, indices)
|
||||
|
||||
free(indices)
|
||||
|
||||
|
||||
cdef void update_center_sparse(
|
||||
int cluster_idx,
|
||||
const floating[::1] X_data, # IN
|
||||
const int[::1] X_indices, # IN
|
||||
const int[::1] X_indptr, # IN
|
||||
const floating[::1] sample_weight, # IN
|
||||
const floating[:, ::1] centers_old, # IN
|
||||
floating[:, ::1] centers_new, # OUT
|
||||
floating[::1] weight_sums, # INOUT
|
||||
const int[::1] labels, # IN
|
||||
int *indices) noexcept nogil: # TMP
|
||||
"""Update of a single center for sparse MinibatchKMeans"""
|
||||
cdef:
|
||||
int n_samples = sample_weight.shape[0]
|
||||
int n_features = centers_old.shape[1]
|
||||
floating alpha
|
||||
int n_indices
|
||||
int k, sample_idx, feature_idx
|
||||
|
||||
floating wsum = 0
|
||||
|
||||
# indices = np.where(labels == cluster_idx)[0]
|
||||
k = 0
|
||||
for sample_idx in range(n_samples):
|
||||
if labels[sample_idx] == cluster_idx:
|
||||
indices[k] = sample_idx
|
||||
wsum += sample_weight[sample_idx]
|
||||
k += 1
|
||||
n_indices = k
|
||||
|
||||
if wsum > 0:
|
||||
# Undo the previous count-based scaling for this cluster center:
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx] * weight_sums[cluster_idx]
|
||||
|
||||
# Update cluster with new point members
|
||||
for k in range(n_indices):
|
||||
sample_idx = indices[k]
|
||||
for feature_idx in range(X_indptr[sample_idx], X_indptr[sample_idx + 1]):
|
||||
centers_new[cluster_idx, X_indices[feature_idx]] += X_data[feature_idx] * sample_weight[sample_idx]
|
||||
|
||||
# Update the count statistics for this center
|
||||
weight_sums[cluster_idx] += wsum
|
||||
|
||||
# Rescale to compute mean of all points (old and new)
|
||||
alpha = 1 / weight_sums[cluster_idx]
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] *= alpha
|
||||
else:
|
||||
# No sample was assigned to this cluster in this batch of data
|
||||
for feature_idx in range(n_features):
|
||||
centers_new[cluster_idx, feature_idx] = centers_old[cluster_idx, feature_idx]
|
||||
|
||||
@@ -1,26 +1,26 @@
|
||||
cluster_extension_metadata = {
|
||||
'_dbscan_inner':
|
||||
{'sources': [cython_gen_cpp.process('_dbscan_inner.pyx')]},
|
||||
'_hierarchical_fast':
|
||||
{'sources': [cython_gen_cpp.process('_hierarchical_fast.pyx'), metrics_cython_tree]},
|
||||
'_k_means_common':
|
||||
{'sources': [cython_gen.process('_k_means_common.pyx')], 'dependencies': [openmp_dep]},
|
||||
'_k_means_lloyd':
|
||||
{'sources': [cython_gen.process('_k_means_lloyd.pyx')], 'dependencies': [openmp_dep]},
|
||||
'_k_means_elkan':
|
||||
{'sources': [cython_gen.process('_k_means_elkan.pyx')], 'dependencies': [openmp_dep]},
|
||||
'_k_means_minibatch':
|
||||
{'sources': [cython_gen.process('_k_means_minibatch.pyx')], 'dependencies': [openmp_dep]},
|
||||
}
|
||||
|
||||
foreach ext_name, ext_dict : cluster_extension_metadata
|
||||
py.extension_module(
|
||||
ext_name,
|
||||
[ext_dict.get('sources'), utils_cython_tree],
|
||||
dependencies: [np_dep] + ext_dict.get('dependencies', []),
|
||||
subdir: 'sklearn/cluster',
|
||||
install: true
|
||||
)
|
||||
endforeach
|
||||
|
||||
subdir('_hdbscan')
|
||||
cluster_extension_metadata = {
|
||||
'_dbscan_inner':
|
||||
{'sources': [cython_gen_cpp.process('_dbscan_inner.pyx')]},
|
||||
'_hierarchical_fast':
|
||||
{'sources': [cython_gen_cpp.process('_hierarchical_fast.pyx'), metrics_cython_tree]},
|
||||
'_k_means_common':
|
||||
{'sources': [cython_gen.process('_k_means_common.pyx')], 'dependencies': [openmp_dep]},
|
||||
'_k_means_lloyd':
|
||||
{'sources': [cython_gen.process('_k_means_lloyd.pyx')], 'dependencies': [openmp_dep]},
|
||||
'_k_means_elkan':
|
||||
{'sources': [cython_gen.process('_k_means_elkan.pyx')], 'dependencies': [openmp_dep]},
|
||||
'_k_means_minibatch':
|
||||
{'sources': [cython_gen.process('_k_means_minibatch.pyx')], 'dependencies': [openmp_dep]},
|
||||
}
|
||||
|
||||
foreach ext_name, ext_dict : cluster_extension_metadata
|
||||
py.extension_module(
|
||||
ext_name,
|
||||
[ext_dict.get('sources'), utils_cython_tree],
|
||||
dependencies: [np_dep] + ext_dict.get('dependencies', []),
|
||||
subdir: 'sklearn/cluster',
|
||||
install: true
|
||||
)
|
||||
endforeach
|
||||
|
||||
subdir('_hdbscan')
|
||||
|
||||
Reference in New Issue
Block a user