function [ D, psin, beta ] = BioDictDel( D, psin, beta, j ) 
% BioDictDel deletes a vector from a basis selected from a given dictionary and
% appropriately modifies biorthogonal functions, orthonormal functions spanning the same 
% space as the basis, and unselected dictionary atoms subtracted by their components in
% the selected basis.  
%  
% The dictionary is then rearranged to have one to one correspondence with  beta (dual 
% functions) and psin.
%  
% Usage:  [ D, psin, beta ] = BioDictDel( D, psin, beta, j );
%
% Inputs:
%   D       dictionary, first kb functions are the selected basis
%   psin    psin(:,1:kb) orthonormal functions spanning the same space as D(:,1:kb)
%           psin(:,kb+1:end)=D(:,kb+1:end) without component from D(:,1:kb)   
%   beta    kb biorthogonal functions to D(:,1:kb), kb=size(beta,2)
%   j       index of function to eliminate  
%  
% Outputs:
%   D       updated (rearranged) dictionary  
%   psin    psin(:,1:kb-1) updated (recalculated) orthonormal functions spanning the same
%           space as the chosen atoms psin(:,kb:end) updated (recalculated) unchosen
%           dictionary atoms subtracted by their components in the selected basis.
%   beta    updated biorthogonal functions to D(:,1:kb-1) 
%  
% References:
%    L. Rebollo-Neira, "Recursive bi-orthogonalisation approach and orthogonal projectors",
%    math-ph/0209026 (2002).  
%  
% See also: NBioDictDel, BioDelete, NBioDelete.

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/
  
name='BioDictDel';
k=size(D,2);
kb=size(beta,2); %number of selected functions (dimension)
%simple tests of input data
if (j<1) | (j>kb)
  fprintf('%s: Required index %d is out of range (%d atoms)\n',name,j,kb);
  return;
end  
b=beta(:,j)/norm(beta(:,j));
beta(:,j)=[];
beta=beta-b*(b'*beta);
%downdating of psin
R=psin(:,1:kb)'*D(:,1:kb);
dummy=R(:,j);R(:,j)=[];R(:,kb)=dummy;
for i=j:kb-1
  p=[i i+1];
  [G,R(p,i)]=planerot(R(p,i)); %find appropriate Givens rotation
  R(p,i+1:kb)=G*R(p,i+1:kb);   %apply Givens rotation to R(p,i+1:k)
  psin(:,p)=psin(:,p)*G';      %apply Givens rotation to psin(:,p)
end
ud=D(:,j);D(:,j)=[];D(:,k)=ud;
ud=psin(:,kb);psin(:,kb)=[];%psin(:,k)=ud;
psin(:,k)=ud*(D(:,k)'*ud);

for i=kb:k-1
  psin(:,i)=psin(:,i)+(ud'*D(:,i))*ud;
end  

%Last modification Laura REBOLLO-NEIRA (2009)
%
%Copyright (C) 2006 Miroslav ANDRLE and Laura REBOLLO-NEIRA
%
%This program is free software; you can redistribute it and/or modify it under the terms 
%of the GNU General Public License as published by the Free Software Foundation; either 
%version 2 of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%See the GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License along with this program;
%if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
%Boston, MA  02110-1301, USA.
