parsimony.utils.start

1 # -*- coding: utf-8 -*- 2 """ 3 Created on Thu Mar 28 15:35:26 2013 4 5 Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved. 6 7 @author: Tommy Löfstedt 8 @email: tommy.loefstedt@cea.fr 9 @license: BSD 3-clause. 10 """ 11 import abc 12 import numpy as np 13 14 from . import maths 15 16 __all__ = ['BaseStartVector', 'IdentityStartVector', 'RandomStartVector', 17 'OnesStartVector', 'ZerosStartVector']

18 19 20 -class BaseStartVector(object):

21 """Base class for start vector generation. 22 23 Parameters 24 ---------- 25 normalise : Bool. Whether or not to normalise the vector that is returned. 26 27 seed : Integer or None. The seed to the pseudo-random number generator. If 28 none, no seed is used. The seed is set at initialisation, so if the 29 RNG is used in between initialisation and utilisation, then the 30 random numbers will change. Default is None. The seed is not used 31 by all implementing classes. 32 """ 33 __metaclass__ = abc.ABCMeta 34

35 - def __init__(self, normalise=True, seed=None):

36 super(BaseStartVector, self).__init__() 37 38 self.normalise = normalise 39 self.seed = seed 40 if seed is not None: 41 np.random.seed(seed)

42 43 @abc.abstractmethod

44 - def get_vector(self, size):

45 46 raise NotImplementedError('Abstract method "get_vector" must be ' 47 'specialised!')

48

49 50 -class IdentityStartVector(BaseStartVector):

51 """A pre-determined start vector. 52 53 Parameters 54 ---------- 55 vector : Numpy array. The predetermined start vector 56 57 Examples 58 -------- 59 >>> from parsimony.utils.start_vectors import IdentityStartVector 60 >>> start_vector = IdentityStartVector(np.array([[0.5], [2.0], [0.3], 61 ... [1.0]])) 62 >>> start_vector.get_vector() 63 array([[ 0.5], 64 [ 2. ], 65 [ 0.3], 66 [ 1. ]]) 67 """

68 - def __init__(self, vector, **kwargs):

69 70 super(IdentityStartVector, self).__init__(**kwargs) 71 72 self.vector = vector

73

74 - def get_vector(self, *args, **kwargs):

75 """Return the predetermined start vector 76 77 Examples 78 -------- 79 >>> from parsimony.utils.start_vectors import IdentityStartVector 80 >>> start_vector = IdentityStartVector(np.array([[0.5], [2.0], [0.3], 81 ... [1.0]])) 82 >>> start_vector.get_vector() 83 array([[ 0.5], 84 [ 2. ], 85 [ 0.3], 86 [ 1. ]]) 87 """ 88 return self.vector

89

90 91 -class RandomStartVector(BaseStartVector):

92 """A start vector of uniformly distributed random values. 93 94 Parameters 95 ---------- 96 normalise : Bool. If True, normalise the randomly created vectors. 97 Default is True. 98 99 seed : Integer or None. The seed to the pseudo-random number generator. If 100 none, no seed is used. The seed is set at initialisation, so if the 101 RNG is used in between initialisation and utilisation, then the 102 random numbers will change. Default is None. 103 104 limits : List or tuple. A list or tuple with two elements, the lower and 105 upper limits of the uniform distribution. If normalise=True, then 106 these limits may not be honoured. Default is (0.0, 1.0). 107 108 Examples 109 -------- 110 >>> from parsimony.utils.start_vectors import RandomStartVector 111 >>> 112 >>> # Without normalization 113 >>> start_vector = RandomStartVector(normalise=False, seed=42) 114 >>> random = start_vector.get_vector(3) 115 >>> print random 116 [[ 0.37454012] 117 [ 0.95071431] 118 [ 0.73199394]] 119 >>> print maths.norm(random) 120 1.25696186254 121 >>> 122 >>> # With normalization 123 >>> start_vector_normalized = RandomStartVector(normalise=True, seed=2) 124 >>> random_normalized = start_vector_normalized.get_vector(3) 125 >>> print random_normalized 126 [[ 0.62101956] 127 [ 0.03692864] 128 [ 0.78292463]] 129 >>> print maths.norm(random_normalized) 130 1.0 131 >>> 132 >>> # With limits 133 >>> start_vector_normalized = RandomStartVector(normalise=True, seed=2, 134 ... limits=(-1, 1)) 135 >>> random_limits = start_vector_normalized.get_vector(3) 136 >>> print random_limits 137 [[-0.1330817 ] 138 [-0.98571123] 139 [ 0.10326001]] 140 >>> print maths.norm(random_limits) 141 1.0 142 143 """

144 - def __init__(self, limits=(0.0, 1.0), **kwargs):

145 146 super(RandomStartVector, self).__init__(**kwargs) 147 148 self.limits = limits

149

150 - def get_vector(self, size):

151 """Return randomly generated vector of given shape. 152 153 Parameters 154 ---------- 155 size : Positive integer. Size of the vector to generate. The shape of 156 the output is (size, 1). 157 158 Examples 159 -------- 160 >>> from parsimony.utils.start_vectors import RandomStartVector 161 >>> start_vector = RandomStartVector(normalise=False, seed=42) 162 >>> random = start_vector.get_vector(3) 163 >>> print random 164 [[ 0.37454012] 165 [ 0.95071431] 166 [ 0.73199394]] 167 >>> 168 >>> start_vector = RandomStartVector(normalise=False, seed=1, 169 ... limits=(-1, 2)) 170 >>> random = start_vector.get_vector(3) 171 >>> print random 172 [[ 0.25106601] 173 [ 1.16097348] 174 [-0.99965688]] 175 """ 176 l = float(self.limits[0]) 177 u = float(self.limits[1]) 178 179 size = int(size) 180 vector = np.random.rand(size, 1) * (u - l) + l # Random vector. 181 182 if self.normalise: 183 return vector * (1.0 / maths.norm(vector)) 184 else: 185 return vector

186

187 188 -class OnesStartVector(BaseStartVector):

189 """A start vector of ones. 190 191 Parameters 192 ---------- 193 normalise : Bool. If True, normalise the randomly created vectors 194 Default is False 195 196 Examples 197 -------- 198 >>> from parsimony.utils.start_vectors import OnesStartVector 199 200 # Without normalization 201 >>> start_vector = OnesStartVector(normalise=False) 202 >>> ones = start_vector.get_vector(3) 203 >>> print ones 204 [[ 1.] 205 [ 1.] 206 [ 1.]] 207 >>> print maths.norm(ones) 208 1.73205080757 209 210 # With normalization 211 >>> start_vector_normalized = OnesStartVector(normalise=True) 212 >>> ones_normalized = start_vector_normalized.get_vector(3) 213 >>> print ones_normalized 214 [[ 0.57735027] 215 [ 0.57735027] 216 [ 0.57735027]] 217 >>> print maths.norm(ones_normalized) 218 1.0 219 """

220 - def __init__(self, normalise=False, **kwargs):

221 222 super(OnesStartVector, self).__init__(normalise=normalise, **kwargs)

223

224 - def get_vector(self, size):

225 """Return vector of ones of chosen shape 226 227 Parameters 228 ---------- 229 size : Positive integer. Size of the vector to generate. The shape of 230 the output is (size, 1). 231 232 Examples 233 -------- 234 >>> from parsimony.utils.start_vectors import OnesStartVector 235 >>> start_vector = OnesStartVector() 236 >>> ones = start_vector.get_vector(3) 237 >>> print ones 238 [[ 1.] 239 [ 1.] 240 [ 1.]] 241 """ 242 size = int(size) 243 vector = np.ones((size, 1)) # Using a vector of ones. 244 245 if self.normalise: 246 return vector * (1.0 / maths.norm(vector)) 247 else: 248 return vector

249

250 251 -class ZerosStartVector(BaseStartVector):

252 """A start vector of zeros. 253 254 Use with care! Be aware that using this in algorithms that are not aware 255 may result in division by zero since the norm of this start vector is 0. 256 257 Examples 258 -------- 259 >>> from parsimony.utils.start_vectors import ZerosStartVector 260 >>> start_vector = ZerosStartVector() 261 >>> zeros = start_vector.get_vector(3) 262 >>> print zeros 263 [[ 0.] 264 [ 0.] 265 [ 0.]] 266 """

267 - def __init__(self, **kwargs):

268 269 kwargs.pop('normalise', False) # We do not care about this argument. 270 271 super(ZerosStartVector, self).__init__(normalise=False, **kwargs)

272

273 - def get_vector(self, size):

274 """Return vector of zeros of chosen shape. 275 276 Parameters 277 ---------- 278 size : Positive integer. Size of the vector to generate. The shape of 279 the output is (size, 1). 280 281 Examples 282 -------- 283 >>> from parsimony.utils.start_vectors import ZerosStartVector 284 >>> start_vector = ZerosStartVector() 285 >>> zeros = start_vector.get_vector(3) 286 >>> print zeros 287 [[ 0.] 288 [ 0.] 289 [ 0.]] 290 """ 291 size = int(size) 292 w = np.zeros((size, 1)) # Using a vector of zeros. 293 294 return w

295 296 297 #class LargestStartVector(BaseStartVector): 298 # 299 # def __init__(self, normalise=True, **kwargs): 300 # 301 # super(LargestStartVector, self).__init__(normalise=normalise, **kwargs) 302 # 303 # def get_vector(self, X, axis=1): 304 # if X == None: 305 # raise ValueError('A matrix X must be must be given.') 306 # 307 # idx = np.argmax(np.sum(X ** 2.0, axis=axis)) 308 # if axis == 0: 309 # w = X[:, [idx]] # Using column with largest sum of squares 310 # else: 311 # w = X[[idx], :].T # Using row with largest sum of squares 312 # 313 # if self.normalise: 314 # return w * (1.0 / norm(w)) 315 # else: 316 # return w 317 318 319 #class GaussianCurveVector(BaseStartVector): 320 # """A start vector with the shape of a Gaussian curve. 321 # 322 # The gaussian is computed with respect to the numbers of dimension in a 323 # supposed image. The output is thus a reshaped vector corresponsing to a 1-, 324 # 2-, 3- or higher-dimensional Gaussian curve. 325 # """ 326 # 327 # def __init__(self, **kwargs): 328 # 329 # super(GaussianCurveVector, self).__init__(**kwargs) 330 # 331 # def get_vector(self, shape=None, size=None, mean=None, cov=None, dims=2): 332 # """ Computes a Gaussian curve-shaped starting vector. 333 # 334 # Parameters: 335 # shape : A tuple. The shape of the start vector. 336 # 337 # size : A tuple. The size of the supposed image. Must have the form (Z, 338 # Y, X). 339 # 340 # mean : A numpy array. The mean vector of the Gaussian. Default is zero. 341 # 342 # cov : A numpy array. The covariance matrix of the Gaussian. Default is 343 # the identity. 344 # 345 # dims : A scalar. The number of dimensions of the output image. Default 346 # is 2. 347 # """ 348 # if size != None: 349 # p = 1 350 # for i in xrange(dims): 351 # p *= size[i] 352 # if axis == 1: 353 # shape = (p, 1) 354 # else: 355 # shape = (1, p) 356 # else: 357 # if X != None: 358 # p = X.shape[axis] 359 # shape = (p, 1) 360 # else: # Assumes shape != None 361 # p = shape[0] * shape[1] 362 # 363 # size = [0] * dims 364 # for i in xrange(dims): # Split in equal-sized hypercube 365 # size[i] = round(float(p) ** (1.0 / float(dims))) 366 # 367 # if mean == None: 368 # mean = [float(s - 1.0) / 2.0 for s in size] 369 # if cov == None: 370 # S = np.diag([s ** (1.0 / dims) for s in size]) 371 # invS = np.linalg.pinv(S) 372 # else: 373 ## S = np.diag(np.diag(cov)) 374 # S = np.asarray(cov) 375 # invS = np.linalg.pinv(S) 376 # 377 # a = np.arange(size[0]) 378 # ans = np.reshape(a, (a.shape[0], 1)).tolist() 379 # for i in xrange(1, dims): 380 # b = np.arange(size[i]).tolist() 381 # ans = [y + [x] for x in b for y in ans] 382 # 383 # X = np.zeros((size)) 384 # for x in ans: 385 # i = tuple(x) 386 # x = np.array([x]) - np.array(mean) 387 # v = np.dot(x, np.dot(invS, x.T)) 388 # X[i] = v[0, 0] 389 # 390 # X = np.exp(-0.5 * X) 391 # X *= (1.0 / np.sum(X)) 392 # 393 ## s = [] 394 ## X = 0 395 ## for i in xrange(dims): 396 ## x = np.arange(size[i]) - mean[i] 397 ## x = np.reshape(x, [size[i]] + s) 398 ## X = X + invS[i, i] * (x ** 2.0) 399 ## s.append(1) 400 # 401 # w = np.reshape(X, (p, 1)) 402 # 403 # if self.normalise: 404 # return w * (1.0 / norm(w)) 405 # else: 406 # return w 407 # 408 # 409 #class GaussianCurveVectors(BaseStartVector): 410 # """A start vector with multibple Gaussian curve shapes. 411 # 412 # The gaussians are in an imagined 1D or 2D image. The output is a reshaped 413 # vector corresponsing to a 1- or 2-dimensional image. 414 # """ 415 # 416 # def __init__(self, num_points=3, normalise=True, **kwargs): 417 # super(GaussianCurveVectors, self).__init__(normalise=normalise, 418 # **kwargs) 419 # 420 # self.num_points = num_points 421 # 422 # def get_vector(self, X=None, axis=1, shape=None, size=None, 423 # mean=None, cov=None, dims=2): 424 # """ Computes a starting vector with set of Gaussian curve-shapes. 425 # 426 # Parameters: 427 # X : The matrix for which we need a start vector. Used in 428 # conjunction with axis to determine the shape of the start 429 # vector. 430 # 431 # axis : The axis along X which the shape is taken. 432 # 433 # shape : The shape of the start vector, may be passed instead of X. 434 # 435 # size : The size of the supposed image. Must have the form (Z, Y, X). 436 # May be passed instead of X or shape. 437 # 438 # means : The mean vectors of the Gaussians. Default is random. 439 # 440 # covs : The covariance matrices of the Gaussians. Default is random. 441 # 442 # dims : The number of dimensions of the output image. Default is 2. 443 # """ 444 # if size != None: 445 # p = 1 446 # for i in xrange(dims): 447 # p *= size[i] 448 # if axis == 1: 449 # shape = (p, 1) 450 # else: 451 # shape = (1, p) 452 # else: 453 # if X != None: 454 # p = X.shape[axis] 455 # shape = (p, 1) 456 # else: # Assumes shape != None 457 # p = shape[0] * shape[1] 458 # 459 # size = [0] * dims 460 # for i in xrange(dims): # Split in equal-sized hypercube 461 # size[i] = round(float(p) ** (1.0 / float(dims))) 462 # 463 # means = np.random.rand(1, 2) 464 # for i in xrange(1, self.num_points): 465 # dist = 0.0 466 # p_best = 0 467 # for j in xrange(20): 468 # p = np.random.rand(1, 2) 469 # dist_curr = np.min(np.sqrt(np.sum((means - p) ** 2.0, axis=1))) 470 # if dist_curr > dist: 471 # p_best = p 472 # dist = dist_curr 473 # if dist_curr > 0.3: 474 # break 475 # means = np.vstack((means, p_best)) 476 # 477 # means[means < 0.05] = 0.05 478 # means[means > 0.95] = 0.95 479 # means[:, 0] *= size[0] 480 # means[:, 1] *= size[1] 481 # means = means.tolist() 482 # 483 # covs = [0] * self.num_points 484 # for i in xrange(self.num_points): 485 # S1 = np.diag((np.abs(np.diag(np.random.rand(2, 2))) * 0.5) + 0.5) 486 # 487 # S2 = np.random.rand(2, 2) 488 # S2 = (((S2 + S2.T) * 0.5) - 0.5) * 0.9 # [0, 0.45] 489 # S2 = S2 - np.diag(np.diag(S2)) 490 # 491 # S = S1 + S2 492 # 493 # S *= 1.0 / np.max(S) 494 # 495 # S *= float(min(size)) 496 # 497 # covs[i] = S.tolist() 498 # 499 # vector = GaussianCurveVector(normalise=False) 500 # 501 # X = np.zeros(shape) 502 # for i in xrange(self.num_points): 503 # X = X + vector.get_vector(size=size, dims=dims, 504 # mean=means[i], cov=covs[i]) 505 # 506 # w = np.reshape(X, size) 507 # 508 # if self.normalise: 509 # return w * (1.0 / norm(w)) 510 # else: 511 # return w 512 513 514 if __name__ == "__main__": 515 import doctest 516 doctest.testmod() 517

Source Code for Module parsimony.utils.start_vectors