% 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] = HilbertFunc(numVars, pointLoc, dualLoc, tol, compMon) 
    % compute the Hilbert function of the points and dual bases
    
    % INPUT
    % numVars = number of variables
    % pointLoc = string of the name of the file containing the points
    % dualLoc = string of the name of the file containing the dual bases
    % 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;

    % setup the dual bases
    [mult, ~, ~, colSize, basis, rV] = readInMultStructs(numPts, dualLoc);

    % check for errors in the dual bases
    if rV == -1
        reg = 0;
        hilbertFunc = 0;
        stdMon = 0;
        return;
    end;
    
    % initialize
    k = 1;
    
    % preallocate h
    h = zeros(10,1);

    % setup A_0
    [Ak] = setupAk(numVars, numPts, Pts, mult, colSize, basis, k - 1);
    h(1) = rank(Ak,tol);

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

        % setup A_{k-1} = [A_{k-2} B_{k-1}]
        Ak = [Ak, setupBk(numVars, numPts, Pts, mult, colSize, basis, k - 1)];

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

        % see if we have the correct rank
        if h(k) == size(Ak,1)
            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 [mult, depth, hilbertFunc, colSize, basis, rV] = readInMultStructs(numPts, fileLocation)
    % read in the dual bases 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!']);
        mult = 0; 
        depth = 0; 
        hilbertFunc = 0; 
        colSize = 0; 
        basis = 0;
        rV = -1;
        return;
    end;

    % setup mult, depth & colSize
    mult = zeros(numPts, 1);
    depth = zeros(numPts, 1);
    colSize = zeros(numPts, 1);
    hilbertFunc = zeros(numPts, 0);
    basis = zeros(numPts, 0, 0);

    % loop through to read in the dual bases
    for k = 1:numPts
        % read in the path number & the length
        [num] = fscanf(RF, '%d', 2);
        % store the length
        len = num(2);
        % setup the depth = len - 1
        depth(k) = len - 1;
        % read in the Hilbert function
        [hFunc] = fscanf(RF, '%d', len);
        % setup hilbertFunc
        hilbertFunc(k,1:len) = hFunc;
        % read in the size of the dual basis (mult == rows)
        [num] = fscanf(RF, '%d', 2);
        mult(k) = num(1);
        colSize(k) = num(2);
        % read in the basis
        for j = 1:mult(k)
            for l = 1:colSize(k)
                [b] = fscanf(RF, '%e', 2);
                basis(k,j,l) = b(1) + 1i*b(2);
            end;
        end;
    end;

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

function [Ak] = setupAk(numVars, numPts, Pts, mult, colSize, basis, k)
    % setup the matrix A_k(0)
    
    % compute the sum of the multiplicity & number of monomials
    m = sum(mult);
    q = nchoosek(numVars + k, k);

    % setup the monomials up to order k
    monomials = setupMonomials(numVars, k);

    % initialize Ak
    Ak = zeros(m, q);

    % setup Ak
    p = 1;
    for j = 1:numPts
        for k = 1:mult(j)     
            for l = 1:q
                % setup the lth column based on the lth monomial, the jth point
                % and its kth dual basis vector
                Ak(p, l) = setupAkEntry(monomials(l,:), Pts(j,:), basis(j,k,1:colSize(j)), colSize(j), monomials);
            end;
            % update p
            p = p + 1;
        end;
    end;

return;

function [Ak] = setupAkEntry(currMon, currPt, currBasis, basisSize, basisMon)
    % setup an entry of A_k(0)
    
    % find the number that we need to do
    num = min(basisSize, size(basisMon, 1));

    % initialize
    Ak = 0;

    for j = 1:num
        % d_a (x^b)(\hat(x)) = combination(b,a) * \hat(x)^(b-a)
        [c, d] = combination(currMon, basisMon(j,:));

        if not(c == 0)
            % compute and add on
            Ak = Ak + currBasis(j) * c * prod(currPt .^ d);        
        end;
    end;
return;

function [c, d] = combination(I, K)
    % compute \binom{I}{K}
    
    % compute I - K
    d = I - K;

    % make sure I >= K
    if min(d) >= 0
        % compute nchoosek(I,K)
        c = 1;
        n = max(size(d));
        for j = 1:n
            c = c * nchoosek(I(j),K(j));
        end;
    else
        c = 0;
    end;
return;
    
function [Bk] = setupBk(numVars, numPts, Pts, mult, colSize, basis, k)
    % setup the matrix Bk
    
    % compute the sum of the multiplicity & number of monomials
    m = sum(mult);
    q = nchoosek(numVars - 1 + k, k);

    % setup the monomials up to order k
    allMon = setupMonomials(numVars, k);

    % setup the monomials up to order k
    monomials = setupMonomialsOrder(numVars, k);

    % initialize Bk
    Bk = zeros(m, q);

    % setup Bk
    p = 1;
    for j = 1:numPts
        for k = 1:mult(j)     
            for l = 1:q
                % setup the lth column based on the lth monomial, the jth point
                % and its kth dual basis vector
                Bk(p, l) = setupAkEntry(monomials(l,:), Pts(j,:), basis(j,k,1:colSize(j)), colSize(j), allMon);
            end;
            % update p
            p = p + 1;
        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;