--restart
needsPackage "NumericalAlgebraicGeometry"
needsPackage "Bertini" ;
--load "JoseAwesome.m2";

fL = "/Users/abrahammcs/Documents/Mathematics/dump2"

Tolerance'Error := 0.000000000001;
--NAGtrace 0;
------------------------
-- createWitness
-----------------------
-- input:
--    F -- a list of polynomial equations --inside CC[vars]
--    d -- dimension
-- output:
--    W -- a set of witness points for V(F)
-- NOTE
--    right now, we make some change of rings
--    but if F is not in CC[vars], then when we obtain 
--    the witness set W, slice(W) will return the 0 ideal
--    
--    to avoid this, we need to input F in CC[vars]
--    so we can assume F\in CC[vars] and we can remove
--    all the ring changes
----------------------
--////////////////////// 
-- Once Bertini Package gets reliable 
-- we construct the witness set from Bertini
-- using bertiniPosDimSolve
--//////////////////////
createWitness = method()
createWitness(List,ZZ) := (F,d) ->(
    RF := ring matrix{F};
    vrs := gens RF;
    RC := CC_53[vrs];
    n:=#vrs; -- need to create the witness points wrt CC_53
    coeffs := random(CC_53^(n+1),CC_53^(d));
    L0 := flatten entries (sub(vars RF|matrix{{1_RF}},RC)*coeffs);
    Pts := solveSystem join(F/(i->sub(i,RC)),L0);
    --Pts := bertiniZeroDimSolve join(F/(i->sub(i,RC)),L0);
    witnessSet(ideal F, ideal L0, Pts)
    )

--------------------------
-- create Slice through P
--------------------------
-- pSlice
--------------------------
-- Input
--    p -- point given as a list of numbers in CC[vars]
--    k -- ZZ representing the dimension of the slice
-- Output
--    L' -- List of Equations of a random real plane
--     	   of dimension k through p
-------------
-- one way to push P to R is
-- P = {0.6_R, 0_R}
--
--- or flatten entries sub(matrix{P},R)
pSlice = method()
pSlice(List,ZZ) := (P,k) ->(
    Pm := matrix{P};
    Rc := ring Pm;
    n := numgens Rc;
    coeffs := sub(random(RR^k,RR^n),Rc);
    B := coeffs*transpose(Pm);
    flatten entries ((coeffs|B)*transpose(vars Rc|matrix{{-1_Rc}}))    
    )
------------------
--    dist
------------------
-- compute the (square of) the euclidean distance
-- of two points in the affine space A^n
------------------
-- Input:
--    v,u    -- two points given as list of coordinates
-- Output:
--    d	   -- a number representing (v-u)*(v-u)^t = ||v - u||^2
dist =method()
dist(List,List) := (v,u) ->(
    d := v-u;
    dm := matrix {d};
    first first entries (dm*transpose(dm))
    )

--------------------
-- gradDescHomot
--------------------
-- given a system F and two points x0 and P
-- computes the points in F which are locally
-- close to P using gradient descent homotopy
--------------------
-- Input
--    F	   -- List of Poly Equations (List)
--    x0    -- witness point of F (List)
--    P	   -- point outside V(F) (List)
-- Output
--    y0    -- point in F satisfying
--    	  arg min {dist(x0,P) : x0 \in V(F)}
---------------------
--     Note (6 March, 2015)
-- We need to modify this function to include the gradient \nabla\OF
-- as part of the input of the function. Right now, it only works for
-- Euclidean distance... but we want more general.
--
-- Also, it is using Antons tracking function but we can use Bertini 
---------------------
gradDescHomot = method()
gradDescHomot(List,List,List) := (F,X,P) ->(
    MF := matrix{F};
    Rc := ring MF;
    n := #F; --number of equations
    N := numgens Rc; -- number of variables
    gradF := diff(transpose vars Rc, MF); --gradient of F
    RH := (coefficientRing Rc)[gens Rc,l_0..l_n,t];--homotopy ring
    -----------------------------------
    -- maybe we need to increase this cleaning to a higher precision
    K:=(X- P)/clean_0.00001;
    K = apply(K, k->sub(k, RH));
    gradF = sub(gradF,RH);
    lambs := transpose matrix{toList(l_0..l_n)};
    h0 := (transpose(sub(matrix{gens Rc - P},RH))|gradF)*lambs;-- equations for the linear span of X-P and the gradient of F
    F1 := transpose(sub(MF,RH));
    -- we need to take a random real affine projection of lambdas
    a := sub(random(RR^1,RR^(n+1)),RH);
    aeq := (matrix{{1_RH}}|a)*(lambs||matrix{{-1_RH}}); -- equations for the random values of lambda
    a0 := last flatten entries a;
		-- we need to change a0 for a random complex number for the gamma trick t*g*K:
    Hmt := F1|| h0-transpose(t*a0*matrix{K})||aeq; -- homotopy system
    --entries Hmt
    Rh := (coefficientRing RH)[drop(gens RH,-1)]; -- homotopy ring for t=0 
    map't'0 := map(Rh, RH,vars Rh | matrix{{0_CC}});
    map't'1 := map(Rh, RH, vars Rh | matrix{{1_CC}});
    Solstart:= X | prepend(a0,apply(n,i->0)); -- initial solution = (x,a0,0...0)
    Start := flatten(entries(map't'1(Hmt))); -- start system
    Target := flatten(entries(map't'0(Hmt)));--target system 
    --apply(Start, T-> sub(T,matrix{Solstart}))/clean_0.0001
    pts := track(Start,Target,{Solstart}, tStepMin=>1e-10);
    pts = refine(Target,pts, ErrorTolerance=>1e-20);
    --<< "this is the type of point we found: "<< pts/status<< "  and is it real? "<< pts/isRealPoint<<endl;
    pts =flatten(pts/coordinates)--;
    --take(pts,N)
    )

--------------------
-- selectRealPts
--------------------
-- given a list of complex points, we will select
-- those whose coordinates have an imaginary part
-- strictly smaller than an Error Tolerance Bound
--------------------
-- Input
--    W	   -- List of Points
-- Output
--    WR    -- sublist of real points
---------------------
selectRealPts = method()
selectRealPts List := W ->(
    n:=#first(W);
    istrue := toList(n:true);
    imgparts := apply(W, w->(
	    w/(i->abs(imaginaryPart(i))<Tolerance'Error)
	    ));
    indx := positions(imgparts, i->i==istrue);
    W_indx
    )


--------------------
-- critPointRandSlice
--------------------
-- given a variety V(F) and point P outside
-- finds all crit points in V(F) that minimizes
-- (locally) the distance to P by taking a random slice
--------------------
-- Input :
--    P -- a point given as a list of real numbers
--    W -- a witness set for V(F)
--    d -- the codimension of V(F)
--------------------
--    Q -- a set of real points that are in V(F) and  
--    	   are local minimal to the distance to P
--    	   taken from a random slice through P
--------------------
critPointRandSlice = method()
critPointRandSlice (List,WitnessSet,ZZ) := (P,W,d) ->(
    L' := pSlice(P,d);
    F := equations W;
    L := slice W;
    wpts := W#Points/coordinates;
    W' := track(join(F,L), join(F,L'), wpts);
    W'Real := selectRealPts((W')/coordinates);
    out:= apply(W'Real, x0->(
	    y0 := gradDescHomot(F,x0,P);
	    if (norm(sub(matrix{F},matrix{y0}))<0.0000001) then(
		y0
		)else(
		<<"This is a point outside the variety "<< y0<<endl;
		)
	    ));
    delete(null, out)
    )


--------------------
-- uniquePoints
--------------------
-- function to collect only
-- the unique solutions of a list
-- of many random samples
-- i.e. it removes the repetitions
-- by comparing if the numerical sols are equal
-------------------
-- Input:
--    Sampl -- list of points (given as lists)
-- Output:
--    Unique -- list of points that are not repeated
-------------------
uniquePoints = method()
uniquePoints List := Pts ->(
    Unique := {first Pts};
    scan(Pts, p->(
	    newpt:=true;
	    scan(Unique, q->(
		    if areEqual(p,q) then(
			newpt=false;
			break
			);
		    ));
	    if newpt then Unique = append(Unique, p);
	    ));
    Unique
    )

-------------------
-- critPointsfromW
-------------------
-- input:    
--    W    -- list of real points in V(F)
--    F    -- list of equations
--    P    -- coordinates of the point outside V
-------------------
--    Crit  -- list of critical points gotten from
--             running gradDescHomot starting from 
--     	       the points in W
-------------------
critPointsFromW = method()
critPointsFromW (List,List,List) := (W,F,P) ->(
    Crit:=apply(W, x0->(
	    y0 := gradDescHomot(F,x0,P);
	    if (norm(sub(matrix{F},matrix{y0}))<0.000001) then y0
	    ));
    Crit
    )

------------------------
--    TEST for pSlice
------------------------
--/// TEST
--R = CC_53[x,y]
--F = {y^2-x^3+3*x^2-2*x}
--P = {0.6,0}/(i->promote(i, R))

--pSlice(P,1)

--R = CC[x,y,z]
--F = {x^4*z-y^5}
--
--P = {1_R,2_R,1_R}
--pSlice(P,2)
--///
------------------------

-- lagrangeEq
-- Input: ??
--
lagrangeEq= (jacH,gradObjectiveFunction,theData,pUnknowns,lambdas)->(
    if class jacH ===Matrix then jacH=entries jacH;
    if #jacH_0=!= #pUnknowns then return error "Incorrect number of pUnknowns or jacH is sized wrong";
    if #jacH=!= #lambdas then return error "Incorrect number of lambdas or jacH is sized wrong";
    if #theData=!= #theData then return error "#theData and #pUnknowns do not agreee";
    nonZeroData:=theData;
    return for j to #pUnknowns-1 list (
	 gradObjectiveFunction_j-(sum for i to #lambdas-1 list lambdas_i*((jacH_i)_j ) ) )
	)