--
-- Config Bertini and load package
--

-- Put path to Bertini in file config
-- or define path to bertini here
bertini_path = get "./config"

needsPackage("Bertini", Configuration=>{"BERTINIexecutable"=>bertini_path});
needsPackage("Graphs");

load "opt_tools.m2";

--
-- General commments
--

-- squaring and filtering, not the way to do it
-- try computing startpoints carefully

--
-- General config
--

eps = 1e-8
nvars = 3;
-- only FF=QQ works for now
FF = QQ;

nvarsm1 = nvars-1;
RX = FF[x_0..x_nvarsm1];
RY = FF[y_0..y_nvarsm1];

--
-- Define input
--

-- random
-- when FF=QQ use I to get complex random
-- f0r = randcompint(RX, 3, 1);
-- f0i = randcompint(RX, 3, 1);
-- g0r = randcompint(RY, 2, 1);
-- g0i = randcompint(RY, 2, 1);

-- pairf = ratnormalcomplex(RX)
-- f0r = pairf#0;
-- f0i = pairf#1;
-- opts = for i from 0 to nvars-1 list ((gens RX)_i) => ((gens RY)_i);
-- g0r = sub(f0r, opts);
-- g0i = sub(f0i, opts);

-- gensf = gens f0r;
-- f1 = gensf_(0,0);
-- f2 = gensf_(0,1);
-- fsq = ideal(f1^2+f2^2-1/10);
-- opts = for i from 0 to nvars-1 list x_i => y_i;
-- g0r = sub(fsq, opts);
-- g0i = sub(fsq, opts);

-- when FF=QQ use I to get complex random
-- f0 = ratnormal(RX, true);
-- g0 = ratnormal(RY, true);

-- ratnormal
pairf = ratnormalcomplex(RX)
f0r = pairf#0;
f0i = pairf#1;
pairg = ratnormalcomplex(RY)
g0r = pairg#0;
g0i = pairg#1;

-- elliptic curve
-- pairf = ellipticver(RX,2);
-- f0r = pairf#0;
-- f0i = pairf#1;
-- pairg = ellipticver(RY,2);
-- g0r = pairg#0;
-- g0i = pairg#1;

--
--
--

numgenX0 = numColumns(gens f0r);
numgenY0 = numColumns(gens g0r);

-- Fill in dimensions
dimX = 1;
dimY = 1;

Cf = random(RX^(nvars-dimX), RX^numgenX0);
Cg = random(RY^(nvars-dimY), RY^numgenY0);

fr = if (nvars-dimX) < numgenX0 then (ideal(Cf*transpose(gens f0r))) else f0r;
gr = if (nvars-dimY) < numgenY0 then (ideal(Cg*transpose(gens g0r))) else g0r;

fi = if (nvars-dimX) < numgenX0 then (ideal(Cf*transpose(gens f0i))) else f0i;
gi = if (nvars-dimY) < numgenY0 then (ideal(Cg*transpose(gens g0i))) else g0i;

numgenX = numColumns(gens fr);
numgenY = numColumns(gens gr);

--
--
--

dfr = jacobian(fr);
dgr = jacobian(gr);

dfi = jacobian(fi);
dgi = jacobian(gi);


R = FF[x_0..x_nvarsm1,nu_0,lambda_1..lambda_numgenX,y_0..y_nvarsm1,nu_1,mu_1..mu_numgenY,s,I]/ideal(I^2+1);

-- do something when FF is not CC or QQ
gamma = if (FF === QQ) then (random(FF)+random(FF)*I) else random(FF);

beta = 1-s + s*gamma;

dfR = sub(dfr,R) + I*sub(dfi,R);
dgR = sub(dgr,R) + I*sub(dgi,R);

f = ideal (for i from 0 to numgenX-1 list (gens(sub(fr,R)))_(0,i) + I*(gens(sub(fi,R)))_(0,i));
g = ideal (for i from 0 to numgenY-1 list (gens(sub(gr,R)))_(0,i) + I*(gens(sub(gi,R)))_(0,i));


varsRX = sub(transpose(vars RX), R);
varsRY = sub(transpose(vars RY), R);

lamblist = toList(lambda_1..lambda_numgenX);
mulist = toList(mu_1..mu_numgenY);
lambdas = transpose(matrix({lamblist}));
mus = transpose(matrix({mulist}));

-- general data point
qpointx = if (FF === QQ) then (sub(random(FF^nvars, FF^1), R)+sub(random(FF^nvars, FF^1), R)*I) else sub(random(FF^nvars, FF^1), R);
qpointy = qpointx;
--qpointy = if (FF === QQ) then (sub(random(FF^nvars, FF^1), R)+sub(random(FF^nvars, FF^1), R)*I) else sub(random(FF^nvars, FF^1), R);

MX = nu_0*(beta*(varsRX-qpointx)-(1-s)*(varsRY-qpointx))+beta*dfR*lambdas;
MY = nu_1*(beta*(varsRY-qpointy)-(1-s)*(varsRX-qpointy))+beta*dgR*mus;

eqnsmat = MX || MY || sub(transpose(gens f), R) || sub(transpose(gens g), R);
eqns = for i from 0 to 2*nvars+numgenX+numgenY-1 list eqnsmat_(i, 0);
eqnssub = for i from 0 to #eqns-1 list sub(sub(eqns#i,R),s=>1);

indexf = toList(0..nvarsm1)|toList(2*nvars..(2*nvars-1+numgenX));
indexg = toList(nvars..(2*nvars-1))|toList(2*nvars+numgenX..(2*nvars+numgenX+numgenY-1));

eqnsf = eqnssub_indexf;
eqnsg = eqnssub_indexg;

RXext = FF[x_0..x_nvarsm1,nu_0,lambda_1..lambda_numgenX,I];

eqnsfRXext = for i from 0 to #eqnsf-1 list sub(eqnsf#i,RXext);

dirX = currentDirectory()|"bertini_outX";
if not isDirectory(dirX) then mkdir dirX;

solnsf = bertiniZeroDimSolve(eqnsfRXext, TopDirectory => dirX, AffVariableGroup=>{toList(x_0..x_nvarsm1)}, HomVariableGroup=>{{nu_0}|toList(lambda_1..lambda_numgenX)});

run ("cat "|dirX|"/bertini_session.log");

-- filtering out "singular" solutions based on Bertini default CONDNUMTHRESHOLD = 1e8;
solnsf = for i from 0 to #solnsf-1 list (if (solnsf#i)#ConditionNumber>1e8 then continue; solnsf#i);

-- define f0c as complex CC-polynomial?

RXc = CC[x_0..x_nvarsm1];

f0c = ideal(sub(gens(f0r), RXc)+matrix({for i from 0 to numColumns(gens f0i)-1 list ii*sub((gens(f0i))_(0,i), RXc)}));

if ((nvars-dimX) < numgenX0) then (solnsf = filter({nvars,numgenX}, solnsf, f0c, eps));

RYext = FF[y_0..y_nvarsm1,nu_1,mu_1..mu_numgenY, I];

eqnsgRYext = for i from 0 to #eqnsg-1 list sub(eqnsg#i,RYext);

dirY = currentDirectory()|"bertini_outY";
if not isDirectory(dirY) then mkdir dirY;

solnsg = bertiniZeroDimSolve(eqnsgRYext, TopDirectory => dirY, AffVariableGroup=>{toList(y_0..y_nvarsm1)}, HomVariableGroup=>{{nu_1}|toList(mu_1..mu_numgenY)});

run ("cat "|dirY|"/bertini_session.log");

-- filtering out "singular" solutions based on Bertini default CONDNUMTHRESHOLD = 1e8;
solnsg = for i from 0 to #solnsg-1 list (if (solnsg#i)#ConditionNumber>1e8 then continue; solnsg#i);

RYc = CC[y_0..y_nvarsm1];

g0c = ideal(sub(gens(g0r), RYc)+matrix({for i from 0 to numColumns(gens g0i)-1 list ii*sub((gens(g0i))_(0,i), RYc)}));

if ((nvars-dimY) < numgenY0) then (solnsg = filter({nvars, numgenY}, solnsg, g0c, eps));

startprel = toList(set(solnsf)**set(solnsg));

-- when X=Y
--Lpairs = subsets(set(solnsf),2);
--#Lpairs
--startprel = for i from 0 to #Lpairs-1 list toSequence(toList(Lpairs#i));

start = for i from 0 to #startprel-1 list (coordinates((startprel#i)#0)|coordinates((startprel#i)#1));

-- not the best solution. Cannot get variable groups in order I want in main homotopy
start = for i from 0 to #start-1 list (start#i)_(toList(0..numgenX))|(start#i)_(toList((numgenX+1+nvars)..(numgenX+1+nvars+numgenY)))|(start#i)_(toList((numgenX+1)..(numgenX+nvars)))|(start#i)_(toList((numgenX+numgenY+2+nvars)..(numgenX+numgenY+1+2*nvars)));

dirmain = currentDirectory()|"bertini_out_main";
if not isDirectory(dirmain) then mkdir dirmain;

-- S0 = bertiniTrackHomotopy(s, eqns, start);

--S0 = bertiniUserHomotopy(t, {s=>t}, eqns, start, TopDirectory => dirmain, AffVariableGroup=>{toList(x_0..x_nvarsm1|lambda_1..lambda_numgenX|y_0..y_nvarsm1|mu_1..mu_numgenY)});
S0 = bertiniUserHomotopy(t, {s=>t}, eqns, start, TopDirectory => dirmain, AffVariableGroup=>{toList(x_0..x_nvarsm1|y_0..y_nvarsm1)}, HomVariableGroup=>{{nu_0}|toList(lambda_1..lambda_numgenX),{nu_1}|toList(mu_1..mu_numgenY)});

-- do something like this instead:
-- S0reg = for i from 0 to #S0-1 when (status S0_i==Regular) list S0_i;
count = 0;

for i from 0 to #S0-1 do {
--print peek(S0_i);
-- Based on CONDNUMTHRESHOLD default value 1e8
if (((S0_i)#ConditionNumber)<1e8) then count=count+1;
};

run ("cat "|dirmain|"/bertini_session.log");

count
#start