

% SSpline   Computes the estimate of a general linear model using the Stable Spline kernel
% [M,iper,Ak]=SSpline(y,U,p,lab,mv,mb,cn,red,LP,LP2,iperstart)
%
% Outputs:
% M: returns the estimated model in IDPOLY format
% iper: returns the estimated hyperparameter vector
% Ak: complexity of the nonparametric model evaluated by Akaike criterion
% (computed only if both LP and LP2 are set to 0)
%
% Inputs:
% y: matrix containing the output data. The number of columns equals the number of system outputs
% U: matrix containing the input data. The number of columns equals the number of system inputs
% p: number of predictor coefficients to estimate
% lab: "yy" = estimate also a noise model, "ny" = output error model, "yn" = time series
% mv:  1 - estimate only one variance $\lambda^2$ common to all the predictor impulse responses
%      0 - estimate one variance $\lambda_k^2$ for each predictor impulse response
% mb:  1 - estimate only one parameter $\beta$ (decay rate of the variance)
%      0 - estimate one $\beta_k$ for each predictor impulse response
% cn:  dimension of the parametric component of the prior that corresponds to the 
% number of poles with negative real part introduced in the model
% red: number of input-output data used to estimate the hyperparameter vector (can not exceed column dimension of y minus p)
% LP:  1 - introduce an exponential hyperprior on the scale factors of the kernels to obtain a sparse solution
%      0 - disabled
% LP2: 1 - introduce an exponential hyperprior on the scale factors of the kernels and compute a suboptimal sparse solution
%      0 - disabled
% iperstart: starting point for the hyperparameter vector (optional)

function [M,iper,Ak,nv]=SSpline(y,U,p,lab,mv,mb,cn,red,LP,LP2,iperstart);
warning off
n=size(y,1)-p;%dati a partire da p2+1
ipn=-0.5*ones(1,cn);%starting point for the bias space (poles with negative real part)
R=0;%here and below: particular options which are disabled
ipp=[];
BL=[];
Ak=[];%AIC


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% m, handling initial conditions, yred and nv%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[A,m,nv,yred,ytrue]=Abuild(y,U,p,lab);%m = size(A,3)
nv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P1, P2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if mv==1;
    P1=1;%common variance to all i.r.
else
    P1=m;%m variances, to be multiplied accounting for the outputs number 
end
if mb==1;
    P2=1;%one beta
else
    P2=m;%m beta 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% P3 e P4 bias space coefficients in BiasN,BiasP and P5,P6 %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(ipn);
   BiasN=[];%no pole with negative real part to estimate
else
   BiasN=poly(ipn);
   BiasN=BiasN(2:end);
end
if isempty(ipp);
   BiasP=[];%no pole with positive real part to estimate
else
   BiasP=poly(ipp);
   BiasP=BiasP(2:end);
end
P3=size(BiasN,2);
P4=size(BiasP,2);
if LP==0;
    P5=[];
else
    P5=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if size(y,2)==1;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MISO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%% HYPERPARAMETERS ESTIMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %[Lamb  beta  BiasN  BiasP  LP]
    %[P1     P2     P3     P4   P5]
    ipstart=[5*ones(1,P1)   0.5*ones(1,P2)  BiasN  BiasP   15*P5];
    if nargin==11;
        ipstart=iperstart;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    %options = optimset('MaxIter',300*(P1+P2+P3+P4+P5+3),'MaxFunEvals',300*(P1+P2+P3+P4+P5+3));
    [iper,objval]=fminsearch('minloglik',ipstart,foptions,A(1:red,:,:),yred(1:red),P1,P2,P3,P4,P5,R,LP,BL,cn,nv);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%% FURTHER SPARSITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (LP==1 & mv~=1);
        for i=1:m;
            iperr=iper;
            iperr(i)=-inf;
            [obj]=minloglik(iperr,A(1:red,:,:),yred(1:red),P1,P2,P3,P4,P5,R,LP,BL,cn,nv);
            if obj<=objval;
                iper=iperr;
                objval=obj;
            end
        end
    end
    if (LP==0 & LP2==1 & mv==1);
        P1=m;
        P5=1;
        iper=[iper(1)*ones(1,m) iper(2:end) iper(1)];%le varianze da 1 a m e la varianza del prior esponenziale sulle radici di exp(iper) alla fine
        [objval]=minloglik(iper,A(1:red,:,:),yred(1:red),P1,P2,P3,P4,P5,R,1,BL,cn,nv);
        for i=1:m;
            iperr=iper;
            iperr(i)=-inf;
            [obj]=minloglik(iperr,A(1:red,:,:),yred(1:red),P1,P2,P3,P4,P5,R,1,BL,cn,nv);
            if obj<=objval;
                iper=iperr;
                objval=obj;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%% AIC for selecting the complexity of the bias space %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    npar=max(size(iper));ndat=n;
    if LP==0 & LP2==0
        Ak=2*objval+2*npar;
    else
        Ak=[];
    end
    %+[2*npar*(npar+1)]/[ndat-npar-1]; 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION ESTIMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    BL=[];%minloglik restituisce le C giuste solo se BL=[]
    [obj,C,Ey]=minloglik(iper,A,y,P1,P2,P3,P4,P5,R,LP,BL,cn,nv);%nuovo il cn
    Ey_inv=inv(Ey);
    for i=1:m;
        eta(:,:,i)=C(:,:,i)*A(:,:,i)'*Ey_inv*yred;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%% THE MODEL TO RETURN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [eta,M,D,S]=creaModel(eta,nv,lab,m);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MIMO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%% HYPERPARAMETERS ESTIMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nout=size(yred,2);
    if mv==1;
        ip1=5;
    else
        ip1=5*ones(1,P1*nout);
    end
    if mb==1;
        ip2=0.5;
    else
        ip2=0.5*ones(1,P2*nout);
    end
    %[      Lamb/ip1                      beta/ip2                 BiasN    BiasP     LP]
    %[P1*nout | 1(=P1) se mv=1      P2*nout | 1(=P2) se mb=1        P3       P4       P5]
    ipstart=[ip1 ip2  BiasN  BiasP   15*P6];
    Res=size(ipstart,2)-size(ip1,2)-size(ip2,2);
    options = optimset('MaxIter',300*(P1+P2+P3+P4+P5+3),'MaxFunEvals',300*(P1+P2+P3+P4+P5+3));
    [iper,objval]=fminsearch('minloglikMIMO',ipstart,foptions,A(1:red,:,:),yred(1:red,:),P1,P2,P3,P4,P5,R,LP,BL,cn,nv,mv,mb,Res);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%FURTHER SPARSITY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if LP==1 & mv~=1;
        for i=1:size(ip1,2);
            iperr=iper;
            iperr(i)=-inf;
            [obj]=minloglikMIMO(iperr,A(1:red,:,:),yred(1:red,:),P1,P2,P3,P4,P5,R,LP,BL,cn,nv,mv,mb,Res);
            if obj<=objval;
                iper=iperr;
                objval=obj;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%% AIC for selecting the complexity of the bias space %%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    npar=size(ipstart,2);ndat=n;
    if LP==0 & LP2==0
        Ak=2*objval+2*npar;
    else
        Ak=[];
    end
    %+[2*npar*(npar+1)]/[ndat-npar-1]; 2 per AIC, log(n) per BIC
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION ESTIMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    BL=[];%minloglik restituisce le C giuste solo se BL=[]
    [obj,C,Ey]=minloglikMIMO(iper,A,y,P1,P2,P3,P4,P5,R,LP,BL,cn,nv,mv,mb,Res);
    for i=1:size(C,4);%=nout
        Ey_inv=inv(Ey(:,:,i));
        for j=1:nout;
            for i=1:m;
                eta(:,:,i,j)=C(:,:,i,j)*A(:,:,i)'*Ey_inv*yred(:,j);
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%% THE MODEL TO RETURN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    M=[];D=[];S=[];%in the current version no MIMO model is returned but just the predictor impulse responses
%     [eta,M,D,S]=creaModel(eta,nv,lab,m);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end