% Copyright 2015 by Jonathan D. Hauenstein

% Deflate family of examples modified from Example 5.1 in [Li & Zhi JCS 2012]. 
% For any s ≥ 2, the system has s variables, s polynomials of degrees at most 3. 
% The system has a root at the origin with multiplicity 2^s and breadth 2, i.e. its Jacobian matrix at the origin has co-rank 2. 
% "Certifying isolated singular points and their multiplicity structure" 
% by J.D. Hauenstein, B. Mourrain, and A. Szanto

function [Fnf,Vnf] = deflate51b(s)
    % Input: d - integer >= 2
    % Output:
    %         Fnf - polynomial system from normal form deflation
    %         Vnf - variables arising from normal form deflation

    % error checking
    d = 2^(s-1);
    Fmm = []; Vmm = []; Fnf = []; Vnf = [];
    if nargin < 1 || d < 2
        disp(['The degree d must be at least 2!']);
        return;
    end;

    % setup the original system
    X = sym('x',[1,s]);
    F = [X(1)^3+X(1)^2-X(2)^2];
    for i = 2:s-1
        F = [F,X(i)^3+X(i)^2-X(i+1)];
    end;
    F = [F, X(s)^2]
    
     
    % setup E -- leading terms
    E = zeros(2*d,s);
    for j = 1:d
        E(2*j-1,1) =  j-1; % x1^(j-1)
        E(2*j,1) = j-1; % x2 * x1^(j-1)
         E(2*j,2) = 1;
    end;

    % setup normal form deflation
    [Fnf,Vnf, M] = deflate(F,X,E);

	% display summary of MM deflation
    disp(['Normal form']);
    disp(['  functions: ',num2str(length(Fnf))]);
    disp(['  variables: ',num2str(length(Vnf))]);

return;
 
% We din't use this function, because the Macaulay matrix becomes too large.   
function [MM] = constructMM(F,X,deg)
    % construct MM up to degree deg ignoring the trivial first column

    numF = length(F);
    numX = length(X);

    if deg <= 0 || numF < 1 || numX < 1
        MM = [];
        return;
    elseif deg == 1
        MM = jacobian(F,X);
        return;
    end;

    % deg >= 2
    numRows = numF*nchoosek(numX+deg-1,deg-1);
    numCols = nchoosek(numX+deg,deg)-1;

    MM = zeros(numRows,numCols)*X(1);

    % setup degrees
    DegC = [];
    for j = 1:(deg-1)
      DegC = [DegC;setupDegree(j,numX)];
    end;
    DegR = [setupDegree(0,numX);DegC];
    DegC = [DegC;setupDegree(deg,numX)];

    % setup MM
    currRow = 0;
    for j = 1:size(DegR,1)
        for k = 1:numF
            currRow = currRow + 1;
            for l = 1:size(DegC,1)
                f = F(k);
                for m = 1:size(DegR,2)
                    f = f*X(m)^DegR(j,m);
                end;
                mult = 1;
                for m = 1:size(DegC,2)
                    f = diff(f,X(m),DegC(l,m));
                    mult = mult*factorial(DegC(l,m));
                end;
                MM(currRow,l) = f/mult;
            end;
        end;
    end;

return;

function [System,Variables,M] = deflate(F,X,E)
% Computes the deflated system using the parametric normal form algorithm.
% It exploits the structure of the basis to minimize the number of new variables introduced.

    mu = size(E,1);

    % setup M & newVariables
    IX = eye(length(X),length(X));
    newVariables = zeros(0,0)*X(1);
    M = zeros(mu,mu,length(X))*X(1);
    for l = 1:length(X)
      for k = 1:mu-1
       % test to see if in E
        inE = 0;
        for j = k+1:mu  
            if norm(E(j,:) - E(k,:) - IX(l,:)) == 0 
            % entry is 1
               M(j,k,l) = 1;
               inE = 1;
            end;
        end;
        if inE == 0
           %test if we can use parameters introduced earlier
          NoNew=0;
          for a=1:k-1
               for b=1:l-1
                   if norm( E(k,:)-E(a,:)-IX(b,:)) == 0 
                      A=a;
                      B=b;
                      NoNew=1;
                      break;
                      break;
                   end;
                end;
           end;
           if NoNew == 1
              M(:,k,l)=M(:,:,B)*M(:,A,l);
           else
           for j = k+1:mu
              % need to create a variable
              c = ['u_',num2str(l),'_',num2str(j),'_',num2str(k)];
              eval(['syms ',c]);
              newVariables = [newVariables, eval(c)];
              M(j,k,l) = eval(c);           
            end;
          end;
        end;
      end;
    end;
    
    % setup commutation relations
    Comm = [];
    for j = 1:length(X)
      for k = j+1:length(X)
        Comm = [Comm; reshape(M(:,:,j)*M(:,:,k) - M(:,:,k)*M(:,:,j),[mu^2,1])]; 
      end; 
    end; 
    Comm = unique(Comm);

    % setup other polynomials
    Poly = [];
    Deg  = [];
    for j = 1:3
      Deg = [Deg;setupDegree(j,length(X))];
    end;

    for j = 1:length(F)
      % setup degree 0
      P = F(j)*[1;zeros(mu-1,1)];
      for k = 1:size(Deg,1)
        Mult = 1/prod(factorial(Deg(k,:)));
        Q = eye(mu,mu);
        p = F(j);
        for l = 1:length(X)
          p = diff(p,X(l),Deg(k,l));
          if p ==0  
            break;
          else Q = Q*M(:,:,l)^Deg(k,l);
          end
        end;
        P = P + Mult*p*(Q*[1;zeros(mu-1,1)]);
      end;
      Poly = [Poly;P];
    end;

    System = [Poly;Comm];
    Variables = transpose([X newVariables]);
    System = System(subs(System,Variables,rand(size(Variables))) ~= 0);

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;