function [HH, H1, HH2, HH3] = TrainVolterra(XX,Y,Memory,order,Sym,LambdaDecay1,LambdaDecay2,LambdaReg1,LambdaReg0)
% Training the Volterra-model 
%
%INPUTS
% XX = Structure from function InpVolterraSeries
% Y = Training data,(Nx1)
%Memory = memory of the system
%order = order of the system
%SYM = Boolean(true/false), choose if we want to reduce the identification problem due to symmetry
%LambdaDecay1 = 1st order Lambda decay regularization factor
%LambdaDecay2 = 2nd order Lambda decay regularization factor
%LambdaReg1 = Lambda, 1st order Tikhonov regularization
%LambdaReg0 = Lambda, 0st order Tikhonov regularization
 
 
%Outputs:
%HH = Full structure of output kernel matrix, se Wu 2014,(1xLength)
%H1 = First kernel(1x(Memory + 1))
%HH2 = Matrix of second kernel((Memory++1)x(Memory+1))
%HH3 = 3dim Matrix of third kernel((Memory+1)x(Memory+1)x(Memory+1))

% turning the array if not column
if isrow(Y)
    Y = Y';
end

%% Train

[LL, LL0,LLDecay1,LLDecay2] = L_reg_order2(Memory,order);

if Sym
    S_red = SelectionMatrixSymmetricReduction(Memory+1,order);
    XX = XX*S_red;
    LL = S_red'*LL*S_red;
    LLDecay1 = S_red'*LLDecay1*S_red;
    LLDecay2 = S_red'*LLDecay2*S_red;
    LL0 = S_red'*LL0*S_red;
end

X_full_precalc = XX.'*XX;

if order ==2
     H_vec_full_LS=(X_full_precalc +  LambdaReg1*LL + LambdaDecay1*LLDecay1 + LambdaDecay2*LLDecay2+LambdaReg0*LL0)\XX.'*Y; 
elseif order ==1
    H_vec_full_LS=(X_full_precalc + LambdaReg1*LL + LambdaDecay1*LLDecay1+LambdaReg0*LL0)\XX.'*Y;
end


%% Splitting HH into different elements
HH = H_vec_full_LS;
if Sym
    HH = S_red*HH;
end

H1 = HH(1:Memory+1);
HH2 = NaN;
HH3 = NaN;

if order >1
    H2 = HH(Memory+2:Memory+1+(Memory+1)^2);
    HH2 = reshape(H2,length(H1),length(H1));
    H2up = triu(HH2,1);
    H2Low = tril(HH2,-1);
    H2Diag = diag(diag(HH2));
    HH2 = (H2up + H2Low + 2.*H2Diag+H2up' + H2Low')./2;
    
end

if order >2
    H3 = HH(Memory+2+(Memory+1)^2:end);
    HH3 = reshape(H3,length(H1),length(H1),length(H1));
end
end

