root/library/utia_legacy/ticket_12/ld2ld.m @ 798

function [LD, D] = ld2ld(LD, strs, strt)
% convert L'DL decomposition described by a source structure strs
%             to another corresponding to a target structure strt
% LD    = ld2ld(LD, strs, strt)
% LD    = ld2ld(L,D)      composes the factors L and D into the matrix LD
% [L,D] = ld2ld(LD)     decomposes the matrix LD into the factors L and D
%
% LD  : L'DL decomposition of a positive definite matrix 
%       LD ... D placed instead of the unit diagonal of L 
% strs: source structure
% strt: target structure
%
% Designed : M. Karny, P. Nedoma
% Updated  : April 2000, MK July 2004
% Project  : post-ProDaCTool

nPsi = length(LD);

% call [L,D] = ld2ld(LD)
if nargin==1, D = diag(diag(LD));
   LD = LD - D + eye(nPsi);
   return;
end

% call LD    = ld2ld(L,D)
if nargin==2, LD = LD + strs - eye(nPsi);
   return;
end                                        

% call LD    = ld2ld(LD, strs, strt)
if isempty(strt), LD=LD(1,1); return; end
dim  = size(strt,2);
lens = size(strs,2);
if lens<length(LD)-1 | lens>length(LD)
   error('incompatible source structure and LD');
end
% add a fictive output channel (1000)
if lens<length(LD)
   strs = [  [1000;0], strs ];
   strt = [  [1000;0], strt ];
   lens = lens + 1;
   dim  = dim  + 1;
end
for i=dim:-1:1
    chn = strt(1,i);
    td  = strt(2,i);
    j   = find(strs(1,:)==chn & strs(2,:)==td);
%chn,td,strs, strt, j, dim, 
    if isempty (j), error('target structure item is not in source'); end
    LD = ldperm(LD,j);                                 % column j to 1th position
    strs = [strs(:,j),strs(:,1:j-1),strs(:,j+1:lens)]; % rebuild strs
%strs, pause
end
% cut and transform L'DL if target is shorter
if lens>dim
   for i=2:dim
       LD = ldperm(LD,dim);
   end
   for i=dim+1:lens
       LD = ldperm(LD,lens);
   end
   LD   = LD(lens-dim+1:lens, lens-dim+1:lens);
   LD   = ldperm(LD, dim);
end

if nargout==2
   [LD, D] = ld2ld(LD);
end