1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.clustering;
23
24 import java.util.ArrayList;
25 import java.util.Collection;
26 import java.util.Collections;
27 import java.util.List;
28
29 import org.hipparchus.clustering.distance.DistanceMeasure;
30 import org.hipparchus.clustering.distance.EuclideanDistance;
31 import org.hipparchus.exception.LocalizedCoreFormats;
32 import org.hipparchus.exception.MathIllegalArgumentException;
33 import org.hipparchus.exception.MathIllegalStateException;
34 import org.hipparchus.linear.MatrixUtils;
35 import org.hipparchus.linear.RealMatrix;
36 import org.hipparchus.random.JDKRandomGenerator;
37 import org.hipparchus.random.RandomGenerator;
38 import org.hipparchus.util.FastMath;
39 import org.hipparchus.util.MathArrays;
40 import org.hipparchus.util.MathUtils;
41
42 /**
43 * Fuzzy K-Means clustering algorithm.
44 * <p>
45 * The Fuzzy K-Means algorithm is a variation of the classical K-Means algorithm, with the
46 * major difference that a single data point is not uniquely assigned to a single cluster.
47 * Instead, each point i has a set of weights u<sub>ij</sub> which indicate the degree of membership
48 * to the cluster j.
49 * <p>The algorithm then tries to minimize the objective function:
50 * \[
51 * J = \sum_{i=1}^C\sum_{k=1]{N} u_{i,k}^m d_{i,k}^2
52 * \]
53 * with \(d_{i,k}\) being the distance between data point i and the cluster center k.
54 * </p>
55 * <p>The algorithm requires two parameters:</p>
56 * <ul>
57 * <li>k: the number of clusters
58 * <li>fuzziness: determines the level of cluster fuzziness, larger values lead to fuzzier clusters
59 * </ul>
60 * <p>Additional, optional parameters:</p>
61 * <ul>
62 * <li>maxIterations: the maximum number of iterations
63 * <li>epsilon: the convergence criteria, default is 1e-3
64 * </ul>
65 * <p>
66 * The fuzzy variant of the K-Means algorithm is more robust with regard to the selection
67 * of the initial cluster centers.
68 * </p>
69 *
70 * @param <T> type of the points to cluster
71 */
72 public class FuzzyKMeansClusterer<T extends Clusterable> extends Clusterer<T> {
73
74 /** The default value for the convergence criteria. */
75 private static final double DEFAULT_EPSILON = 1e-3;
76
77 /** The number of clusters. */
78 private final int k;
79
80 /** The maximum number of iterations. */
81 private final int maxIterations;
82
83 /** The fuzziness factor. */
84 private final double fuzziness;
85
86 /** The convergence criteria. */
87 private final double epsilon;
88
89 /** Random generator for choosing initial centers. */
90 private final RandomGenerator random;
91
92 /** The membership matrix. */
93 private double[][] membershipMatrix;
94
95 /** The list of points used in the last call to {@link #cluster(Collection)}. */
96 private List<T> points;
97
98 /** The list of clusters resulting from the last call to {@link #cluster(Collection)}. */
99 private List<CentroidCluster<T>> clusters;
100
101 /**
102 * Creates a new instance of a FuzzyKMeansClusterer.
103 * <p>
104 * The euclidean distance will be used as default distance measure.
105 *
106 * @param k the number of clusters to split the data into
107 * @param fuzziness the fuzziness factor, must be > 1.0
108 * @throws MathIllegalArgumentException if {@code fuzziness <= 1.0}
109 */
110 public FuzzyKMeansClusterer(final int k, final double fuzziness) throws MathIllegalArgumentException {
111 this(k, fuzziness, -1, new EuclideanDistance());
112 }
113
114 /**
115 * Creates a new instance of a FuzzyKMeansClusterer.
116 *
117 * @param k the number of clusters to split the data into
118 * @param fuzziness the fuzziness factor, must be > 1.0
119 * @param maxIterations the maximum number of iterations to run the algorithm for.
120 * If negative, no maximum will be used.
121 * @param measure the distance measure to use
122 * @throws MathIllegalArgumentException if {@code fuzziness <= 1.0}
123 */
124 public FuzzyKMeansClusterer(final int k, final double fuzziness,
125 final int maxIterations, final DistanceMeasure measure)
126 throws MathIllegalArgumentException {
127 this(k, fuzziness, maxIterations, measure, DEFAULT_EPSILON, new JDKRandomGenerator());
128 }
129
130 /**
131 * Creates a new instance of a FuzzyKMeansClusterer.
132 *
133 * @param k the number of clusters to split the data into
134 * @param fuzziness the fuzziness factor, must be > 1.0
135 * @param maxIterations the maximum number of iterations to run the algorithm for.
136 * If negative, no maximum will be used.
137 * @param measure the distance measure to use
138 * @param epsilon the convergence criteria (default is 1e-3)
139 * @param random random generator to use for choosing initial centers
140 * @throws MathIllegalArgumentException if {@code fuzziness <= 1.0}
141 */
142 public FuzzyKMeansClusterer(final int k, final double fuzziness,
143 final int maxIterations, final DistanceMeasure measure,
144 final double epsilon, final RandomGenerator random)
145 throws MathIllegalArgumentException {
146
147 super(measure);
148
149 if (fuzziness <= 1.0d) {
150 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
151 fuzziness, 1.0);
152 }
153 this.k = k;
154 this.fuzziness = fuzziness;
155 this.maxIterations = maxIterations;
156 this.epsilon = epsilon;
157 this.random = random;
158
159 this.membershipMatrix = null;
160 this.points = null;
161 this.clusters = null;
162 }
163
164 /**
165 * Return the number of clusters this instance will use.
166 * @return the number of clusters
167 */
168 public int getK() {
169 return k;
170 }
171
172 /**
173 * Returns the fuzziness factor used by this instance.
174 * @return the fuzziness factor
175 */
176 public double getFuzziness() {
177 return fuzziness;
178 }
179
180 /**
181 * Returns the maximum number of iterations this instance will use.
182 * @return the maximum number of iterations, or -1 if no maximum is set
183 */
184 public int getMaxIterations() {
185 return maxIterations;
186 }
187
188 /**
189 * Returns the convergence criteria used by this instance.
190 * @return the convergence criteria
191 */
192 public double getEpsilon() {
193 return epsilon;
194 }
195
196 /**
197 * Returns the random generator this instance will use.
198 * @return the random generator
199 */
200 public RandomGenerator getRandomGenerator() {
201 return random;
202 }
203
204 /**
205 * Returns the {@code nxk} membership matrix, where {@code n} is the number
206 * of data points and {@code k} the number of clusters.
207 * <p>
208 * The element U<sub>i,j</sub> represents the membership value for data point {@code i}
209 * to cluster {@code j}.
210 *
211 * @return the membership matrix
212 * @throws MathIllegalStateException if {@link #cluster(Collection)} has not been called before
213 */
214 public RealMatrix getMembershipMatrix() {
215 if (membershipMatrix == null) {
216 throw new MathIllegalStateException(LocalizedCoreFormats.ILLEGAL_STATE);
217 }
218 return MatrixUtils.createRealMatrix(membershipMatrix);
219 }
220
221 /**
222 * Returns an unmodifiable list of the data points used in the last
223 * call to {@link #cluster(Collection)}.
224 * @return the list of data points, or {@code null} if {@link #cluster(Collection)} has
225 * not been called before.
226 */
227 public List<T> getDataPoints() {
228 return points;
229 }
230
231 /**
232 * Returns the list of clusters resulting from the last call to {@link #cluster(Collection)}.
233 * @return the list of clusters, or {@code null} if {@link #cluster(Collection)} has
234 * not been called before.
235 */
236 public List<CentroidCluster<T>> getClusters() {
237 return clusters;
238 }
239
240 /**
241 * Get the value of the objective function.
242 * @return the objective function evaluation as double value
243 * @throws MathIllegalStateException if {@link #cluster(Collection)} has not been called before
244 */
245 public double getObjectiveFunctionValue() {
246 if (points == null || clusters == null) {
247 throw new MathIllegalStateException(LocalizedCoreFormats.ILLEGAL_STATE);
248 }
249
250 int i = 0;
251 double objFunction = 0.0;
252 for (final T point : points) {
253 int j = 0;
254 for (final CentroidCluster<T> cluster : clusters) {
255 final double dist = distance(point, cluster.getCenter());
256 objFunction += (dist * dist) * FastMath.pow(membershipMatrix[i][j], fuzziness);
257 j++;
258 }
259 i++;
260 }
261 return objFunction;
262 }
263
264 /**
265 * Performs Fuzzy K-Means cluster analysis.
266 *
267 * @param dataPoints the points to cluster
268 * @return the list of clusters
269 * @throws MathIllegalArgumentException if the data points are null or the number
270 * of clusters is larger than the number of data points
271 */
272 @Override
273 public List<CentroidCluster<T>> cluster(final Collection<T> dataPoints)
274 throws MathIllegalArgumentException {
275
276 // sanity checks
277 MathUtils.checkNotNull(dataPoints);
278
279 final int size = dataPoints.size();
280
281 // number of clusters has to be smaller or equal the number of data points
282 if (size < k) {
283 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
284 size, k);
285 }
286
287 // copy the input collection to an unmodifiable list with indexed access
288 points = Collections.unmodifiableList(new ArrayList<>(dataPoints));
289 clusters = new ArrayList<>();
290 membershipMatrix = new double[size][k];
291 final double[][] oldMatrix = new double[size][k];
292
293 // if no points are provided, return an empty list of clusters
294 if (size == 0) {
295 return clusters;
296 }
297
298 initializeMembershipMatrix();
299
300 // there is at least one point
301 final int pointDimension = points.get(0).getPoint().length;
302 for (int i = 0; i < k; i++) {
303 clusters.add(new CentroidCluster<T>(new DoublePoint(new double[pointDimension])));
304 }
305
306 int iteration = 0;
307 final int max = (maxIterations < 0) ? Integer.MAX_VALUE : maxIterations;
308 double difference;
309
310 do {
311 saveMembershipMatrix(oldMatrix);
312 updateClusterCenters();
313 updateMembershipMatrix();
314 difference = calculateMaxMembershipChange(oldMatrix);
315 } while (difference > epsilon && ++iteration < max);
316
317 return clusters;
318 }
319
320 /**
321 * Update the cluster centers.
322 */
323 private void updateClusterCenters() {
324 int j = 0;
325 final List<CentroidCluster<T>> newClusters = new ArrayList<>(k);
326 for (final CentroidCluster<T> cluster : clusters) {
327 final Clusterable center = cluster.getCenter();
328 int i = 0;
329 double[] arr = new double[center.getPoint().length];
330 double sum = 0.0;
331 for (final T point : points) {
332 final double u = FastMath.pow(membershipMatrix[i][j], fuzziness);
333 final double[] pointArr = point.getPoint();
334 for (int idx = 0; idx < arr.length; idx++) {
335 arr[idx] += u * pointArr[idx];
336 }
337 sum += u;
338 i++;
339 }
340 MathArrays.scaleInPlace(1.0 / sum, arr);
341 newClusters.add(new CentroidCluster<T>(new DoublePoint(arr)));
342 j++;
343 }
344 clusters.clear();
345 clusters = newClusters;
346 }
347
348 /**
349 * Updates the membership matrix and assigns the points to the cluster with
350 * the highest membership.
351 */
352 private void updateMembershipMatrix() {
353 for (int i = 0; i < points.size(); i++) {
354 final T point = points.get(i);
355 double maxMembership = Double.MIN_VALUE;
356 int newCluster = -1;
357 for (int j = 0; j < clusters.size(); j++) {
358 double sum = 0.0;
359 final double distA = FastMath.abs(distance(point, clusters.get(j).getCenter()));
360
361 if (distA != 0.0) {
362 for (final CentroidCluster<T> c : clusters) {
363 final double distB = FastMath.abs(distance(point, c.getCenter()));
364 if (distB == 0.0) {
365 sum = Double.POSITIVE_INFINITY;
366 break;
367 }
368 sum += FastMath.pow(distA / distB, 2.0 / (fuzziness - 1.0));
369 }
370 }
371
372 double membership;
373 if (sum == 0.0) {
374 membership = 1.0;
375 } else if (sum == Double.POSITIVE_INFINITY) {
376 membership = 0.0;
377 } else {
378 membership = 1.0 / sum;
379 }
380 membershipMatrix[i][j] = membership;
381
382 if (membershipMatrix[i][j] > maxMembership) {
383 maxMembership = membershipMatrix[i][j];
384 newCluster = j;
385 }
386 }
387 clusters.get(newCluster).addPoint(point);
388 }
389 }
390
391 /**
392 * Initialize the membership matrix with random values.
393 */
394 private void initializeMembershipMatrix() {
395 for (int i = 0; i < points.size(); i++) {
396 for (int j = 0; j < k; j++) {
397 membershipMatrix[i][j] = random.nextDouble();
398 }
399 membershipMatrix[i] = MathArrays.normalizeArray(membershipMatrix[i], 1.0);
400 }
401 }
402
403 /**
404 * Calculate the maximum element-by-element change of the membership matrix
405 * for the current iteration.
406 *
407 * @param matrix the membership matrix of the previous iteration
408 * @return the maximum membership matrix change
409 */
410 private double calculateMaxMembershipChange(final double[][] matrix) {
411 double maxMembership = 0.0;
412 for (int i = 0; i < points.size(); i++) {
413 for (int j = 0; j < clusters.size(); j++) {
414 double v = FastMath.abs(membershipMatrix[i][j] - matrix[i][j]);
415 maxMembership = FastMath.max(v, maxMembership);
416 }
417 }
418 return maxMembership;
419 }
420
421 /**
422 * Copy the membership matrix into the provided matrix.
423 *
424 * @param matrix the place to store the membership matrix
425 */
426 private void saveMembershipMatrix(final double[][] matrix) {
427 for (int i = 0; i < points.size(); i++) {
428 System.arraycopy(membershipMatrix[i], 0, matrix[i], 0, clusters.size());
429 }
430 }
431
432 }