close all
%---------- Generate training and test set ----------
clc

%%%% -------------------*******------------------
load sunspot2003
[x,xmin,xmax]=scale(sunspot2003,0,1);  

plot(x);
  xlabel('Year');
  ylabel('Sunspot activity');
  title('Sunspot benchmark data');
%=====================------------------------------
nlag=11;
[x0,train,test]=tsgenf_hzh(x,304,244,1,nlag);

%------------******-----------------
%  W1    : Input-to-hidden layer weights. The matrix dimension is
%          dim(W1) = [(# of hidden units) * (inputs + 1)] (the 1 is due to the bias)
%  W2    : hidden-to-output layer weights.
%           dim(W2) = [(outputs)  *  (# of hidden units + 1)]
rand('seed',0);
  W1         = rand(6,nlag+1)-0.5;  % Weights to hidden layer 權值加減-0。5 +0。5有助于預測的準確性
  W2         = rand(1,7)-0.5;  % Weights to output
  NetDef = ['HHHHHH'       
            'L-----'];     
        
    maxiter = 1000;
    stop_crit = 1e-12;
    lambda=1;
    D=1;%加動量項有助于快速收斂
    trparms=[maxiter stop_crit lambda D];
           
    [W1,W2,PI_vector,iter,lambda]=marq(NetDef,W1,W2,train(:,1:nlag)',train(:,nlag+1)',trparms);
        
    % -----------  Validate Network  -----------
    [Y_sim,E0,PI]=nneval(NetDef,W1,W2,train(:,1:nlag)',train(:,nlag+1)');
    
    %[Y_sim,E,PI,tr,E1,PI1] = nneval_hzh(NetDef,W1,W2,test(:,1:nlag)',test(:,nlag+1)');%增加了遞歸方式預測輸出
    [Y_sim,E] = nneval(NetDef,W1,W2,test(:,1:nlag)',test(:,nlag+1)');
    % -----------  Plot Cost function  -----------
    figure
    semilogy(PI_vector)
    title('Criterion evaluated after each iteration')
    xlabel('Iteration (epoch)')
    ylabel('Criterion')
    grid
    
    sum(E0.*E0)/length(train)
    sum(E.*E)/length(test)
    %sum(E1.*E1)/length(test)