% Copyright 2015 by Jonathan D. Hauenstein

% Deflate Example 2 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_2()
    % 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_2_A();
    
    % display summary of method A
    disp(['Method A']);
    disp(['  functions: ',num2str(length(FA))]);
    disp(['  variables: ',num2str(length(VA))]);    
    
    [FB,VB] = deflate53_2_B();
    
    % display summary of method B
    disp(['Method B']);
    disp(['  functions: ',num2str(length(FB))]);
    disp(['  variables: ',num2str(length(VB))]);    
    
    [FC,VC] = deflate53_2_C();
    
    % display summary of method C
    disp(['Method C']);
    disp(['  functions: ',num2str(length(FC))]);
    disp(['  variables: ',num2str(length(VC))]);    

	[FD,VD] = deflate53_2_D();
    
    % display summary of method D
    disp(['Method D']);
    disp(['  functions: ',num2str(length(FD))]);
    disp(['  variables: ',num2str(length(VD))]);    

return;

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

    syms x y z;
    X = [x y z];

    F = [x^4;
         x^2*y + y^4;
         z + z^2 - 7*x^3 - 8*x^2];
     
    % first iteration -- Jacobian has 2 dim'l null space at Pt = [0 0 -1]
    J = jacobian(F,X);

    % deflate
    syms u;
    X = [X u];
    Fd = [F;J*(randi([-10,10],3,1)*u + randi([-10,10],3,1))];

    % second iteration -- Jacobian still has 2 dim'l null space
    Jd = jacobian(Fd,X);
    
    % deflate
    syms v1 v2;
    X = [X v1 v2];
    Fdd = [Fd;Jd*(randi([-10,10],4,2)*[v1;v2] + randi([-10,10],4,1))];

    % third iteration -- Jacobian has 1 dim'l null space
    Jdd = jacobian(Fdd,X);
    
    % deflate
    syms w1 w2 w3 w4 w5;
    X = [X w1 w2 w3 w4 w5];
    Fddd = [Fdd;Jdd*(randi([-10,10],6,5)*[w1;w2;w3;w4;w5] + randi([-10,10],6,1))];

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

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

    syms x y z;
    X = [x y z];

    F = [x^4;
         x^2*y + y^4;
         z + z^2 - 7*x^3 - 8*x^2];

    % first iteration -- Jacobian has 2 dim'l null space at Pt = [0 0 -1]
    J = jacobian(F,X);

    % deflate -- add 2x2 minors containing (3,3)
    Fd = F;
    for j = 1:2
        for k = 1:2
            Fd = [Fd;det(J([j,3],[k,3]))];
        end;
    end;

    % second iteration -- Jacobian still has 2 dim'l null space
    Jd = jacobian(Fd,X);
    
    % deflate -- add 2x2 minors containing (3,3)
    Fdd = Fd;
    for j = [1:2,4:7]
        for k = 1:2
            Fdd = [Fdd;det(Jd([j,3],[k,3]))];
        end;
    end;

    % 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_2_C()
    % deflate using method C -- "kerneling" deflation

    syms x y z;
    X = [x y z];

    F = [x^4;
         x^2*y + y^4;
         z + z^2 - 7*x^3 - 8*x^2];

    % first iteration -- Jacobian has 2 dim'l null space at Pt = [0 0 -1]
    J = jacobian(F,X);
    A = J(3,3);
    B = J(3,1:2);
    C = J(1:2,3);
    D = J(1:2,1:2);
    N = det(A)*(D - C*inv(A)*B);

    % deflate
    Fd = [F;reshape(N,[numel(N),1])];
    
    % second iteration -- Jacobian has 2 dim'l null space with same
    % nonzero location
    Jd = jacobian(Fd,X);
    A = Jd(3,3);
    B = Jd(3,1:2);
    C = Jd([1:2,4:end],3);
    D = Jd([1:2,4:end],1:2);
    N = det(A)*(D - C*inv(A)*B);
    
    % deflate
    Fdd = [Fd;reshape(N,[numel(N),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_2_D()
    % deflate using method D -- |i|=1 deflation

    syms x y z;
    X = [x y z];

    F = [x^4;
         x^2*y + y^4;
         z + z^2 - 7*x^3 - 8*x^2];

    % first iteration -- Jacobian has 2 dim'l null space at Pt = [0 0 -1]
    J = jacobian(F,X);
    A = J(3,3);
    B = J(3,1:2);
    N = det(A)*[eye(2,2);-inv(A)*B]*randi([-10,10],2,1);

    % deflate
    Fd = [F;J(1:2,:)*N];
    
    % second iteration -- Jacobian has 2 dim'l null space with same
    % nonzero location
    Jd = jacobian(Fd,X);
    
    % deflate
    Fdd = [Fd;Jd([1:2,4:end],:)*N];

    % third iteration -- Jacobian has 1 dim'l null space
    Jdd = jacobian(Fdd,X);    
    A = Jdd([3,9],[2,3]);
    B = Jdd([3,9],1);
    N = det(A)*[1;-inv(A)*B];
    
    % deflate
    Fddd = [Fdd;Jdd([1:2,4:8],:)*N];
    
    % jacobian is full rank -- find distinct nonzero polynomials
    FD = unique(Fddd);
    FD = FD(subs(FD,X,rand(size(X))) ~= 0);
    VD = X;
    
return;
