1
|
|
|
"""!
|
2
|
|
|
|
3
|
|
|
@brief Cluster analysis algorithm: Expectation-Maximization Algorithm (EMA).
|
4
|
|
|
@details Implementation based on article:
|
5
|
|
|
-
|
6
|
|
|
|
7
|
|
|
@authors Andrei Novikov ([email protected])
|
8
|
|
|
@date 2014-2017
|
9
|
|
|
@copyright GNU Public License
|
10
|
|
|
|
11
|
|
|
@cond GNU_PUBLIC_LICENSE
|
12
|
|
|
PyClustering is free software: you can redistribute it and/or modify
|
13
|
|
|
it under the terms of the GNU General Public License as published by
|
14
|
|
|
the Free Software Foundation, either version 3 of the License, or
|
15
|
|
|
(at your option) any later version.
|
16
|
|
|
|
17
|
|
|
PyClustering is distributed in the hope that it will be useful,
|
18
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
19
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
20
|
|
|
GNU General Public License for more details.
|
21
|
|
|
|
22
|
|
|
You should have received a copy of the GNU General Public License
|
23
|
|
|
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
24
|
|
|
@endcond
|
25
|
|
|
|
26
|
|
|
"""
|
27
|
|
|
|
28
|
|
|
|
29
|
|
|
import numpy;
|
|
|
|
|
30
|
|
|
|
31
|
|
|
from pyclustering.cluster import cluster_visualizer;
|
32
|
|
|
from pyclustering.cluster.center_initializer import kmeans_plusplus_initializer;
|
33
|
|
|
from pyclustering.cluster.kmeans import kmeans;
|
34
|
|
|
|
35
|
|
|
from pyclustering.utils import pi, calculate_ellipse_description;
|
36
|
|
|
|
37
|
|
|
import matplotlib.pyplot as plt;
|
|
|
|
|
38
|
|
|
from matplotlib import patches;
|
|
|
|
|
39
|
|
|
|
40
|
|
|
|
41
|
|
|
|
42
|
|
|
def gaussian(data, mean, covariance):
|
43
|
|
|
dimension = float(len(data[0]));
|
44
|
|
|
|
45
|
|
|
if (dimension != 1.0):
|
46
|
|
|
inv_variance = numpy.linalg.pinv(covariance);
|
47
|
|
|
else:
|
48
|
|
|
inv_variance = 1.0 / covariance;
|
49
|
|
|
|
50
|
|
|
divider = (pi * 2.0) ** (dimension / 2.0) * numpy.sqrt(numpy.linalg.norm(covariance));
|
51
|
|
|
right_const = 1.0 / divider;
|
52
|
|
|
|
53
|
|
|
result = [];
|
54
|
|
|
|
55
|
|
|
for point in data:
|
56
|
|
|
mean_delta = point - mean;
|
57
|
|
|
point_gaussian = right_const * numpy.exp( -0.5 * mean_delta.dot(inv_variance).dot(numpy.transpose(mean_delta)) );
|
58
|
|
|
result.append(point_gaussian);
|
59
|
|
|
|
60
|
|
|
return result;
|
61
|
|
|
|
62
|
|
|
|
63
|
|
|
|
64
|
|
|
class ema_observer:
|
65
|
|
|
def __init__(self):
|
66
|
|
|
self.__means_evolution = [];
|
67
|
|
|
self.__covariances_evolution = [];
|
68
|
|
|
self.__clusters_evolution = [];
|
69
|
|
|
|
70
|
|
|
|
71
|
|
|
def get_iterations(self):
|
72
|
|
|
return len(self.__means);
|
|
|
|
|
73
|
|
|
|
74
|
|
|
|
75
|
|
|
def get_means(self):
|
76
|
|
|
return self.__means_evolution;
|
77
|
|
|
|
78
|
|
|
|
79
|
|
|
def get_covariances(self):
|
80
|
|
|
return self.__covariances_evolution;
|
81
|
|
|
|
82
|
|
|
|
83
|
|
|
def notify(self, means, covariances, clusters):
|
84
|
|
|
self.__means_evolution.append(means);
|
85
|
|
|
self.__covariances_evolution.append(covariances);
|
86
|
|
|
self.__clusters_evolution.append(clusters);
|
87
|
|
|
|
88
|
|
|
|
89
|
|
|
|
90
|
|
|
class ema_visualizer:
|
91
|
|
|
@staticmethod
|
92
|
|
|
def show_clusters(clusters, sample, covariances, means, display = True):
|
93
|
|
|
visualizer = cluster_visualizer();
|
94
|
|
|
visualizer.append_clusters(clusters, sample);
|
95
|
|
|
figure = visualizer.show(display = False);
|
96
|
|
|
|
97
|
|
|
if (len(sample[0]) == 2):
|
98
|
|
|
ema_visualizer.__draw_ellipses(figure, visualizer, clusters, covariances, means);
|
99
|
|
|
|
100
|
|
|
if (display is True):
|
101
|
|
|
plt.show();
|
102
|
|
|
|
103
|
|
|
return figure;
|
104
|
|
|
|
105
|
|
|
|
106
|
|
|
@staticmethod
|
107
|
|
|
def __draw_ellipses(figure, visualizer, clusters, covariances, means):
|
108
|
|
|
print(len(clusters));
|
109
|
|
|
print([len(cluster) for cluster in clusters]);
|
110
|
|
|
print(clusters);
|
111
|
|
|
|
112
|
|
|
ax = figure.get_axes()[0];
|
113
|
|
|
|
114
|
|
|
for index in range(len(clusters)):
|
115
|
|
|
angle, width, height = calculate_ellipse_description(covariances[index]);
|
116
|
|
|
color = visualizer.get_cluster_color(index, 0);
|
117
|
|
|
|
118
|
|
|
ema_visualizer.__draw_ellipse(ax, means[index][0], means[index][1], angle, width, height, color);
|
119
|
|
|
|
120
|
|
|
|
121
|
|
|
@staticmethod
|
122
|
|
|
def __draw_ellipse(ax, x, y, angle, width, height, color):
|
123
|
|
|
ellipse = patches.Ellipse((x, y), width, height, alpha=0.2, angle=angle, linewidth=2, fill=True, zorder=2, color=color);
|
124
|
|
|
ax.add_patch(ellipse);
|
125
|
|
|
|
126
|
|
|
|
127
|
|
|
class ema:
|
128
|
|
|
def __init__(self, data, amount_clusters, means = None, variances = None, observer = None, tolerance = 0.00001):
|
129
|
|
|
self.__data = numpy.array(data);
|
130
|
|
|
self.__amount_clusters = amount_clusters;
|
131
|
|
|
self.__tolerance = tolerance;
|
132
|
|
|
self.__observer = observer;
|
133
|
|
|
|
134
|
|
|
self.__means = means;
|
135
|
|
|
if ((means is None) or (variances is None)):
|
136
|
|
|
self.__means, self.__variances = self.__get_initial_parameters(data, amount_clusters);
|
137
|
|
|
|
138
|
|
|
self.__rc = [ [0.0] * len(self.__data) for _ in range(amount_clusters) ];
|
139
|
|
|
self.__pic = [1.0] * amount_clusters;
|
140
|
|
|
self.__clusters = [];
|
141
|
|
|
self.__gaussians = [ [] for _ in range(amount_clusters) ];
|
142
|
|
|
self.__stop = False;
|
143
|
|
|
|
144
|
|
|
|
145
|
|
|
def process(self):
|
146
|
|
|
self.__clusters = None;
|
147
|
|
|
|
148
|
|
|
previous_likelihood = -200000;
|
149
|
|
|
current_likelihood = -100000;
|
150
|
|
|
|
151
|
|
|
while( (self.__stop is False) and (abs(previous_likelihood - current_likelihood) > self.__tolerance) ):
|
152
|
|
|
self.__expectation_step();
|
153
|
|
|
self.__maximization_step();
|
154
|
|
|
|
155
|
|
|
previous_likelihood = current_likelihood;
|
156
|
|
|
current_likelihood = self.__log_likelihood();
|
157
|
|
|
self.__stop = self.__get_stop_condition();
|
158
|
|
|
|
159
|
|
|
self.__clusters = self.__extract_clusters();
|
160
|
|
|
|
161
|
|
|
|
162
|
|
|
def get_clusters(self):
|
163
|
|
|
return self.__clusters;
|
164
|
|
|
|
165
|
|
|
|
166
|
|
|
def get_centers(self):
|
167
|
|
|
return self.__means;
|
168
|
|
|
|
169
|
|
|
|
170
|
|
|
def get_covariances(self):
|
171
|
|
|
return self.__variances;
|
172
|
|
|
|
173
|
|
|
|
174
|
|
|
def __notify(self):
|
175
|
|
|
if (self.__observer is not None):
|
176
|
|
|
clusters = self.__extract_clusters();
|
177
|
|
|
self.__notify(self.__means, self.__variances, clusters);
|
|
|
|
|
178
|
|
|
|
179
|
|
|
|
180
|
|
|
def __extract_clusters(self):
|
181
|
|
|
clusters = [ [] for _ in range(self.__amount_clusters) ];
|
182
|
|
|
for index_point in range(len(self.__data)):
|
183
|
|
|
candidates = [];
|
184
|
|
|
for index_cluster in range(self.__amount_clusters):
|
185
|
|
|
candidates.append((index_cluster, self.__rc[index_cluster][index_point]));
|
186
|
|
|
|
187
|
|
|
index_winner = max(candidates, key = lambda candidate : candidate[1])[0];
|
188
|
|
|
clusters[index_winner].append(index_point);
|
189
|
|
|
|
190
|
|
|
clusters = [ cluster for cluster in clusters if len(cluster) > 0 ];
|
191
|
|
|
return clusters;
|
192
|
|
|
|
193
|
|
|
|
194
|
|
|
def __log_likelihood(self):
|
195
|
|
|
likelihood = 0.0;
|
196
|
|
|
|
197
|
|
|
for index_point in range(len(self.__data)):
|
198
|
|
|
particle = 0.0;
|
199
|
|
|
for index_cluster in range(self.__amount_clusters):
|
200
|
|
|
particle += self.__pic[index_cluster] * self.__gaussians[index_cluster][index_point];
|
201
|
|
|
|
202
|
|
|
likelihood += numpy.log(particle);
|
203
|
|
|
|
204
|
|
|
return likelihood;
|
205
|
|
|
|
206
|
|
|
|
207
|
|
|
def __probabilities(self, index_cluster, index_point):
|
208
|
|
|
divider = 0.0;
|
209
|
|
|
for i in range(self.__amount_clusters):
|
210
|
|
|
divider += self.__pic[i] * self.__gaussians[i][index_point];
|
211
|
|
|
|
212
|
|
|
rc = self.__pic[index_cluster] * self.__gaussians[index_cluster][index_point] / divider;
|
213
|
|
|
return rc;
|
214
|
|
|
|
215
|
|
|
|
216
|
|
|
def __expectation_step(self):
|
217
|
|
|
for index in range(self.__amount_clusters):
|
218
|
|
|
self.__gaussians[index] = gaussian(self.__data, self.__means[index], self.__variances[index]);
|
219
|
|
|
|
220
|
|
|
for index_cluster in range(self.__amount_clusters):
|
221
|
|
|
for index_point in range(len(self.__data)):
|
222
|
|
|
self.__rc[index_cluster][index_point] = self.__probabilities(index_cluster, index_point);
|
223
|
|
|
|
224
|
|
|
|
225
|
|
|
def __maximization_step(self):
|
226
|
|
|
self.__pic = [];
|
227
|
|
|
self.__means = [];
|
228
|
|
|
self.__variances = [];
|
229
|
|
|
|
230
|
|
|
amount_impossible_clusters = 0;
|
231
|
|
|
|
232
|
|
|
for index_cluster in range(self.__amount_clusters):
|
233
|
|
|
mc = numpy.sum(self.__rc[index_cluster]);
|
234
|
|
|
|
235
|
|
|
if (mc == 0.0):
|
236
|
|
|
amount_impossible_clusters += 1;
|
237
|
|
|
continue;
|
238
|
|
|
|
239
|
|
|
self.__pic.append( mc / len(self.__data) );
|
240
|
|
|
self.__means.append( self.__update_mean(self.__rc[index_cluster], mc) );
|
241
|
|
|
self.__variances.append( self.__update_covariance(self.__means[-1], self.__rc[index_cluster], mc) );
|
242
|
|
|
|
243
|
|
|
self.__amount_clusters -= amount_impossible_clusters;
|
244
|
|
|
|
245
|
|
|
|
246
|
|
|
def __get_stop_condition(self):
|
247
|
|
|
for covariance in self.__variances:
|
248
|
|
|
if (numpy.linalg.norm(covariance) == 0.0):
|
249
|
|
|
return True;
|
250
|
|
|
|
251
|
|
|
return False;
|
252
|
|
|
|
253
|
|
|
|
254
|
|
|
def __update_covariance(self, means, rc, mc):
|
255
|
|
|
covariance = 0.0;
|
256
|
|
|
for index_point in range(len(self.__data)):
|
257
|
|
|
deviation = numpy.array( [ self.__data[index_point] - means ]);
|
258
|
|
|
covariance += rc[index_point] * deviation.T.dot(deviation);
|
259
|
|
|
|
260
|
|
|
covariance = covariance / mc;
|
261
|
|
|
return covariance;
|
262
|
|
|
|
263
|
|
|
|
264
|
|
|
def __update_mean(self, rc, mc):
|
265
|
|
|
mean = 0.0;
|
266
|
|
|
for index_point in range(len(self.__data)):
|
267
|
|
|
mean += rc[index_point] * self.__data[index_point];
|
268
|
|
|
|
269
|
|
|
mean = mean / mc;
|
270
|
|
|
return mean;
|
271
|
|
|
|
272
|
|
|
|
273
|
|
|
def __get_initial_parameters(self, sample, amount_clusters):
|
274
|
|
|
initial_centers = kmeans_plusplus_initializer(sample, amount_clusters).initialize();
|
275
|
|
|
kmeans_instance = kmeans(sample, initial_centers, ccore = True);
|
276
|
|
|
kmeans_instance.process();
|
277
|
|
|
|
278
|
|
|
means = kmeans_instance.get_centers();
|
279
|
|
|
|
280
|
|
|
covariances = [];
|
281
|
|
|
initial_clusters = kmeans_instance.get_clusters();
|
282
|
|
|
for initial_cluster in initial_clusters:
|
283
|
|
|
cluster_sample = [ sample[index_point] for index_point in initial_cluster ];
|
284
|
|
|
covariances.append(numpy.cov(cluster_sample, rowvar = False));
|
285
|
|
|
|
286
|
|
|
return means, covariances; |
This can be caused by one of the following:
1. Missing Dependencies
This error could indicate a configuration issue of Pylint. Make sure that your libraries are available by adding the necessary commands.
2. Missing __init__.py files
This error could also result from missing
__init__.py
files in your module folders. Make sure that you place one file in each sub-folder.