Package mvpa :: Package misc :: Module data_generators

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.misc.data_generators

      1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
      2  # vi: set ft=python sts=4 ts=4 sw=4 et:
      3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      4  #
      5  #   See COPYING file distributed along with the PyMVPA package for the
      6  #   copyright and license terms.
      7  #
      8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      9  """Miscelaneous data generators for unittests and demos"""
     10
     11  __docformat__ = 'restructuredtext'
     12
     13  import numpy as N
     14
     15  from mvpa.datasets import Dataset
     16
     17  if __debug__:
     18      from mvpa.base import debug
     19
     20 -def multipleChunks(func, n_chunks, *args, **kwargs):
     21      """Replicate datasets multiple times raising different chunks
     22
     23      Given some randomized (noisy) generator of a dataset with a single
     24      chunk call generator multiple times and place results into a
     25      distinct chunks
     26      """
     27      for chunk in xrange(n_chunks):
     28          dataset_ = func(*args, **kwargs)
     29          dataset_.chunks[:] = chunk + 1
     30          if chunk == 0:
     31              dataset = dataset_
     32          else:
     33              dataset += dataset_
     34
     35      return dataset
     36
     37
     38 -def dumbFeatureDataset():
     39      """Create a very simple dataset with 2 features and 3 labels
     40      """
     41      data = [[1, 0], [1, 1], [2, 0], [2, 1], [3, 0], [3, 1], [4, 0], [4, 1],
     42              [5, 0], [5, 1], [6, 0], [6, 1], [7, 0], [7, 1], [8, 0], [8, 1],
     43              [9, 0], [9, 1], [10, 0], [10, 1], [11, 0], [11, 1], [12, 0],
     44              [12, 1]]
     45      regs = ([1] * 8) + ([2] * 8) + ([3] * 8)
     46
     47      return Dataset(samples=data, labels=regs)
     48
     49
     50 -def dumbFeatureBinaryDataset():
     51      """Very simple binary (2 labels) dataset
     52      """
     53      data = [[1, 0], [1, 1], [2, 0], [2, 1], [3, 0], [3, 1], [4, 0], [4, 1],
     54              [5, 0], [5, 1], [6, 0], [6, 1], [7, 0], [7, 1], [8, 0], [8, 1],
     55              [9, 0], [9, 1], [10, 0], [10, 1], [11, 0], [11, 1], [12, 0],
     56              [12, 1]]
     57      regs = ([0] * 12) + ([1] * 12)
     58
     59      return Dataset(samples=data, labels=regs)
     60
     61
     62
     63 -def normalFeatureDataset(perlabel=50, nlabels=2, nfeatures=4, nchunks=5,
     64                           means=None, nonbogus_features=None, snr=3.0):
     65      """Generate a univariate dataset with normal noise and specified means.
     66
     67      :Keywords:
     68        perlabel : int
     69           Number of samples per each label
     70        nlabels : int
     71           Number of labels in the dataset
     72        nfeatures : int
     73           Total number of features (including bogus features which carry
     74           no label-related signal)
     75        nchunks : int
     76           Number of chunks (perlabel should be multiple of nchunks)
     77        means : None or list of float or ndarray
     78           Specified means for each of features among nfeatures.
     79        nonbogus_features : None or list of int
     80           Indexes of non-bogus features (1 per label)
     81        snr : float
     82           Signal-to-noise ration assuming that signal has std 1.0 so we
     83           just divide random normal noise by snr
     84
     85      Probably it is a generalization of pureMultivariateSignal where
     86      means=[ [0,1], [1,0] ]
     87
     88      Specify either means or nonbogus_features so means get assigned
     89      accordingly
     90      """
     91
     92      data = N.random.standard_normal((perlabel*nlabels, nfeatures))/N.sqrt(snr)
     93      if (means is None) and (not nonbogus_features is None):
     94          if len(nonbogus_features) > nlabels:
     95              raise ValueError, "Can't assign simply a feature to a " + \
     96                    "class: more nonbogus_features than labels"
     97          means = N.zeros((len(nonbogus_features), nfeatures))
     98          # pure multivariate -- single bit per feature
     99          for i in xrange(len(nonbogus_features)):
    100              means[i, nonbogus_features[i]] = 1.0
    101      if not means is None:
    102          # add mean
    103          data += N.repeat(N.array(means, ndmin=2), perlabel, axis=0)
    104      # bring it 'under 1', since otherwise some classifiers have difficulties
    105      # during optimization
    106      data = 1.0/(N.max(N.abs(data))) * data
    107      labels = N.concatenate([N.repeat('L%d' % i, perlabel)
    108                                  for i in range(nlabels)])
    109      chunks = N.concatenate([N.repeat(range(nchunks),
    110                                       perlabel/nchunks) for i in range(nlabels)])
    111      ds = Dataset(samples=data, labels=labels, chunks=chunks, labels_map=True)
    112      ds.nonbogus_features = nonbogus_features
    113      return ds
    114
    115 -def pureMultivariateSignal(patterns, signal2noise = 1.5, chunks=None):
    116      """ Create a 2d dataset with a clear multivariate signal, but no
    117      univariate information.
    118
    119      ::
    120
    121        %%%%%%%%%
    122        % O % X %
    123        %%%%%%%%%
    124        % X % O %
    125        %%%%%%%%%
    126      """
    127
    128      # start with noise
    129      data = N.random.normal(size=(4*patterns, 2))
    130
    131      # add signal
    132      data[:2*patterns, 1] += signal2noise
    133
    134      data[2*patterns:4*patterns, 1] -= signal2noise
    135      data[:patterns, 0] -= signal2noise
    136      data[2*patterns:3*patterns, 0] -= signal2noise
    137      data[patterns:2*patterns, 0] += signal2noise
    138      data[3*patterns:4*patterns, 0] += signal2noise
    139
    140      # two conditions
    141      regs = N.array(([0] * patterns) + ([1] * 2 * patterns) + ([0] * patterns))
    142
    143      return Dataset(samples=data, labels=regs, chunks=chunks)
    144
    145
    146 -def normalFeatureDataset__(dataset=None, labels=None, nchunks=None,
    147                             perlabel=50, activation_probability_steps=1,
    148                             randomseed=None, randomvoxels=False):
    149      """ NOT FINISHED """
    150      raise NotImplementedError
    151
    152      if dataset is None and labels is None:
    153          raise ValueError, \
    154                "Provide at least labels or a background dataset"
    155
    156      if dataset is None:
    157          nlabels = len(labels)
    158      else:
    159          nchunks = len(dataset.uniquechunks)
    160
    161      N.random.seed(randomseed)
    162
    163      # Create a sequence of indexes from which to select voxels to be used
    164      # for features
    165      if randomvoxels:
    166          indexes = N.random.permutation(dataset.nfeatures)
    167      else:
    168          indexes = N.arange(dataset.nfeatures)
    169
    170      allind, maps = [], []
    171      if __debug__:
    172          debug('DG', "Ugly creation of the copy of background")
    173
    174      dtype = dataset.samples.dtype
    175      if not N.issubdtype(dtype, N.float):
    176          dtype = N.float
    177      totalsignal = N.zeros(dataset.samples.shape, dtype=dtype)
    178
    179      for l in xrange(len(labels)):
    180          label = labels[l]
    181          if __debug__:
    182              debug('DG', "Simulating independent labels for %s" % label)
    183
    184          # What sample ids belong to this label
    185          labelids = dataset.idsbylabels(label)
    186
    187
    188          # What features belong here and what is left over
    189          nfeatures = perlabel * activation_probability_steps
    190          ind, indexes = indexes[0:nfeatures], \
    191                         indexes[nfeatures+1:]
    192          allind += list(ind)              # store what indexes we used
    193
    194          # Create a dataset only for 'selected' features NB there is
    195          # sideeffect that selectFeatures will sort those ind provided
    196          ds = dataset.selectFeatures(ind)
    197          ds.samples[:] = 0.0             # zero them out
    198
    199          # assign data
    200          prob = [1.0 - x*1.0/activation_probability_steps
    201                  for x in xrange(activation_probability_steps)]
    202
    203          # repeat so each feature gets itw own
    204          probabilities = N.repeat(prob, perlabel)
    205          if __debug__:
    206              debug('DG', 'For prob=%s probabilities=%s' % (prob, probabilities))
    207
    208          for chunk in ds.uniquechunks:
    209              chunkids = ds.idsbychunks(chunk) # samples in this chunk
    210              ids = list(set(chunkids).intersection(set(labelids)))
    211              chunkvalue = N.random.uniform() # random number to decide either
    212                                          # to 'activate' the voxel
    213              for id_ in ids:
    214                  ds.samples[id_, :] = \
    215                                  (chunkvalue <= probabilities).astype('float')
    216
    217          maps.append(N.array(probabilities, copy=True))
    218
    219          signal = ds.map2Nifti(ds.samples)
    220          totalsignal[:, ind] += ds.samples
    221
    222      # figure out average variance across all 'working' features
    223      wfeatures = dataset.samples[:, allind]
    224      meanstd = N.mean(N.std(wfeatures, 1))
    225      if __debug__:
    226          debug('DG', "Mean deviation is %f" % meanstd)
    227
    228      totalsignal *= meanstd * options.snr
    229      # add signal on top of background
    230      dataset.samples += totalsignal
    231
    232      return dataset
    233
    234 -def getMVPattern(s2n):
    235      """Simple multivariate dataset"""
    236      return multipleChunks(pureMultivariateSignal, 6,
    237                            5, s2n, 1)
    238
    239
    240 -def wr1996(size=200):
    241      """Generate '6d robot arm' dataset (Williams and Rasmussen 1996)
    242
    243      Was originally created in order to test the correctness of the
    244      implementation of kernel ARD.  For full details see:
    245      http://www.gaussianprocess.org/gpml/code/matlab/doc/regression.html#ard
    246
    247      x_1 picked randomly in [-1.932, -0.453]
    248      x_2 picked randomly in [0.534, 3.142]
    249      r_1 = 2.0
    250      r_2 = 1.3
    251      f(x_1,x_2) = r_1 cos (x_1) + r_2 cos(x_1 + x_2) + N(0,0.0025)
    252      etc.
    253
    254      Expected relevances:
    255      ell_1      1.804377
    256      ell_2      1.963956
    257      ell_3      8.884361
    258      ell_4     34.417657
    259      ell_5   1081.610451
    260      ell_6    375.445823
    261      sigma_f    2.379139
    262      sigma_n    0.050835
    263      """
    264      intervals = N.array([[-1.932, -0.453], [0.534, 3.142]])
    265      r = N.array([2.0, 1.3])
    266      x = N.random.rand(size, 2)
    267      x *= N.array(intervals[:, 1]-intervals[:, 0])
    268      x += N.array(intervals[:, 0])
    269      if __debug__:
    270          for i in xrange(2):
    271              debug('DG', '%d columnt Min: %g Max: %g' %
    272                    (i, x[:, i].min(), x[:, i].max()))
    273      y = r[0]*N.cos(x[:, 0] + r[1]*N.cos(x.sum(1))) + \
    274          N.random.randn(size)*N.sqrt(0.0025)
    275      y -= y.mean()
    276      x34 = x + N.random.randn(size, 2)*0.02
    277      x56 = N.random.randn(size, 2)
    278      x = N.hstack([x, x34, x56])
    279      return Dataset(samples=x, labels=y)
    280
    281
    282 -def sinModulated(n_instances, n_features,
    283                    flat=False, noise=0.4):
    284      """ Generate a (quite) complex multidimensional non-linear dataset
    285
    286      Used for regression testing. In the data label is a sin of a x^2 +
    287      uniform noise
    288      """
    289      if flat:
    290          data = (N.arange(0.0, 1.0, 1.0/n_instances)*N.pi)
    291          data.resize(n_instances, n_features)
    292      else:
    293          data = N.random.rand(n_instances, n_features)*N.pi
    294      label = N.sin((data**2).sum(1)).round()
    295      label += N.random.rand(label.size)*noise
    296      return Dataset(samples=data, labels=label)
    297
    298 -def chirpLinear(n_instances, n_features=4, n_nonbogus_features=2,
    299                  data_noise=0.4, noise=0.1):
    300      """ Generates simple dataset for linear regressions
    301
    302      Generates chirp signal, populates n_nonbogus_features out of
    303      n_features with it with different noise level and then provides
    304      signal itself with additional noise as labels
    305      """
    306      x = N.linspace(0, 1, n_instances)
    307      y = N.sin((10 * N.pi * x **2))
    308
    309      data = N.random.normal(size=(n_instances, n_features ))*data_noise
    310      for i in xrange(n_nonbogus_features):
    311          data[:, i] += y[:]
    312
    313      labels = y + N.random.normal(size=(n_instances,))*noise
    314
    315      return Dataset(samples=data, labels=labels)
    316
    317
    318 -def linear_awgn(size=10, intercept=0.0, slope=0.4, noise_std=0.01, flat=False):
    319      """Generate a dataset from a linear function with AWGN
    320      (Added White Gaussian Noise).
    321
    322      It can be multidimensional if 'slope' is a vector. If flat is True
    323      (in 1 dimesion) generate equally spaces samples instead of random
    324      ones. This is useful for the test phase.
    325      """
    326      dimensions = 1
    327      if isinstance(slope, N.ndarray):
    328          dimensions = slope.size
    329
    330      if flat and dimensions == 1:
    331          x = N.linspace(0, 1, size)[:, N.newaxis]
    332      else:
    333          x = N.random.rand(size, dimensions)
    334
    335      y = N.dot(x, slope)[:, N.newaxis] \
    336          + (N.random.randn(*(x.shape[0], 1)) * noise_std) + intercept
    337
    338      return Dataset(samples=x, labels=y)
    339
    340
    341 -def noisy_2d_fx(size_per_fx, dfx, sfx, center, noise_std=1):
    342      """
    343      """
    344      x = []
    345      y = []
    346      labels = []
    347      for fx in sfx:
    348          nx = N.random.normal(size=size_per_fx)
    349          ny = fx(nx) + N.random.normal(size=nx.shape, scale=noise_std)
    350          x.append(nx)
    351          y.append(ny)
    352
    353          # whenever larger than first function value
    354          labels.append(N.array(ny < dfx(nx), dtype='int'))
    355
    356      samples = N.array((N.hstack(x), N.hstack(y))).squeeze().T
    357      labels = N.hstack(labels).squeeze().T
    358
    359      samples += N.array(center)
    360
    361      return Dataset(samples=samples, labels=labels)
    362
    363
    364 -def linear1d_gaussian_noise(size=100, slope=0.5, intercept=1.0,
    365                              x_min=-2.0, x_max=3.0, sigma=0.2):
    366      """A straight line with some Gaussian noise.
    367      """
    368      x = N.linspace(start=x_min, stop=x_max, num=size)
    369      noise = N.random.randn(size)*sigma
    370      y = x * slope + intercept + noise
    371      return Dataset(samples=x.reshape((-1, 1)), labels=y)
    372