% Copyright 2015 by Jonathan D. Hauenstein

% Deflate Caprasse system from Section 5.2 of 
% "Certifying isolated singular points and their multiplicity structure" 
% by J.D. Hauenstein, B. Mourrain, and A. Szanto

function [Fmm,Vmm,Fnf,Vnf] = deflate52_structure()
    % Output: Fmm - polynomial system from Macaulay matrix deflation
    %         Vmm - variables arising from Macaulay matrix deflation
    %         Fnf - polynomial system from normal form deflation
    %         Vnf - variables arising from normal form deflation

    % setup the original system
    syms x1 x2 x3 x4;
    X = [x1 x2 x3 x4];
    F = [x1^3*x3-4*x1^2*x2*x4-4*x1*x2^2*x3-2*x2^3*x4-4*x1^2-4*x1*x3+10*x2^2+10*x2*x4-2; 
    	 x1*x3^3-4*x1*x3*x4^2-4*x2*x3^2*x4-2*x2*x4^3-4*x1*x3+10*x2*x4-4*x3^2+10*x4^2-2; 
         2*x1*x2*x4+x2^2*x3-2*x1-x3; 
         x1*x4^2+2*x2*x3*x4-x1-2*x3];
    
    % setup the Macaulay matrix up to degree 2
    MM = constructMM(F,X,2);

    % setup the null space -- dimension 3
    NullCols = 3; % # of columns of null space
    Cols = size(MM,2); % # of columns of MM = # of rows of null space
    
    A = zeros(Cols-NullCols,NullCols)*x1;
    for j = (NullCols+1):Cols
        for k = 1:NullCols
            c = ['n_',num2str(j),'_',num2str(k)];
            eval(['syms ',c]);
            A(j-NullCols,k) = eval(c);
        end;
    end;
    Vmm = [X symvar(A)];
    N = randi([-10,10],Cols,Cols)*[eye(NullCols,NullCols);A];
    
    % setup Fmm
    Q = MM*N;
    Fmm = F;
    for j = 1:NullCols
        Fmm = [Fmm;Q(:,j)];
    end;
        
    % display summary of MM deflation
    disp(['Macaulay matrix']);
    disp(['  functions: ',num2str(length(Fmm))]);
    disp(['  variables: ',num2str(length(Vmm))]);
    
    % setup E -- leading terms
    E = [0 0 0 0; 1 0 0 0; 0 1 0 0; 1 1 0 0]; % 1, x1, x2, x1*x2
    
    % setup normal form deflation
    [Fnf,Vnf] = deflate(F,X,E);

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

return;
    
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] = deflate(F,X,E)

    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 j = 2:mu
        for k = 1:j-1  
          if norm(E(j,:) - E(k,:) - IX(l,:)) == 0 
            % entry is 1
            M(j,k,l) = 1;
          else
            % test to see if in E
            inE = 0;
            for m = 1:mu
              if norm(E(m,:) - E(k,:) - IX(l,:)) == 0
                inE = 1;
              end;
            end;

            if inE == 0
              % 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:(mu-1)
      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)
          Q = Q*M(:,:,l)^Deg(k,l);
          p = diff(p,X(l),Deg(k,l));
        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;