-----------------------------------------------------
-- Examples for 
-- "Determinantal representations of hyperbolic curves 
--  via polynomial homotopy continuation"
-- by Anton Leykin and Daniel Plaumann 
-- (May 2016)
-----------------------------------------------------
restart
load "HyperbolicPolynomialHomotopy.m2"
setRandomSeed 0
(XYT,PhiDR,Phi) = determinantalRepresentation(4,CC);
DR = matrix{toList(1..10)|{1,2,3,4_CC}} -- some real (D,R)
getDandR(4,DR)

-- direct (real) N-path, i.e., s\in [0,1]
time DRend'Nuij = trackNuijHomotopy(DR,XYT,PhiDR,Phi,Randomized=>false)

-- take a "random" piecewise linear path for s in CC
random's = 1+ii
time DRsegment1 = trackNuijHomotopy(DR,XYT,PhiDR,Phi,TargetS=>random's,Randomized=>false)
time DRsegment2 = trackNuijHomotopy(DRsegment1,XYT,PhiDR,Phi,StartS=>random's,NuijPathEndDR=>DR,Randomized=>false)

-- check the results
CC'XYT = CC[gens XYT]
p = sub( sub(PhiDR, vars XYT|DR), CC'XYT )
p'Nuij = makeNuijHomotopy(p,Randomized=>false);

sXYT = ring p'Nuij;
at's = s -> map(CC'XYT,sXYT,matrix{{toCC s}}|vars CC'XYT); -- specialize value of s

p'end'Nuij = sub( sub(PhiDR, vars XYT|DRend'Nuij), CC'XYT )
assert( clean_0.0001 (
    	p'end'Nuij - (at's 0) p'Nuij -- this should be close to zero
    	) == 0 )
p'end'Nuij'random = sub( sub(PhiDR, vars XYT|DRsegment2), CC'XYT )
assert( clean_0.0001 (
	p'end'Nuij - (at's 0) p'Nuij -- this should be close to zero
    	) == 0 )

------------------------------------------------------------------------
-- Nuij path via SLPs (expanding determinants takes a long time for d>5)
------------------------------------------------------------------------
restart
load "HyperbolicPolynomialHomotopy.m2"
setRandomSeed 1; 
d = 6 -- we could successfully compute with up to d = 10
{*
d = 8
d = 10
*}
bc = binomial(d+1,2)
DR = matrix{toList ((1..bc)/(a->a%3+1.)|(-d//2..d-1-d//2))} -- "random" (D,R)
{* 
-- for large d, the pair (D,R) above could be too close to the singular locus; take another random choice, e.g.:
DR = matrix{toList ((1..bc)/(a->random RR)|(-d//2..d-1-d//2))}
*}
{* 
-- for large d, higher precision may be necessary, e.g.: 
DR = sub(DR,RR_100)
*}
getDandR(d,DR)

NAGtrace 10

(homotopyRing, PhiDR'minus'nuij'homotopy) = makeNuijPathSLP(d,DR,1,0,Randomized=>false) 
elapsedTime solution = trackHomotopy ((homotopyRing,PhiDR'minus'nuij'homotopy), {transpose DR}, 
    Software=>M2); -- takes 28 steps for d = 6
                   -- 1810.05 seconds 	  
	 --  d = 8 -- 18781.8 seconds (Randomized, 13847.5 on RedHat)
    	    	     
-- ALTERNATIVELY: use randomize N-path: 
(homotopyRing, PhiDR'minus'nuij'homotopy) = makeNuijPathSLP(d,DR,1,0,Randomized=>true) 
elapsedTime solution = trackHomotopy ((homotopyRing,PhiDR'minus'nuij'homotopy), {transpose DR}, 
    Software=>M2); -- (usually finishes faster than nonrandomized) 
fixed'polyDR = matrix first solution
(D,R) = getDandR(d,fixed'polyDR)/clean_0.0001

----------------------------------------------------------------------------------------------------
{* latex and external string:
tex D	    
tex R
print toString fixed'polyDR
-- RESTART here to track N-path from the fixed point with the det. rep. computed above:
restart
load "HyperbolicPolynomialHomotopy.m2"
d = 6
fixed'polyDR = matrix {{6, 2.51352, 1.19571, 4.04309, 1.42786, -1.98597, 6, 3.08656, .468873+2e-7*ii, 2.38468, 1.05948, 6, .785785, 4.66027, 2.29433, 6, 1.6226, .933245, 6, 3.50198, 6, .222847, 1.18893, 2.99274, 5.77514, 9.83747, 15.9829}}
*}

CC[x,y,t];
target'poly = -36*x^6-157*x^4*y^2-20*x^3*y^3-109*x^2*y^4+246*x*y^5-92*y^6-12*x^3*y^2*t+90*x^2*y^3*t+10*x*y^4*t+76*y^5*t+49*x^4*t^2+156*x^2*y^2*t^2-16*x*y^3*t^2+132*y^4*t^2+12*x*y^2*t^3-14*y^3*t^3-14*x^2*t^4-27*y^2*t^4+t^6;
(homotopyRing, PhiDR'minus'nuij'homotopy) = makeNuijPathSLP(d,target'poly,0,1,Randomized=>false)
NAGtrace 10; gbTrace = 0;
target'polyDR = matrix first trackHomotopy ((homotopyRing,PhiDR'minus'nuij'homotopy), {transpose fixed'polyDR}, 
    Software=>M2) -- takes 15 steps
(D,R) = getDandR(d,target'polyDR)/clean_0.0001

{* latex and external string:
tex D	    
tex R
print toString clean_0.0001 target'polyDR
-- CHECK:
restart
load "HyperbolicPolynomialHomotopy.m2"
d=6
target'polyDR = matrix {{0, .596508, -1.43241, 2.00316, 1.10471, -.725394, 0, .739773, 1.79407, .0604427, -1.60948, 0, 1.56816, 1.66137, -.165953, 0, .839374, 2.00885, 0, 1.57679, 0, -3, -2, -1, 1, 2, 3}}
(D,R) = getDandR(d,target'polyDR)/clean_0.0001
hPoly = hyperbolicPolynomial (D,R)
use ring hPoly 
(x,y,t) = toSequence gens ring hPoly;
target'poly = -36*x^6-157*x^4*y^2-20*x^3*y^3-109*x^2*y^4+246*x*y^5-92*y^6-12*x^3*y^2*t+90*x^2*y^3*t+10*x*y^4*t+76*y^5*t+49*x^4*t^2+156*x^2*y^2*t^2-16*x*y^3*t^2+132*y^4*t^2+12*x*y^2*t^3-14*y^3*t^3-14*x^2*t^4-27*y^2*t^4+t^6;
assert(clean_0.01(hPoly - target'poly)  == 0)
*}