MATLAB Tutorial – CCN Course 2012 Malte J. Rasch July 17, 2012

Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
MATLAB Tutorial – CCN Course 2012
How to code a neural network simulation
Malte J. Rasch
National Key Laboratory of Cognitive Neuroscience and Learning
Beijing Normal University
China
July 17, 2012
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Overview
Basic introduction to Matlab
Learn to code a neural network simulation
Further exercises with solutions
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Why MATLAB ?
Pro:
Matrix-like numerical operations very fast and easy to use
Good plotting capabilities
Script language for programming small to medium sized
problems in applied mathematics (rapid prototyping)
Widely used in the neuroscience community for data analysis
as well as computational projects
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Why MATLAB ?
Pro:
Matrix-like numerical operations very fast and easy to use
Good plotting capabilities
Script language for programming small to medium sized
problems in applied mathematics (rapid prototyping)
Widely used in the neuroscience community for data analysis
as well as computational projects
Contra:
Support of symbolic/analytic expressions less advanced
Mathematica, Maple
Often not flexible/ fast enough for big projects
e.g. Python much more versatile
specialized software for detailed/large neural network
simulations
(NEURON, PCSIM, NEST, GENESIS, ...)
Very expensive, especially when using on a cluster
free alternative: (scientific) python
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Starting MATLAB
Easy: Just click on MATLAB symbol...
Step 2
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Scalar expressions
Binary operations: work as expected, use = + - * / ^.
Example (compute y =
a2 x
2+a
+ b)
>> a = 2;
>> b = 1;
>> x = 0.5;
>> y = a^2*x/(2+a) + b;
>> y
y =
1.500
Unary operations: called as functions with (), eg.
sin cos tan atan sqrt log gamma
Example ( compute y =
√
sin x
ln x )
>> x = 0.5;
>> y = sqrt(sin(x))/log(x);
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
MATLAB == MATrix LABoratory
Not suprisingly, in MATLAB everything is about matrices.
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
MATLAB == MATrix LABoratory
Not suprisingly, in MATLAB everything is about matrices.
However, the matrix-like datastructure in MATLAB is better called
a n-dimensional array, because it can be manipulated in
non-algebraic ways.
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
MATLAB == MATrix LABoratory
Not suprisingly, in MATLAB everything is about matrices.
However, the matrix-like datastructure in MATLAB is better called
a n-dimensional array, because it can be manipulated in
non-algebraic ways.
(Almost) all functions will work on arrays as well
usually element-wise
Many MATLAB functions will produce arrays as output
Array operations much faster than for-looped element-wise
operation
Introduction
Arrays
Coding
Plotting
Network model:
Arrays: Overview
How to initialize arrays
Indexing
Calculating with arrays
Step 1
Step 2
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Calculating with arrays is
straight-forward
however, carefully check
the size of matrices
if element-wise or matrix-like operations are intended
which matrix dimension to operate on
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Calculating with arrays is
straight-forward
however, carefully check
the size of matrices
if element-wise or matrix-like operations are intended
which matrix dimension to operate on
Example (compute yi = Axi + b with b, xi ∈ R2 )
>> A = [1,0.2;0.4,1];
>> b = [1;2] + 0.1;
>> x = 2*randn(2,1);
>> y = A * x + b
y =
4.5535
1.4856
>> N = 5;
% same as b = [b,b,b,b,b]
>> bi = repmat(b,[1,N]);
>> xi = 2*randn(2,N);
>> xi(:,1) = x;
>> yi = A * xi + bi
yi =
4.5535 2.3126 [..] -0.9021
1.4856 6.8091 [..] -0.1080
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Syntax for initializing arrays
1
Implicitly, using function returning an array
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Syntax for initializing arrays
1
2
Implicitly, using function returning an array
By explicit concatenation
Concatenation of columns of arrays
[ arr1, arr2, ... , arrn]
Concatenation of rows of arrays
[ arr1; arr2; ... ; arrn]
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Syntax for initializing arrays
1
2
Implicitly, using function returning an array
By explicit concatenation
Concatenation of columns of arrays
[ arr1, arr2, ... , arrn]
Concatenation of rows of arrays
[ arr1; arr2; ... ; arrn]
Note: an scalar is also regarded as an array (of size [1,1]).
Note 2: arrays must have matching sizes.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Syntax for initializing arrays
1
2
Implicitly, using function returning an array
By explicit concatenation
Concatenation of columns of arrays
[ arr1, arr2, ... , arrn]
Concatenation of rows of arrays
[ arr1; arr2; ... ; arrn]
Note: an scalar is also regarded as an array (of size [1,1]).
Note 2: arrays must have matching sizes.
Example (Concatenation)
>> A = [1,2;3,4]
A =
1
2
3
4
>> B = [A,A]
B =
1
2
1
2
3
4
3
4
>> C = [A;A]
C =
1
2
3
4
1
2
3
4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Syntax for initializing arrays: implicitly
Functions that pre-allocate memory and set each array element to
an initial value:
zeros – all zero ND-arrays
ones – all one ND-arrays
rand – random ND-arrays (equally in [0, 1])
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Syntax for initializing arrays: implicitly
Functions that pre-allocate memory and set each array element to
an initial value:
zeros – all zero ND-arrays
ones – all one ND-arrays
rand – random ND-arrays (equally in [0, 1])
colon (:) – for linear sequences, syntax:
istart:[step:]iend
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Syntax for initializing arrays: implicitly
Functions that pre-allocate memory and set each array element to
an initial value:
zeros – all zero ND-arrays
ones – all one ND-arrays
rand – random ND-arrays (equally in [0, 1])
colon (:) – for linear sequences, syntax:
istart:[step:]iend
– many others, for details type: help command
randn, linspace, logspace, ndgrid, repmat,...
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Syntax for initializing arrays: implicitly
Functions that pre-allocate memory and set each array element to
an initial value:
zeros – all zero ND-arrays
ones – all one ND-arrays
rand – random ND-arrays (equally in [0, 1])
colon (:) – for linear sequences, syntax:
istart:[step:]iend
– many others, for details type: help command
randn, linspace, logspace, ndgrid, repmat,...
Example (Functions initializing arrays)
>> A = zeros(3,3);
>> A = ones(4,4,4);
>> size(A)
ans =
4
4
4
>> x =
x =
3.0
>> A =
1
1
3:-0.5:1
2.5 2.0
ones(2)
1
1
1.5
1.0
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Indexing arrays
1
Subscript of a matrix:
access the (i, j)-th element of a 2D-matrix A of dimension (m, n)
>> A(i,j)
Note: The first index is always 1 (not 0 as in most other languages)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Indexing arrays
1
Subscript of a matrix:
access the (i, j)-th element of a 2D-matrix A of dimension (m, n)
>> A(i,j)
Note: The first index is always 1 (not 0 as in most other languages)
2
Linear index of a matrix:
access the (i, j)th element of the 2D-matrix A of dimension (m, n)
>> linearidx = i + (j-1)*m;
>> A(linearidx)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Indexing arrays
1
Subscript of a matrix:
access the (i, j)-th element of a 2D-matrix A of dimension (m, n)
>> A(i,j)
Note: The first index is always 1 (not 0 as in most other languages)
2
Linear index of a matrix:
access the (i, j)th element of the 2D-matrix A of dimension (m, n)
>> linearidx = i + (j-1)*m;
>> A(linearidx)
3
“Slice” indexing with “:”
access the i-th row and jth column in A, respectively
>> A(i,:)
>> A(:,j)
get all elements as a cocatenated column vector
>> A(:)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Indexing arrays (2)
4
Multiple indices
vectors of linear indices can be used
>> A([1, 4, 5 ,6])
access the 1st to 4th rows of the 2D-matrix A of dimension (m, n)
>> A(1:4,:)
access the 2nd (m,n)-slice of a 3D-matrix B of dimension (m, n, p)
>> B(:,:,2)
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Indexing arrays (2)
4
Multiple indices
vectors of linear indices can be used
>> A([1, 4, 5 ,6])
access the 1st to 4th rows of the 2D-matrix A of dimension (m, n)
>> A(1:4,:)
access the 2nd (m,n)-slice of a 3D-matrix B of dimension (m, n, p)
>> B(:,:,2)
5
Logical indexing
logical matrices of the same size as A can be used as indices
>> A(A>0)
>> A(find(A>0))
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Calculating with arrays
1
Element-wise interpretation
For instance, sin cos log etc.
Reserved symbols, .* ./ .^
2
“true” matrix interpretation (with dot product)
3
Operations on one specified dimensions of the matrix
Symbols * / ^ etc.
For instance, sum mean max etc.
4
Array manipulations
eg. reshape repmat permute circshift tril
Example (element-wise product and dot product)
>> A = ones(2,2);
>> A.*A
ans =
1
1
1
1
>> A*A
ans =
2
2
2
2
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Useful operations on specified dimensions
Often used functions on array elements include
size – number of elements in a dimension
sum – sum of elements along an array dimension
prod – product of elements
all, any – logical AND or OR, respectively
mean – mean of elements
max – maximal element of an array dimension
– and many more, eg. min var std median diff
Functions are usually called as
res = func(arr,dim);
Note: max and min are exceptions! Here: res = max(arr,[],dim);
For details type: help func
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
MATLAB is a script language
Scripts are blocks of code which can be called within MATLAB or
within another script.
They are text-files with extensions “.m”.
They should contain all commands associated with a scientific
project.
(at least to easily reproduce and the calculations)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
MATLAB is a script language
Scripts are blocks of code which can be called within MATLAB or
within another script.
They are text-files with extensions “.m”.
They should contain all commands associated with a scientific
project.
(at least to easily reproduce and the calculations)
There are basically two types of m-files
1 m-file script
A squential list of MATLAB commands
The variable scope is that of the caller.
2
m-file function
m-file functions accept input parameters and deliver outputs.
The variable scope is different from that of the caller.
Introduction
Arrays
Coding
Plotting
Network model:
m-file scripts
How to write an m-file script?
Step 1
Step 2
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file scripts
How to write an m-file script?
1
open a text-editor of your choice or use the editor provided
with MATLAB
>> edit myscript
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file scripts
How to write an m-file script?
1
open a text-editor of your choice or use the editor provided
with MATLAB
>> edit myscript
2
Write all your calculations in a text-file with extension “.m”
(here myscript.m)
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file scripts
How to write an m-file script?
1
open a text-editor of your choice or use the editor provided
with MATLAB
>> edit myscript
2
Write all your calculations in a text-file with extension “.m”
(here myscript.m)
3
Save file in the current working directory
(or addpath to search path )
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file scripts
How to write an m-file script?
1
open a text-editor of your choice or use the editor provided
with MATLAB
>> edit myscript
2
Write all your calculations in a text-file with extension “.m”
(here myscript.m)
3
Save file in the current working directory
(or addpath to search path )
4
Call your script by calling it from the “Command Window”
>> myscript;
(no compilation needed)
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
m-file scripts
How to write an m-file script?
1
open a text-editor of your choice or use the editor provided
with MATLAB
>> edit myscript
2
Write all your calculations in a text-file with extension “.m”
(here myscript.m)
3
Save file in the current working directory
(or addpath to search path )
4
Call your script by calling it from the “Command Window”
>> myscript;
(no compilation needed)
5
Note: The variable scope is that of the caller. This means that
the variables of the script are now present in the workspace
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
m-file scripts
Excercise #1 (How to write an m-file script)
Write an m-file computing the Stirling approximation
n n
√
n! ≈ 2πn
e
Hint: π is defined in MATLAB as pi, and e x as exp(x).
Solution #1
Reminder:
1 Open a text-editor
>> edit myscript
2 Write calculations in text-file with extension “.m”
3 Save file into current working directory
4 Call your script
>> myscript;
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
m-file functions
How to write an m-file function?
Step 1
Step 2
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
m-file functions
How to write an m-file function?
Identical to writing an m-file script, except:
Step 2
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file functions
How to write an m-file function?
Identical to writing an m-file script, except:
In the first line input and output arguments are defined:
function outarg = myfun(inarg1,inarg2,...);
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file functions
How to write an m-file function?
Identical to writing an m-file script, except:
In the first line input and output arguments are defined:
function outarg = myfun(inarg1,inarg2,...);
Call the function with
>> result = myfun(p1,p2,...);
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file functions
How to write an m-file function?
Identical to writing an m-file script, except:
In the first line input and output arguments are defined:
function outarg = myfun(inarg1,inarg2,...);
Call the function with
>> result = myfun(p1,p2,...);
Note: The variable scope is different from that of the caller. That
means, that variables defined in the function body are NOT
present in the workspace after calling the function
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file functions
How to write an m-file function?
Identical to writing an m-file script, except:
In the first line input and output arguments are defined:
function outarg = myfun(inarg1,inarg2,...);
Call the function with
>> result = myfun(p1,p2,...);
Note: The variable scope is different from that of the caller. That
means, that variables defined in the function body are NOT
present in the workspace after calling the function
Note 2: Input arguments are always referenced by value
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file functions
How to write an m-file function?
Identical to writing an m-file script, except:
In the first line input and output arguments are defined:
function outarg = myfun(inarg1,inarg2,...);
Call the function with
>> result = myfun(p1,p2,...);
Note: The variable scope is different from that of the caller. That
means, that variables defined in the function body are NOT
present in the workspace after calling the function
Note 2: Input arguments are always referenced by value
Note 3: When called the m-file file name is used.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
m-file functions
Excercise #2 (writing a m-file function)
Write an m-file function computing the Stirling approximation
n n
√
n! ≈ 2πn
e
for a given n
Solution #2
Reminder:
In the first line input and output arguments are defined:
function outarg = myfun(inarg1,inarg2,...);
Call the function with
>> result = myfun(p1,p2,...);
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Basic syntax
The basic syntax is similar to other script languages, eg. “python”.
MATLAB has the usual flow control structures, namely
if ..
else for decisions
for-loop
while-loop
switch ..
try ..
case for multiple if-else decisions
catch for handling errors
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Basic syntax: flow control (1)
if-else block syntax:
1
2
3
4
5
if scalar condition
expressions
else
expressions
end
Relational operators, eg.: == (equals), || (or), && (and), ~ (not)
for details type: help relop
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Basic syntax: flow control (1)
if-else block syntax:
1
2
3
4
5
if scalar condition
expressions
else
expressions
end
Relational operators, eg.: == (equals), || (or), && (and), ~ (not)
for details type: help relop
Example (if-else)
a = rand ( 1 ) ;
i f a == 0 . 5
f p r i n t f ( ’ you a r e v e r y l u c k y ! ’ ) ;
end
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Programming: basic flow control (2)
for-loop block syntax:
1
2
3
4
for i = array
% i==a r r a y ( j ) i n t h e j −t h l o o p
expressions
end
(one can also use break and continue keywords)
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Programming: basic flow control (2)
for-loop block syntax:
1
2
3
4
for i = array
% i==a r r a y ( j ) i n t h e j −t h l o o p
expressions
end
(one can also use break and continue keywords)
Example (for loop)
a =0;
for i = 1:100
a = a+i ;
end
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Plotting with MATLAB
Very flexible plotting tools
Many functions for a variety of plot types and graphics
drawing
Simple cell response
Response profile
2
〈 |⋅|2 〉
RF filter
K
x,y,t
θ,ρ,ω
1.5
1
θ
2
y
y
(+ phase invariance)
0.5
θ3
vary θ
t
t
x
Low speed
E
High speed
D
Low speed
90°
180°
High speed
Horz. RF
Vert. RF
10
5
TF [Hz]
→ Movement dir. →
C
→ Movement dir. →
Static
B
0
0°
:
x
Avg. response
Stimulus
A
0
Response
−5
−10
2
°
0
°
90
°
°
180 0
°
90
°
°
180 0
°
90
°
180
0
SF [cyc/deg]
−2
−2
0
2
SF [cyc/deg]
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Overview
Just two most useful plotting commands
plot – plots all kinds of lines, dots, markers in 2D
imagesc – plots an 2-D pseudo image
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Overview
Just two most useful plotting commands
plot – plots all kinds of lines, dots, markers in 2D
imagesc – plots an 2-D pseudo image
Note: For an overview of 2D plotting commands: help graph2d
Note 2: For fancy graphs see: help specgraph
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
How to plot
How to plot:
1
Open a figure (window) with
>> figure;
2
Create an axes object with
>> axes;
3
Type the plotting command of your choice
>> plot(x,y,’b--’,’LineWidth’,2);
Note: graphics commands draw into the current axes (gca) on the current
figure (gcf). They automatically create a new figure and axes if necessary.
4
Make the plot nicer by adding labels and setting limits, eg:
>>
>>
>>
>>
xlabel(’Time [s]’);
ylabel(’Response [spks/sec]’);
xlim([-1,1]); ylim([-2,2]);
title(’Simulation’)
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
plot command
Basic syntax:
handle = plot(X,Y,linespec,optname1,val1,...);
X,Y – x- and y-values to plot. If Y is a 2-D array, all columns are
plotted as different lines
linespec – a string with a short hand for color and line style and
marker type. See help plot for an overview. Eg,
linespec = ’:ko’
plots dotted (:) black line (k) with a circle at each given
coordinate (xi , yi )
optname1 – a string specifing the option name, eg. ’LineWidth’
val1 – a corresponding value.
handle – graphics handle for later reference, eg. with
>> set(handle,optname1,val1)
Tip: To get an overview over all possible options, try get(handle)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Plotting example: empirical PDF of Gaussian
1
figure ;
2
3
4
5
6
% p l o t G a u s s i a n PDF
x = −4:0.01:4;
y = exp ( −0.5∗ x . ˆ 2 ) / s q r t ( 2 ∗ p i ) ;
p l o t ( x , y , ’−k ’ , ’ L i n e w i d t h ’ , 2 ) ;
0.45
0.35
9
10
11
12
% e m p i r i c a l pdf
% h i s t o g r a m o f 1000 random v a l u e s
[ N, z ]= h i s t ( r a n d n ( 1 0 0 0 , 1 ) ) ;
%n o r m a l i z a t i o n
p = N/sum (N) / ( z (2) − z ( 1 ) ) ;
13
14
15
16
17
20
21
0.3
0.25
0.2
0.15
0.1
0.05
h o l d on ; %a d d s a l l f o l l o w i n g g r a p h s
%t o t h e c u r r e n t a x e s
plot ( z , p , ’ rx ’ , ’ MarkerSize ’ ,10);
hold o f f ;
18
19
Probability density
7
8
true PDF
empirical PDF
0.4
x l a b e l ( ’ Random v a r i a b l e x ’ )
ylabel ( ’ Probability density ’ )
l e g e n d ( ’ t r u e PDF ’ , ’ e m p i r i c a l PDF ’ )
0
−4
−3
−2
−1
0
1
Random variable x
2
3
4
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Goal of tutorial
We will program a neural network simulation together.
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
We will practice on the way:
Writing scripts
Usage of array notation
How to integrate ODEs
How to plot results
How to simulate neurons and synapses
How to program a quite realistic network simulation
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Network model simulation
What has to be done in principle:
n neurons, excitatory and inhibitory, are inter-connected with
synapses.
The each neuron and synapse follows a particular dynamic
over time.
The simulation calculates the dynamic and spiking activity of
each neuron for each t.
The network gets some input and reacts in a certain and the
response of the network is plotted.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
We will proceed in 4 successive steps
1
2
3
4
Simulate
Simulate
Simulate
Simulate
a single neuron with current step input
a single neuron with Poisson input
1000 neurons (no recurrent connections)
complete network (conductance based)
Network activity
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Which neuron model to use?
Biophysical model (i.e. Hodgkin-Huxley model)
Cm
X
dVm
1
=−
(Vm − VL ) −
gi (t)(Vm − Ei ) + Isyn + Iapp
dt
Rm
i
Including non-linear dynamics of many channels in gi (t)
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Which neuron model to use?
Biophysical model (i.e. Hodgkin-Huxley model)
Cm
X
dVm
1
=−
(Vm − VL ) −
gi (t)(Vm − Ei ) + Isyn + Iapp
dt
Rm
i
Including non-linear dynamics of many channels in gi (t)
Mathematical simplification (Izhikevich, book chapter 8)
X
v˙ = (0.04v + 5) v + 150 − u −
wj gj (v − Ej ) + Iapp
j
u˙ = a (b v − u)
v ← c,
u ← u + d,
if
v ≥ 35
with g˙j = −gj /τg and gj ← gj + 1 if vj ≥ 35. It is b = 0.2, c = −65, and
d = 8, a = 0.02 for exc. neurons and d = 2, a = 0.1 for inh. neurons.
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Neuron model
peak 30 mV
v(t)
de
reset d
cay
with r
ate a
parameter b
reset c
if v = 30 mV,
then v c, u u + d
RZ
0.25
0.2
LTS,TC
FS
RS,IB,CH
sensitivity b
0 0.02
0.1
parameter a
RS
IB
4
2
0.05
u(t)
regular spiking (RS)
8
parameter d
v'= 0.04v 2+5v +140 - u + I
u'= a(bv - u)
FS,LTS,RZ
CH
TC
-65
-55
-50
parameter c
intrinsically bursting (IB)
chattering (CH)
fast spiking (FS)
thalamo-cortical (TC)
resonator (RZ)
low-threshold spiking (LTS)
v(t)
I(t)
thalamo-cortical (TC)
20 mV
40 ms
-63 mV
-87 mV
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 1: Simulate a single neuron with injected current
Exercise 1
Simulate the neuron model for 1000ms and plot the resulting
voltage trace. Apply a current step (Iapp = 7pA) between time
200ms and 700ms.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 1: Simulate a single neuron with injected current
Exercise 1
Simulate the neuron model for 1000ms and plot the resulting
voltage trace. Apply a current step (Iapp = 7pA) between time
200ms and 700ms.
Neuron model:
v˙
u˙
=
=
v
(0.04v + 5) v + 140 − u + Iapp
a (b v − u)
← c, u ← u + d,
if v ≥ 35
with d = 8, a = 0.02, b = 0.2,
c = −65 for an excitatory neuron.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
2
Reserve memory for voltage trace v and u (of length T = 1000/∆t)
HowTo
and set first element to −70 and −14, respectively.
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
2
3
Reserve memory for voltage trace v and u (of length T = 1000/∆t)
HowTo
and set first element to −70 and −14, respectively.
Loop over T − 1 time steps and do for each step t
HowTo
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
2
3
Reserve memory for voltage trace v and u (of length T = 1000/∆t)
HowTo
and set first element to −70 and −14, respectively.
Loop over T − 1 time steps and do for each step t HowTo
1 set Iapp ← 7 if t∆t is between 200 and 700 (otherwise 0)
HowTo
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
2
3
Reserve memory for voltage trace v and u (of length T = 1000/∆t)
HowTo
and set first element to −70 and −14, respectively.
Loop over T − 1 time steps and do for each step t HowTo
1 set Iapp ← 7 if t∆t is between 200 and 700 (otherwise 0)
HowTo
2
if vt < 35: update element t + 1
according to HowTo
vt+1
ut+1
←
←
HowTo
of v and u
vt + ∆t {(0.04 vt + 5) vt − ut + 140 + Iapp }
ut + ∆t a (b vt − ut )
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
2
3
Reserve memory for voltage trace v and u (of length T = 1000/∆t)
HowTo
and set first element to −70 and −14, respectively.
Loop over T − 1 time steps and do for each step t HowTo
1 set Iapp ← 7 if t∆t is between 200 and 700 (otherwise 0)
HowTo
2
if vt < 35: update element t + 1
according to HowTo
vt+1
ut+1
3
←
←
HowTo
of v and u
vt + ∆t {(0.04 vt + 5) vt − ut + 140 + Iapp }
ut + ∆t a (b vt − ut )
if vt ≥ 35: set vt+1 ← c and ut+1 ← ut + d (and optional
vt ← 35)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 1 in detail:
Open MATLAB
HowTo
HowTo
and create a new file (script) that will simulate the neuron
. If you do not know how to use a command get some help
HowTo
Proceed as follows:
1
Initialize parameter values (∆t = 0.5ms, a = 0.02, d = 8, · · · )
HowTo
2
3
Reserve memory for voltage trace v and u (of length T = 1000/∆t)
HowTo
and set first element to −70 and −14, respectively.
Loop over T − 1 time steps and do for each step t HowTo
1 set Iapp ← 7 if t∆t is between 200 and 700 (otherwise 0)
HowTo
2
if vt < 35: update element t + 1
according to HowTo
vt+1
ut+1
3
←
←
HowTo
of v and u
vt + ∆t {(0.04 vt + 5) vt − ut + 140 + Iapp }
ut + ∆t a (b vt − ut )
if vt ≥ 35: set vt+1 ← c and ut+1 ← ut + d (and optional
vt ← 35)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 2: Single neuron with synaptic input
Exercise 2
Simulate the neuron model for 1000ms and plot the resulting
voltage trace. Assume that 100 synapses are attached to the
neuron, with each pre-synaptic neuron firing with a Poisson
process of rate frate = 2 Hz between time 200ms and 700ms.
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 2: Single neuron with synaptic input
Exercise 2
Simulate the neuron model for 1000ms and plot the resulting
voltage trace. Assume that 100 synapses are attached to the
neuron, with each pre-synaptic neuron firing with a Poisson
process of rate frate = 2 Hz between time 200ms and 700ms.
Synaptic input model:
Isyn
=
X
g˙ jin
=
gjin /τg
gjin
←
gjin + 1,
wjin gjin (t)(Ejin − v (t))
j
if rj (t) < frate ∆t
with τg = 10ms, weights wjin = 0.07,
Ej = 0, rj (t) ∈ [0, 1] uniform
random numbers drawn for each
step t, and j = 1 . . . 100.
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 2 in detail:
Use the last script, save it under a new file name, and add the necessary lines.
Proceed as follows:
1
Initialize new parameter values (τg = 10, frate = 0.002ms−1 )
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 2 in detail:
Use the last script, save it under a new file name, and add the necessary lines.
Proceed as follows:
1
Initialize new parameter values (τg = 10, frate = 0.002ms−1 )
2
Reserve memory and initialize gin = (gjin ), win = (wjin ), and
E = (Ej ) (vectors of length nin = 100) with constant elements
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 2 in detail:
Use the last script, save it under a new file name, and add the necessary lines.
Proceed as follows:
1
Initialize new parameter values (τg = 10, frate = 0.002ms−1 )
2
Reserve memory and initialize gin = (gjin ), win = (wjin ), and
E = (Ej ) (vectors of length nin = 100) with constant elements
3
Inside the for-loop change/add the following:
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 2 in detail:
Use the last script, save it under a new file name, and add the necessary lines.
Proceed as follows:
1
Initialize new parameter values (τg = 10, frate = 0.002ms−1 )
2
Reserve memory and initialize gin = (gjin ), win = (wjin ), and
E = (Ej ) (vectors of length nin = 100) with constant elements
3
Inside the for-loop change/add the following:
1
set pj = 1 if r ≤ frate ∆t (otherwise 0) in the case when i∆t is
between 200 and 700 (otherwise 0). r = (rj ) is a vector of
uniform random numbers of length nin
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 2 in detail:
Use the last script, save it under a new file name, and add the necessary lines.
Proceed as follows:
1
Initialize new parameter values (τg = 10, frate = 0.002ms−1 )
2
Reserve memory and initialize gin = (gjin ), win = (wjin ), and
E = (Ej ) (vectors of length nin = 100) with constant elements
3
Inside the for-loop change/add the following:
1
set pj = 1 if r ≤ frate ∆t (otherwise 0) in the case when i∆t is
between 200 and 700 (otherwise 0). r = (rj ) is a vector of
uniform random numbers of length nin
2
before the vt update: implement the conductance dynamics g
and set Iapp (using array notation HowTo carefully HowTo ):
gjin
Iapp
gjin
←
←
←
gjin + pj
win · gin ⊙ Ein − win · gin ⊙ vt
(1 − ∆t/τg ) gjin
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3: Simulate 1000 neurons (not inter-connected)
Exercise 3
Simulate 1000 neurons for 1000 ms and plot the resulting spikes.
Assume that each neuron receives (random) 10% of the 100
Poisson spike trains of rate frate = 2 Hz between time 200 ms and
700 ms. Note that the neurons are not yet inter-connected.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3: Simulate 1000 neurons (not inter-connected)
Exercise 3
Simulate 1000 neurons for 1000 ms and plot the resulting spikes.
Assume that each neuron receives (random) 10% of the 100
Poisson spike trains of rate frate = 2 Hz between time 200 ms and
700 ms. Note that the neurons are not yet inter-connected.
Excitatory and inhibitory neurons:
A neuron is, with probability pI = 0.2,
a (fast-spiking) inhibitory neuron
(a = 0.1, d = 2), others are (regular
spiking) excitatory neurons (a = 0.02
and d = 8). Weights of the input
synapse j to inhibitory neuron i is
w in = 0.07 if connected (otherwise 0).
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 3 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize new parameter values (n = 1000)
Step 2
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize new parameter values (n = 1000)
2
Initialize 2 logical vectors (for indexing HowTo ) kinh and kexc of
length n, where kinh has a 1 in element i with probability p = 0.2
(marking an inhibitory neuron) and 0 otherwise. kexc = ¬kinh .
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize new parameter values (n = 1000)
2
Initialize 2 logical vectors (for indexing HowTo ) kinh and kexc of
length n, where kinh has a 1 in element i with probability p = 0.2
(marking an inhibitory neuron) and 0 otherwise. kexc = ¬kinh .
3
Reserve memory and initialize vij , uij (now being n × T matrices).
Set parameters ai and di according to kexc and kinh .
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize new parameter values (n = 1000)
2
Initialize 2 logical vectors (for indexing HowTo ) kinh and kexc of
length n, where kinh has a 1 in element i with probability p = 0.2
(marking an inhibitory neuron) and 0 otherwise. kexc = ¬kinh .
3
Reserve memory and initialize vij , uij (now being n × T matrices).
Set parameters ai and di according to kexc and kinh .
4
The weights wijin = 0.07 now form a n × nin matrix. Set 90 %
random elements to 0 to account for the connection probability.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize new parameter values (n = 1000)
2
Initialize 2 logical vectors (for indexing HowTo ) kinh and kexc of
length n, where kinh has a 1 in element i with probability p = 0.2
(marking an inhibitory neuron) and 0 otherwise. kexc = ¬kinh .
3
Reserve memory and initialize vij , uij (now being n × T matrices).
Set parameters ai and di according to kexc and kinh .
4
The weights wijin = 0.07 now form a n × nin matrix. Set 90 %
random elements to 0 to account for the connection probability.
5
Inside the for-loop change/add the following:
1
Same update equations (for vi,t+1 and ui,t+1 ) but use array
notation to update all i neurons simultaneously.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Step 3 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize new parameter values (n = 1000)
2
Initialize 2 logical vectors (for indexing HowTo ) kinh and kexc of
length n, where kinh has a 1 in element i with probability p = 0.2
(marking an inhibitory neuron) and 0 otherwise. kexc = ¬kinh .
3
Reserve memory and initialize vij , uij (now being n × T matrices).
Set parameters ai and di according to kexc and kinh .
4
The weights wijin = 0.07 now form a n × nin matrix. Set 90 %
random elements to 0 to account for the connection probability.
5
Inside the for-loop change/add the following:
1
6
Same update equations (for vi,t+1 and ui,t+1 ) but use array
notation to update all i neurons simultaneously.
Plot the spike raster. Plot black dots at {(t, i)|vit ≥ 35} for
excitatory neuron i. Use red dots for inhibitory neurons.
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4: Simulate recurrent network
Exercise 4
Simulate 1000 neurons as before but with added recurrent
connections.
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4: Simulate recurrent network
Exercise 4
Simulate 1000 neurons as before but with added recurrent
connections.
Recurrent synaptic activations
A neuron i is sparsely (with probability
prc = 0.1) connected to a neuron j.
Thus neuron i receives an additional
current Iisyn of the form:
Iisyn
=
n
X
wij gj (t) (Ej − vi (t))
j=1
Weights are Gamma distributed (scale
0.003, shape 2). Inh. to exc.
connections are twice as strong.
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 4 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize and allocate memory for the new variables ( g = (gj ), Ej ).
Set Ej = −85 if j is an inhibitory neuron (otherwise 0).
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 4 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
2
Initialize and allocate memory for the new variables ( g = (gj ), Ej ).
Set Ej = −85 if j is an inhibitory neuron (otherwise 0).
Reserve memory and initialize weights W = (wij ) to zero.
Randomly choose 10% of the matrix elements.
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 4 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize and allocate memory for the new variables ( g = (gj ), Ej ).
Set Ej = −85 if j is an inhibitory neuron (otherwise 0).
2
Reserve memory and initialize weights W = (wij ) to zero.
Randomly choose 10% of the matrix elements.
3
Set the chosen weight matrix elements to values drawn from a
Gamma distribution of scale 0.003 and shape 2 HowTo
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 4 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize and allocate memory for the new variables ( g = (gj ), Ej ).
Set Ej = −85 if j is an inhibitory neuron (otherwise 0).
2
Reserve memory and initialize weights W = (wij ) to zero.
Randomly choose 10% of the matrix elements.
3
Set the chosen weight matrix elements to values drawn from a
Gamma distribution of scale 0.003 and shape 2 HowTo
4
Make the weight matrix “sparse”
HowTo
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 4 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize and allocate memory for the new variables ( g = (gj ), Ej ).
Set Ej = −85 if j is an inhibitory neuron (otherwise 0).
2
Reserve memory and initialize weights W = (wij ) to zero.
Randomly choose 10% of the matrix elements.
3
Set the chosen weight matrix elements to values drawn from a
Gamma distribution of scale 0.003 and shape 2 HowTo
4
Make the weight matrix “sparse”
5
Scale weights from inh. to exc. neurons by the factor of 2
HowTo
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Step 4 in detail:
Modify the last script (after saving it under new name).
Proceed as follows:
1
Initialize and allocate memory for the new variables ( g = (gj ), Ej ).
Set Ej = −85 if j is an inhibitory neuron (otherwise 0).
2
Reserve memory and initialize weights W = (wij ) to zero.
Randomly choose 10% of the matrix elements.
3
Set the chosen weight matrix elements to values drawn from a
Gamma distribution of scale 0.003 and shape 2 HowTo
4
Make the weight matrix “sparse”
5
Scale weights from inh. to exc. neurons by the factor of 2
6
Inside the for-loop change/add the following:
1 add the equations for recurrent conductances gj
I
gj
syn
gj
←
←
←
HowTo
gj + 1, if vj (t − 1) ≥ 35
W · (g ⊙ E) − (W · g) ⊙ v
(1 − ∆t/τg ) gj
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Congratulation !
You have just coded and simulated a
quite realistic network model !
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Optional exercises
Step 3
Step 4
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Exercises
Exercise (p-series)
Calculate the p-series (generalization of the Harmonic Series) for a
given p up to a given m
µm =
m
X
1
np
i=1
Use array notation.
Advise: Use array notation and avoid for-loops whereever you can!
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Exercises
Exercise (blob movie)
1
Generate a two arrays with, x and y , ranging from -2 to 2
(and about 100 elements)
2
Generate two 100 by 100 grid-matrices, X and Y using
meshgrid with x and y as input (look at X and Y to
understand what meshgrid is doing).
3
calculate a matrix Z of the same size as X and Y where each
2
2
element is given by zi = e −(xi +yi ) .
write a loop with t ranging from 0 to 1000 where you
4
1
2
3
plot the matrix Z (using imagesc)
circular shift the matrix Z by 1 element (using circshift)
force the plot to be drawn (using drawnow)
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Exercise (Logical indexing and basic plotting)
1
Generate a 100 by 2 matrix M of Gaussian random values
2
Plot a point cloud (first vs. second column of M) using blue
circles
3
Calculate how many values of M are positive
4
Erase all rows in M which have at least one negative value
5
Plot the points of the modified array M using red crosses in
the same graph as above.
6
Set both axes limits to −3 to 3 and plot dotted gray lines on
the coordinate axes (where y = 0 or x = 0).
7
Label the axes and write a legend
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Exercises
Exercise (Poisson spike trains)
Write a function that generates a homogeneous Poisson spike train
of rate λ having exactly N spikes. Use array notations (cumsum).
Hint: Poisson spike intervals t are exponentially distributed. They
can be generated by inverse transform sampling: Given uniform
random variables u in 0 to 1 valid inter-spike intervals can be
calculated as
log u
t=−
λ
Advise: Use array notation and avoid for-loops whereever you can!
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Exercise (More on spike trains)
Generate a long Poisson spike train (with the function from the
last exercise). Compute the mean and standard deviation of the
inter-spike interval distribution.
Further, write a function that counts the number of spike times
falling into time-bins of length ∆t.
Hint: Use diff to get the intervals from spike times.
Hint 2: Use histc to bin the spike-times
Exercises
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Exercises
Exercise (Plot 2-D Gaussian point cloud)
Write a function which plots a point cloud of n random
samples drawn from a 2D Normal distribution (with 0
mean) and variance Σ = (RD)T (RD), where the
rotation matrix is defined as
cos θ − sin θ
R=
sin θ
cos θ
5
0
-5
-5
0
5
and D is a diagonal matrix of the standard deviation in
the principal directions
σ1 0
D=
0 σ2
The function should accept parameters σ1 , σ2 , and θ.
Hint: use randn to produce independent Normal-distributed random vectors xi
and transform them according to yi = RDxi .
Introduction
Arrays
Coding
Plotting
Network model:
Step 1
Step 2
Step 3
Step 4
Exercises
Advanced exercise
Exercise (Optimizing Matlab code)
Optimize “bad” Matlab code for generating self-organizing maps.
See provided code and description in som_exercise.zip.
Howtos
How to start Matlab
Easy: Just click on Matlab symbol...
back
Howtos
How to get help
Each Matlab function has
a header which describes its
usage
Just type:
>> help command
or
>> doc command
Alternatively:
http://www.mathworks.com/access/helpdesk/help/techdoc/
back
Howtos
How to plot in general
How to plot:
1
Open a figure (window) with
>> figure;
2
Issue plotting command of your choice
>> plot(x,y,’b--’,’LineWidth’,2);
HowTo
Note: graphics commands draw into the current axes (gca) on the
current figure (gcf). They automatically create a new figure and axes if
necessary.
3
Make the plot nicer by adding labels and setting limits, eg:
>>
>>
>>
>>
back
xlabel(’Time [s]’);
ylabel(’Response [spks/sec]’);
xlim([-1,1]); ylim([-2,2]);
title(’Simulation’)
Howtos
How to use the plot command
Basic syntax:
handle = plot(X,Y,linespec,optname1,val1,...);
X,Y – x- and y-values to plot. If Y is a 2-D array, all columns are
plotted as different lines
linespec – a string with a short hand for color and line style and
marker type. See help plot for an overview. Eg,
linespec = ’:ko’
plots dotted (:) black line (k) with a circle at each given
coordinate (xi , yi )
optname1 – a string specifing the option name, eg. ’LineWidth’
val1 – a corresponding value.
handle – graphics handle for later reference, eg. with
>> set(handle,optname1,val1)
Tip: To get an overview over all possible options, try get(handle)
back
Howtos
How to initialize parameters
Just use the syntax
>> parname = value;
Example
>> a = 2;
>> vreset = 0;
>> tau = 0.02
tau =
0.02
back
Howtos
How to use scalar expressions
Binary operations: work as expected, use = + - * / ^
Example (compute y =
a2 x
2+a
+ b)
>> a = 2;
>> b = 1;
>> x = 0.5;
>> y = a^2*x/(2+a) + b;
>> y
y =
1.500
back
Howtos
How to initialize arrays
1
2
Implicitly, using function returning an array
By explicit concatenation
Concatenation of columns of arrays
[ arr1, arr2, ... , arrn]
Concatenation of rows of arrays
[ arr1; arr2; ... ; arrn]
Note: an scalar is also regarded as an array (of size [1,1]).
Note 2: arrays must have matching sizes.
Example (Concatenation)
>> A = [1,2;3,4]
A =
1
2
3
4
>> B = [A,A]
B =
1
2
1
2
3
4
3
4
>> C = [A;A]
C =
1
2
3
4
1
2
3
4
Howtos
How to pre-allocate memory
Functions for pre-allocating memory include:
colon (:) – for linear sequences
zeros – all zero array of given size
ones – all one array of given size
rand – random array of given size (equally in [0, 1])
To improve performance arrays should always be pre-allocated!
Example (Functions initializing arrays)
>> A = zeros(3,3);
>> A = ones(4,4,4);
>> size(A)
ans =
4
4
4
back
>> x =
x =
3.0
>> A =
1
1
3:-0.5:1
2.5 2.0
ones(2)
1
1
1.5
1.0
Howtos
How to write an m-file script?
1
2
3
open a text-editor of your choice or use the editor provided
with Matlab
>> edit myscript
Write all your calculations in a text-file with extension “.m”
(here myscript.m)
Save file in the current working directory
(or addpath to search path )
4
Call your script by calling it from the “Command Window”
>> myscript;
Example (myscript.m)
1
2
3
%t h i s i s my f i r s t s c r i p t . I t d i s p l a y s a random number
random number = rand ( 1 ) ;
f p r i n t f ( ’A random number : %1.4 f : \ n ’ , random number ) ;
back
Howtos
How to use basic syntax: if-clause
if-else block syntax:
if scalar condition
expressions
else
expressions
end
1
2
3
4
5
Relational operators, eg.: == (equals), || (or), && (and), ~ (not)
for details type: help relop
Example (if-else)
a = rand ( 1 ) ;
i f a == 0 . 5
f p r i n t f ( ’ you a r e v e r y l u c k y ! ’ ) ;
end
back
Howtos
How to use a for-loop
for-loop block syntax:
for i = array
% i==a r r a y ( j ) i n t h e j −t h l o o p
expressions
end
1
2
3
4
(one can also use break and continue keywords)
Example (for loop)
a =0;
for i = 1:100
a = a+i ;
end
back
Howtos
How to index arrays
1
Subscript of a matrix:
access the (i, j)-th element of a 2D-matrix W of dimension (m, n)
>> W(i,j) = 1
Note: The first index is always 1 (not 0 as in most other languages)
2
Linear index of a matrix:
access the (i, j)th element of the 2D-matrix W of dimension (m, n)
>> linearidx = i + (j-1)*m;
>> W(linearidx)
3
“Slice” indexing with “:”
access the i-th row and jth column in W , respectively
>> wi = W(i,:)
>> wj = W(:,j)
get all elements as a concatenated column vector
>> W(:)
back
Howtos
How to index arrays (2)
4
Multiple indices
vectors of linear indices can be used
>> W([1, 4, 5 ,6])
access the 1st to 4th rows of the 2D-matrix W of dimension (m, n)
>> W(1:4,:)
access the 2nd (m,n)-slice of a 3D-matrix V of dimension (m, n, p)
>> W(:,:,2)
5
Logical indexing
logical matrices of the same size as W can be used as index (very
useful)
>> W(W>0)
>> W(find(W>0)) = 1
back
Howtos
Calculating with arrays
1
Element-wise interpretation
For instance, sin cos log etc.
Reserved symbols, .* ./ .^
2
“true” matrix interpretation (with dot product)
3
Operations on one specified dimensions of the matrix
4
Array manipulations
Symbols * / ^ etc.
For instance, sum mean max etc.
eg. reshape repmat permute circshift tril
Example (element-wise product and dot product)
>> A = ones(2,2);
>> A.*A
ans =
1
1
1
1
back
>> A*A
ans =
2
2
2
2
Howtos
Calculating with arrays is
straightforward
however, carefully check
the size of matrices
if element-wise or matrix-like operations are intended
which matrix dimension to operate on
Example (compute yi = W xi + b with b, xi ∈ R2 )
>> W = [1,0.2;0.4,1];
>> b = [1;2] + 0.1;
>> x = 2*randn(2,1);
>> y = W * x + b
y =
4.5535
1.4856
back
>> N = 5;
% same as b = [b,b,b,b,b]
>> bi = repmat(b,[1,N]);
>> xi = 2*randn(2,N);
>> xi(:,1) = x;
>> yi = W * xi + bi
yi =
4.5535 2.3126 [..] -0.9021
1.4856 6.8091 [..] -0.1080
Howtos
Gamma distribution
The Gamma probability density function is defined as
p(x|k, θ) =
x
1
x k−1 e − θ
k
Γ(k)θ
with shape k and scale θ and x ∈ [0, ∞).
Example (Gamma random numbers)
>> shape = 2;
>> scale = 0.003;
>> n = 1e7;
>> r = gamrnd(shape,scale,n,1);
>> size(r)
ans =
1000000
1
>> hist(r,1000)
back
Howtos
How to use sparse matrices
In Matlab it is often more efficient (faster) to use sparse matrices
instead of regular matrices if the majority of matrix elements are 0.
To generate a sparse matrix:
>> W = sparse(W);
In general, the syntax for using sparse matrices is the same as for
regular matrices.
Example (Sparse matrix)
>> W = double(rand(100,100)>0.9);
>> Ws = sparse(W);
>> y = Ws*x; % the same as y=W*x but faster
back
Howtos
Solution to Excercise #1
1
%compute t h e S t i r l i n g Formula
2
3
n = 50;
4
5
n f a c t o r i a l = s q r t ( 2 ∗ p i ∗n ) ∗ ( n/ exp ( 1 ) ) ˆ n
back to text
Howtos
Solution to Excercise #2
1
2
function n f a c t o r i a l = s t i r f a c (n ) ;
%compute t h e S t i r l i n g Formula
3
4
n f a c t o r i a l = s q r t ( 2 ∗ p i ∗n ) ∗ ( n/ exp ( 1 ) ) ˆ n ;
back to text