% Copyright 2015 by Jonathan D. Hauenstein

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

function [FA,VA,FB,VB,FC,VC,FD,VD] = deflate53_1()
    % Output: FA - polynomial system from method A
    %         VA - variables arising from method A
    %         FB - polynomial system from method B
    %         VB - variables arising from method B
    %         FC - polynomial system from method C
    %         VC - variables arising from method C
    %         FD - polynomial system from method D
    %         VD - variables arising from method D
    
    [FA,VA] = deflate53_1_A();
    
    % display summary of method A
    disp(['Method A']);
    disp(['  functions: ',num2str(length(FA))]);
    disp(['  variables: ',num2str(length(VA))]);    
    
    [FB,VB] = deflate53_1_B();
    
    % display summary of method B
    disp(['Method B']);
    disp(['  functions: ',num2str(length(FB))]);
    disp(['  variables: ',num2str(length(VB))]);    
    
    [FC,VC] = deflate53_1_C();
    
    % display summary of method C
    disp(['Method C']);
    disp(['  functions: ',num2str(length(FC))]);
    disp(['  variables: ',num2str(length(VC))]);    
    
	[FD,VD] = deflate53_1_D();
    
    % display summary of method D
    disp(['Method D']);
    disp(['  functions: ',num2str(length(FD))]);
    disp(['  variables: ',num2str(length(VD))]);    

return;

function [FA,VA] = deflate53_1_A()
    % deflate using method A -- intrinsic slicing deflation

    syms x1 x2 x3 x4;
    X = [x1 x2 x3 x4];

    F = [x1^4 - x2*x3*x4;
         x2^4 - x1*x3*x4;
         x3^4 - x1*x2*x4;
         x4^4 - x1*x2*x3];

    % first iteration -- Jacobian is identically zero at Pt = [0 0 0 0]
    J = jacobian(F,X);

    % deflate -- use random vector
    Fd = [F;J*randi([-10,10],4,1)];

    % second iteration -- Jacobian is identically zero at Pt
    Jd = jacobian(Fd,X);
    
    % deflate -- use random vector
    Fdd = [Fd;Jd*randi([-10,10],4,1)];

    % jacobian is full rank -- find distinct nonzero polynomials
    FA = unique(Fdd);
    FA = FA(subs(FA,X,rand(size(X))) ~= 0);
    VA = X;
    
return;

function [FB,VB] = deflate53_1_B()
    % deflate using method B -- isosingular defl. with full rank submatrix

    syms x1 x2 x3 x4;
    X = [x1 x2 x3 x4];

    F = [x1^4 - x2*x3*x4;
         x2^4 - x1*x3*x4;
         x3^4 - x1*x2*x4;
         x4^4 - x1*x2*x3];

    % first iteration -- Jacobian is identically zero at Pt = [0 0 0 0]
    J = jacobian(F,X);

    % deflate -- add all entries of Jacobian
    Fd = [F;reshape(J,16,1)];

    % second iteration -- Jacobian is identically zero at Pt
    Jd = jacobian(Fd,X);

    % deflate -- add all entries of Jacobian
    Fdd = [Fd;reshape(Jd,80,1)];

    % jacobian is full rank -- find distinct nonzero polynomials
    FB = unique(Fdd);
    FB = FB(subs(FB,X,rand(size(X))) ~= 0);
    VB = X;
    
return;

function [FC,VC] = deflate53_1_C()
    % deflate using method C -- "kerneling" deflation

    syms x1 x2 x3 x4;
    X = [x1 x2 x3 x4];

    F = [x1^4 - x2*x3*x4;
         x2^4 - x1*x3*x4;
         x3^4 - x1*x2*x4;
         x4^4 - x1*x2*x3];

    % first iteration -- Jacobian is identically zero at Pt = [0 0 0 0]
    J = jacobian(F,X);

    % deflate -- add all partial derivatives
    Fd = [F;reshape(J,[numel(J),1])];

    % second iteration -- Jacobian is identically zero at Pt
    Jd = jacobian(Fd,X);
    
    % deflate -- add all partial derivatives
    Fdd = [Fd;reshape(Jd,[numel(Jd),1])];

    % jacobian is full rank -- find distinct nonzero polynomials
    FC = unique(Fdd);
    FC = FC(subs(FC,X,rand(size(X))) ~= 0);
    VC = X;
    
return;

function [FD,VD] = deflate53_1_D()
    % deflate using method D -- |i|=1 deflation

    syms x1 x2 x3 x4;
    X = [x1 x2 x3 x4];

    F = [x1^4 - x2*x3*x4;
         x2^4 - x1*x3*x4;
         x3^4 - x1*x2*x4;
         x4^4 - x1*x2*x3];

    % first iteration -- Jacobian is identically zero at Pt = [0 0 0 0]
    J = jacobian(F,X);

    % deflate -- use random vector
    Fd = [F;J*randi([-10,10],4,1)];

    % second iteration -- Jacobian is identically zero at Pt
    Jd = jacobian(Fd,X);
    
    % deflate -- use random vector
    Fdd = [Fd;Jd*randi([-10,10],4,1)];

    % jacobian is full rank -- find distinct nonzero polynomials
    FD = unique(Fdd);
    FD = FD(subs(FD,X,rand(size(X))) ~= 0);
    VD = X;
    
return;
