%一個最基本的rbf算法，學習算法采用偽逆函數
%
%
%
%
%
%
%
%
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
clc;
%定義用到的變量
r=1;                  %領域半徑，人為設定
nodenum=0;            %隱層節點數，由實際數據與r共同決定
inputdata=[];         %輸入矩陣
inputpath=' ';        %存儲原始輸入數據的路徑
nodeout=[];           %隱層輸入陣　
netout=[];            %網絡輸出陣
weight=[];            %輸出權值陣，也是唯一的權值
inputnum=0;           %輸入維數
outputnum=0;          %輸出維數
center=[];            %聚類中心
numtrain=0;           %學習樣本的個數
row=0;                %學習樣本的個數
simrow=0;             %全部樣本的個數
numtest=0;            %泛化樣本的個數
strength=1;           %歸一化處理時用到的范圍，一般為1
yout=[];              %輸出的期望值

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%變量的初始化
r=1;
inputnum=9;
inputpath='e:\yinjia\data\yearsun.dat';
outputnum=1;
simrow=290;
source=load(inputpath);
%保存一個source的副本，為反歸一化做準備
copysource=source;
%歸一化處理
source=normalize(source,strength);
yout=source(inputnum+1:simrow+inputnum);
inputdata=phasespace(source,inputnum,simrow);
row=250;
numtrain=row;
numtest=simrow-row;
%把第一人輸入變量定為初始的中心
%中心的位置在初始化后，不會再變化
center=inputdata(1,:)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%開始進行簡單的聚類
%其基本思想是：給定一個固定的鄰域半徑，在中心點鄰域半徑內的點就算作
%是該鄰域的點；半徑外的點，取第一個不是半徑內的點做為新的中心點
%設計一個函數iscenter,來計算一個點是否為中心點。
for step=2:row
    if iscenter(inputdata(step,:)',center,r)
        center=[center inputdata(step,:)'];
        nodenum=nodenum+1;
    end % this end for iscenter(inputdata(step,:))
end % this end for step=2:row
%聚類完成.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%準備進行神經網絡的前向計算
[centerrow nodenum]=size(center);
%網絡的初始化
nodeout=zeros(row,nodenum);
netout=zeros(row,outputnum);
weight=zeros(nodenum,outputnum);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%網絡的計算
for step=1:row
    for step1=1:nodenum
        nodeout(step,step1)=GaussRadialBasisFunction(inputdata(step,:)',...
            center(:,step1),r);
    end % this end for step1=1:nodenum
end %this end for step=1:row
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用偽逆來計算權值
weight=pinv(nodeout)*yout(1:250);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%進行仿真
for step=1:simrow
    for step1=1:nodenum
        nodeout(step,step1)=GaussRadialBasisFunction(inputdata(step,:)',...
            center(:,step1),r);
    end % this end for step1=1:nodenum
end %this end for step=1:simrow

mydata=nodeout*weight;
%反歸一化
mydata=unnomal(mydata,copysource,strength);
rldata=copysource(inputnum+1:inputnum+simrow);
plot(rldata);hold on;plot(mydata,'r');

%計算一個評價函數
rmsetest=(norm(rldata(numtrain+1:simrow)-mydata(numtrain+1:simrow))...
    ^2/(numtest-1))^0.5