###############################################################################
#
#  The procedure varder calculates the variational derivative of
#  a functional L with respect to u.
#
###############################################################################

varder:=proc(L,u,x)
option `Authors: B. Deconinck`;
local N,k,lag,v,diffvec,m;
N:=hdifforder(L,u,x);
# N is the order of the highest derivative of u w.r.t. x occuring in L
lag:=L;
for k from N by -1 to 1 do
 lag:=subs(diff(u,[seq(x,j=1..k)])=v[k],lag)
od;
lag:=subs(u=v[0],lag);
diffvec:=[seq(diff(lag,v[k]),k=0..N)];
diffvec:=map2(subs,v[0]=u,diffvec);
diffvec:=map2(subs,[seq(v[k]=diff(u,[seq(x,m=1..k)]),k=1..N)],diffvec);
# The last result serves as the output of the procedure
(expand(diffvec[1]+add((-1)^k*diff(diffvec[k+1],[seq(x,m=1..k)]),k=1..N)))
end:

###############################################################################
#
# The procedure hdifforder finds the order of the highest
# derivative of u w.r.t. x occuring in L.
#
###############################################################################

hdifforder:=proc(L,u,x)
option `Authors: B. Deconinck`;
local k,lijst,flag;
k:=0; flag:=1; lijst:=[];
while flag=1 do
 if has(L,diff(u,lijst)) then
  lijst:=[op(lijst),x];
  k:=k+1
 else flag:=0
 fi
od;
# the last result serves as the output of the procedure
if k=0 then 0
else k-1
fi;
end:

###############################################################################
#
# The procedure varder calculates the higher-order
# variational derivative of a functional L with respect to u.
# This procedure cannot be used to compute the variational
# derivative. We require n to be at least 1.
#
###############################################################################

vardern:=proc(L,u,n,x)
option `Authors: B. Deconinck`;
local N,k,lag,v,diffvec,m;
N:=hdifforder(L,u,x);
if N<n then return(0)
fi;
# N is the order of the highest derivative of u w.r.t. x occuring in L
lag:=L;
for k from N by -1 to n do
 lag:=subs(diff(u,[seq(x,j=1..k)])=v[k],lag)
od;
diffvec:=[seq(diff(lag,v[k]),k=n..N)];
diffvec:=map2(subs,[seq(v[k]=diff(u,[seq(x,m=1..k)]),k=n..N)],diffvec);
# The last result serves as the output of the procedure
expand(diffvec[1]+add((-1)^(k-n)*binomial(k,n)*diff(diffvec[k-n+1],[seq(x,m=n..k-1)]),k=n+1..N))
end:

################################################################################################
#
# This is a modified homotopy integral, a list of terms is returned
# Here u is a vector of unknown functions
# that depend on an independent variable x.
#
################################################################################################

vechomotopyintegral:=proc(LL,u,x)
option `Authors: B. Deconinck, M. Nivala`;
local sh,result,newans,newflux,aa,limset,ans,i,r,m,lambda,lambda0,flux,scaledflux,summedflux,N,n,L,h,a,s;
L:=expand(eval(LL));
n:=nops(u); # number of unknown functions
# N is the vector of differential orders of u in L
N:=[seq(hdifforder(L,u[r],x),r=1..n)];
# flux is the integrand before scaling of the arguments
# to obtain the maximum number of terms, we take this as a list and a sum
flux:=expand([seq(op([u[r]*vardern(L,u[r],1,x),seq(diff(u[r]*vardern(L,u[r],i+1,x),[seq(x,m=1..i)]),i=1..N[r]-1)]),r=1..n)]);
scaledflux:=expand(subs({seq(u[r]=lambda*u[r],r=1..n)},flux));
sh:=normal(convert(scaledflux,`+`));
ans:=expand({seq(int(expand(scaledflux[i]/lambda),lambda),i=1..nops(scaledflux)),int(expand(sh/lambda),lambda)});
#remove expressions that were not integrated
for s in ans do
 if has(s,int) then
  ans:=ans minus {s}
 fi;
od;
#get as many terms as possible
newans:=[NULL];
for i from 1 to nops(ans) do
 if whattype(ans[i])=`+` then
  newans:=[op(newans),op(ans[i])];
 else
  newans:=[op(newans),ans[i]];
 fi;
od;
#normalize the solution and remove constants
result:={NULL};
for i from 1 to nops(newans) do
 if has(newans[i],x) then
  result:={op(result),SPLIT(newans[i],x)[1]};
 fi;
od;
[op(subs(lambda=1,result))];
end:

###############################################################################
#
# whatfun picks out undefined functions in an expression
#
###############################################################################

whatfun:=proc(E)
option `Authors: M. Nivala`;
local n,s,sid;
s:={NULL};
sid:=[op(indets(E,`function`))];
for n from 1 to nops(sid) do
 if not(has(sid[n],`diff`)) and not(type(op(0,sid[n]),`procedure`))
  then s:={op(s),sid[n]};
 fi;
od;
[op(s)];
end:

##############################################################################
#
# orderset orders a set from large to small
#
##############################################################################

orderset:=proc(s)
option `Authors: B. Deconinck, M. Nivala`;
    local ss,leftover,k;
    ss:=NULL;
    leftover:=s;
    for k to nops(s) do
        ss:=ss,max(op(leftover));
        leftover:=leftover minus {max(op(leftover))}
    od;
    [ss]
end:

#################################################################################
#
# funlexorderconverts a list of functional expressions to strings, orders them
# according to lexorder, and returns the funtional expressions
#
#################################################################################

funlexorder:=proc(LIST)
option `Authors: M. Nivala`;
local Newlist,i,subrules,s;
subrules:=[seq(convert(LIST[i],`string`)=LIST[i],i=1..nops(LIST))];
Newlist:=[seq(op(1,subrules[i]),i=1..nops(LIST))];
s:=sort(Newlist,`lexorder`);
subs(subrules,s);
end:

#######################################################################################
#
# orders a set according to order of derivatives, from highest to lowest,
# grouping those of the same order
#
#######################################################################################

dorderlist:=proc(L,u,x)
option `Authors: M. Nivala`;
local finlist,j,diffset,dlist,i,MAX,difflist,ord;
difflist:=[seq(hdifforder(L[i],u,x),i=1..nops(L))];
diffset:=orderset({op(difflist)});
finlist:=[NULL];
for i in diffset do
 dlist:=[NULL];
 for j from 1 to nops(L) do
  if i=difflist[j] then
   dlist:=[op(dlist),L[j]];
  fi;
 od;
finlist:=[op(finlist),dlist];
od;
finlist;
end:

###################################################################################
#
# torderset orders a set according to number of products in a term
#
###################################################################################

torderset:=proc(E)
option `Authors: M. Nivala`;
local w,fin,n,s;
s:={seq(nops(op(1,E[n]))+(n/(nops(E)+1)),n=1..nops(E))};
s:=orderset(s);
w:={seq(nops(op(1,E[n]))+(n/(nops(E)+1))=E[n],n=1..nops(E))};
fin:=subs(w,s);
end:

###################################################################################
#
# SPLIT splits the product T into an x-dependent part and a constant part
#
###################################################################################

SPLIT:=proc(T,x)
option `Authors: M. Nivala`;
local xpart,coeffpart,s,opT,xset,coeffset;
xset:={NULL};
coeffset:={NULL};
if whattype(T)=`*` then
 opT:=[op(T)];
else
 opT:=[T];
fi;
for s in opT do
 if has(s,x) then
  xset:={op(xset),s};
 else
  coeffset:={op(coeffset),s};
 fi;
od;
xpart:=convert(xset,`*`);
coeffpart:=convert(coeffset,`*`);
return([xpart,coeffpart]);
end:

#################################################################################
#
# SUBS substitutes ssubrules into EE
#
#################################################################################

SUBS:=proc(ssubrules,EE,x)
option `Authors: M. Nivala`;
local E,Elist,subrules;
subrules:=torderset(ssubrules);
E:=expand(EE);
if not(whattype(E)=`+`) then
 Elist:=[E];
else
 Elist:=convert(E,`list`);
fi;
Elist:=[seq(SPLIT(Elist[i],x),i=1..nops(Elist))];
Elist:=subs(subrules,Elist);
E:=add(convert(Elist[i],`*`),i=1..nops(Elist));
end:

##################################################################################
#
# GROUP goups terms with the same x-dependence in E
#
##################################################################################

GROUP:=proc(E,x)
option `Authors: M. Nivala`;
local EL,subrules,EE,Esubs,newsubrules,i,b;
if not(whattype(E)=`+`) then
  EL:=[E];
else
 EL:=convert(E,`list`);
fi;
subrules:=[op({seq(SPLIT(EL[i],x)[1]=b[i],i=1..nops(EL))})];
subrules:=torderset(subrules);
Esubs:=SUBS(subrules,E,x);
Esubs:=collect(Esubs,{seq(b[i],i=1..nops(subrules))});
newsubrules:=[seq(op(2,subrules[i])=op(1,subrules[i]),i=1..nops(subrules))];
EE:=subs(newsubrules,Esubs);
end:

################################################################################################
#
# RRED performs gaussian elimination with a minimal number of row switching operations
# it row-reduces A while preserving the order of the rows
# returns the first maxsstep-1 linearly independent rows
#
################################################################################################

RRED:=proc(A,sstep,B,maxsstep)
option `Authors: M. Nivala`;
local step,i,j,newA,flag1,rowdims,newB;
rowdims:=linalg[rowdim](A);
if sstep=maxsstep then
 return(B);
fi;
step:=sstep;
newA:=A;
flag1:=0;
for i from 1 to rowdims do
  if newA[i,step]<>0 and flag1=0 then
   flag1:=i;
  fi;
od;
if flag1<>0 then
 newA:=linalg[pivot](newA,flag1,step);
 if step=1 then
  newB:=linalg[matrix]([linalg[row](newA,flag1)]);
 else
  newB:=linalg[stackmatrix](B,linalg[row](newA,flag1));
 fi;
 if rowdims=1 then
  return(newB);
 else
  newA:=linalg[delrows](newA,flag1..flag1);
  RRED(newA,step+1,newB,maxsstep);
 fi;
else
 RRED(newA,step+1,newB,maxsstep);
fi;
end:

#################################################################################################################################
#
# INT is used to compute integrals of expressions that depend on unknown functions that depend on an independent variable x.
# All terms that are total derivatives are integrated.
# The order of the highest derivative in the remaining integral is minimized.
#
#################################################################################################################################

INT:=proc(EE,x)
option `Authors: M. Nivala`;
local dcheck,limset,eqnset2,RE,param, eqnset,v,slnset,novan,Flist,intnogood,nogood,subDFL,intxset,sln,Esub,Fsub,DFLsub,termlist,FL,xset,E,u,goo,A,Ar,Ag,bob,syst,i,j,m,DFL,s,n,fresult,result,a,E1,DF1,R,subrules,b,EL,F1,F,DF,k;

E:=expand(eval(EE));

#determine the unknown functions
u:=whatfun(E);
if u=[NULL] then return(int(E,x))
fi;

#check if E is exact
#dcheck:=[seq(varder(E,u[i],x),i=1..nops(u))];
#if dcheck=[seq(0,i=1..nops(u))] then
 #return(expand(int(E,x)));
#fi;

#convert E to a list
if not(whattype(E)=`+`) then
 EL:=[E];
else
 EL:=convert(E,`list`);
fi;

#remove terms which have no functional or x dependence to be integrated directly
xset:={NULL};
for k in EL do
if not(has(k,u)) or not(has(k,x)) then
 xset:={op(xset),k};
fi;
od;
xset:=expand(convert(xset,`+`));
intxset:=int(xset,x);
E:=expand(eval(E-xset));
if not(whattype(E)=`+`) then
 EL:=[E];
else
 EL:=convert(E,`list`);
fi;

#FL is a list of terms from the homotopy integral of E
FL:=vechomotopyintegral(E,u,x); #print(FFFFFFF,FL);
if FL=[NULL] then
 return(int(E,x)+intxset);
fi;

#DFL is the derivative of FL with respect to x
DFL:=[seq(expand(diff(FL[n],x)),n=1..nops(FL))];

#form substitution rules
#subDFL is for substitution rules
subDFL:={NULL};
for k in DFL do
 if not(whattype(k)=`+`) then
   subDFL:={op(subDFL),k};
 else
  subDFL:={op(subDFL),op(k)};
 fi;
od;
termlist:=[op({op(EL),op(FL),op(subDFL)})];
#remove constant coefficients
termlist:=[op({seq(SPLIT(termlist[i],x)[1],i=1..nops(termlist))})];
#orders termlist according to the order of the highest derivative of each element, from highest to lowest, groouping those of the
#same order
termlist:=dorderlist(termlist,u,x);#print(dord,termlist);
#for consistency, implement an ordering of funtions and their derivatives according to lexorder for those terms which are of the
#same differential order
termlist:=[seq(op(funlexorder(termlist[i])),i=1..nops(termlist))];#print(orderedsub,termlist);

#perform the substitutions
subrules:=(seq(termlist[n]=b[n],n=1..nops(termlist)));
#order subrules according to number of products in each term
subrules:=torderset({subrules});#print(SUBRULES,subrules);
Esub:=SUBS(subrules,E,x);
DFLsub:=[seq(SUBS(subrules,DFL[i],x),i=1..nops(DFL))];

#R is what remains after integrating E
#R is optimized, removing higest order derivaites first
R:=expand(Esub-sum(a[m]*DFLsub[m],m=1..nops(DFL)));#print(E-sum(a[m]*DFL[m],m=1..nops(DFL)));
R:=collect(R,{seq(b[m],m=1..nops(termlist))});#print(RRR,R);
#form the linear system to be solved
syst:={NULL};
for m from 1 to nops(termlist) do
 if has(coeff(R,b[m]),[seq(a[n],n=1..nops(DFL))]) then
  syst:=[op(syst),coeff(R,b[m])=0];
 fi;
od;

#A is the overdetermined system for the coefficients, a
#Ar is the first nops(a) linearly independent rows of A
A:=linalg[genmatrix](syst,[seq(a[j],j=1..nops(DFLsub))],flag);#print(AAAA,A);
Ar:=RRED(A,1,A,nops(DFLsub)+1);#print(ArAr,Ar);

slnset:=linalg[backsub](Ar,`false`,`v`);#print(slnset);
sln:=seq(a[m]=slnset[m],m=1..nops(DFLsub));#print(solution,sln);
sln:=eval(sln,v=0);#print(solution,sln);
assign(sln);

#all terms which were not integrated using homotopy, attempt to integrate each independently without homotopy
nogood:=expand(E-add(a[i]*DFL[i],i=1..nops(DFL)));#print(nogood);
if not(whattype(nogood)=`+`) then
  nogood:=[nogood];
else
 nogood:=convert(nogood,`list`);
fi;
intnogood:=expand(add(int(nogood[i],x),i=1..nops(nogood)));#print(intnogood);
#print(expand(E-add(a[i]*DFL[i],i=1..nops(DFL))+xset+novan+RE),"Was Integrated By Maple");

fresult:=add(a[j]*FL[j],j=1..nops(FL))+combine(intnogood)+intxset;#print(fdfdfdfd,fresult);
unassign(a);

#check if correct
if simplify(expand(EE-diff(fresult,x)))<>0 then
 print("No integration done.");
 return(int(E,x)+intxset)
else
 return(fresult);
fi;
end:

########################################################################################################
#
# ouputs the expression G modulo total derivatives
#
########################################################################################################

DerivativeMod:=proc(G,x)
local IG,n,intpart,Glist;
IG:=INT(G,x);
if not(whattype(IG)=`+`) then
 Glist:=[IG];
else
 Glist:=convert(IG,`list`);
fi;
intpart:=0;
for n from 1 to nops(Glist) do
 if has(Glist[n],`int`) then intpart:=diff(Glist[n],x);
 fi;
od;
intpart;
end:

#########################################
#
# checks if INT is correct
#
#########################################

INTcheck:=proc(E,F,x)
normal(expand(E-diff(F,x)));
end:

##############################################
#
#
#
##############################################

chi:=proc(n,x)
option remember;
local result;
if n=0 then result:=-u(x);
 else if n=1 then result:=diff(u(x),x);
         else result:=-diff(chi(n-1,x),x)-add(chi(n-k-2,x)*chi(k,x),k=0..n-2);
      fi;
fi;
result;
end: