2013/den13: SKIM.m

File SKIM.m, 3.9 KB (added by jtapson, 5 years ago)

Basic SKIM network script

Line 
1% SKIM algorithm to find spatio-temporal pattern using synaptic kernels
2% assumes an input array X with inputs in columns, and target values in last
3% column
4% Written (if you could call it that) by Jonathan Tapson at Telluride
5% This version 2 July 2013
6%
7% Assumes a data array X with row format (input neuron 1, input neuron 2,..
8%                                      ...output neuron)
9%
10% Only one output in this code; trivial to extend it
11% Time series is column-wise, each row is one time sample
12%
13
14% First separate X into inputs Xi and output Y
15
16sizex=size(X);              % size of data array
17Y = X(:,sizex(2));          % output array (last column of X)
18Xi=X(:,1:(sizex(2)-1));     % input array
19
20[l,w]=size(Y);              % size of output
21[l,W]=size(Xi);             % size of input
22Yo = zeros(l,1);            % continuous output series 
23Stn=1;                      % function is stationary (set to 0 for adaptive)
24LearningRate=1;             % if nonstationary, you may set learning rate
25                            % Number of neurons to which to spread input:
26HiddenLayerSize=400;        % hidden layer size
27H=zeros(HiddenLayerSize,1);                         % initial value of neurons
28HiddenLayerWeights=rand(W,HiddenLayerSize)-0.5;     % initialise with
29                                                    % random weights, zero mean
30
31max_tc=200;                 % maximum time constant                                                     
32max_alpha=200;              % maximum time length of synaptic function
33Alpha_ts=max_tc*rand(HiddenLayerSize,1);           % randomised synaptic time constants
34
35
36% NB hidden layer weights must be scaled so that function pushes hidden
37% outputs well into nonlinear region.
38
39Weights=zeros(HiddenLayerSize,w);   % initial value of linear weights
40P=eye(HiddenLayerSize);             % initialise correlation matrix
41E=zeros(l,w);                       % Error record
42
43Hlin = zeros(l,HiddenLayerSize);    % array to store linear portion of synaptic inputs
44Ha = zeros(l,HiddenLayerSize);      % array to store alpha outputs at synapses
45
46scale=1;                            % may need to scale inputs to push into nonlinearity
47Hlin = scale*Xi*HiddenLayerWeights;     % calculate linear inputs (should be spikes)
48Hlin = 1 ./ (1 + exp(-5.0*Hlin))-0.5;   % sigmoid activation, zero mean
49for i=1:(l-max_alpha)               % now calculate synaptic function contributions
50    for j=1:HiddenLayerSize
51        if (Hlin(i,j)~=0)           % if spike
52            Ha(i:(i+max_alpha-1),j)=Ha(i:(i+max_alpha-1),j)...
53                +Hlin(i,j).*alpha(Alpha_ts(j),max_alpha);       % alpha function
54                         %.*del_gauss(Alpha_ts(j),max_alpha);  % delay plus gaussian
55                         %.*del_alpha(Alpha_ts(j),max_alpha);  % delay plus alpha
56                         %.*res_den(Alpha_ts(j),max_alpha);  % resonant
57                         %synapse
58               
59        end
60    end
61end
62
63% Now caculate error
64
65for i=1:l
66
67    H=Ha(i,:)';
68    y=Weights'*H;                   % output value
69    Yo(i,1)=y;                      % output is saved as the function output
70    E(i,1)=Y(i,1)-y;                % calculate error
71           
72   % Adapted Greville's Method (not as fast as possible, see AvS
73   
74    psi=P*H;                        % measure of hidden layer activity
75    nrm=real(sqrt(1+H'*P*H));       % normalisation of activity
76    e=E(i,1)*nrm;                   % normalise error
77    Weights=Weights+(e/nrm)*(psi/nrm); % update weights
78    P=P-(psi/nrm)*(psi'/nrm);       % update correlations
79    if (Stn==0)                     % if non-stationary data
80        P=P+(1/(HiddenLayerSize/LearningRate)...
81            *((1-exp(-abs(e))).*eye(HiddenLayerSize)));  % if data nonstationary,
82                                         % adapt correlation matrix
83                                         % according to error
84    end
85end