% Numerical computation of the Hilbert function of a zero-scheme
% Z. Griffin, J. Hauenstein, A. Sommese, and C. Peterson
% Copyright July 2011

function [reg, hilbertFunc, stdMon] = HilbertFuncVeronese(numVars, pointLoc, tol, compMon) 
    % compute the Hilbert function of the points
    
    % INPUT
    % numVars = number of variables
    % pointLoc = string of the name of the file containing the points
    % tol = tolerance for rank determination
    % compMon = boolean which determines whether to compute the standard monomials
    
    % OUTPUT
    % reg = index of regularity
    % hilbertFunc = Hilbert function up to index of regularity
    % stdMon = standard monomials
    
    % setup the points
    [numPts, Pts, rV] = readInMultPoints(numVars, pointLoc);

    % check for errors in the points
    if rV == -1
        reg = 0;
        hilbertFunc = 0;
        stdMon = 0;
        return;
    end;

    % initialize the homogenizing coordinate to 1
    Pts = [ones(numPts, 1) Pts];

    % preallocate h
    h = zeros(10,1);

    % initialize
    k = 0;

    while (1)
        % increment k
        k = k + 1;

        % compute the Veronese embedding of the Pts
        Ak = VeroneseMatrix(Pts, k-1, 1);

        % compute the rank
        h(k) = rank(Ak,tol);

        % see if we have the correct rank
        if h(k) == numPts
            break;
        end;
    end;

    % copy h to hilbertFunc
    hilbertFunc = h(1:k);

    % setup reg
    reg = k - 1;

    if compMon
        % setup the monomials
        M = setupMonomials(numVars, reg);

        % compute the pivot columns
        [~, jb] = rref(Ak, tol);

        % setup the monomial support
        stdMon = M(jb,:);
    else
        stdMon = 0;
    end;

return;

function [numPts, Pts, rV] = readInMultPoints(numVars, fileLocation)
    % read in the points from a file
    
    % open file
    RF = fopen(fileLocation, 'r');

    % make sure the file exists
    if RF == -1
        % file does not exist
        disp(['The file ''', fileLocation, ''' does not exist!']);
        numPts = 0;
        Pts = 0;
        rV = -1;
        return;
    end;

    % read in the number of points
    numPts = fscanf(RF, '%d', 1);

    % setup Pts
    Pts = zeros(numPts, numVars);

    % loop through to read in the points
    for k = 1:numPts
        % read in the point
        for j = 1:numVars
            tempR = fscanf(RF, '%e', 1);
            tempI = fscanf(RF, '%e', 1);
            Pts(k,j) = tempR + 1i * tempI;
        end;
    end;

    % close file
    fclose(RF);
    
    % setup return value
    rV = 0;
return;

function [VM] = VeroneseMatrix(Pts, deg, normalize)
    % compute the Veronese embedding of 'Pts' in degree 'deg'

    % default
    if nargin < 3
        normalize = 1;
    end;

    N = size(Pts,1); % number of points
    numVars = size(Pts,2); % number of variables

    % initialize VM
    VM = zeros(N,nchoosek(numVars + deg - 1, deg));

    % setup A
    A = setupDegree(deg, numVars);
    
    % loop to setup VM
    for j = 1:N
        % normalize Pts
        if normalize == 0
            P = Pts(j,:);
        else
            P = Pts(j,:) / norm(Pts(j,:));
        end;
        
        % compute embedding
        if normalize == 0
            VM(j,:) = veronese(P,A);
        else
            v = veronese(P,A);
            VM(j,:) = v/norm(v);
        end;
    end;
return;

function [M] = veronese(X, A)
    % compute the Veronese embedding based on A

    N = size(A,1);
    M = zeros(1,N);
    for j = 1:N
        M(1,j) = prod(X.^A(j,:));
    end;
return;

function [A] = setupDegree(deg, n)
    % setup degrees for Veronese embedding

	numRows = nchoosek(n+deg-1,deg);
    A = zeros(numRows,n);
    A(1,1) = deg;

    for j = 1:numRows-1
        % copy
        A(j+1,:) = A(j,:);
        % find target entry
        target = -1;
        for k = n-1:-1:1
            if A(j+1,k) ~= 0
                target = k;
                break;
            end;
        end;
        % move entry
        A(j+1,target) = A(j+1,target) - 1;
        A(j+1,target+1) = A(j+1,target+1) + 1;
        % move all of final entry to target+1 if not the same!
        if target ~= n-1
            A(j+1,target+1) = A(j+1,target+1) + A(j+1,n);
            A(j+1,n) = 0;
        end;
    end;
return;

function [monomials] = setupMonomials(numVars, order)
    % setup the monomials up to the given order

    q = nchoosek(numVars + order, order);
    monomials = zeros(q,numVars);

    % setup the monomials for each order
    loc = 1;
    for j = 0:order
        % compute the number of monomials of order j
        M = setupMonomialsOrder(numVars, j);
        rows = size(M,1);
        monomials(loc:loc+rows-1,1:numVars) = M;
        % update loc
        loc = loc + rows;
    end;
return;
    
function [M] = setupMonomialsOrder(numVars, order)
    % setup the monomials of the given order

    if (order < 0)
        M = 0;
    elseif (order == 0)
        M = zeros(1, numVars);
    else % order >= 1
        q = nchoosek(numVars - 1 + order, order);
        M = zeros(q, numVars);

        % setup the first one
        M(1,1) = order;
        % loop through the other ones
        for j = 2:q
            % copy over M(j-1,:)
            M(j,:) = M(j-1,:);
        
            % find where we can increment
            target = numVars - 1;
            while (target >= 0)
                if not(M(j,target) == 0)
                    break;
                else
                    target = target - 1;
                end;
            end;
        
            % increment
            M(j,target) = M(j,target) - 1;
            M(j,target+1) = M(j,target+1) + 1;
        
            if not(target + 1 == numVars)
                M(j,target+1) = M(j,target+1) + M(j,numVars);
                M(j,numVars) = 0;
            end;
        end;
    end;
return;