mvpa.mappers.som

21 """Mapper using a self-organizing map (SOM) for dimensionality reduction. 22 23 This mapper provides a simple, but pretty fast implementation of a 24 self-organizing map using an unsupervised training algorithm. It performs a 25 ND -> 2D mapping, which can for, example, be used for visualization of 26 high-dimensional data. 27 28 This SOM implementation uses squared Euclidean distance to determine 29 the best matching Kohonen unit and a Gaussian neighborhood influence 30 kernel. 31 """

32 - def __init__(self, kshape, niter, learning_rate=0.005, 33 iradius=None):

34 """ 35 :Parameters: 36 kshape: (int, int) 37 Shape of the internal Kohonen layer. Currently, only 2D Kohonen 38 layers are supported, although the length of an axis might be set 39 to 1. 40 niter: int 41 Number of iteration during network training. 42 learning_rate: float 43 Initial learning rate, which will continuously decreased during 44 network training. 45 iradius: float | None 46 Initial radius of the Gaussian neighborhood kernel radius, which 47 will continuously decreased during network training. If `None` 48 (default) the radius is set equal to the longest edge of the 49 Kohonen layer. 50 """ 51 # init base class 52 Mapper.__init__(self) 53 54 self.kshape = N.array(kshape, dtype='int') 55 56 if iradius is None: 57 self.radius = self.kshape.max() 58 else: 59 self.radius = iradius 60 61 # learning rate 62 self.lrate = learning_rate 63 64 # number of training iterations 65 self.niter = niter 66 67 # precompute whatever can be done 68 # scalar for decay of learning rate and radius across all iterations 69 self.iter_scale = self.niter / N.log(self.radius) 70 71 # the internal kohonen layer 72 self._K = None

73 74

75 - def train(self, ds):

76 """Perform network training. 77 78 :Parameter: 79 ds: Dataset 80 All samples in the dataset will be used for unsupervised training of 81 the SOM. 82 """ 83 # XXX initialize with clever default, e.g. plain of first two PCA 84 # components 85 self._K = N.random.standard_normal(tuple(self.kshape) + (ds.nfeatures,)) 86 87 # units weight vector deltas for batch training 88 # (height x width x #features) 89 unit_deltas = N.zeros(self._K.shape, dtype='float') 90 91 # precompute distance kernel between elements in the Kohonen layer 92 # that will remain constant throughout the training 93 # (just compute one quadrant, as the distances are symmetric) 94 # XXX maybe do other than squared Euclidean? 95 dqd = N.fromfunction(lambda x, y: (x**2 + y**2)**0.5, 96 self.kshape, dtype='float') 97 98 # for all iterations 99 for it in xrange(1, self.niter + 1): 100 # compute the neighborhood impact kernel for this iteration 101 # has to be recomputed since kernel shrinks over time 102 k = self._computeInfluenceKernel(it, dqd) 103 104 # for all training vectors 105 for s in ds.samples: 106 # determine closest unit (as element coordinate) 107 b = self._getBMU(s) 108 109 # train all units at once by unfolding the kernel (from the 110 # single quadrant that is precomputed), cutting it to the 111 # right shape and simply multiply it to the difference of target 112 # and all unit weights.... 113 infl = N.vstack(( 114 N.hstack(( 115 # upper left 116 k[b[0]:0:-1, b[1]:0:-1], 117 # upper right 118 k[b[0]:0:-1, :self.kshape[1] - b[1]])), 119 N.hstack(( 120 # lower left 121 k[:self.kshape[0] - b[0], b[1]:0:-1], 122 # lower right 123 k[:self.kshape[0] - b[0], :self.kshape[1] - b[1]])) 124 )) 125 unit_deltas += infl[:,:,N.newaxis] * (s - self._K) 126 127 # apply cumulative unit deltas 128 self._K += unit_deltas 129 130 if __debug__: 131 debug("SOM", "Iteration %d/%d done: ||unit_deltas||=%g" % 132 (it, self.niter, N.sqrt(N.sum(unit_deltas **2)))) 133 134 # reset unit deltas 135 unit_deltas.fill(0.)

136 137

138 - def _computeInfluenceKernel(self, iter, dqd):

139 """Compute the neighborhood kernel for some iteration. 140 141 :Parameters: 142 iter: int 143 The iteration for which to compute the kernel. 144 dqd: array (nrows x ncolumns) 145 This is one quadrant of Euclidean distances between Kohonen unit 146 locations. 147 """ 148 # compute radius decay for this iteration 149 curr_max_radius = self.radius * N.exp(-1.0 * iter / self.iter_scale) 150 151 # same for learning rate 152 curr_lrate = self.lrate * N.exp(-1.0 * iter / self.iter_scale) 153 154 # compute Gaussian influence kernel 155 infl = N.exp((-1.0 * dqd) / (2 * curr_max_radius * iter)) 156 infl *= curr_lrate 157 158 # hard-limit kernel to max radius 159 # XXX is this really necessary? 160 infl[dqd > curr_max_radius] = 0. 161 162 return infl

163 164

165 - def _getBMU(self, sample):

166 """Returns the ID of the best matching unit. 167 168 'best' is determined as minimal squared Euclidean distance between 169 any units weight vector and some given target `sample` 170 171 :Parameters: 172 sample: array 173 Target sample. 174 175 :Returns: 176 tuple: (row, column) 177 """ 178 # TODO expose distance function as parameter 179 loc = N.argmin(((self.K - sample) ** 2).sum(axis=2)) 180 181 # assumes 2D Kohonen layer 182 return (N.divide(loc, self.kshape[1]), loc % self.kshape[1])

183 184

185 - def forward(self, data):

186 """Map data from the IN dataspace into OUT space. 187 188 Mapping is performs by simple determining the best matching Kohonen 189 unit for each data sample. 190 """ 191 return N.array([self._getBMU(d) for d in data])

192 193

194 - def reverse(self, data):

195 """Reverse map data from OUT space into the IN space. 196 """ 197 # simple transform into appropriate array slicing and 198 # return the associated Kohonen unit weights 199 return self.K[tuple(N.transpose(data))]

200 201

202 - def getInSize(self):

203 """Returns the size of the entity in input space""" 204 return self.K.shape[-1]

205 206

207 - def getOutSize(self):

208 """Returns the size of the entity in output space""" 209 return self.K.shape[:-1]

210 211

212 - def selectOut(self, outIds):

213 """Limit the OUT space to a certain set of features. 214 215 This is currently not implemented. Moreover, although it is technically 216 possible to implement this functionality, it is unsure whether it is 217 meaningful in the context of SOMs. 218 """ 219 raise NotImplementedError

220 221

222 - def getInId(self, outId):

223 """Translate a feature id into a coordinate/index in input space. 224 225 This is not meaningful in the context of SOMs. 226 """ 227 raise NotImplementedError

228 229

230 - def isValidOutId(self, outId):

231 """Validate feature id in OUT space. 232 """ 233 return (outId >= 0).all() and (outId < self.kshape).all()

234 235

236 - def __repr__(self):

237 s = Mapper.__repr__(self).rstrip(' )') 238 # beautify 239 if not s[-1] == '(': 240 s += ' ' 241 s += 'kshape=%s, niter=%i, learning_rate=%f, iradius=%f)' \ 242 % (str(tuple(self.kshape)), self.niter, self.lrate, 243 self.radius) 244 return s

245 246

247 - def _accessKohonen(self):

248 """Provide access to the Kohonen layer. 249 250 With some care. 251 """ 252 if self._K is None: 253 raise RuntimeError, \ 254 'The SOM needs to be trained before access to the Kohonen ' \ 255 'layer is possible.' 256 257 return self._K

258 259 260 K = property(fget=_accessKohonen)

Source Code for Module mvpa.mappers.som