%
% converting a proper cone
%     A : face vectors (column vectors)
%     U : span vectors (column vectors)
%
function U = PCC_convert_proper_cone ( A )
    A = PCC_simplify(A);
    n = size(A,1); % dimension
    m = size(A,2); % number of face vectors
    % note n <= m

    % calculating n span vectors from selected n face vectors
    U = [];
    Uset = [];
    for k=1:n
        s = [1:k-1,k+1:n];
        N = null(A(:,s)');
        N = N(:,1);
        if A(:,k)'*N < -eps
            U = [U, N];
            Uset = [Uset, s'];
        elseif A(:,k)'*N > eps
            U = [U, -N];
            Uset = [Uset, s'];
        end  
    end
    
    % calculating m - n span vectors one-by-one
    for k=n+1:m
        ns = size(U,2); % number of span vectors
        s = A(:,k)'*U;
        index = find(s > eps); % span vector that inner product with ak > 0
        if length(index) > 1
            disp('calculation fails');
            index
            pause;
        end
        jp = index(1);
        if (1 < jp && jp < ns) % jp = 2, 3, ..., ns-1
            sjmp = union(intersect(Uset(:,jp-1),Uset(:,jp)),k);
            Njmp = null(A(:,sjmp)');
            Njmp = Njmp(:,1);        
            sjpm = union(intersect(Uset(:,jp),Uset(:,jp+1)),k);
            Njpm = null(A(:,sjpm)');
            Njpm = Njpm(:,1);
            index = [1:jp-1,jp+1:ns];
            Njmp = negative_vector(Njmp, A, index);
            Njpm = negative_vector(Njpm, A, index);
            U = [ U(:,[1:jp-1]), Njmp, Njpm, U(:,[jp+1:ns]) ];
            Uset = [ Uset(:,[1:jp-1]), sjmp', sjpm', Uset(:,[jp+1:ns]) ];
        elseif (jp == 1) % jp = 1
            sjmp = union(intersect(Uset(:,ns),Uset(:,jp)),k);
            Njmp = null(A(:,sjmp)');
            Njmp = Njmp(:,1);        
            sjpm = union(intersect(Uset(:,jp),Uset(:,jp+1)),k);
            Njpm = null(A(:,sjpm)');
            Njpm = Njpm(:,1);
            index = [jp+1:ns];
            Njmp = negative_vector(Njmp, A, index);
            Njpm = negative_vector(Njpm, A, index);
            U = [ Njmp, Njpm, U(:,[jp+1:ns]) ];
            Uset = [ sjmp', sjpm', Uset(:,[jp+1:ns]) ];
        elseif (jp == ns) % jp = ns
            sjmp = union(intersect(Uset(:,jp-1),Uset(:,jp)),k);
            Njmp = null(A(:,sjmp)');
            Njmp = Njmp(:,1);        
            sjpm = union(intersect(Uset(:,jp),Uset(:,1)),k);
            Njpm = null(A(:,sjpm)');
            Njpm = Njpm(:,1);
            index = [1:jp-1];
            Njmp = negative_vector(Njmp, A, index);
            Njpm = negative_vector(Njpm, A, index);
            U = [ U(:,[1:jp-1]), Njmp, Njpm ];
            Uset = [ Uset(:,[1:jp-1]), sjmp', sjpm' ];
        else
            disp('calculation fails');
            fprintf('[1, %d]   + %d\n', ns, jp);
            pause;
        end
    end
end

function vout = negative_vector (v, A, index)
    s = A(:,index)'*v;
    if (prod(s <= eps)) % all negative
        vout = v;
    elseif (prod(s >= -eps)) % all positive
        vout = -v;
    else
        disp('calculation fails');
        v
        index
        s
        s <= 0
        s >= 0
        pause;
    end
end
