close all
%---------- Generate training and test set ----------
%clear all
clc
load sunspot2003;
org_data=scale(sunspot2003,0,1); % load the normalised originate data sunspots
%%%% -------------------*******------------------
load emd_sunspot
x=imf(7,:);  
[x,xmin,xmax]=scale(x,0,1);
plot(x);
  xlabel('Year');
  ylabel('Sunspot activity');
  title('Sunspot benchmark data');
%=====================------------------------------
nlag=3;
[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(8,nlag+1)-0.5;  % Weights to hidden layer 權(quán)值加減-0。5 +0。5有助于預(yù)測的準(zhǔn)確性
  W2         = rand(1,9)-0.5;  % Weights to output
  NetDef = ['HHHHHHHH'       
            'L-------'];     
        
    maxiter = 1000;
    stop_crit = 1e-12;
    lambda=1;
    D=1;%加動(dòng)量項(xiàng)有助于快速收斂
    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,E,PI] = nneval(NetDef,W1,W2,train(:,1:nlag)',train(:,nlag+1)');
    [Y_sim1,E,PI,tr,E1,PI1] = nneval_hzh(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(E.*E)/length(test)
    sum(E1.*E1)/length(test)