%construct the phase space
load st1.dat;
s1=st1(8001:8800,3);
[c,l]=wavedec(s1,5,'db1');
[thr,sorh,keepapp]=ddencmp('den','wv',s1);
s=wdencmp('gbl',c,l,'db1',5,thr,sorh,keepapp);
s=wden(s1,'heursure','s','one',3,'sym8');
%First construct the input samples,namely:RPS
num=length(s);%total data length
K=1;%time delay
%c=zeros(5,16);
m=45;%embedding dimension
n=num-(m-1)*K;%Number of phase space state vectors
Sample=zeros(n,m);
for i=1:n
    for j=1:m
        Sample(i,j)=s(i+(j-1)*K);%reconstruction with time delay
    end
end
%Then construct target set
target=Sample(:,end);
%train SVM
[nsv, beta, bias,H,svi]  = svr(Sample(1:250,:),target(2:251),'rbf',10000,'einsensitive',0.0001);
%check the effect
%let's see the approximation
%reconstruct the identification result
TST_target= H(:,svi)*beta(svi) +bias;
nnt=1:250;
figure;
plot(nnt(svi),TST_target(svi),'r+',nnt,target(1:250),'bs-',nnt,TST_target,'ko--');
legend('support vectors','original data','approximation result');
xlabel('Time Series');
ylabel('Pressure');
%Now it's forecasting applications
%svm+one step
%construct a sample
%%  single step forcasting
for k=251:307
    XT=Sample(k,:);
    svm(k-250) = svroutput(Sample(1:250,:),XT,'rbf',beta,bias);
end

%Now rvm test
apha=0.5;
X_rvm=Sample(1:250,:);
Y_rvm=target(2:251);
[n_rvm,m_rvm]=size(X_rvm);
h_rvm=zeros(n_rvm,n_rvm+1);
h_rvm(1:n_rvm,1)=1;
for i=1:n_rvm
    for j=2:n_rvm+1
        dbtmp=norm(X_rvm(i,1:m_rvm)-X_rvm(j-1,1:m_rvm));
        h_rvm(i,j)=exp(-0.5*dbtmp*dbtmp/(apha*apha));
    end
end    

sigma=1e-6;
beta_r=inv(sigma*eye(n_rvm+1)+h_rvm'*h_rvm)*h_rvm'*Y_rvm;
pbeta_r=beta_r;
delta=1e-3;
stime=cputime;
for i=1:10000
    htmp=norm(Y_rvm-h_rvm*beta_r);
    sigma=htmp*htmp/n_rvm;
    ss=diag(abs(beta_r));
    beta_r=ss*inv(sigma*eye(n_rvm+1)+ss*h_rvm'*h_rvm*ss)*ss*h_rvm'*Y_rvm;
    if i>1
        if norm(beta_r-pbeta_r)/norm(pbeta_r)<delta
        %if sigma<delta
            break
        end    
    end 
    pbeta_r=beta_r;
end    
fprintf('Execution time: %4.1f seconds\n',cputime - stime);

ind = find((abs(beta_r)>1e-4));
k=size(ind);

for i=251:307
    XT_rvm= Sample(i,:);
    dbtmp=0;
    for j=1:k
        mm=ind(j);
        if mm==1
            dbtmp=dbtmp+beta_r(mm);
        else
            tmpp=norm(XT_rvm-X_rvm(mm-1,1:m_rvm));
            dbtmp=dbtmp+beta_r(mm)*exp(-0.5*tmpp*tmpp/(apha*apha));
        end    
    end
    rvm_single(i-250) = dbtmp;
end    
svm_error=mse(svm'-target(251:307))
rvm_error=mse(rvm_single'-target(251:307))



%Now check the effect
figure;
plot(251:307,svm,'bs-',251:307,target(251:307),'ro--',251:307,rvm_single,'^g-')
legend('one step forecasting result','original data')
xlabel('Time Series');
ylabel('Pressure');
%Now see multi step forecasting effect
%first construct the samples
step=2;
for i=251:307
    XT=Sample(i-step+1,:);% time window move backward
    for j=1:step
        y_SVM=svroutput(Sample(1:250,:),XT,'rbf',beta,bias);
        XT=[y_SVM Sample(i-step+1,2:end)];
    end
    Multi_SVM(i-250)=y_SVM;
end

figure;
    plot(251:307,Multi_SVM,'bs-',251:307,target(251:307),'ro:')
    xlabel('Forecasting step')
    ylabel('Pressure')
    legend('forecasting result','original data')