function V = constructV(varargin);
% [ V ] = constructV( Om, [k] )
% Constructs a basis for the space V_k(Om), as described in Charina, Mejstrik, 2018.
%
% Input: 
%   Om          The set for which V_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 V_k, where dim V_k>0 are computed
%
% Options:
%   '01'        give input as a 0/1 matrix
%
% Output:
%   V           matrix/or cell array of matrices, each column is one basis vector for the space V_i i \in k
%
% E.g.: vdisp(constructVbar([0 1 2 3]))
%       Om=constructOmega('2_butterfly'); constructV(Om,1)
%       constructV([1 0 1; 0 1 1; 1 0 1],1,'01')
%
% See also: constructVbar, constructU
%
% Reference:
%   Set V_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

% XX Verhalten angleichen mit constructU. Dort ist 'k' ein obligator parameter
% XX test if transition matrices or subdivision scheme is given instead of set Omega


% Parse Input
Om=varargin{1};
[zeroone,varargin]=parsem('01',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);
    varargin(1)=[];
    kflag=0; %flag if 'k' is given or not
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}=constructV_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 [ V] = constructV_worker( varargin)
%for Omega a matrix of vectors
Om=varargin{1};
k=varargin{2};

if(k>8); cprintf('err','dimVk: Function may return false value, since k>8.\n'); end;

sizeOm=size(Om,2);
dim=size(Om,1);
    
    
    mu=mixvector(0:k,dim);  %all possible multiindices which sum is less then k
    [~,osum]=sort(sum(mu,1));
    mu=mu(:,osum);
    mu=mu(:,sum(mu,1)<=k & sum(mu,1)>=0);
    mu=fliplr(mu);
    MV=zeros(size(mu,2),sizeOm);
    for j=1:size(mu,2);
        for i=1:sizeOm;
            MV(j,i)=prod(Om(:,i).^mu(:,j));
        end
    end

V=null(MV,'r'); 
%V=V';

end

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

