-
Notifications
You must be signed in to change notification settings - Fork 1
/
som_batchtrain.py
405 lines (351 loc) · 14.9 KB
/
som_batchtrain.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 25 09:19:39 2024
@author: sdd380
"""
import numpy as np
import inspect
from som_set import som_set
from som_train_struct import som_train_struct
from som_unit_dists import som_unit_dists
import time
import math
import matplotlib.pyplot as plt
from scipy.sparse import coo_matrix
from som_set import som_set
from datetime import datetime
def som_batchtrain(sMap, D, *args, **kwargs):
############# DETAILED DESCRIPTION ########################################
#
# som_batchtrain
#
# PURPOSE
#
# Trains a Self-Organizing Map using the batch algorithm.
# Trains the given SOM (sM or M above) with the given training data
# (sD or D) using batch training algorithm. If no optional arguments
# (argID, value) are given, a default training is done. Using optional
# arguments the training parameters can be specified. Returns the
# trained and updated SOM and a train struct which contains
# information on the training.
#
# REFERENCES
#
# Kohonen, T., "Self-Organizing Map", 2nd ed., Springer-Verlag,
# Berlin, 1995, pp. 127-128.
# Kohonen, T., "Things you haven't heard about the Self-Organizing
# Map", In proceedings of International Conference
# on Neural Networks (ICNN), San Francisco, 1993, pp. 1147-1156.
# REQUIRED INPUT ARGUMENTS
#
# sM The map to be trained.
# (struct) map struct
# (matrix) codebook matrix (field .data of map struct)
# Size is either [munits dim], in which case the map grid
# dimensions (msize) should be specified with optional arguments,
# or [msize(1) ... msize(k) dim] in which case the map
# grid dimensions are taken from the size of the matrix.
# Lattice, by default, is 'rect' and shape 'sheet'.
# D Training data.
# (struct) data struct
# (matrix) data matrix, size [dlen dim]
# check arguments
num_args = len(args)
if num_args < 2:
raise ValueError(f"Too few arguments: expected at least 2, got {num_args}")
if num_args > float('inf'):
raise ValueError(f"Too many arguments: expected at most inf, got {num_args}")
# map
struct_mode = isinstance(sMap, dict)
if struct_mode == True:
sTopol = sMap['topol']
else:
dlens, ndims = sMap.shape
if ndims(sMap) > 2:
si= sMap.shape
dim = si[0, -1]
msize = si[0:-1]
munits, dim = sMap['codebook'].shape
# data
if isinstance(D, dict):
data_name = D['name']
D = D['data']
else:
sig = inspect.signature(som_batchtrain)
params = sig.parameters
param_list = list(params.keys())
data_name = param_list[1]
nonempty = np.where(np.sum(np.isnan(D), axis=1) < dim)[0]
D = D[nonempty,:]
dlen, ddim = D.shape
if dim != ddim:
raise ValueError('Map and data input space dimensions disagree.')
# args
sTrain = som_set('som_train', *['algorithm', 'neigh', 'mask', 'data_name'],
**{'algorithm': 'batch', 'neigh': sMap['neigh'], 'mask': sMap['mask'], 'data_name': data_name})
radius = []
tracking = 1
weights = np.ones(1)
i=0
while i<= len(args)-1:
argok = 1
if isinstance(args[i], str):
if args[i]== 'msize':
sTopol['msize'] = kwargs[args[i]]
elif args[i] == 'lattice':
sTopol['lattice']= kwargs[args[i]]
elif args[i] == 'shape':
sTopol['shape']= kwargs[args[i]]
elif args[i] == 'mask':
sTrain['mask']= kwargs[args[i]]
elif args[i] == 'neigh':
sTrain['neigh']= kwargs[args[i]]
elif args[i] == 'trainlen':
sTrain['trainlen'] = kwargs[args[i]]
elif args[i] == 'tracking':
tracking = kwargs[args[i]]
elif args[i] == 'weights':
weights = kwargs[args[i]]
elif args[i] == 'radius_ini':
sTrain['radius_ini'] = np.array(kwargs[args[i]])
elif args[i] == 'radius_fin':
sTrain['radius_fin'] = np.array(kwargs[args[i]])
elif args[i] == 'radius':
l = len(kwargs[args[i]])
if l==1:
sTrain['radius_ini']=np.array(kwargs[args[i]])
else:
sTrain['radius_ini']= np.array(kwargs[args[i][1]])
sTrain['radius_fin']= np.array(kwargs[args[i][-1]])
if l>2:
radius = np.array(kwargs[args[i]])
elif args[i] in {'sTrain','train','som_train'}:
sTrain = kwargs[args[i]]
elif args[i] in {'topol','sTopol','som_topol'}:
sTopol = kwargs[args[i]]
if np.prod(sTopol['msize']) != munits:
raise ValueError('Given map grid size does not match the codebook size.')
# unambiguous values
elif args[i] in {'hexa','rect'}:
sTopol['lattice'] = kwargs[args[i]]
elif args[i] in {'sheet','cyl','toroid'}:
sTopol['shape'] = kwargs[args[i]]
elif args[i] in {'gaussian','cutgauss','ep','bubble'}:
sTrain['neigh'] = kwargs[args[i]]
else:
argok=0
elif isinstance(args[i], dict) and 'type' in kwargs[args[i]]:
if kwargs[args[i]]['type'] == 'som_topol':
sTopol= kwargs[args[i]]
if np.prod(sTopol['msize']) != munits:
raise ValueError('Given map grid size does not match the codebook size.')
elif kwargs[args[i]]['type'] == 'som_train':
sTrain = kwargs[args[i]]
else:
argok =0
else:
argok=0
if argok == 0:
print("(som_batchtrain) Ignoring invalid argument #" + str(i))
i=i+1
# take only weights of non-empty vectors
if len(weights) > dlen:
weights = weights[nonempty]
# trainlen
if bool(radius):
sTrain['trainlen'] = len(radius)
# check topology
if bool(struct_mode):
if sTopol['lattice'] != sMap['topol']['lattice'] or sTopol['shape'] != sMap['topol']['shape'] or np.any(sTopol['msize'] != sMap['topol']['msize']):
print("Changing the original map topology")
sMap['topol'] = sTopol
# complement the training struct
sTrain = som_train_struct(*['sTrain','sMap', 'dlen'],**{'sTrain': sTrain,'sMap': sMap, 'dlen': dlen})
if not bool(np.any(sTrain['mask'])):
sTrain['mask'] = np.ones((dim,1))
#### INITIALIZE
M = sMap['codebook']
mask = sTrain['mask']
trainlen = sTrain['trainlen']
# neigbourhood radius
if trainlen == 1:
radius = np.array([sTrain['radius_ini']])
elif len(radius)<=2:
r0 = sTrain['radius_ini']
r1 = sTrain['radius_fin']
radius = r1 + (np.flip(np.arange(trainlen))/ (trainlen - 1)) * (r0 - r1)
# distance between map units in the output space
# Since in the case of gaussian and ep neighborhood functions, the
# equations utilize squares of the unit distances and in bubble case
# it doesn't matter which is used, the unitdistances and neighborhood
# radiuses are squared.
Ud = som_unit_dists(sTopol)
Ud = Ud**2
radius = radius **2
# zero neighborhood radius may cause div-by-zero error
eps = np.finfo(float).eps
if np.isscalar(radius):
if radius==0:
radius=eps
else:
radius = np.array(radius, dtype=np.float64)
if len(radius) == 1:
if radius[0] == 0:
radius[0] = eps
else:
# Replace all zeros with eps
radius[radius == 0] = eps
# The training algorithm involves calculating weighted Euclidian distances
# to all map units for each data vector. Basically this is done as
# for i=1:dlen,
# for j=1:munits,
# for k=1:dim
# Dist(j,i) = Dist(j,i) + mask(k) * (D(i,k) - M(j,k))^2;
# end
# end
# end
# where mask is the weighting vector for distance calculation. However, taking
# into account that distance between vectors m and v can be expressed as
# |m - v|^2 = sum_i ((m_i - v_i)^2) = sum_i (m_i^2 + v_i^2 - 2*m_i*v_i)
# this can be made much faster by transforming it to a matrix operation:
# Dist = (M.^2)*mask*ones(1,d) + ones(m,1)*mask'*(D'.^2) - 2*M*diag(mask)*D'
# Of the involved matrices, several are constant, as the mask and data do
# not change during training. Therefore they are calculated beforehand.
# For the case where there are unknown components in the data, each data
# vector will have an individual mask vector so that for that unit, the
# unknown components are not taken into account in distance calculation.
# In addition all NaN's are changed to zeros so that they don't screw up
# the matrix multiplications and behave correctly in updating step.
# D[147, :] = np.nan
Known = (~np.isnan(D)).astype(int)
D[~Known.astype(bool)] = 0
W1 = (mask * np.ones((1, dlen)) * Known.T)
# constant matrices
WD = 2 * np.dot(np.diag(mask.flatten()), D.T)
dconst = np.dot((D**2),mask).T
start = time.time()
qe = np.zeros((trainlen,1))
##### Action
# With the 'blen' parameter you can control the memory consumption
# of the algorithm, which is in practive directly proportional
# to munits*blen. If you're having problems with memory, try to
# set the value of blen lower.
blen = min((munits,dlen))
# reserve some space
bmus = np.zeros((1,dlen))
ddists = np.zeros((1,dlen))
for t in range(0, trainlen):
# batchy train - this is done a block of data (inds) at a time
# rather than in a single sweep to save memory consumption.
# The 'Dist' and 'Hw' matrices have size munits*blen
# which - if you have a lot of data - would be HUGE if you
# calculated it all at once. A single-sweep version would
# look like this:
# Dist = (M.^2)*W1 - M*WD; #+ W2*D2
# [ddists, bmus] = min(Dist);
# (notice that the W2*D2 term can be ignored since it is constant)
# This "batchy" version is the same as single-sweep if blen=dlen.
ddists = np.zeros(dlen)
bmus = np.zeros(dlen, dtype=int)
i0 = 0
while i0 + 1 <= dlen:
inds = np.arange(i0, min(dlen, i0 + blen))
i0 += blen
# Calculate Dist
Dist = np.square(M) @ W1[:, inds] - M @ WD[:, inds]
ddists[inds] = np.min(Dist, axis=0)
bmus[inds] = np.argmin(Dist, axis=0)
# tracking
if tracking >0:
ddists = ddists + dconst
ddists[ddists<0] = 0
qe[t] = np.mean(np.sqrt(ddists))
trackplot(M,D,tracking,start,t,qe)
H = np.zeros_like(Ud)
if sTrain['neigh']=='bubble':
H[Ud <= radius[t]] = 1
elif sTrain['neigh']=='gaussian':
H = np.exp(-Ud / (2 * radius[t]))
elif sTrain['neigh']=='cutgauss':
H = np.exp(-Ud / (2 * radius[t])) * (Ud <= radius[t])
elif sTrain['neigh']=='ep':
H = (1 - Ud / radius[t]) * (Ud <= radius[t])
else:
raise ValueError(f"Unknown neighborhood type: {sTrain['neigh']}")
# update
# In principle the updating step goes like this: replace each map unit
# by the average of the data vectors that were in its neighborhood.
# The contribution, or activation, of data vectors in the mean can
# be varied with the neighborhood function. This activation is given
# by matrix H. So, for each map unit the new weight vector is
#
# m = sum_i (h_i * d_i) / sum_i (h_i),
#
# where i denotes the index of data vector. Since the values of
# neighborhood function h_i are the same for all data vectors belonging to
# the Voronoi set of the same map unit, the calculation is actually done
# by first calculating a partition matrix P with elements p_ij=1 if the
# BMU of data vector j is i.
# Ensure bmus, weights, and np.arange(len(bmus)) have the same length
weights1=np.full(len(bmus), weights)
P = coo_matrix((weights1, (bmus, np.arange(len(bmus)))), shape=(munits, dlen))
S = np.dot(H , P @ D)
A = np.dot(H, P @ Known)
# only update units for which the "activation" is nonzero
nonzero = np.where(A>0)
M = np.zeros_like(A)
M[nonzero[0]] = S[nonzero[0]]/A[nonzero[0]]
# Build / clean up the return arguments
if tracking > 0:
print('\n')
# update structures
sTrain = som_set(sTrain, *['time'], **{'time':datetime.now().strftime('%Y-%m-%d %H:%M:%S')})
if struct_mode == 1:
sMap = som_set(sMap, *['codebook', 'mask', 'neigh'], **{'codebook': M, 'mask': sTrain['mask'], 'neigh': sTrain['neigh']})
sMap['trainhist'].append(sTrain)
else:
sMap = np.reshape(M, orig_size)
return sMap
def trackplot(M,D,tracking,start,n,qe):
l=len(qe)
elap_t = time.time() - start
# Total time
tot_t = elap_t * l / (n+1)
# Print training progress
print(f'\rTraining: {elap_t:.0f} / {tot_t:.0f} s', end='')
# Switch statement equivalent
tracking = 2 # Example value for tracking
if tracking == 1:
pass # Case 1: No action
elif tracking == 2:
# Plot quantization error after each epoch
plt.figure(figsize=(8, 6))
plt.plot(range(1, n + 1), qe[:n], 'b-', label='First part')
plt.plot(range(n + 1, l + 1), qe[n:l], 'r-', label='Second part')
plt.title('Quantization error after each epoch')
plt.legend()
plt.grid(True)
plt.draw()
plt.pause(0.001) # Needed for interactive plotting
else:
# Default case: subplot visualization
plt.figure(figsize=(10, 8))
plt.subplot(2, 1, 1)
plt.plot(range(1, n + 1), qe[:n], 'b-', label='First part')
plt.plot(range(n + 1, l + 1), qe[n:l], 'r-', label='Second part')
plt.title('Quantization error after each epoch')
plt.legend()
plt.grid(True)
plt.subplot(2, 1, 2)
# Example data M and D (replace with your actual data)
M = np.random.rand(10, 2)
D = np.random.rand(10, 2)
plt.plot(M[:, 0], M[:, 1], 'ro', label='Map units (o)')
plt.plot(D[:, 0], D[:, 1], 'b+', label='Data vectors (+)')
plt.title('First two components of map units (o) and data vectors (+)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.draw()
plt.pause(0.001) # Needed for interactive plotting
plt.show()