脉冲神经网络之Tempotron(二)

                        脉冲神经网络之Tempotron代码示例

        上一篇从原理的角度大致介绍了脉冲神经网络的神经元模型以及Tempotron监督学习方法,这一章记录了Tempotron的代码实现。这份代码是使用matlab编写,用脉冲神经网络实现对26个字母分类,我们从细节去解释该代码:

function TempotronClassify()
% Tempotron: a neuron that learns spike timing-based decisions
% Rober Gutig 2006 Nature Neuroscience
clear; clc;
NumImages=26;
for i=1:NumImages
    ImageName=strcat('Icon16X16\Letter-',char('A'+i-1),'-black-icon');% 从icon16X16文件夹中读取所有图片
    ImageMatrix=imread(ImageName,'bmp');% 读取图片为灰度图,保存在矩阵中,取反
    ImageMatrix=~ImageMatrix;  % make the white pixel be 0, and black be 1;
    TrainPtns(:,i)=image2ptn(ImageMatrix);
end
上面代码片段就是从Icon16X16\文件夹中读取图片文件,图片从A到Z,这个代码的图片数据很简单,只有0和1,由黑白两色构成字母,读取的数据以16*16矩阵形式保存,之后调用image2ptn()方法,将矩阵转换。这个方法实现的就是对外界刺激进行编码,也就是将图片数据转换成脉冲序列的形式。打开该方法查看代码:

%% convert a image to a spike train saved in a vector
RandParts=1;

spikeTrain=zeros(32,1);
if RandParts==1
    loadR=1;
    AR=A(:);
    if loadR==0
        R = randperm(size(A,1)*size(A,2));
        save('RandIndex','R');
    else
        load('RandIndex','R');
    end
    numRandNeus=32;
    for i=1:numRandNeus
        IndexR=R((1+(i-1)*8):(8+(i-1)*8));
        Bits=AR(IndexR);
        Bits=Bits';
        spikeTime=bi2de(Bits,'left-msb');    % 二进制转十进制
        if spikeTime==0
            spikeTime=2^8;  % put 0 to the end of the interval
        end
        spikeTrain(i)=spikeTime;
    end
end
这里的rand函数目的只是打乱16*16矩阵中的行,我是不太理解这个=-=,AR是个256*1的数组,是将A(图片矩阵16*16)展平(flatten)的结果,Bits=AR(IndexR);这句简单的说就是从AR中随机抽取8个数字,构成bits,然后将bits转换成spikeTime(十进制数字),代表脉冲发放的时间。综上原先的图片16*16矩阵被转换成32*8矩阵,8代表8位二进制编码,范围在0~255之间,也就是说时间窗口大小为256毫秒,然后将8为二进制数转为十进制,代表脉冲发放的时间,该时间在(0~255)之间。总结起来,就是对于一张16*16的黑白图片,转换成了一个32*1的数组,数组中有32个时间窗口,每个元素代表脉冲发送时间。
TrainPtns=TrainPtns*1e-3;  % scale to ms
nAfferents = size(TrainPtns,1);
nPtns = NumImages;
%nOutputs = 1;    %%%%%%%%%%   1
nOutputs = 5;

TrainPtns是32*26的矩阵,也就是将所有转换后的spikeTrain拼接一起。nAfferents代表了输入脉冲数量,这里就是32,由这32个输入确定一个输出(字母)。输出就定义为5,也就是5个二进制数表示一个字母的编号。

loadData=0;% 是否载入已保存的模型

V_thr = 1; V_rest = 0;
T = 256e-3;         % pattern duration ms
dt = 1e-3;
tau_m = 20e-3; % tau_m = 15e-3;???
tau_s = tau_m/4;
% K(t?ti)=V0(exp[?(t?ti)/τ]–exp[?(t?ti)/τs])
aa = exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s);

V0 = 1/max(exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s));
lmd = 2e-2;%1e-2/V0;   % optimal performance lmd=3e-3*T/(tau_m*nAfferents*V0)  1e-4/V0;
maxEpoch = 200;
mu = 0.99;  % momentum factor
上面定义了阈值电压V_thr,这是按经验来取,这里取1,复位电压则取0.定义了时间窗口大小为256毫秒,单位时间是1毫秒,aa则是为了核函数K的计算简化定义。即预先定义公式:
K(t-ti)=V0(exp[-(t-ti)/τ]–exp[-(t-ti)/τs])
V0也被预先定义。最大训练次数为200次(maxEpoch),lmd类似于学习率,而mu则是动量,该代码用了动量优化模型。

if loadData ==0 %初始化网络
    weights = 1e-2*randn(nAfferents,nOutputs);  % 1e-3*randn(nAfferents,1);
    save('weights0','weights');
else
    load('weights0','weights');
end
随机初始化权重,或者读取保存的权重文件。
Class = de2bi(1:26,'left-msb'); Class=Class';

correctRate=zeros(1,maxEpoch);
dw_Past=zeros(nAfferents,nPtns,nOutputs);  % momentum for accelerating learning.上一个权重的更新,用于动量计算
for epoch=1:maxEpoch    
    
    Class_Tr = false(nOutputs,nPtns);  % actual outputs of training
    for pp=1:nPtns                 
        for neuron=1:nOutputs
            Vmax=0; tmax=0;
            fired=false;        
            Vm1=zeros(1,256); indx1= 1; % trace pattern 1
            for t=dt:dt:T
                Vm = 0; 
                if fired==false
                    Tsyn=find(TrainPtns(:,pp)<=t+0.1*dt);    % no cut window
                else
                    Tsyn=find(TrainPtns(:,pp)<=t_fire+0.1*dt); % shut down inputs
                end
                if ~isempty(Tsyn)                    
                    A1=TrainPtns(:,pp);
                    A2=A1(Tsyn);
                    K =V0*(exp(-(t-A2)/tau_m)-exp(-(t-A2)/tau_s)); % the kernel value for each fired afferent
                    A1=weights(:,neuron);
                    firedWeights=A1(Tsyn);
                    Vm = Vm + firedWeights'*K ;
                end

                Vm = Vm + V_rest;
                if Vm>=V_thr && fired==false % fire
                    fired=true;
                    t_fire=t;
                    Class_Tr(neuron,pp)=true;
                end
                if Vm>Vmax
                    Vmax=Vm; tmax=t;
                end

                %if pp==1
                    Vm1(indx1)=Vm;
                    indx1=indx1+1;
                %end
            end

            %if pp==1
                figure(1); plot(dt:dt:T,Vm1);
                title(strcat('Image ',char('A'+pp-1),'; neuron: ',num2str(neuron))); drawnow;
            %end
            if Vmax<=0
                tmax=max(TrainPtns(:,pp));
            end
            
            if Class_Tr(neuron,pp)~=Class(neuron,pp) % error
                
                Tsyn=find(TrainPtns(:,pp)<=tmax+0.1*dt); 

                if ~isempty(Tsyn)                    
                    A1=TrainPtns(:,pp);
                    A2=A1(Tsyn);
                    K =V0*(exp(-(tmax-A2)/tau_m)-exp(-(tmax-A2)/tau_s)); % the kernel value for each fired afferent
                    A1=weights(:,neuron);
                    dwPst=dw_Past(:,pp,neuron);
                    if fired==false    % LTP
                        Dw=lmd*K;
                    else           % LTD
                        Dw=-1.1*lmd*K;
                    end
                    A1(Tsyn) = A1(Tsyn) + Dw + mu*dwPst(Tsyn);
                    weights(:,neuron)=A1;
                    dwPst(Tsyn)=Dw;
                    dw_Past(:,pp,neuron) = dwPst;
                end                
            end            
            
        end  % end of one neuron computation
        
   end % end of one image
   %CC=isequal(Class,Class_Tr);
   %correctRate(epoch)=sum(Class==Class_Tr)/length(Class);
   CC = bi2de(Class_Tr','left-msb');
end
save('TrainedWt','weights');
figure(2); plot(1:maxEpoch,correctRate,'-b.');
end
class代表最后的类别,由5位二进制编码表示。最外层循环是最大训练次数,第二层nPtns代表了26个样本,最后一层循环nOutputs为5,代表对每个输出神经元进行一次计算。之后将Vmax,Tmax(代表最大膜电位出现时间)初始化为0,fired表示是否发出脉冲,这里初始化为false,Vm1记录了此输出神经元膜电位(在一个时间窗口内)的变化,之后进入循环,循环一个时间窗口长度,对一个时间窗口内的每个单位时间,若该神经元已经发送过脉冲,则shut down所有脉冲输入;若该神经元还未发出过脉冲(fired = false), 就执行Tsyn=find(TrainPtns(:,pp)<=t+0.1*dt);也就是在当前样本(pp)的脉冲序列中,找到比时间t+0.1*dt小的脉冲发送时间,也就是在t这个时间点之前有脉冲输入发生,此时A2就保存了该脉冲发送的时间(也就是ti),之后就可调用K(t - ti)来计算此时的核函数;下一步就是找到发送这个脉冲的输入神经元和该输出神经元之间连接的权重数值,用K*Weights就代表了该脉冲输入对膜电位的影响,最后将该贡献值累加到Vm上,得到此刻膜电位大小。

对整个时间窗口计算完后,将Vm加上复位电压,如果该值超过了阈值电压,则发放脉冲;若该值超过了已有的Vmax,则更新Vmax。之后进入判断语句:

if Class_Tr(neuron,pp)~=Class(neuron,pp) 神经元的输出和实际输出不相符时,分为两个情况:

(1)实际输出为0,但发放了脉冲。

此时应该抑制脉冲发放,于是权重应该减小:

Dw=-1.1*lmd*K;

(2)实际输出为1,但没有发放脉冲:

此时应该增大神经元的刺激,权重应该增加:

Dw=lmd*K;
最后保存此次循环的权重。于是单层的脉冲神经网路监督学习就这样完成了。下面附上完整代码:

function TempotronClassify()
% Tempotron: a neuron that learns spike timing-based decisions
% Rober Gutig 2006 Nature Neuroscience
clear; clc;
NumImages=26;
for i=1:NumImages
    ImageName=strcat('Icon16X16\Letter-',char('A'+i-1),'-black-icon');% 从icon16X16文件夹中读取所有图片
    ImageMatrix=imread(ImageName,'bmp');% 读取图片为灰度图,保存在矩阵中,取反
    ImageMatrix=~ImageMatrix;  % make the white pixel be 0, and black be 1;
    TrainPtns(:,i)=image2ptn(ImageMatrix);
end
TrainPtns=TrainPtns*1e-3;  % scale to ms
nAfferents = size(TrainPtns,1);
nPtns = NumImages;
%nOutputs = 1;    %%%%%%%%%%   1
nOutputs = 5;

loadData=0;% 是否载入已保存的模型

V_thr = 1; V_rest = 0;
T = 256e-3;         % pattern duration ms
dt = 1e-3;
tau_m = 20e-3; % tau_m = 15e-3;???
tau_s = tau_m/4;
% K(t?ti)=V0(exp[?(t?ti)/τ]–exp[?(t?ti)/τs])
aa = exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s);

V0 = 1/max(exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s));
lmd = 2e-2;%1e-2/V0;   % optimal performance lmd=3e-3*T/(tau_m*nAfferents*V0)  1e-4/V0;
maxEpoch = 200;
mu = 0.99;  % momentum factor
% generate patterns (each pattern consists one spik-e per afferent)

if loadData ==0 %初始化网络
    weights = 1e-2*randn(nAfferents,nOutputs);  % 1e-3*randn(nAfferents,1);
    save('weights0','weights');
else
    load('weights0','weights');
end
%Class = logical(eye(nOutputs));     % desired class label for each pattern
%Class = false(1,26); Class(26)=true;
Class = de2bi(1:26,'left-msb'); Class=Class';

correctRate=zeros(1,maxEpoch);
dw_Past=zeros(nAfferents,nPtns,nOutputs);  % momentum for accelerating learning.上一个权重的更新,用于动量计算
for epoch=1:maxEpoch    
    
    Class_Tr = false(nOutputs,nPtns);  % actual outputs of training
    for pp=1:nPtns 
 %       Class_Tr = false(nOutputs,1);  % actual outputs of training
                
        for neuron=1:nOutputs
            Vmax=0; tmax=0;
            fired=false;        
            Vm1=zeros(1,256); indx1= 1; % trace pattern 1
            for t=dt:dt:T
                Vm = 0; 
                if fired==false
                    Tsyn=find(TrainPtns(:,pp)<=t+0.1*dt);    % no cut window
                else
                    Tsyn=find(TrainPtns(:,pp)<=t_fire+0.1*dt); % shut down inputs
                end
                if ~isempty(Tsyn)                    
                    A1=TrainPtns(:,pp);
                    A2=A1(Tsyn);
                    K =V0*(exp(-(t-A2)/tau_m)-exp(-(t-A2)/tau_s)); % the kernel value for each fired afferent
                    A1=weights(:,neuron);
                    firedWeights=A1(Tsyn);
                    Vm = Vm + firedWeights'*K ;
                end

                Vm = Vm + V_rest;
                if Vm>=V_thr && fired==false % fire
                    fired=true;
                    t_fire=t;
                    Class_Tr(neuron,pp)=true;
                end
                if Vm>Vmax
                    Vmax=Vm; tmax=t;
                end

                %if pp==1
                    Vm1(indx1)=Vm;
                    indx1=indx1+1;
                %end
            end

            %if pp==1
                figure(1); plot(dt:dt:T,Vm1);
                title(strcat('Image ',char('A'+pp-1),'; neuron: ',num2str(neuron))); drawnow;
            %end
            if Vmax<=0
                tmax=max(TrainPtns(:,pp));
            end
            
            if Class_Tr(neuron,pp)~=Class(neuron,pp) % error
                
                Tsyn=find(TrainPtns(:,pp)<=tmax+0.1*dt); 

                if ~isempty(Tsyn)                    
                    A1=TrainPtns(:,pp);
                    A2=A1(Tsyn);
                    K =V0*(exp(-(tmax-A2)/tau_m)-exp(-(tmax-A2)/tau_s)); % the kernel value for each fired afferent
                    A1=weights(:,neuron);
                    dwPst=dw_Past(:,pp,neuron);
                    if fired==false    % LTP
                        Dw=lmd*K;
                    else           % LTD
                        Dw=-1.1*lmd*K;
                    end
                    A1(Tsyn) = A1(Tsyn) + Dw + mu*dwPst(Tsyn);
                    weights(:,neuron)=A1;
                    dwPst(Tsyn)=Dw;
                    dw_Past(:,pp,neuron) = dwPst;
                end                
            end            
            
        end  % end of one neuron computation
        
   end % end of one image
   %CC=isequal(Class,Class_Tr);
   %correctRate(epoch)=sum(Class==Class_Tr)/length(Class);
   CC = bi2de(Class_Tr','left-msb');
end
save('TrainedWt','weights');
figure(2); plot(1:maxEpoch,correctRate,'-b.');
end

%%将图片编码为脉冲序列并保存在向量中
function spikeTrain=image2ptn(A)  
%% convert a image to a spike train saved in a vector
RandParts=1;
% A1=A';
% B=[A1(:);A(:)];
% numPixels=length(B);
% numInputNeurons=numPixels/8; % 64 neurons
% spikeTrain=zeros(numInputNeurons,1);
% for i=1:numInputNeurons
%     Bits=B((1+(i-1)*8):(8+(i-1)*8));
%     Bits=Bits';
%     spikeTime=bi2de(Bits,'left-msb');    
%     if spikeTime==0
%         spikeTime=2^8;  % put 0 to the end of the interval
%     end
%     spikeTrain(i)=spikeTime;
% end

spikeTrain=zeros(32,1);
if RandParts==1
    loadR=1;
    AR=A(:);
    if loadR==0
        R = randperm(size(A,1)*size(A,2));
        save('RandIndex','R');
    else
        load('RandIndex','R');
    end
    numRandNeus=32;
    for i=1:numRandNeus
        IndexR=R((1+(i-1)*8):(8+(i-1)*8));
        Bits=AR(IndexR);
        Bits=Bits';
        spikeTime=bi2de(Bits,'left-msb');    % 二进制转十进制
        if spikeTime==0
            spikeTime=2^8;  % put 0 to the end of the interval
        end
        spikeTrain(i)=spikeTime;
    end
end
end



发布了8 篇原创文章 · 获赞 41 · 访问量 4万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览