# This file contains the implementation of our algorithms to solve
# NLP and Mp-NLP problems. This file contains also an
#implementation from the following web-site to compute
#comprehensive Groebner systems:
#https://amirhashemi.iut.ac.ir/softwares

with(LinearAlgebra):
with(Groebner):
with(Algebraic):
with(PolynomialIdeals):

############################### KKT conditions
KKT := proc (J, G, H, U, X) 
local N; 
N := [seq(diff(J, U[i])+add(mu[j]*(diff(G[j], U[i])), j = 1 .. nops(G))+add(lambda[j]*(diff(H[j], U[i])), j = 1 .. nops(H)), i = 1 .. nops(U)), seq(mu[j]*G[j], j = 1 .. nops(G))]; 
RETURN([op(N), op(H)]); 
end proc:

###### NLP solve
NLP := proc (J, G, H, U) 
local F, T, Grob, ns, rv, M, E, MU, flag, e, S, Sat, u, g, Elist; 
F := KKT(J, G, H, U, []); 
T := tdeg(seq(mu[i], i = 1 .. nops(G)), op(U)); 
Grob := Basis(F, T); 
if not IsZeroDimensional(`<,>`(op(Grob))) then 
    error "The KKT ideal is not zero-dimensional."; 
end if; 
ns, rv := NormalSet(Grob, T); 
M := MultiplicationMatrix(J, ns, rv, Grob, T); 
E := Eigenvalues(M); 
Elist := sort(convert(E, list), proc (p, q) options operator, arrow; evalb(p < q) end proc); 
MU := [seq(mu[i], i = 1 .. nops(G))]; 
flag := false; 
for e in Elist while not flag do 
    S := solve([op(Grob), J-e]); 
    Sat := true; 
    for u in MU while Sat do 
        if eval(u, S) < 0 then 
            Sat := false; 
        end if; 
    end do; 
    if Sat then 
        for g in G while Sat do 
            if 0 < eval(g, S) then 
                Sat := false; 
            end if; 
        end do; 
    end if; 
    if Sat then 
        RETURN(seq(u = eval(u, S), u = U), Min(J) = eval(J, S)); 
    end if; 
end do;
end proc:
#################################### Example
print("Example for NLP");
printf("%-1s %1s %1s%1s    : %3a \n",The, objective, function, is, 2*x^2-5*y):
 printf("%-1s %1s %1s   : %3a \n",The,constraints, are, {100<=x, x<=200, 80<=y, y<=170, x+y>=200}):
NLP(2*x^2-5*y, [100-x, x-200, 80-y, y-170, 200-x-y], [], [x, y]);





############################Select Algorithm###################### 
sff := proc (x) 
evalb(type(x, 'constant')): 
end proc:
################
lt := proc (f, T) 
RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f, T)[2]): 
end proc: 
#############Normalize Algorithm################
nrm := proc (F) 
local KK, h, i, NM; 
#option trace; 
if F = [] then 
RETURN(F); 
end if; 
KK := F; 
if type(denom(KK[1]), monomial) then 
h := LeadingMonomial(denom(KK[1]), plex(op(indets(denom(KK[1]))))); 
else 
h := denom(KK[1]); 
end if; 
for i from 2 to nops(KK) do 
if type(denom(KK[i]), monomial) then 
h := lcm(h, LeadingMonomial(denom(KK[i]), plex(op(indets(denom(KK[i])))))); 
else 
h := lcm(h, denom(KK[i])); 
end if; 
end do; 
NM := simplify(expand(h*KK)); 
RETURN(NM); 
end proc:
###############convert a list of polynomial to the set "FACVAR"##########
fac := proc (L, T) 
local A, AA, AAA, N, P, p, B, C, i; 
#option trace;
A := L; 
AA := [seq(factor(i), `in`(i, A))]; 
AAA := NULL; 
B := []; 
for i to nops(AA) do 
    p := AA[i]; 
    if irreduc(p) = true then 
        AAA := AAA, p: 
    else 
        AAA := AAA, op(convert(p, list)): 
    end if 
end do; 
P := NULL; 
for i to nops([AAA]) do 
    if `not`(type(AAA[i], 'constant')) then 
        P := P, [AAA][i]: 
    end if 
end do; 
N := NULL; 
for i to nops([P]) do 
    if irreduc([P][i]) = false and type([P][i], 'constant') <> true then 
        N := N, op([P][i])[1]: 
    else 
        N := N, [P][i]: 
    end if 
end do; 
RETURN({N}): 
end proc:

################## New Condition #################
new := proc (N, W, f, R, T) 
local ff, test, NN, WW, cd, KK, ww, i, cc, ccc, k, M, www,fff; 
#option trace:
ff := simplify(expand(f)); 
test := true; 
NN := N; 
KK := []; 
while test and NN <> [1] and ff <> 0 do 
    while type(LeadingCoefficient(ff, T), 'constant') = true and ff <> 0 do 
        KK := [op(KK), lt(ff, T)]; 
        ff := simplify(expand(ff-simplify(expand(lt(ff, T))))) 
    end do; 
    if RadicalMembership(LeadingCoefficient(simplify(ff), T), simplify(`<,>`(NN))) = true then 
        NN := [op({op(NN), LeadingCoefficient(ff, T)})]; 
        ff := simplify(ff-simplify(expand(lt(ff, T)))): 
    else 
        test := false:
    end if :
end do; 
NN := Basis(NN, R); 
WW := [op({seq(NormalForm(W[i], NN, R), i = 1 .. nops(W))})]; 
ff := op(nrm([NormalForm(ff, NN, R)])); 
if WW <> [] then 
    #`WW&Assign;nrm`(WW);
     WW:=nrm(WW):
    ww := NULL; 
    for i to nops(WW) do 
        if not type(WW[i], 'constant') = true then 
            ww := ww, WW[i] 
        end if 
    end do; 
    WW := [ww] 
end if; 
fff := ff; 
while type(LeadingCoefficient(fff, T), 'constant') = true and fff <> 0 do 
    fff := simplify(expand(fff-simplify(expand(lt(fff, T))))) 
end do; 
cc := fac([LeadingCoefficient(fff, T)]); 
ccc := {seq(op(nrm([k/LeadingCoefficient(k, R)])), `in`(k, cc))};
cd := `minus`({seq(expand(kk), `in`(kk, ccc))}, {-op(WW), op(WW)}); 
k := add(KK[j], j = 1 .. nops(KK)); 
M := [op(fac(WW, R))]; 
WW := [op({seq(ii/LeadingCoefficient(ii, R), `in`(ii, M))})]; 
www := [seq(nrm([ee]), `in`(ee, WW))]; 
if www <> [] then 
    www := [op(www)]; 
www := [seq(op(jj), `in`(jj, www))] 
end if; 
if nops(cd)<>0 then
     RETURN([{op([nrm([op(cd minus select(sff,cd))][1])])}, ff+k, NN, www]); 
else
     RETURN([{}, ff+k, NN, www]);
fi:
end proc:
##################
sf := proc (x) 
evalb(x <> 0): 
end proc: 

##############################
##       G2V Algorithm       ##
##############################
LM:=proc(f)
global tord;
if f<>0 then
    RETURN(LeadingMonomial(op(nrm([f])),tord));
fi:
RETURN(0);
end:
LC:=proc(f)
global tord;
RETURN(LeadingCoefficient(op(nrm([f])),tord));
end:
#############################
refine:=proc(LL,Jsig)
local u,f:
u:=Jsig:
f:=proc(Q)
global u:
if not divide(LM(Q[1]),u) then
    RETURN(true);
fi:
RETURN(false);
end:
RETURN(select(f,LL)):
end:
#############################
JPair2:=proc(P,Q,LMH)
#option trace;
global tord,Polys,Deg,NumOfFirstBuch,m,DDD;
local t,t1,t2,f,L,T:
####
f:=proc(r)
if r<>0 then 
    RETURN(true);
fi:
RETURN(false);
end:
####
T:=tord;
L:=select(f,LMH):
t:=lcm(P[3],Q[3]):
DDD:=max(DDD,degree(t));
t1,t2:=simplify(t/P[3]) , simplify(t/Q[3]):
if TestOrder(t1*LM(P[1]) , t2*LM(Q[1]),T) and NormalForm(t2*LM(Q[1]),L,T)<>0  then
    RETURN(expand([LeadingCoefficient(Polys[P[2]],T)*t2*Q[1],Q[2],t2*Q[3],t2*Q[4],P[2]]));
elif TestOrder(t2*LM(Q[1]) , t1*LM(P[1]),T) and NormalForm(t1*LM(P[1]),L,T)<>0 then 
    RETURN(expand([LeadingCoefficient(Polys[Q[2]],T)*t1*P[1],P[2],t1*P[3],t1*P[4],Q[2]]));
fi:
RETURN(NULL);
end:
#############################
tidy:=proc(P,Q)
global tord;
local T:
T:=tord;
if TestOrder(LM(P[1]),LM(Q[1]),T) then
    RETURN(true);
fi:
RETURN(false);
end:
#############################
RegRed:=proc(P)
#option trace:
local u1,v1,i,flag,Q,u2,v2,t,c,cmpt,T:
global tord,r,Grob,Polys,NumOfSup,R,GG,tp:
T:=tord;
u1,v1:=P[1],expand(P[4]*Polys[P[2]]):
i:=1:
cmpt:=0;
while i<=nops(R) and v1<>0 do
    flag:=false:
    Q:=R[i]:
    u2,v2:=Q[1],expand(Q[4]*Polys[Q[2]]):
    t:=LM(v1)/LM(v2):
    if divide(LM(v1),LM(v2)) and TestOrder(LM(expand(t*u2)) , LM(u1),T) and (P[2]<>Q[2] or cmpt=1) then# and Q[1]<>t*P[1] then
        cmpt:=1;
        c:=LeadingCoefficient(op(nrm([v1])),T)/LeadingCoefficient(op(nrm([v2])),T):
        if LM(u1)=t*LM(u2) then 
	       RETURN(1); 
	    fi;
        u1:=op(nrm([simplify(u1-c*t*u2)])): 
        v1:=SPolynomial(v1, v2, T);
        if LM(expand(t*u2)) = LM(u1) and c<>1 then
            u1:=op(nrm([u1/(1-c)])):
            v1:=op(nrm([v1/(1-c)])):
        fi:
        u1:=op(nrm([NormalForm(u1,Grob,T)])):
        v1:=op(nrm([NormalForm(v1,Grob,T)])):
        i:=1:
        flag:=true:
    else
        i:=i+1:
    fi:
od:
RETURN([op(nrm([u1])),op(nrm([v1]))]);
end:
############################
g2v:=proc(fff,GGG,N,W,RR,T,CPP,RRR)
#option trace:
local n,U,V,LMH,v,u,JP,NF,UV,Pm,g,flag,nr,m,P,k,Y,c,S,test,CP,f,fun,i,funi,G;
global tp,Grob,L,Polys,R,NumOfZero,NumOfF5,H,tord,null,notnull,fflag,iii,DDD,GG,rdz:
f:=fff[2];
CP:=CPP:
GG:=NULL:
tp:=RR;
R:=NULL:
fun:=proc (p, i) if i in {p[2],p[-1]} then RETURN(true) else RETURN(false) end if end proc;
for i from 1 to nops(GGG) do
    g:=NormalForm(GGG[i],N,RR):
    if g=0 then 
    CP:=remove(fun,CP,i);
    funi:=proc(p,j) if p<j then RETURN(p); else RETURN(p-1); fi;end:  
	CP:={seq([p[1],funi(p[2],i),p[3],p[4],funi(p[5],i)], p in CP)}:
    else
	   GG:=GG,g;
	   R:=R,[RRR[i][1],RRR[i][2],LM(g),RRR[i][4]];
    fi;
od;
G:=[GG];if N=[1] then RETURN(true,G,N,W,0,{},{}); ;fi;
R:=[R];
tord:=T;
null, notnull:=N,W:
Polys:=G:
m:=nops(Polys);
LMH:=[seq(LM(Grob[j]), j=1..nops(Grob))]:
g:=expand(NormalForm(f,N,RR)):
v:=op(nrm([NormalForm(g,Grob,T)])):
if v=0 and CP={} then 
    RETURN(true,G,N,W,0,{},{}):
elif v=0 then
    JP:=CP:
else
    
    u:=fff[1];
    Y := new(null, notnull, op(nrm([v])), RR, T);
                    c[dec] := Y[1];
                    S:=Y[2]: 
                    null := Y[3]; 
                    notnull := Y[4]; 
                    if S=0 then RETURN(true,G,N,W,0,{},{}): fi:
                    if c[dec] = {} then 
                        if S <> 0 then 
			                 v:=simplify(op(nrm([S])));
			                 nr:=nops(R):
			                 R:=[op(R),[u,nr+1,LM(v),1]]:
			                 Polys:=[op(Polys),v]:
			                 JP:=[op({op(CP),seq(JPair2(R[-1],R[j],LMH),j=1..nr)})]:
			                 JP:=sort(JP,tidy);
			             end if: 
                    elif S <> 0 then 
                        RETURN(false,Polys,null,notnull,[u,S],CP,R);
                    end if: 
fi:
while nops(JP)<>0 do #print(nops(JP));
        P:=JP[1]:
        JP:=JP[2..-1]:
        NF:=RegRed(P):
        if NF<>1 then
            u:=NF[1]:
            v:=NF[2]:
            if v=0 then
                #H:=[op(H),u]:
                LMH:=[op(LMH),LM(u)]:
                JP:=refine(JP,LM(u)):
                rdz:=rdz+1;
            else
                    v:=NormalForm(v,N,RR);
                    ##print(null, notnull, op(nrm([v])), RR, T);
                    Y := new(null, notnull, op(nrm([v])), RR, T);
                    c[dec] := Y[1];
                    S:=Y[2]: 
                    null := Y[3]; 
                    notnull := Y[4]; 
                    if c[dec] = {} then 
                        if S <> 0 then 
			                 v:=simplify(op(nrm([S])));
			                 nr:=nops(R):
			                 R:=[op(R),[u,nr+1,LM(v),1]]:
			                 Polys:=[op(Polys),v]:
			                 JP:=[op({op(JP),seq(JPair2(R[-1],R[j],LMH),j=1..nr)})]:
			                 JP:=sort(JP,tidy);
			             end if: 
                    elif S <> 0 then 
                        RETURN(false,Polys,null,notnull,[u,S],JP,R);
                    end if: 
            fi:
        fi:
od:
RETURN(true,Polys,null,notnull,0,{},{});
end:
#######################################
#######################################
canspec := proc (LL, M, R) 
local N, W, WW, WWW, NN, NNN, test, i,h, t, facN, facW, x, NNNN, L, w, flag, ww, www, MM,NnN,n; 
#option trace;
N := Basis(LL,R); 
W := M; 
WW := seq(NormalForm(W[i], N, R), i = 1 .. nops(W));
WWW := seq(factor([WW][i]), i = 1 .. nops([WW])); 
test := true; 
t:= true; 
NN := N; 
h := product(W[i], i = 1 .. nops(W)); 
if RadicalMembership(h, `<,>`(NN)) = true then
    test := false;
    NN := {1};
    RETURN([test, NN,[WW]]):
end if; 
while t do
    t := false;
    facN := [seq([op(i)], `in`(i, [seq(fac([NN[i]]), i = 1 .. nops(NN))]))];
    facW := [op(fac([WWW], R))];
    NNN := [];
    L := NULL;
    for i from 1 to nops(facN) do
        for x in facN[i] do
           flag := true;
           for w in facW while flag do
             if divide(w, x) then
                flag := false;
                L := L, x:
             end if:
           end do:
           
        end do;
        NnN[i]:= [op(`minus`({op(facN[i])}, {L}))];
        n[i]:=simplify(expand(product(NnN[i][j], j = 1 .. nops(NnN[i]))));
     end do;   
     NNNN :=[seq(n[i],i=1..nops(facN))];
    if {op(NNNN)} <> {op(NN)} then 
        t := true;
        NN := Basis(NNNN, R);
        WW := seq(NormalForm([WWW][i], NN, R), i = 1 .. nops([WWW]));
        WWW := seq(factor([WW][i]), i = 1 .. nops([WW]))
    end if:
end do; 
ww := NULL; 
for i to nops([WWW]) do
    if type(WWW[i], 'constant') = true then
    else
        ww := ww, [WWW][i]:
    end if
end do; 
WWW := ww; 
MM := [op(fac([WWW], R))]; 
WWW :=[op({seq(i/LeadingTerm(i, R)[1], `in`(i, MM))})]; 
www :=[seq(nrm([ee]), `in`(ee, WWW))]; 
if www <> [] then
    www := [op(www)];
    www := [seq(op(jj), `in`(jj, www))]:
end if; 
RETURN([test, NN, www]): 
end proc:

##########################################
branch := proc (v, ffff,BBB, NNN, WWW, R, T,CPP,RR) 
local cd,l,ll,i,ff,Y,X,BB,pivot,test,TTT,Bb,g,TT,vv,LL,gg,B,CN,N,W,NN,WW,f,CP,RRRR;
global LIST, Grob,termord1, termord2, termord3, iii, DDD,ind,rdz; 
#option trace:
f:=ffff[2];
CN := canspec(NNN, WWW, R); 
B:=BBB:
N:=CN[2];
W:=CN[3];
cd := []; 
Bb:=B: 
Y := new(N, W, f, R, T);
cd := [op(Y[1])];
ff := simplify(expand(Y[2])); 
NN := Y[3]; 
WW := Y[4];  
if v = [] then 
    vv := NULL: 
else 
    vv := op(v): 
end if;  
TT[vv] := [v, Bb, NN, WW]; 
if cd = [] then 
    X:= g2v([ffff[1],ff],Bb, NN, WW, R, T,CPP,RR); 
    test := X[1]; 
    BB := X[2]; 
    NN := X[3]; 
    WW := X[4]; 
    ff:=X[5];
    CP:=X[6];
    RRRR:=X[7]:
    if test  then 
        #TT[vv] := [[vv], nrm(InterReduce(BB,T)), NN, WW]; 
        TT[vv] := [[vv], nrm(BB), NN, WW]; 
	if CP={} then
	    LIST := LIST, TT[vv]:
        
	fi;
    else 
        branch([vv],ff, BB, NN, WW, R, T,CP,RRRR); 
    end if: 
else 
    branch([0,vv], [ffff[1],ff],Bb, [op(N),op(cd)], W, R, T,CPP,RR); 
    branch([1,vv], [ffff[1],ff],Bb, N, [op(W),op(cd)], R, T,CPP,RR);
end if; 
end proc:
################## Incremental Montes Algorithm ##############
IncMontes:=proc(f,Grobsys,R,T)
local L,i,G,N,W,ff,CPP;
global LIST, DDD,rdz,Grob;
#option trace;
L:=NULL:
for i from 1 to nops(Grobsys) do 
    CPP:={};
    LIST:=NULL:
    G:=Grobsys[i];
    Grob:=G[2]:
    N:=G[3]:
    W:=G[4]:
    ff:=NormalForm(f,N,R);
    if ff=0 then
        LIST:=LIST,[G[1],Grob,G[3],G[4]]: 
    else 
	 branch([], [1,ff],Grob, N, W, R, T, {},[seq([0,i,LM(Grob[i]),1],i=1..nops(Grob))]);
    fi:
    L:=L,LIST;
od:
LIST:=L;
RETURN(LIST);
end:

#######################
selpoly1 := proc (f, g) 
local ZPf, ZPg;
global termord1, termord2; 
ZPf := LeadingCoefficient(f, termord2); 
ZPg := LeadingCoefficient(g, termord2); 
RETURN(TestOrder(ZPg, ZPf, termord1)): 
end proc:
#######################
################## F5Montes Algorithm ##############
Montes:=proc(FF,R,T)
local f,t1,t2,b1,b2,F;
global LIST,DDD,ind,rdz,Grob,termord1, termord2;
local IndPolys,POLYS;
#option trace;
t1,b1:=kernelopts(cputime,bytesused);
termord1, termord2:=R,T:
IndPolys:=[seq(i,i=1..nops(FF))];
F:=InterReduce(FF,prod(T,R));
POLYS:=FF;
rdz:=0;
DDD:=0;
LIST:=NULL;
f:=F[1];
Grob:=[];
ind:=nops(F)-1;
branch([], [1,f],[], [], [], R, T,{},[]);
for f in F[2..-1] do
    ind:=ind-1;#print(ind);
    IncMontes(f,[LIST],R,T); 
od;
t2,b2:=kernelopts(cputime,bytesused);
#printf("%-1s %1s %1s%1s    : %3a %3a\n",The, cpu, time, is,t2-t1,(sec)):
# printf("%-1s %1s %1s   : %3a %3a\n",The,used,memory,b2-b1,(bytes)):
# printf("%-1s %1s %1s        : %3g\n",Num,of,Rds,rdz):
# printf("%-1s %1s %1s        : %3a\n",Num,of,sys,nops([LIST])):
# printf("%-1s %1s %1s %1s    : %3a\n",The, max, deg, is,DDD):
RETURN(LIST);
end:


################################
MPNLP := proc (J, G, H, U, X) 
local F, T, n, N, i, ns, rv, M, E, GS; 
F := KKT(J, G, H, U, X); 
T := tdeg(seq(mu[i], i = 1 .. nops(G)), seq(lambda[i], i = 1 .. nops(H)), op(U)); 
GS := [Montes(F, tdeg(op(X)), T)]; 
n := nops(GS); 
N := NULL; 
for i to n do 
    if IsZeroDimensional(`<,>`(op(LeadingMonomial(GS[i][2], T)))) then 
        ns, rv := NormalSet(GS[i][2], T); 
        M := MultiplicationMatrix(J, ns, rv, GS[i][2], T); 
        E := Eigenvalues(M); 
        N := N, [convert(E, set), op(GS[i][3 .. 4])]: 
    end if: 
end do; 
RETURN(N): 
end proc:

############ Example MPNLP(Objective, non-equality constraints, equality constraints, variables, parameters)
print();
print("Example for Mp-NLP");
printf("%-1s %1s %1s%1s    : %3a \n",The, objective, function, is, 2*x[1]^2+3*x[2]):
 printf("%-1s %1s %1s   : %3a \n",The,constraints, are, {0<=x[1], 0<=x[2],x[1]+x[2]-(130*(1/100))*p[1]<=0, p[1]-x[1]-x[2]<=0, 230*p[1]-480*x[1]-85*x[2]<=0, 6*p[1]-5*x[1]-3*x[2]<=0 }):
MPNLP(2*x[1]^2+3*x[2], [-x[1], -x[2], x[1]+x[2]-(130*(1/100))*p[1], p[1]-x[1]-x[2], 230*p[1]-480*x[1]-85*x[2], 6*p[1]-5*x[1]-3*x[2]], [], [x[1], x[2]], [p[1]]);