function [Fit_err_norm,LH_norm,LambdaInterval,Fit_error_normFull,LambdaNormFullDecay1,LambdaNormFullDecay2,LambdaDecayInt21,LambdaDecayInt22] = FindLambda(XX,Y,Memory,order,LambdaIntervalReg1,LambdaDecayInterval1,LambdaDecayInterval2,Sym,RegType)
% Finding the Lambda values for L-curves and double L curves
%
%INPUTS
% XX = Structure from function InpVolterraSeries()
% Y = Training data,(Nx1)
%Memory = memory of the system
%order = order of the system
%LambdaIntReg = The Interval Lambda is looped over, for Tikhonov regularization(1st and 0th)
%LambdeDecayInt1 = The Interval Lambda decay is looped over, axis 1
%LambdeDecayInt2 = The Interval Lambda decay is looped over, axis 2
%Sym = true/false. If one uses symmetry reduction(true false). 
%RegType = 'Reg1', 'Reg0' or 'decay' for 1st order Tikhonov, 0th order tikhonov or decay type regularization
%
%Outputs:
% 
%For Single L-curves: In other words Tikhonov regularization
%Fit_err_norm            %||r*H-F1|||     (1xlength(LambdaInterval))
%LH_norm                 %||L*H||            (1xlength(LambdaInterval))
%LambdaInterval        %%LambdaInterval vector that fits to the different points, Lambda^2   (1xlength(LambdaInterval)) = LambdaIntervalReg1

%For double L-curves:
%FittError_normFull         %||r*H-F1|||   ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))   
%LambdaNormFullDecay1       %Lambda Decay Interval in direction 1
%LambdaNormFullDecay2       %Lambda Decay Interval in direction 2
%LambdaDecayInt21            %||L1*H||, ((length(LambdaDecayInt21)xlength(LambdaDecayInt22)) 
%LambdaDecayInt22            %||L2*H||, ((length(LambdaDecayInt21)xlength(LambdaDecayInt22))  

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

%% Train

%Extract different L-matrixes dependent on type of regularization
[LL, L0,LLDecay1,LLDecay2] = L_reg_order2(Memory,order);

%Exeption of RegType = Reg0, then change the LL matrix. 
if strcmp(RegType,'Reg0')
    LL = L0;
end

%Reduce the matrixes based of symmetry
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;
end


%Calculate some of the mean square error matrixes before the loop
X_full_precalc = XX.'*XX;

%declare empty variables
Fit_err_norm = NaN;
LH_norm= NaN;
Fit_error_normFull= NaN;
LambdaNormFullDecay1= NaN;
LambdaNormFullDecay2= NaN;
LambdaInterval= NaN;
LambdaDecayInt21=NaN;
LambdaDecayInt22=NaN;


if order == 2
    if strcmp(RegType,'Reg1') || strcmp(RegType,'Reg0')
        LambdaInterval = LambdaIntervalReg1;
        for j = 1:length(LambdaInterval)
            Lambda=LambdaInterval(j);
            H_vec_full_LS=(X_full_precalc+  Lambda*LL)\XX.'*Y;  %The actual lest squares equation
            F_est = XX*H_vec_full_LS;  
            Fit_err_norm(j) = norm(Y - F_est,2);
            LH_norm(j) = norm((LL)*H_vec_full_LS,2);
        end   
        plotLcurve(Fit_err_norm,LH_norm,LambdaInterval);

    
    elseif strcmp(RegType,'decay')
        for j = 1:1:length(LambdaDecayInterval1)
            for i = 1:1:length(LambdaDecayInterval2)    
                LambdaDecayInt21(j,i) = LambdaDecayInterval1(j);
                LambdaDecayInt22(j,i) = LambdaDecayInterval2(i);        
                H_vec_full_LS=(X_full_precalc + LambdaDecayInterval2(i) * LLDecay2 +  LambdaDecayInterval1(j) * LLDecay1)\XX.'*Y;
                F_est = XX*H_vec_full_LS;
                Fit_error_normFull(j,i) = norm(Y - F_est,2);
                LambdaNormFullDecay1(j,i) = norm((LLDecay1)*H_vec_full_LS,2); 
                LambdaNormFullDecay2(j,i) = norm((LLDecay2)*H_vec_full_LS,2); 
            end
        end
    end
    
    
elseif order==1 && ~strcmp(RegType,'decay')
        LambdaInterval = LambdaDecayInterval1;
        for j = 1:length(LambdaInterval)
            Lambda=LambdaInterval(j);
            H_vec_full_LS=(X_full_precalc+Lambda*LL )\XX.'*Y;
            F_est = XX*H_vec_full_LS;
            Fit_err_norm(j) = norm(Y - F_est,2);
            LH_norm(j) = norm((LL)*H_vec_full_LS,2);
        end    
elseif order==1 && strcmp(RegType,'decay')
        LambdaInterval = LambdaDecayInterVall;
        for j = 1:length(LambdaInterval)
            Lambda=LambdaInterval(j);
            H_vec_full_LS=(X_full_precalc + Lambda*LLDecay1)\XX.'*Y;
            F_est = XX*H_vec_full_LS;
            Fit_err_norm(j) = norm(Y - F_est,2);
            LH_norm(j) = norm((LLDecay1)*H_vec_full_LS,2);
        plotLcurve(Fit_err_norm,LH_norm,LambdaInterval);
        end      
end 
toc
end