function M=lsdtrpm(DATA,ORDERS,Kr)

%   LSDTRPM computes LS-estimates of Discrete-Time MIMO RPM models 
%   (RPM : reinitialized partial moment)
%   M=lsdtrpm(DATA,ORDERS,Kr)
% 
%   M: returned with the estimated parameters of the RPM model in an IDPOLY
%   object format
%   A(q) y(t) = B(q) u(t-nk)
%   See HELP IDPOLY for all details.
%  
%   DATA: The estimation data as an IDDATA object. 
%   Data must be regularly sampled.
%   See  HELP IDDATA.
%     
%   ORDERS = [na nb nk], the orders and delays of the above model.
%   For multi-output systems, ORDERS has as many rows as there are outputs.
%   na is then an ny|1 column vector whose j:th row gives the order of the
%   polynomials (deniminators of the row transfer functions) relating 
%   the j:th output. Similarly nb and nk are ny|nu matrices. (ny:# of 
%   outputs, nu:# of inputs).
% 
%	Optional input:
%	Kr     : interval of reinitialization
%            Default : Kr = 20
%
%	Author  : Régis Ouvrard, University of Poitiers
%	e-mail  : regis.ouvrard@univ-poitiers.fr
%	Last Revised : 30 March 2010
%   Descritpion  : Rapport Interne N° 20100225OR
%                  http://laii.univ-poitiers.fr/ouvrard/CTRPM/

%   This routine has been tested with MATLAB Version 7.3 and updates 

%% Preliminaries 
[N,ny]=size(DATA.y); 
[N,nu]=size(DATA.u);

[nr,nc]=size(ORDERS);
if nc<(nu+2)    % If the delays are not specified
    ORDERS(1:ny,nu+2:2*nu+1)=zeros(ny,nc-1);
end
na=ORDERS(1:ny,1);
nb=ORDERS(1:ny,2:nu+1);
nk=ORDERS(1:ny,nu+2:2*nu+1);

if ( nargin<3 )
	Kr = [];
end;

if ( isempty(Kr) )
	Kr = 20*ones(1,ny);
end;

%% Loop to estimate ny MISO models 
for counter1=1:ny

    %% Inputs with delays 
    for counter2=1:nu,
        u(:,counter2) = [zeros(nk(counter1,counter2),1);DATA.u(1:N-nk(counter1,counter2),counter2)];
    end

    %% Extract the current output 
    y=DATA.y(:,counter1);

    %% Compute the DT RPM filter sequence mi 
    den_constant=factorial(na(counter1)-1)*factorial(Kr(counter1))/factorial(Kr(counter1)-na(counter1));
    for i=0:Kr(counter1)-na(counter1)
        m(i+1)=factorial(i+na(counter1)-1)/factorial(i)*factorial(Kr(counter1)-i)/factorial(Kr(counter1)-i-na(counter1))/den_constant;
    end

    %% Compute the functions alpha(k) (and gamma(k))
    % alpha_0
    clear alpha
    alpha(:,:)=zeros(N,na(counter1)+1);
    for i=2:N
        if i<=(Kr(counter1)-na(counter1))   % To remove the transient effect of the DT RPM FIR filter
            temp=[flipud(y(1:i-1));zeros(Kr(counter1)-na(counter1)-i+1,1)];
        else
            temp=flipud(y(i-Kr(counter1)+na(counter1):i-1));
        end
        alpha(i,1)=-m(2:Kr(counter1)-na(counter1)+1)*temp;
    end
    %alpha_1
    alpha(2:N,2)=alpha(1:N-1,1)-y(1:N-1);
    %alpha_n
    for i=3:na(counter1)+1
        alpha(2:N,i)=alpha(1:N-1,i-1);
    end

    %% Compute the functions beta(k) 
    clear beta
    beta(:,:)=zeros(N,sum(nb(counter1,:)));
    for counter2=1:nu,
        % beta_0
        for i=1:N
            if i<=(Kr(counter1)-na(counter1))      % To remove the transient effect of the DT RPM FIR filter
                temp=[flipud(u(1:i,counter2));zeros(Kr(counter1)-na(counter1)-i+1,1)];
            else
                temp=flipud(u(i-Kr(counter1)+na(counter1):i,counter2));
            end
            beta0(i,1)=m(1:Kr(counter1)-na(counter1)+1)*temp;
        end
        %beta_n
        beta(:,sum(nb(counter1,1:counter2-1))+1)=beta0;
        for i=2:nb(counter1,counter2)
            beta(2:N,sum(nb(counter1,1:counter2-1))+i)=beta(1:N-1,sum(nb(counter1,1:counter2-1))+i-1);
        end
    end

    %% Generate regressor 
    phi=[alpha(:,2:na(counter1)+1) beta]';

    %% Get LS-estimate 
    theta=inv(phi*phi')*phi*(y-alpha(:,1));

    %% Extract numerator and denominator and build the MISO IDPOLY model
    num=zeros(nu,max(nb(counter1,:)));
    for counter2=1:nu
        num(counter2,1:nb(counter1,counter2))=[theta(na(counter1)+sum(nb(counter1,1:counter2-1))+1:na(counter1)+sum(nb(counter1,1:counter2-1))+nb(counter1,counter2))'];
    end
    M1=idpoly([1 theta(1:na(counter1))'],num,[],[],[],0,DATA.Ts);
    for counter2=1:nu
        M1.InputDelay(counter2)=nk(counter1,counter2);
    end

    %% Build the MIMO model in tf form 
    M(counter1,1:nu)=tf(M1);
    set(M,'OutputName',[M.OutputName(1:counter1-1,1);'y' num2str(counter1)])

end
