% Copyright 2015 by Jonathan D. Hauenstein

% Deflate Example 4 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_4()
    % 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_4_A();
    
    % display summary of method A
    disp(['Method A']);
    disp(['  functions: ',num2str(length(FA))]);
    disp(['  variables: ',num2str(length(VA))]);    
    
    [FB,VB] = deflate53_4_B();
    
    % display summary of method B
    disp(['Method B']);
    disp(['  functions: ',num2str(length(FB))]);
    disp(['  variables: ',num2str(length(VB))]);    
    
    [FC,VC] = deflate53_4_C();
    
    % display summary of method C
    disp(['Method C']);
    disp(['  functions: ',num2str(length(FC))]);
    disp(['  variables: ',num2str(length(VC))]);    

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

return;

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

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

    F = [2*x1 + 2*x1^2 + 2*x2 + 2*x2^2 + x3^2 - 1; 
    	 (x1 + x2 - x3 - 1)^3 - x1^3; 
    	 (2*x1^3 + 2*x2^2 + 10*x3 + 5*x3^2 + 5)^3 - 1000*x1^5];

    % 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))];

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);

    % deflate
    syms r1 r2 r3 r4 r5 r6 r7 r8 r9 r10;
    X = [X r1 r2 r3 r4 r5 r6 r7 r8 r9 r10];
    Fdddd = [Fddd;Jddd*(randi([-10,10],11,10)*[r1;r2;r3;r4;r5;r6;r7;r8;r9;r10] + randi([-10,10],11,1))];

    % fifth iteration -- Jacobian still has 1 dim'l null space
    Jdddd = jacobian(Fdddd,X);

    % deflate
    syms s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20;
    X = [X s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20];
    Fddddd = [Fdddd;Jdddd*(randi([-10,10],21,20)*[s1;s2;s3;s4;s5;s6;s7;s8;s9;s10;s11;s12;s13;s14;s15;s16;s17;s18;s19;s20] + randi([-10,10],21,1))];
    
    % jacobian is full rank -- find distinct nonzero polynomials
    FA = unique(Fddddd);
    FA = FA(subs(FA,X,rand(size(X))) ~= 0);
    VA = X;
    
return;

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

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

    F = [2*x1 + 2*x1^2 + 2*x2 + 2*x2^2 + x3^2 - 1; 
    	 (x1 + x2 - x3 - 1)^3 - x1^3; 
    	 (2*x1^3 + 2*x2^2 + 10*x3 + 5*x3^2 + 5)^3 - 1000*x1^5];

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

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

    % third iteration -- Jacobian has 1 dim'l null space
    Jdd = jacobian(Fdd,X);

    % deflate -- add 3x3 minors containing [1,12] X [1,2]
    Fddd = Fdd;
    for j = [2:11,13:19]
        Fddd = [Fddd;det(Jdd([1,12,j],:))];
    end;

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);

    % deflate -- add 3x3 minors containing [1,12] X [1,2]
    Fdddd = Fddd;
    for j = [2:11,13:36]
        Fdddd = [Fdddd;det(Jddd([1,12,j],:))];
    end;

    % fifth iteration -- Jacobian still has 1 dim'l null space
    Jdddd = jacobian(Fdddd,X);

	% deflate -- add 3x3 minors containing [1,12] X [1,2]
    Fddddd = Fdddd;
    for j = [2:11,13:70]
        Fddddd = [Fddddd;det(Jdddd([1,12,j],:))];
    end;

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

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

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

    F = [2*x1 + 2*x1^2 + 2*x2 + 2*x2^2 + x3^2 - 1; 
    	 (x1 + x2 - x3 - 1)^3 - x1^3; 
    	 (2*x1^3 + 2*x2^2 + 10*x3 + 5*x3^2 + 5)^3 - 1000*x1^5];

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

    % deflate
    Fd = [F;reshape(N,[numel(N),1])];

    % second iteration -- Jacobian still has 2 dim'l null space
    Jd = jacobian(Fd,X);
	A = Jd(1,1);
    B = Jd(1,2:3);
    C = Jd(2:end,1);
    D = Jd(2:end,2:3);
    N = det(A)*(D - C*inv(A)*B);

    % deflate
    Fdd = [Fd;reshape(N,[numel(N),1])];

    % third iteration -- Jacobian has 1 dim'l null space
    Jdd = jacobian(Fdd,X);
    A = Jdd([1,10],1:2);
    B = Jdd([1,10],3);
    C = Jdd([2:9,11:end],1:2);
    D = Jdd([2:9,11:end],3);
    N = det(A)*(D - C*inv(A)*B);
    
    % deflate
    Fddd = [Fdd;reshape(N,[numel(N),1])];

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);
    A = Jddd([1,10],1:2);
    B = Jddd([1,10],3);
    C = Jddd([2:9,11:end],1:2);
    D = Jddd([2:9,11:end],3);
    N = det(A)*(D - C*inv(A)*B);

    % deflate
    Fdddd = [Fddd;reshape(N,[numel(N),1])];

    % fifth iteration -- Jacobian still has 1 dim'l null space
    Jdddd = jacobian(Fdddd,X);
    A = Jdddd([1,10],1:2);
    B = Jdddd([1,10],3);
    C = Jdddd([2:9,11:end],1:2);
    D = Jdddd([2:9,11:end],3);
    N = det(A)*(D - C*inv(A)*B);

    % deflate
    Fddddd = [Fdddd;reshape(N,[numel(N),1])];    

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

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

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

    F = [2*x1 + 2*x1^2 + 2*x2 + 2*x2^2 + x3^2 - 1; 
    	 (x1 + x2 - x3 - 1)^3 - x1^3; 
    	 (2*x1^3 + 2*x2^2 + 10*x3 + 5*x3^2 + 5)^3 - 1000*x1^5];

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

    % deflate
    Fd = [F;J(2:3,:)*N];

    % second iteration -- Jacobian still has 2 dim'l null space
    Jd = jacobian(Fd,X);
    
    % deflate
    Fdd = [Fd;Jd(2:end,:)*N];

    % third iteration -- Jacobian has 1 dim'l null space
    Jdd = jacobian(Fdd,X);
    A = Jdd([1,8],[1,2]);
    B = Jdd([1,8],3);
    N = det(A)*[-inv(A)*B;1];
    
    % deflate
    Fddd = [Fdd;Jdd([2:7,9:end],:)*N];

    % fourth iteration -- Jacobian still has 1 dim'l null space
    Jddd = jacobian(Fddd,X);

    % deflate
    Fdddd = [Fddd;Jddd([2:7,9:end],:)*N];

    % fifth iteration -- Jacobian still has 1 dim'l null space
    Jdddd = jacobian(Fdddd,X);

    % deflate
    Fddddd = [Fdddd;Jdddd([2:7,9:end],:)*N];    

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

