Introduction

Modular toolkit for Data Processing (MDP) is a data processing framework written in Python.

From the user's perspective, MDP consists of a collection of trainable supervised and unsupervised algorithms or other data processing units (nodes) that can be combined into data processing flows and more complex feed-forward network architectures. Given a sequence of input data, MDP takes care of successively training or executing all nodes in the network. This structure allows to specify complex algorithms as a sequence of simpler data processing steps in a natural way. Training can be performed using small chunks of input data, so that the use of very large data sets becomes possible while reducing the memory requirements. Memory usage can also be minimized by defining the internals of the nodes to be single precision.

The base of readily available algorithms includes Principal Component Analysis (PCA and NIPALS), four flavors of Independent Component Analysis (CuBICA, FastICA, TDSEP, and JADE), Slow Feature Analysis, Independent Slow Feature Analysis, Gaussian Classifiers, Growing Neural Gas, Fisher Discriminant Analysis, Factor Analysis, Restricted Boltzmann Machine, and many more. The full list of implemented nodes can be found in the Node List section.

From the developer's perspective, MDP is a framework to make the implementation of new supervised and unsupervised algorithms easier. The basic class Node takes care of tedious tasks like numerical type and dimensionality checking, leaving the developer free to concentrate on the implementation of the training and execution phases. The node then automatically integrates with the rest of the library and can be used in a flow together with other nodes. A node can have multiple training phases and even an undetermined number of phases. This allows for example the implementation of algorithms that need to collect some statistics on the whole input before proceeding with the actual training, or others that need to iterate over a training phase until a convergence criterion is satisfied. The ability to train each phase using chunks of input data is maintained if the chunks are generated with iterators. Moreover, crash recovery is optionally available: in case of failure, the current state of the flow is saved for later inspection.

MDP has been written in the context of theoretical research in neuroscience, but it has been designed to be helpful in any context where trainable data processing algorithms are used. Its simplicity on the user side together with the reusability of the implemented nodes make it also a valid educational tool.

As its user and contributor base is steadily increasing, MDP appears as a good candidate for becoming a common repository of user-supplied, freely available, Python implemented data processing algorithms.

Quick Start

Using MDP is as easy as:

>>> import mdp
>>> # perform pca on some data x
...
>>> y = mdp.pca(x)
>>> # perform ica on some data x using single precision
...
>>> y = mdp.fastica(x, dtype='float32')

A complete list of all short-cut functions like pca or fastica can be obtained as follows:

>>> dir(mdp.helper_funcs)
['__builtins__', '__doc__', '__file__', '__name__',
'cubica', 'factor_analysis', 'fastica', 'get_eta',
'isfa', 'mdp', 'pca', 'sfa', 'sfa2', 'whitening']

MDP is of course much more than this: it allows to combine different algorithms and other data processing elements (nodes) into data processing sequences (flows), and more general feed-forward architectures (with the new hinet subpackage). Moreover, it provides a framework that makes the implementation of new algorithms easy and intuitive.

MDP requires the numerical Python extensions numpy or scipy. In its namespace MDP offers references to the main modules numpy or scipy, and the subpackages linalg, random, and fft as mdp.numx, mdp.numx_linalg, mdp.numx_rand, and mdp.numx_fft. This is done to possibly support additional numerical extensions in the future. At import time MDP will select scipy if available, otherwise numpy will be loaded. You can force the use of a numerical extension by setting the environment variable MDPNUMX=numpy or MDPNUMX=scipy.

Nodes

A node is the basic unit in MDP and it represents a data processing element, like for example a learning algorithm, a filter, a visualization step, etc. Each node can have one or more training phases, during which the internal structures are learned from training data (e.g. the weights of a neural network are adapted or the covariance matrix is estimated) and an execution phase, where new data can be processed forwards (by processing the data through the node) or backwards (by applying the inverse of the transformation computed by the node if defined). The Node class is designed to make the implementation of new algorithms easy and intuitive, for example by setting automatically input and output dimension and by casting the data to match the numerical type (e.g. float or double) of the internal structures. Node was designed to be applied to arbitrarily long sets of data: the internal structures can be updated incrementally by sending chunks of the input data (this is equivalent to online learning if the chunks consists of single observations, or to batch learning if the whole data is sent in a single chunk). A Node can be copied or saved using the corresponding copy and save methods.

Flows

A flow consists in an acyclic graph of nodes (currently only node sequences are implemented). The data is sent to an input node and is successively processed by the following nodes on the graph. The general flow implementation automatizes the training, execution, and inverse execution (if defined) of the whole graph. Training can be supervised and can consist of multiple phases. Crash recovery is optionally available: in case of failure the current state of the flow is saved for later inspection. A subclass of the basic flow class (CheckpointFlow) allows user-supplied checkpoint functions to be executed at the end of each phase, for example to save the internal structures of a node for later analysis. Flow objects are Python containers. Most of the builtin list methods are available. A Flow can be saved or copied using the corresponding save and copy methods.

>>> inp = mdp.numx_rand.random((1000, 20))

>>> inp = (inp - mdp.numx.mean(inp, 0))/mdp.numx.std(inp, 0)

>>> inp[:,5:] /= 10.0

>>> x = mdp.utils.mult(inp,mdp.numx_rand.random((20, 20)))

>>> inp_test = mdp.numx_rand.random((1000, 20))
>>> inp_test = (inp_test - mdp.numx.mean(inp_test, 0))/mdp.numx.std(inp_test, 0)
>>> inp_test[:,5:] /= 10.0
>>> x_test = mdp.utils.mult(inp_test, mdp.numx_rand.random((20, 20)))

>>> rec = flow.inverse(out)

>>> cov = mdp.numx.amax(abs(mdp.utils.cov2(x/mdp.numx.std(x,axis=0),
...                                        rec/mdp.numx.std(rec,axis=0))))
>>> print cov
[ 0.99839606  0.99744461  0.99616208  0.99772863  0.99690947
  0.99864056  0.99734378  0.98722502  0.98118101  0.99407939
  0.99683096  0.99756988  0.99664384  0.99723419  0.9985529
  0.99829763  0.9982712   0.99721741  0.99682906  0.98858858]

>>> for node in flow:
...     print repr(node)
...
PCANode(input_dim=20, output_dim=5, dtype='float64')
CuBICANode(input_dim=5, output_dim=5, dtype='float64')
HitParadeNode(input_dim=5, output_dim=5, dtype='float64')
>>>

>>> len(flow)
3
>>> print flow[::2]
[PCANode, HitParadeNode]
>>> nodetoberemoved = flow.pop(-1)
>>> nodetoberemoved
HitParadeNode(input_dim=5, output_dim=5, dtype='float64')
>>> len(flow)
2

>>> dummyflow = flow[1:].copy()
>>> longflow = flow + dummyflow
>>> len(longflow)
3

>>> class BogusExceptNode(mdp.Node):
...    def train(self,x):
...        self.bogus_attr = 1
...        raise Exception, "Bogus Exception"
...    def execute(self,x):
...        raise Exception, "Bogus Exception"
...
>>> flow = mdp.Flow([BogusExceptNode()])

>>> flow.set_crash_recovery(1)

>>> flow.train(x)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  [...]
mdp.linear_flows.FlowExceptionCR:
----------------------------------------
! Exception in node #0 (BogusExceptNode):
Node Traceback:
Traceback (most recent call last):
  [...]
Exception: Bogus Exception
----------------------------------------
A crash dump is available on: "/tmp/MDPcrash_LmISO_.pic"

>>> flow.set_crash_recovery('/home/myself/mydumps/MDPdump.pic')

Iterators

Python allows user-defined classes to support iteration, as described in the Python docs. A convenient implementation of the iterator protocol is provided by generators: see this article for an introduction, and the official PEP for a complete description.

Let us define two bogus node classes to be used as examples of nodes:

>>> class BogusNode(mdp.Node):
...     """This node does nothing."""
...     def _train(self, x):
...         pass
...
>>> class BogusNode2(mdp.Node):
...     """This node does nothing. But it's not trainable nor invertible.
...     """
...     def is_trainable(self): return False
...     def is_invertible(self): return False
...
>>>

This generator generates blocks input blocks to be used as training set. In this example one block is a 2-dimensional time-series. The first variable is [2,4,6,....,1000] and the second one [0,1,3,5,...,999]. All blocks are equal, this of course would not be the case in a real-life example.

In this example we use a progress bar to get progress information.

>>> def gen_data(blocks):
...     for i in mdp.utils.progressinfo(xrange(blocks)):
...         block_x = mdp.numx.atleast_2d(mdp.numx.arange(2,1001,2))
...         block_y = mdp.numx.atleast_2d(mdp.numx.arange(1,1001,2))
...         # put variables on columns and observations on rows
...         block = mdp.numx.transpose(mdp.numx.concatenate([block_x,block_y]))
...         yield block
...
>>>

Let's define a bogus flow consisting of 2 BogusNode:

>>> flow = mdp.Flow([BogusNode(),BogusNode()], verbose=1)

Train the first node with 5000 blocks and the second node with 3000 blocks. Note that the only allowed argument to train is a sequence (list or tuple) of iterators. In case you don't want or need to use incremental learning and want to do a one-shot training, you can use as argument to train a single array of data:

block-mode training

>>> flow.train([gen_data(5000),gen_data(3000)])
Training node #0 (BogusNode)
[===================================100%==================================>]

Training finished
Training node #1 (BogusNode)
[===================================100%==================================>]

Training finished
Close the training phase of the last node

one-shot training using one single set of data for both nodes

>>> flow = mdp.Flow([BogusNode(),BogusNode()])
>>> block_x = mdp.numx.atleast_2d(mdp.numx.arange(2,1001,2))
>>> block_y = mdp.numx.atleast_2d(mdp.numx.arange(1,1001,2))
>>> single_block = mdp.numx.transpose(mdp.numx.concatenate([block_x,block_y]))
>>> flow.train(single_block)

If your flow contains non-trainable nodes, you must specify a None iterator for the non-trainable nodes:

>>> flow = mdp.Flow([BogusNode2(),BogusNode()], verbose=1)
>>> flow.train([None, gen_data(5000)])
Training node #0 (BogusNode2)
Training finished
Training node #1 (BogusNode)
[===================================100%==================================>]

Training finished
Close the training phase of the last node

You can use the one-shot training:

>>> flow = mdp.Flow([BogusNode2(),BogusNode()], verbose=1)
>>> flow.train(single_block)
Training node #0 (BogusNode2)
Training finished
Training node #1 (BogusNode)
Training finished
Close the training phase of the last node

Iterators can be used also for execution (and inversion):

>>> flow = mdp.Flow([BogusNode(),BogusNode()], verbose=1)
>>> flow.train([gen_data(1), gen_data(1)])
Training node #0 (BogusNode2)
Training finished
Training node #1 (IdentityNode)
[===================================100%==================================>]

Training finished
Close the training phase of the last node
>>> output = flow.execute(gen_data(1000))
[===================================100%==================================>]
>>> output = flow.inverse(gen_data(1000))
[===================================100%==================================>]

Execution and inversion can be done in one-shot mode also. Note that since training is finished you are not going to get a warning

>>> output = flow.execute(single_block)
>>> output = flow.inverse(single_block)

If a node requires multiple training phases (e.g., GaussianClassifierNode), Flow automatically takes care of reusing the iterator multiple times. In this case generators are not allowed, since they expire after yielding the last data block. If you try to restart them, they raise a StopIteration exception. General iterators, instead, can always be restarted. For example, you can loop over a list as many times as you need.

However, it is fairly easy to wrap a generator in a simple iterator if you need to:

>>> class SimpleIterator(object):
...     def __init__(self, blocks):
...         self.blocks = blocks
...     def __iter__(self):
...         # this is a generator
...         for i in range(self.blocks):
...             yield generate_some_data()
>>>

Note that if you use random numbers within the iterator, you usually would like to reset the random number generator to produce the same sequence every time:

>>> class RandomIterator(object):
...     def __init__(self):
...         self.state = None
...     def __iter__(self):
...         if self.state is None:
...             self.state = mdp.numx_rand.get_state()
...         else:
...             mdp.numx_rand.set_state(self.state)
...         for i in range(2):
...             yield mdp.numx_rand.random((1,4))
>>> iterator = RandomIterator()
>>> for x in iterator: print x
...
[[ 0.99586495  0.53463386  0.6306412   0.09679571]]
[[ 0.51117469  0.46647448  0.95089738  0.94837122]]
>>> for x in iterator: print x
...
[[ 0.99586495  0.53463386  0.6306412   0.09679571]]
[[ 0.51117469  0.46647448  0.95089738  0.94837122]]

Checkpoints

It can sometimes be useful to execute arbitrary functions at the end of the training or execution phase, for example to save the internal structures of a node for later analysis. This can easily be done by defining a CheckpointFlow. As an example imagine the following situation: you want to perform Principal Component Analysis (PCA) on your data to reduce the dimensionality. After this you want to expand the signals into a nonlinear space and then perform Slow Feature Analysis to extract slowly varying signals. As the expansion will increase the number of components, you don't want to run out of memory, but at the same time you want to keep as much information as possible after the dimensionality reduction. You could do that by specifying the percentage of the total input variance that has to be conserved in the dimensionality reduction. As the number of output components of the PCA node now can become as large as the that of the input components, you want to check, after training the PCA node, that this number is below a certain threshold. If this is not the case you want to abort the execution and maybe start again requesting less variance to be kept.

Let start defining a generator to be used through the whole example:

>>> def gen_data(blocks,dims):
...     mat = mdp.numx_rand.random((dims,dims))-0.5
...     for i in xrange(blocks):
...         # put variables on columns and observations on rows
...         block = mdp.utils.mult(mdp.numx_rand.random((1000,dims)), mat)
...         yield block
...
>>>

Define a PCANode which reduces dimensionality of the input, a PolynomialExpansionNode to expand the signals in the space of polynomials of degree 2 and a SFANode to perform SFA:

>>> pca = mdp.nodes.PCANode(output_dim=0.9)
>>> exp = mdp.nodes.PolynomialExpansionNode(2)
>>> sfa = mdp.nodes.SFANode()

As you see we have set the output dimension of the PCANode to be 0.9. This means that we want to keep at least 90% of the variance of the original signal. We define a PCADimensionExceededException that has to be thrown when the number of output components exceeds a certain threshold:

>>> class PCADimensionExceededException(Exception):
...     """Exception base class for PCA exceeded dimensions case."""
...     pass
...
>>>

Then, write a CheckpointFunction that checks the number of output dimensions of the PCANode and aborts if this number is larger than max_dim:

>>> class CheckPCA(mdp.CheckpointFunction):
...     def __init__(self,max_dim):
...         self.max_dim = max_dim
...     def __call__(self,node):
...         node.stop_training()
...         act_dim = node.get_output_dim()
...         if act_dim > self.max_dim:
...             errstr = 'PCA output dimensions exceeded maximum '+\
...                      '(%d > %d)'%(act_dim,self.max_dim)
...             raise PCADimensionExceededException, errstr
...         else:
...             print 'PCA output dimensions = %d'%(act_dim)
...
>>>

Define the CheckpointFlow:

>>> flow = mdp.CheckpointFlow([pca, exp, sfa])

To train it we have to supply 3 generators and 3 checkpoint functions:

>>> flow.train([gen_data(10, 50), None, gen_data(10, 50)],
...            [CheckPCA(10), None, None])
Traceback (most recent call last):
  File "<stdin>", line 2, in ?
  [...]
__main__.PCADimensionExceededException: PCA output dimensions exceeded maximum (25 > 10)

The training fails with a PCADimensionExceededException. If we only had 12 input dimensions instead of 50 we would have passed the checkpoint:

>>> flow[0] = mdp.nodes.PCANode(output_dim=0.9)
>>> flow.train([gen_data(10, 12), None, gen_data(10, 12)],
...            [CheckPCA(10), None, None])
PCA output dimensions = 6

We could use the built-in CheckpoinSaveFunction to save the SFANode and analyze the results later :

>>> pca = mdp.nodes.PCANode(output_dim=0.9)
>>> exp = mdp.nodes.PolynomialExpansionNode(2)
>>> sfa = mdp.nodes.SFANode()
>>> flow = mdp.CheckpointFlow([pca, exp, sfa])
>>> flow.train([gen_data(10, 12), None, gen_data(10, 12)],
...            [CheckPCA(10),
...             None,
...             mdp.CheckpointSaveFunction('dummy.pic',
...                                        stop_training = 1,
...                                        protocol = 0)])
...
PCA output dimensions = 7

We can now reload and analyze the SFANode:

>>> fl = file('dummy.pic')
>>> import cPickle
>>> sfa_reloaded = cPickle.load(fl)
>>> sfa_reloaded
SFANode(input_dim=35, output_dim=35, dtype='d')

Don't forget to clean the rubbish:

>>> fl.close()
>>> import os
>>> os.remove('dummy.pic')

Hierarchical Networks

The hinet package makes it possible to construct graph-like Node structures, especially hierarchical networks.

>>> file = open("test.html")
>>> file.write('<html>\n<head>\n<title>HiNetHTML Test</title>\n</head>\n<body>\n')
>>> hinet_html = mdp.hinet.HiNetHTML(file)
>>> flow = mdp.Flow([layer])
>>> hinet_html.parse_flow(flow)
>>> file.write('</body>\n</html>')
>>> file.close()

>>> switchboard = mdp.hinet.Rectangular2dSwitchboard(x_in_channels=50,
...                                                  y_in_channels=50,
...                                                  x_field_channels=10,
...                                                  y_field_channels=10,
...                                                  x_field_spacing=5,
...                                                  y_field_spacing=5,
...                                                  in_channel_dim=3)
>>> sfa_dim = 48
>>> sfa_node = mdp.nodes.SFANode(input_dim=switchboard.out_channel_dim,
...                              output_dim=sfa_dim)
>>> sfa2_dim = 32
>>> sfa2_node = mdp.nodes.SFA2Node(input_dim=sfa_dim,
...                                output_dim=sfa2_dim)
>>> flownode = mdp.hinet.FlowNode(mdp.Flow([sfa_node, sfa2_node]))
>>> sfa_layer = mdp.hinet.CloneLayer(flownode,
...                                  n_nodes=switchboard.output_channels)
>>> flow = mdp.Flow([switchboard, sfa_layer])

A real life example (Logistic maps)

We show an application of Slow Feature Analysis to the analysis of non-stationary time series. We consider a chaotic time series generated by the logistic map based on the logistic equation (a demographic model of the population biomass of species in the presence of limiting factors such as food supply or disease), and extract the slowly varying parameter that is hidden behind the time series. This example reproduces some of the results reported in: Laurenz Wiskott, Estimating Driving Forces of Nonstationary Time Series with Slow Feature Analysis. arXiv.org e-Print archive, http://arxiv.org/abs/cond-mat/0312317

Generate the slowly varying driving force, a combination of three sine waves (freqs: 5, 11, 13 Hz), and define a function to generate the logistic map

>>> p2 = mdp.numx.pi*2
>>> t = mdp.numx.linspace(0,1,10000,endpoint=0) # time axis 1s, samplerate 10KHz
>>> dforce = mdp.numx.sin(p2*5*t) + mdp.numx.sin(p2*11*t) + mdp.numx.sin(p2*13*t)
>>> def logistic_map(x,r):
...     return r*x*(1-x)
...
>>>

Note that we define series to be a two-dimensional array. Inputs to MDP must be two-dimensional arrays with variables on columns and observations on rows. In this case we have only one variable:

>>> series = mdp.numx.zeros((10000,1),'d')

Fix the initial condition:

>>> series[0] = 0.6

Generate the time-series using the logistic equation the driving force modifies the logistic equation parameter r:

>>> for i in range(1,10000):
...     series[i] = logistic_map(series[i-1],3.6+0.13*dforce[i])
...
>>>

If you have a plotting package series should look like this:

Define a flow to perform SFA in the space of polynomials of degree 3. We need a node that embeds the time-series in a 10 dimensional space, where different variables correspond to time-delayed copies of the original time-series: the TimeFramesNode(10). Then we need a node that expands the new signal in the space of polynomials of degree 3: the PolynomialExpansionNode(3). Finally we perform SFA onto the expanded signal and keep the slowest feature: SFANode(output_dim=1). We also measure the slowness of the input time-series and of the slow feature obtained by SFA. Therefore we put at the beginning and at the end of the sequence an analysis node that computes the eta-value (a measure of slowness) of its input (see docs for the definition of eta-value): the EtaComputerNode():

>>> sequence = [mdp.nodes.EtaComputerNode(),
...             mdp.nodes.TimeFramesNode(10),
...             mdp.nodes.PolynomialExpansionNode(3),
...             mdp.nodes.SFANode(output_dim=1),
...             mdp.nodes.EtaComputerNode()]
...
>>>
>>> flow = mdp.Flow(sequence, verbose=1)

Since the time-series is short enough to be kept in memory we don't need to define generators and we can feed the flow directly with the whole signal:

>>> flow.train(series)

Since the second and the third nodes are not trainable we are going to get two warnings (Training Interrupted). We can safely ignore them. Execute the flow to get the slow feature

>>> slow = flow(series)

The slow feautre should match the driving force up to a scaling factor, a constant offset and the sign. To allow a comparison we rescale the driving force to have zero mean and unit variance:

>>> resc_dforce = (dforce - mdp.numx.mean(dforce,0))/mdp.numx.std(dforce,0)

Print covariance between the rescaled driving force and the slow feature. Note that embedding the time-series with 10 time frames leads to a time-series with 9 observations less:

>>> mdp.utils.cov2(resc_dforce[:-9],slow)
0.99992501533859179

Print the eta-values of the chaotic time-series and of the slow feature

>>> print 'Eta value (time-series): ', flow[0].get_eta(t=10000)
Eta value (time-series):  [ 3002.53380245]
>>> print 'Eta value (slow feature): ', flow[-1].get_eta(t=9996)
Eta value (slow feature):  [ 10.2185087]

If you have a plotting package you could plot resc_dforce together with slow and see that they match perfectly:

Another real life example (Growing neural gas)

We generate uniformly distributed random data points confined on different 2-D geometrical objects. The Growing Neural Gas Node builds a graph with the same topological structure.

Fix the random seed to obtain reproducible results:

>>> mdp.numx_rand.seed(1266090063)

Some functions to generate uniform probability distributions on different geometrical objects:

>>> def uniform(min_, max_, dims):
...     """Return a random number between min_ and max_ ."""
...     return mdp.numx_rand.random(dims)*(max_-min_)+min_
...
>>> def circumference_distr(center, radius, n):
...     """Return n random points uniformly distributed on a circumference."""
...     phi = uniform(0, 2*mdp.numx.pi, (n,1))
...     x = radius*mdp.numx.cos(phi)+center[0]
...     y = radius*mdp.numx.sin(phi)+center[1]
...     return mdp.numx.concatenate((x,y), axis=1)
...
>>> def circle_distr(center, radius, n):
...     """Return n random points uniformly distributed on a circle."""
...     phi = uniform(0, 2*mdp.numx.pi, (n,1))
...     sqrt_r = mdp.numx.sqrt(uniform(0, radius*radius, (n,1)))
...     x = sqrt_r*mdp.numx.cos(phi)+center[0]
...     y = sqrt_r*mdp.numx.sin(phi)+center[1]
...     return mdp.numx.concatenate((x,y), axis=1)
...
>>> def rectangle_distr(center, w, h, n):
...     """Return n random points uniformly distributed on a rectangle."""
...     x = uniform(-w/2., w/2., (n,1))+center[0]
...     y = uniform(-h/2., h/2., (n,1))+center[1]
...     return mdp.numx.concatenate((x,y), axis=1)
...
>>> N = 2000

Explicitly collect random points from some distributions:

• Circumferences:

  >>> cf1 = circumference_distr([6,-0.5], 2, N)
  >>> cf2 = circumference_distr([3,-2], 0.3, N)
  

• Circles:

  >>> cl1 = circle_distr([-5,3], 0.5, N/2)
  >>> cl2 = circle_distr([3.5,2.5], 0.7, N)
  

• Rectangles:

  >>> r1 = rectangle_distr([-1.5,0], 1, 4, N)
  >>> r2 = rectangle_distr([+1.5,0], 1, 4, N)
  >>> r3 = rectangle_distr([0,+1.5], 2, 1, N/2)
  >>> r4 = rectangle_distr([0,-1.5], 2, 1, N/2)
  

Shuffle the points to make the statistics stationary

>>> x = mdp.numx.concatenate([cf1, cf2, cl1, cl2, r1,r2,r3,r4], axis=0)
>>> x = mdp.numx.take(x,mdp.numx_rand.permutation(x.shape[0]), axis=0)

If you have a plotting package x should look like this:

Create a GrowingNeuralGasNode and train it:

>>> gng = mdp.nodes.GrowingNeuralGasNode(max_nodes=75)

The initial distribution of nodes is randomly chosen:

The training is performed in small chunks in order to visualize the evolution of the graph:

>>> STEP = 500
>>> for i in range(0,x.shape[0],STEP):
...     gng.train(x[i:i+STEP])
...     # [...] plotting instructions
...
>>> gng.stop_training()

See here the animation of training.

Visualizing the neural gas network, we'll see that it is adapted to the topological structure of the data distribution:

Calculate the number of connected components:

>>> n_obj = len(gng.graph.connected_components())
>>> print n_obj
5

Node List

Here is the complete list of implemented nodes. Refer to the API for the full documentation and interface description.

• CuBICANode

  Perform Independent Component Analysis using the CuBICA algorithm.

  Reference: Blaschke, T. and Wiskott, L. (2003). CuBICA: Independent Component Analysis by Simultaneous Third- and Fourth-Order Cumulant Diagonalization. IEEE Transactions on Signal Processing, 52(5), pp. 1250-1256. More information about ICA can be found among others in Hyvarinen A., Karhunen J., Oja E. (2001). Independent Component Analysis, Wiley.

• EtaComputerNode

  Compute the eta values of the normalized training data. The delta value of a signal is a measure of its temporal variation, and is defined as the mean of the derivative squared, i.e. delta(x) = mean(dx/dt(t)^2). delta(x) is zero if 'x' is a constant signal, and increases if the temporal variation of the signal is bigger. The eta value is a more intuitive measure of temporal variation, defined as eta(x) = T/(2*pi) * sqrt(delta(x)). If 'x' is a signal of length 'T' which consists of a sine function that accomplishes exactly 'N' oscillations, then eta(x) = N.

  Reference: Wiskott, L. and Sejnowski, T.J. (2002). Slow Feature Analysis: Unsupervised Learning of Invariances, Neural Computation, 14(4):715-770.

• FANode

  Perform Factor Analysis. The current implementation should be most efficient for long data sets: the sufficient statistics are collected in the training phase, and all EM-cycles are performed at its end. More information about Factor Analysis can be found in Max Welling's classnotes in the chapter "Linear Models".

• FastICANode

  Perform Independent Component Analysis using the FastICA algorithm.

  Reference: Aapo Hyvarinen (1999). Fast and Robust Fixed-Point Algorithms for Independent Component Analysis, IEEE Transactions on Neural Networks, 10(3):626-634. More information about ICA can be found among others in Hyvarinen A., Karhunen J., Oja E. (2001). Independent Component Analysis, Wiley.

• FDANode

  Perform a (generalized) Fisher Discriminant Analysis of its input. It is a supervised node that implements FDA using a generalized eigenvalue approach.

  More information on Fisher Discriminant Analysis can be found for example in C. Bishop, Neural Networks for Pattern Recognition, Oxford Press, pp. 105-112.

• GaussianClassifierNode

  Perform a supervised Gaussian classification. Given a set of labelled data, the node fits a gaussian distribution to each class.

• GrowingNeuralGasNode

  Learn the topological structure of the input data by building a corresponding graph approximation.

  More information about the Growing Neural Gas algorithm can be found in B. Fritzke, A Growing Neural Gas Network Learns Topologies, in G. Tesauro, D. S. Touretzky, and T. K. Leen (editors), Advances in Neural Information Processing Systems 7, pages 625-632. MIT Press, Cambridge MA, 1995.

• HitParadeNode

  Collect the first 'n' local maxima and minima of the training signal which are separated by a minimum gap 'd'.

• ISFANode

  Perform Independent Slow Feature Analysis on the input data.

  More information about ISFA can be found in: Blaschke, T. , Zito, T., and Wiskott, L. Independent Slow Feature Analysis and Nonlinear Blind Source Separation. Neural Computation 19(4):994-1021 (2007).

• JADENode

  Original code contributed by Gabriel Beckers.

  Perform Independent Component Analysis using the JADE algorithm.

  References: Cardoso, J.-F, and Souloumiac, A. Blind beamforming for non Gaussian signals. Radar and Signal Processing, IEE Proceedings F, 140(6): 362-370 (1993), and Cardoso, J.-F. High-order contrasts for independent component analysis. Neural Computation, 11(1): 157-192 (1999). More information about ICA can be found among others in Hyvarinen A., Karhunen J., Oja E. (2001). Independent Component Analysis, Wiley.

• NIPALSNode

  Original code contributed by Michael Schmuker, Susanne Lezius, and Farzad Farkhooi.

  Perform Principal Component Analysis using the NIPALS algorithm. This algorithm is particularyl useful if you have more variable than observations, or in general when the number of variables is huge and calculating a full covariance matrix may be unfeasable. It's also more efficient of the standard PCANode if you expect the number of significant principal components to be a small. In this case setting output_dim to be a certain fraction of the total variance, say 90%, may be of some help.

  Reference for NIPALS (Nonlinear Iterative Partial Least Squares): Wold, H. Nonlinear estimation by iterative least squares procedures. in David, F. (Editor), Research Papers in Statistics, Wiley, New York, pp 411-444 (1966).

  More information about Principal Component Analysis, a.k.a. discrete Karhunen-Loeve transform can be found among others in I.T. Jolliffe, Principal Component Analysis, Springer-Verlag (1986).

• NoiseNode

  Original code contributed by Mathias Franzius.

  Inject multiplicative or additive noise into the input data.

• PCANode

  Filter the input data throug the most significatives of its principal components.

  More information about Principal Component Analysis, a.k.a. discrete Karhunen-Loeve transform can be found among others in I.T. Jolliffe, Principal Component Analysis, Springer-Verlag (1986).

• PolynomialExpansionNode

  Perform expansion in a polynomial space.

• QuadraticExpansionNode

  Perform expansion in the space formed by all linear and quadratic monomials

• RBMNode

  Implementation of a Restricted Boltzmann Machine.

  For more information on RBMs, see Geoffrey E. Hinton (2007) Boltzmann machine. Scholarpedia, 2(5):1668

• RBMWithLabelsNode

  Implementation of a Restricted Boltzmann Machine with softmax labels.

  For more information on RBMs, see Geoffrey E. Hinton (2007) Boltzmann machine Scholarpedia, 2(5):1668

  Hinton, G. E, Osindero, S., and Teh, Y. W. A fast learning algorithm for deep belief nets, Neural Computation, 18:1527-1554 (2006).

• SFANode

  Extract the slowly varying components from the input data.

  More information about Slow Feature Analysis can be found in Wiskott, L. and Sejnowski, T.J., Slow Feature Analysis: Unsupervised Learning of Invariances, Neural Computation, 14(4):715-770 (2002).

• SFA2Node

  Get an input signal, expand it in the space of inhomogeneous polynomials of degree 2 and extract its slowly varying components. The get_quadratic_form method returns the input-output function of one of the learned unit as a mdp.utils.QuadraticForm object.

  More information about Slow Feature Analysis can be found in Wiskott, L. and Sejnowski, T.J., Slow Feature Analysis: Unsupervised Learning of Invariances, Neural Computation, 14(4):715-770 (2002).

• TDSEPNode

  Perform Independent Component Analysis using the TDSEP algorithm. Note that TDSEP, as implemented in this Node, is an online algorithm, i.e. it is suited to be trained on huge data sets, provided that the training is done sending small chunks of data for each time.

  Reference: Ziehe, Andreas and Muller, Klaus-Robert (1998). TDSEP an efficient algorithm for blind separation using time structure. in Niklasson, L, Boden, M, and Ziemke, T (Editors), Proc. 8th Int. Conf. Artificial Neural Networks (ICANN 1998).

• TimeFramesNode

  Copy delayed version of the input signal on the space dimensions.

  For example, for time_frames=3 and gap=2:
  
  [ X(1) Y(1)        [ X(1) Y(1) X(3) Y(3) X(5) Y(5)
    X(2) Y(2)          X(2) Y(2) X(4) Y(4) X(6) Y(6)
    X(3) Y(3)   -->    X(3) Y(3) X(5) Y(5) X(7) Y(7)
    X(4) Y(4)          X(4) Y(4) X(6) Y(6) X(8) Y(8)
    X(5) Y(5)          ...  ...  ...  ...  ...  ... ]
    X(6) Y(6)
    X(7) Y(7)
    X(8) Y(8)
    ...  ...  ]
  

• WhiteningNode

  'Whiten' the input data by filtering it through the most significatives of its principal components. All output signals have zero mean, unit variance and are decorrelated.

Additional utilities

MDP offers some additional utilities of general interest in the mdp.utils module. Refer to the API for the full documentation and interface description.

CovarianceMatrix

This class stores an empirical covariance matrix that can be updated incrementally. A call to the fix method returns the current state of the covariance matrix, the average and the number of observations, and resets the internal data.

Note that the internal sum is a standard __add__ operation. We are not using any of the fancy sum algorithms to avoid round off errors when adding many numbers. If you want to contribute a CovarianceMatrix class that uses such algorithms we would be happy to include it in MDP. For a start see the Python recipe by Raymond Hettinger. For a review about floating point arithmetic and its pitfalls see this interesting article.

DelayCovarianceMatrix

This class stores an empirical covariance matrix between the signal and time delayed signal that can be updated incrementally.

MultipleCovarianceMatrices

Container class for multiple covariance matrices to easily execute operations on all matrices at the same time.

dig_node(node)

Crawl recursively an MDP Node looking for arrays. Return (dictionary, string), where the dictionary is: { attribute_name: (size_in_bytes, array_reference)} and string is a nice string representation of it.

get_node_size(node)

Get 'node' total byte-size using cPickle with protocol=2. (The byte-size is related the memory needed by the node).

progressinfo(sequence, length, style, custom)

A fully configurable text-mode progress info box. To get a progress info box for your loops use it like this:

>>> for i in progressinfo(sequence):
...     do_something(i)

You can also use it with generators, files or any other iterable object, but in this case you have to specify the total length of the sequence:

>>> for line in progressinfo(open_file, nlines):
...     do_something(line)

A few examples of the available layouts:

[===================================73%==============>...................]

Progress:  67%[======================================>                   ]

23% [02:01:28] - [00:12:37]

QuadraticForm

Define an inhomogeneous quadratic form as 1/2 x'Hx + f'x + c. This class implements the quadratic form analysis methods presented in: Berkes, P. and Wiskott, L. On the analysis and interpretation of inhomogeneous quadratic forms as receptive fields. Neural Computation, 18(8): 1868-1895. (2006).

refcast(array, dtype)

Cast the array to 'dtype' only if necessary, otherwise return a reference.

rotate(mat, angle, columns, units)

Rotate in-place a NxM data matrix in the plane defined by the 'columns' when observation are stored on rows. Observations are rotated counterclockwise. This corresponds to the following matrix-multiplication for each data-point (unchanged elements omitted):

[  cos(angle) -sin(angle)     [ x_i ]
   sin(angle)  cos(angle) ] * [ x_j ]

random_rot(dim, dtype)

Return a random rotation matrix, drawn from the Haar distribution (the only uniform distribution on SO(n)). The algorithm is described in the paper Stewart, G.W., The efficient generation of random orthogonal matrices with an application to condition estimators, SIAM Journal on Numerical Analysis, 17(3), pp. 403-409, 1980. For more information see this Wikipedia entry.

symrand(dim_or_eigv, dtype)

Return a random symmetric (Hermitian) matrix with eigenvalues uniformly distributed on (0,1].

To Do

In this section we want to give you an overview about our plans for the development of MDP:

• Add more data processing algorithms.
• Advanced usage of the hinet package will be possible only in presence of an easy and intuitive GUI :)
• Wait for a good guy who wants to contribute a CovarianceMatrix class that uses some of the fancy sum algorithms to avoid round off errors when adding many numbers.

Contributors

In this final section we want to thank all users who have contributed code to the MDP project. Strictly in alphabetical order:

• Gabriel Beckers
• Farzad Farkhooi
• Mathias Franzius
• Susanne Lezius
• Michael Schmuker