function [ V, Om ] = constructVbar(varargin);
% [ V, Om ] = constructVbar( Om, [k], [options] )
% Constructs a basis for the space Vbar_k(Om) as described in Charina, Mejstrik, 2018.
%
% Input: 
%   Om          The set for which Vbar_k shall be constructed as an array of column vectors
%   k           index k, optional. Either an integer or a vector of integers (greater equal zero)
%               If k is empty, then all spaces Vbar_k, where dim Vbar_k>0 are computed
%
% Options:
%   '01'        give Omega as a 0/1 matrix
%
% Output:
%   V           matrix/or cell array of matrices, each column is one basis vector for the space Vbar_i, i \in k
%   Om          The set for which V_k is constructed as an array of column vectors.
%
% E.g.: vdisp(constructVbar([0 1 2 3]))
%       Om=constructOmega('2_butterfly'); constructVbar(Om,1)
%       constructVbar([1 0 1; 0 1 1; 1 0 1],1,'01')
%
% See also: constructV, constructU
%
% Reference:
%   Set Vbar_k is described in:
%       Maria Charina, Thomas Mejstrik,
%       Multiple multivariate subdivision schemes: Matrix and operator approaches, 
%       Journal of Computational and Applied Mathematics, 2018
%
% Written by: tommsch, 2016

Om=varargin{1};
if(isempty(Om)); V=[]; return; end;
[zeroone,varargin]=parsem('01',varargin);
[verbose,varargin]=parsem('verbose',varargin);

if(zeroone); 
    dim=length(size(Om));
    Om=supp(Om,dim,zeros(dim,1));
end;

if(size(varargin,2)>=2 && ~isempty(varargin{2}) && isnumeric(varargin{2}) && ~ischar(varargin{2})); %test if second argument is a valid number for k
    k=varargin{2};
    kflag=1;%flag if 'k' is given or not
    varargin(1:2)=[];
else;
    k=0:size(Om,2);
    kflag=0;%flag if 'k' is given or not
    varargin(1)=[];
end
if(~isempty(varargin)); vprintf('Unkown options: \n %v',varargin,'cpr','err','imp',[0 verbose]); end;

j=0; %idx
V=cell(1,size(k,2)); %allocate space big enough

while(true);
    j=j+1;
    V{j}=constructVbar_worker(Om,k(j));
    if( (kflag && j==size(k,2)) || (~kflag && isempty(V{j})) );        break;     end; %test stopping criterion
end
if(~kflag); V=V(~cellfun('isempty',V)); end;
if(size(k,2)==1); V=V{1}; end;


end

function Vbar = constructVbar_worker( Om, k )
dim=size(Om,1);
Om=characteristic(Om);
Vbar=zeros(0,numel(Om)); % is getting transposed at the end

mu=mixvector(0:k+1,dim);
mu=mu(:,sum(mu,1)==k+1);

for i=1:size(mu,2);

    coriginal=diffsequence(1,mu(:,i)');
    sizecoriginal=size(coriginal);
    sizeOmega=size(Om);
    if(length(sizecoriginal)<length(sizeOmega));
        sizecoriginal(length(sizeOmega))=1;
    end;
    sizecoriginal(sizecoriginal==0)=1;
    Loriginal=sizeOmega-sizecoriginal;
    if(any(Loriginal<0)); continue; end;
    L=num2cell(Loriginal);
    for j=1:size(L,2);
        L{j}=0:L{j};
    end;
    Lstart=mixvector(L);
    Lend=setplus(Loriginal',-Lstart,'stable');

    for j=1:size(Lstart,2);
        c=coriginal;
        c=padarraym(c,Lstart(:,j)','pre');
        c=padarraym(c,Lend(:,j)','post');
        size(c);

        if (isequal((Om~=0).*(c~=0),c~=0));
            Vbar=[Vbar; reshape(c,numel(c),1)']; %#ok<AGROW>
        end;
    end

end

if(isempty(Vbar)); return; end;

Vbar(:,find(~Om))=[]; %#ok<FNDSB>

RANK=rank(Vbar);
for i=1:size(Vbar,1);
    if(rank(Vbar([1:i-1 i+1:end],:))==RANK); 
        Vbar(i,:)=0;  
    end;
end   
Vbar=removezero(Vbar,1);
Vbar=Vbar';
end

function dummy; end %#ok<DEFNU>  %Generates an error, if the 'end' of a function is missing.