(* Last updated: Monday, August 20, 2007 at 16:30 at home at Christchurch  *)

(* Previously updated: Saturday, August 14, 2004 at 2:40am in Christchurch *)

(* HEREMAN changes made in Christchurch *)
(* look for 08/14/2004, 08/15/2004 or CHRISTCHURCH *)

(* Used cross testing between Holly's code DDEDensityFlux.m *)
(* and diffdens code *)
(* In particular, for Taha/Herbst KdV, mKdV, and combined KdV-mKdV cases *)

(* From ugoktas@lie  Mon Apr 27 10:03:46 1998 *)
(* Date: Mon, 27 Apr 1998 09:38:59 -0600 (MDT) *)
(* Subject: diffdens.m *)
(* Added split function, written by Hereman, May 3, 1998 *)
(* Weights are computed now before the rank of rho is entered, May 4th, 1998 *)

(*****************************************************************************)
(*                                                                           *)
(*          *** M A T H E M A T I C A   P R O G R A M ***                    *)
(*                                                                           *)
(*      SYMBOLIC COMPUTATION of CONSERVED DENSITIES for SYSTEMS of           *)
(*         EVOLUTIONARY TYPE DIFFERENTIAL-DIFFERENCE EQUATIONS               *)
(*                                                                           *)
(* program name: diffdens.m                                                  *)
(*                                                                           *)
(* purpose: computation of the conserved densities (no fluxes)               *)
(* with possible compatibility conditions                                    *)
(*                                                                           *)
(* input to diffdens.m: system of evolution type differential difference     *)
(*                      equations of any order, any degree, polynomial       *)
(*                      type, only constant parameters, variable             *)
(*                      coefficients are NOT allowed as parameters           *)
(*                                                                           *)
(* output: density of desired rank (if it exists),                           *)
(*         and compatibility conditions for the parameters (if applicable)   *)
(*                                                                           *)
(* tested on: IBM RISC 6000, IBM Compatible PC 486, SGI Indigo2 XL           *)
(*                                                                           *)
(* language: Mathematica 3.0                                                 *)
(*                                                                           *)
(* authors: Unal Goktas and Willy Hereman                                    *)
(*          Department of Mathematical and Computer Sciences                 *)
(*          Colorado School of Mines                                         *)
(*          Golden, CO 80401-1887, USA                                       *)
(*                                                                           *)
(* Last updated: August 20, 2007                                             *)
(* Previous version: May 4, 1998                                             *)
(*                                                                           *)
(* Copyright 1998-2007                                                       *)
(*                                                                           *)
(*****************************************************************************)

(* 08/04/2007 WH print statement added *)
Print["Loading package diffdens2007.m. Last updated: August 20, 2007."];
Clear["Global`*"];

debugwillyCH = True;

(* ------------------------------------------------------------------------- *)

(*****************************************************************************)
(* commentinter[]: prints a welcome message                                  *)
(*****************************************************************************)

commentinter[] := Block[{},
Print["*********************************************************"];
Print["        WELCOME TO THE MATHEMATICA PROGRAM               "];
Print["         by UNAL GOKTAS and WILLY HEREMAN                "];
Print["    FOR THE COMPUTATION OF CONSERVED DENSITIES.          "];
Print["    (NO FLUXES) AND POSSIBLE COMPATIBILITY CONDITIONS.   "];
Print["            Last updated: August 20, 2007                "];
Print["            Previous version: May 4, 1998                "];
Print["                 Copyright 1998-2007                     "];
Print["*********************************************************"]
]; (* end commentinter *)

(*****************************************************************************)
(* cls: clears the screen                                                    *)
(*****************************************************************************)

cls := Print["\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"];

(*****************************************************************************)
(* printpage[]: a subroutine for menu                                       *)
(*****************************************************************************)

printpage[n_,page_] := Module[
{i,choice,lenpage},
(* *)
(* 08/14/2004 HEREMAN CHRISTCHURCH cls taken out *)
(* cls; *)
lenpage = Length[page[n]];
Print[" "];
Print["  *** MENU INTERFACE ***  (page: ",n,")"];
Print["-------------------------------------------"];
For[i=1,i <= lenpage,i++,
   Print[Part[Part[page[n],i],1]]];
Print[" nn) Next Page"];
Print[" tt) Your System"];
Print[" qq) Exit the Program"];
Print["------------------------------------------"];
choice = Input["ENTER YOUR CHOICE: "];
Return[choice]; ]; (* end printpage *)

(*****************************************************************************)
(* menu: creates the menu                                                    *)
(*****************************************************************************)

(* 08/14/2004 HEREMAN CHR added to the menu *)
(* 08/19/2007 WH reordered the menu *) 
menu := Module[
{counterpage = 1,menulist,numpages,page,
choice1,control,choice2,choice3,lenmenulist,i},
menulist = {
{"  1) Kac-van Moerbeke Equation (d_kdv.m)"},
{"  2) Modified KdV Equation (d_mkdv.m)"},
{"  3) Toda Lattice (d_toda.m)"},
{"  4) Parameterized Toda Lattice (d_ptoda.m)"},
{"  5) Generalized Toda Lattice-1 (d_gtoda1.m)"},
{"  6) Generalized Toda Lattice-2 (d_gtoda2.m)"}, 
{"  7) NLS Equation-Standard Discretization (d_stdnls.m)"},
{"  8) NLS Equation-Ablowitz Ladik Discretization (d_ablnls.m)"},
{"  9) Self-Dual Network Equations (d_dual.m)"},
{"  10) Discretization KdV equation (Taha/Herbst) (d_diskdv.m)"},
{"  11) Discretization mKdV equation (Taha/Herbst) (d_dimkdv.m)"},
{"  12) Discretization combined KdV-mKdV equation (Taha/Herbst) (d_comkdv.m)"}
}; (* closes menulist *)

lenmenulist = Length[menulist];
numpages = Ceiling[lenmenulist/10];
For[i = 1,i <= numpages,i++,
   page[i] = If[lenmenulist >= (i*10),
               menulist[[Table[k,{k,(i-1)*10+1,i*10}]]],
               menulist[[Table[k,{k,(i-1)*10+1,lenmenulist}]]]]];

choice1 = printpage[counterpage,page];

control := (
Switch[choice1,
tt,Print["Make sure that you have prepared the data file for the system"];
   Print["you want to test (similar to the data files we supplied)."];
   choice2 = Input["If your file is ready, press 1, else 2: "];
   If[choice2 === 1,
    choice3 = InputString["Enter the name of your data file: "];
    Get[choice3],
    Print["Aborting the computations!"];
    Print["Prepare your data file first, then start over."];
    Abort[]
     ],
nn,If[counterpage < numpages,counterpage++;
     choice1 = printpage[counterpage,page]; control,
     counterpage = 1; choice1 = printpage[1,page]; control],
qq,Print["Aborting the computations!"];Abort[],
(* *)
(* 08/14/2004, HEREMAN in CHRISTCHURCH added cases to menu *)
(* 08/19/2007 WH reordered the files *)
1,<<d_kdv.m,
2,<<d_mkdv.m,
3,<<d_toda.m,
4,<<d_ptoda.m,
5,<<d_gtoda1.m,
6,<<d_gtoda2.m,
7,<<d_stdnls.m,
8,<<d_ablnls.m,
9,<<d_dual.m,
10,<<d_diskdv.m,
11,<<d_dimkdv.m,
12,<<d_comkdv.m,
_,Print["Aborting the computations!"];Abort[]
]; (* closes Switch *)
); (* control *)

control;
If[Not[StringQ[myfile]],
   myfile = InputString["Enter the name of the output file: "]]; 
]; (* end menu *)

(*****************************************************************************)
(* subscriptform[]: takes an expression and prints it with subscripted form  *)
(*                  and removes the functional dependence                    *)
(*****************************************************************************)

subscriptform[expr_] := expr /. {
  u[i_][n_][x__] :>
          SequenceForm[u,Subscript[i],Subscript[","],Subscript[n]]
  }; (* end subscriptform *)

(*****************************************************************************)
(* split[]: if free coefficients c[i] remain in the density, the density     *)
(*          is further split in independent pieces corresponding to each     *)
(*          coefficient c[i]                                                  *)
(*****************************************************************************)

split[expr_] := Module[
  {listofc,lenlistofc,ci,rhoi,rest}, 
  rest=Expand[expr];
  If[FreeQ[rest,c[__]], 
    Print[" "];
    Print["The density has no free coefficients."], 
    listofc=Union[Cases[rest,c[__],12]];
    lenlistofc=Length[listofc];
    Print[" "];
    Print["There is/are ", lenlistofc, " free coefficient(s) in the density."];
    Print["These free coefficients are: ", listofc, "."];
    Print["Splitting the density in independent pieces."];
     Do[
       ci=Part[listofc,i];
       rhoi=Expand[Coefficient[rest,ci]]; 
       rest=Expand[rest-ci*rhoi];
       Print[" "];
       Print["Part of the density with coefficient ", ci, ":"];
       Print[" "];
       Print[rhoi],{i,1,lenlistofc}]; (* end do *)
     Print[" "];
     Print["Part of the density that had no free coefficient:"];
     Print[" "];
     Print[rest];
     Print[" "] ] (* end if *) ]; (* end split *)

(*****************************************************************************)
(* stripper[]: cancels the numerical factors of each term of the argument    *)
(*             and puts these into a list after cancellation                 *)
(*****************************************************************************)

stripper[givenexpr_,flag1_,flag2_] := Module[
{lenexpr1,list = {},expr1,expr2,expr3,expr3s,coefficient,expr4,
lenexpr3,iexpr,k,s},
expr1 = Expand[givenexpr];
If[Head[expr1] === Plus,lenexpr1 = Length[expr1],lenexpr1 = 1];
For[k = 1,k <= lenexpr1,k++,
   If[lenexpr1 == 1,expr2 = expr1,expr2 = Part[expr1,k]];
   coefficient = 1;
   expr3 = FactorList[expr2];
   lenexpr3 = Length[expr3];
   For[s = 1,s <= lenexpr3,s++,
      expr3s = Part[expr3,s];
      If[FreeQ[expr3s,u],
        If[flag1,
          If[FreeQ[weightpars,Part[expr3s,1]],
            If[flag2 == 3,
              coefficient = coefficient*Part[expr3s,1]^Part[expr3s,2]];
            expr2 = expr2/Part[expr3s,1]^Part[expr3s,2];
            ],
          If[flag2 == 3,
            coefficient = coefficient*Part[expr3s,1]^Part[expr3s,2]];
          expr2 = expr2/Part[expr3s,1]^Part[expr3s,2];
          ],
        expr4 = Part[expr3s,1];
        s = lenexpr3 + 1;
        ];
      ];
   Switch[flag2,
         1,list = Append[list,expr2],
         2,iexpr = Part[Head[expr4],1];
           expr2 = expr2 /. n->(2*n-iexpr);
           list = Union[list,{expr2}],
         3,iexpr = Part[Head[expr4],1];
           expr2 = expr2 /. n->(2*n-iexpr);
           list = Union[list,{{coefficient,expr2}}];
         ];
   ];
Return[list]; ]; (* end stripper *)

(*****************************************************************************)
(* setuniformrank[]: a subroutine for the scaling function                   *)
(*****************************************************************************)

setuniformrank[list_] := Module[
{i,lengthlist,syslist = {}},
lengthlist = Length[list];
For[i = 1,i < lengthlist,i++,
   syslist = Append[syslist,Part[list,i] == Part[list,lengthlist]];
   ];
Return[syslist]; ]; (* end setuniformrank *)

(*****************************************************************************)
(* scaling[]: determines the scaling properties, if the equations are        *)
(*            incompatible it will print the appropriate messages            *)
(*****************************************************************************)

scaling[eqlist_] := Module[
{pointsym,expr,list,syslist,msyslist = {},i,j,k,scalesol,
lenmsyslist,tempscalesol = {},tempmsyslist,trouble,troubleleft,
troubleright,posleft,posright,troublelist},

If[weightpars =!= {},
  pointsym = Union[Table[weightpars[[i]] -> E^weight[weightpars[[i]]],
                     {i,1,Length[weightpars]}],
                   {u[i_Integer][_][t] :> E^weightu[i], t -> E^(-1)}],
  pointsym = {u[i_Integer][_][t] :> E^weightu[i], t -> E^(-1)};
  ];
For[i = 1, i <= noeqs,i++,
   expr[i] = Part[eqlist,i];
   list[i] = stripper[expr[i],True,1];
   list[i] = list[i] /. pointsym;
   list[i] = PowerExpand[Log[list[i]]];
   list[i] = Prepend[list[i],weightu[i]+1];
   syslist[i] = setuniformrank[list[i]];
   msyslist = Union[msyslist,syslist[i]];
   ];
scalesol = Flatten[Solve[msyslist]];
If[MemberQ[msyslist,False],
  Print["Fatal Error! Weights are incorrect."];
  Print["Check your choice(s) for the weights in the data file."];
  Print["Aborting the computations!"];
  CloseLog[]; Abort[],
  If[scalesol === {} && msyslist =!= {True},
    Print["In the given system there is at least one equation with terms"];
    Print["of unequal rank. Scaling properties can not be determined for"];
    Print["this system. The program will try to find the conflict, and,"];
    Print["if successful, provide suggestions to help resolve the conflict."];
    lenmsyslist = Length[msyslist];
    i = 1;
    While[tempscalesol === {} && i <= lenmsyslist,
         tempmsyslist = Complement[msyslist,{Part[msyslist,i]}];
         tempscalesol = Flatten[Solve[tempmsyslist]];
         i++;
         ];
    i--;
    If[tempscalesol =!= {},
      trouble = Part[msyslist,i];
      troubleleft = Part[trouble,1];
      troubleright = Part[trouble,2];
      For[j = 1, j <= noeqs, j++,
         list[j] = Drop[list[j],1];
         posleft = Flatten[Position[list[j],troubleleft,{1}]];
         posright = Flatten[Position[list[j],troubleright,{1}]];
         If[posleft =!= {} && posright =!= {},
           troublelist = Union[
             Table[Part[Expand[u[j][n]'[t]],posright[[k]]],
                     {k,1,Length[posright]}],
             Table[Part[Expand[u[j][n]'[t]],posleft[[k]]],
                     {k,1,Length[posleft]}]
                ];
           Print["The terms"];
           Print[troublelist];
           Print["in equation ",j," are incompatible. Try to introduce an"];
           Print["auxiliary parameter with weight as coefficient of some of"];
           Print["these terms. Aborting the computations!"];
           ];
         ],
    Print["Try to introduce auxiliary parameters with weight as coefficients"];
    Print["into the system. Aborting the computations! "];
    ];
    CloseLog[];
    Abort[],
    Return[scalesol];
    ];
  ]; ]; (* end scaling *)

(*****************************************************************************)
(* subroutine5[]: a subroutine for constructformrho function                 *)
(*****************************************************************************)

subroutine5[nodims_,givenlist_] := Module[
{m,boundary,scale,var,len,tempwei,temppar,list,pair,ctm = 0,k,s},
list[0] = {};
For[m = 1,m <= nodims,m++,
   pair = Part[givenlist,m];
   var = Part[pair,1];
   scale = Part[pair,2];
   If[scale =!= 0,
     If[m == 1,
       boundary = Floor[rhorank/scale];
       list[m] = Table[{i*scale,var^i},{i,0,boundary}];
       ctm = m,
       len = Length[list[m-1]];
       list[m] = {};
       For[k = 1,k <= len,k++,
          tempwei = Part[Part[list[m-1],k],1];
          temppar = Part[Part[list[m-1],k],2];
          boundary = Floor[(rhorank-tempwei)/scale];
          For[s = 0,s <= boundary,s++,
             list[m] = Union[list[m],{{tempwei+s*scale,temppar*var^s}}];
             ctm = m;
             ];
          ];
       ];
     ];
   ];
Return[list[ctm]]; ]; (* end subroutine5 *)

(*****************************************************************************)
(* constructformrho[]: returns the form of rho                               *)
(*****************************************************************************)
 
constructformrho[nodims_,varscalelist_] := Module[
{list,len,tempwei,temppar,inttest,difflist,lendifflist,temp,
formrholist = {},i,j},
Print[" "];
Print["Starting the construction of the form of the density."];
list = subroutine5[nodims,varscalelist];
len = Length[list];
For[i = 1,i <= len,i++,
   tempwei = Part[Part[list,i],1];
   temppar = Part[Part[list,i],2];
   inttest = rhorank-tempwei;
   If[IntegerQ[inttest] && Not[FreeQ[temppar,u]],
     difflist = D[temppar,{t,inttest}];
     difflist = stripper[difflist,True,2];
     lendifflist = Length[difflist];
     For[j = 1,j <= lendifflist,j++,
        temp = Part[difflist,j];
        If[Intersection[formrholist,{temp}] === {},
          formrholist = Union[formrholist,{temp}]; 
          ];
        ];
     ];
   ];

lenformrho = Length[formrholist];
For[i = 1,i <= lenformrho,i++,
   formrho = formrho + c[i]*Part[formrholist,i];
   ];
Print[" "];
Print["For RANK = ",rhorank,", this is the form of the density rho:"];
Print[" "];
Print[subscriptform[formrho]];
Print[" "];

(* begin taken out 
(* 08/14/2004, HEREMAN CHRISTCHURCH *)
(* Form rho saved in input notation for later testing *)
formrhosaved = formrho;
Print["*** Density saved for later testing ***, formrhosaved: "];
Print[formrhosaved];
end taken out *)

If[formrho === 0,
  Print["The only density is the trivial density rho = 0."];
  Print[" "];
  Print["Aborting!"]; 
  CloseLog[]; Abort[] ] ]; (* end constructformrho *)               

(*****************************************************************************)
(* subroutine4[]: a subroutine for evaluate function                         *)
(*****************************************************************************)

subroutine4[givenexpr_] := Module[
{expr1,expr2,freerule = {},part1,lenexpr2,i},
expr1 = Expand[givenexpr];
expr2 = Numerator[Factor[expr1]];
expr2 = FactorList[expr2];
lenexpr2 = Length[expr2];
For[i = 1,i<= lenexpr2,i++,
   part1 = Part[Part[expr2,i],1];
   If[MemberQ[unknownlist,part1],
     freerule = Union[freerule,{part1 -> 1}]];
   ];
If[freerule =!= {},
  Print[" "];
  Print["Caution! There is a common (free) factor in the density."];
  Print[" "];
  Print["Setting ",freerule,"."];
  ];
Return[freerule]; ]; (* end subroutine4 *)

(*****************************************************************************)
(* evaluate[]: computes the density and flux with verification               *)
(*****************************************************************************)

evaluate[expr1_,expr2_,solrules_] := Module[
{newformrho,formfn,formjn = 0,expr,test,factortest,freerule},
newformrho = expr1 /. solrules;

(*
   %%%%%%%%%%%%%%%%%%%%% 
formfn = expr2 /. solrules;
formfn = Expand[formfn];
While[formfn =!= 0,
     expr = Part[formfn,1];
     formjn = formjn + expr;
     formfn = formfn - expr + (expr /. n->n+1);
     formfn = Expand[formfn];
     ]; 
Print[" "];
   %%%%%%%%%%%%%%%%%%%%% 
*)

If[newformrho =!= 0,
  freerule = subroutine4[newformrho];
  rho = newformrho /. freerule;
  rho = Expand[rho];

(* 
   %%%%%%%%%%%%%%%%%%%%%
  flux = formjn /. freerule;
  test = (D[rho,t]-flux + (flux /. n->n+1)) /. solrules;
  factortest = Expand[test];
  If[factortest =!= 0,factortest = Factor[factortest]];
  If[factortest =!= 0,
    Print[" "];
    Print["Automatic verification of the result FAILED!"]];
  Print[" "];
   %%%%%%%%%%%%%%%%%%%%%
*)

  Print["***************************************************************"];
  Print[" "];
  Print["This is the density rho:"];
  Print[" "];
  Print[subscriptform[rho]];
  Print[" "];
  split[subscriptform[rho]];
  Print["***************************************************************"];

(*
   %%%%%%%%%%%%%%%%%%%%%
  Print[" "];
  Print["The corresponding flux J is:"];
  Print[" "];
  Print[subscriptform[flux]];
  Print[" "];
  Print["***************************************************************"];
  Print["Result of explicit verification: D_t (rho) - J_n + J_{n+1} = ",
         factortest];
  Print["***************************************************************"],
   %%%%%%%%%%%%%%%%%%%%% 
*),

  Print["***************************************************************"];
  Print["There is NO DENSITY of this form!"];
  Print["***************************************************************"];
  ]; ]; (* end evaluate *)
 
(*****************************************************************************)
(* constructeqlist[]: forms the system of equations for the c[i]             *)
(*****************************************************************************)

constructeqlist[expr_] := Module[
{expr1,eqlist = {},expr2,templist1,templist2,lentemplist,s,
selectsame,coefficient},
expr1 = Expand[expr];
coefficient = expr1 //. u[i_][n_][t] :> 0;
eqlist = Union[eqlist,{coefficient}];
expr2 = expr1-coefficient;
If[expr2 =!= 0,
  templist1 = stripper[expr2,False,3];
  templist2 = Table[Part[Part[templist1,j],2],{j,1,Length[templist1]}];
  templist2 = Union[templist2];
  lentemplist = Length[templist2];
  For[s = 1,s <= lentemplist,s++,
     selectsame = Select[templist1,SameQ[#[[2]],Part[templist2,s]]&];
     coefficient = Sum[Part[Part[selectsame,j],1],{j,1,Length[selectsame]}];
     eqlist = Union[eqlist,{coefficient}];
     ];
  ];
Return[eqlist]; ]; (* end constructeqlist *)

(*****************************************************************************)
(* mysimplify[]: simplifies the system of equations for the c[i]             *)
(*****************************************************************************)

mysimplify[list_] := Module[
{newlist,lennewlist,k,fterm,simplelist = {}},
Print["Starting the simplification of the system."];
Print[" "];
newlist = Complement[list,{0}];
lennewlist = Length[newlist];
For[k = 1 ,k <= lennewlist,k++,
   fterm = Factor[newlist[[k]]];
   If[Part[fterm,0] === Times && NumberQ[Part[fterm,1]],
     fterm = fterm/Part[fterm,1] ];
   If[Intersection[Expand[simplelist],Expand[{fterm*(-1)}]] === {},
     simplelist = Union[simplelist,{fterm}]];
   ];
Return[simplelist]; ]; (* end mysimplify *)

(*****************************************************************************)
(* analyzer[]: determines the coefficients c[i] that must be zero, used in   *)
(*             the search for compatibility conditions                       *)
(*****************************************************************************)

analyzer[list_,len_] := Module[
{i,partlist,analyzelist = {}},
Print[" "];
Print["Starting the analysis of the system."];
Print[" "];
For[i = 1,i <= len,i++,
   partlist = Part[list,i];
   If[Length[partlist] === 1 && MemberQ[unknownlist,partlist],
     analyzelist = Union[analyzelist,{partlist}],
     If[Length[partlist] === 2 && MemberQ[unknownlist,partlist[[2]]] &&
      (MemberQ[parameters,partlist[[1]]] || MemberQ[weightpars,partlist[[1]]]),
      analyzelist = Union[analyzelist,{partlist[[2]]}] ] ];
   ];
Return[analyzelist]; ]; (* end analyzer *)

(*****************************************************************************)
(* picklhs[]: makes a system, setting LHS == 0                               *)
(*****************************************************************************)

picklhs[list_,k_] := Part[list,k] == 0;

(*****************************************************************************)
(* coefmat[]: creates the coefficient matrix of the system for the c[i]      *)
(*****************************************************************************)

coefmat[system_List,unknowns_List] := Module[
{i,mat = {},lensys},
If[system === {} || unknowns === {}, Print["Fatal Error"],
lensys = Length[system];
For[i = 1,i <= lensys,i++,
   mat = Append[mat,Table[Coefficient[Expand[Part[Part[system,i],1]],
   Part[unknowns,k]],{k,1,Length[unknowns]}]]]];
Return[mat]; ]; (* end coefmat *)

(*****************************************************************************)
(* main[] :                                                                  *)
(*****************************************************************************)

main[] := Module[
{rhot,myeqlist,seqlist,lenseqlist,maineqlist,analyzelist,lengthpar,
complexanalysis,syscondlist = {},i,j,k,syscond,inputlist,
inputpart,inputval,inputrule,comcond,comcondfac,sol,lengthsol,rules,
parameter,newmaineqlist,solc,solrules,comcondfactab},
rhot = D[formrho,t];
rhot = Expand[rhot];

(* 08/14/2004 HEREMAN, CHRISTCHURCH *)
(*
If[debugwillyCH, 
   Print["Time derivative of rho, rhot: "];
   Print[rhot]
  ]; 
*)

(* 08/14/2004 HEREMAN CHRISTCHURCH *)
(* Q: Where does it evaluate the expression on equations ? *)
(* A: Done through the equations because of the special format of *)
(* the left hand side as u[i][n_]'[t] := rhs *)
(* *)
Print[" "];
Print["Starting the construction of the system for the c[i]."];
Print[" "];
myeqlist = constructeqlist[rhot];
seqlist = mysimplify[myeqlist];
lenseqlist = Length[seqlist];
maineqlist = Flatten[Table[picklhs[seqlist,i],{i,1,lenseqlist}]];
Print["This is system for the c[i] (to be solved):"];
Print[" "];
Print[MapAll[Factor,maineqlist]];
Print[" "];
Print["List of unknown coefficients c[i]:"];
Print[" "];
Print[unknownlist];
Print[" "];
Print["Starting the solution of this system with ",
         Length[maineqlist], " equations."];
If[parameters =!= {}, (* then start for parameter test *)
  analyzelist = analyzer[seqlist,lenseqlist];
  lengthpar = Length[parameters];
  Print["The system for the coefficients c[i] has ",
         lengthpar," parameter(s)."];
  complexanalysis=lengthpar+Length[maineqlist]+Length[unknownlist];
  If[complexanalysis > 20,
    system = MapAll[Factor,maineqlist];
    coefmatrix = coefmat[system,unknownlist];
    Print["Rank and form of rho, the system for c[i], its"];
    Print["coefficient matrix, and the lists of unknowns and"];
    Print["parameters, will all be saved in the file worklog.m for"];
    Print["further analysis. For analysis use the Mathematica functions"];
    Print["Reduce, Solve, etc. Load the file worklog.m with the command"];
    Print["<<worklog.m ."];
    Save["worklog.m",myfile,rhorank,formrho,system,coefmatrix,
    unknownlist,parameters];
    Clear[coefmatrix,system];
    ];

  For[i = 1,i <= lengthpar,i++,
     syscondlist = Union[syscondlist,{parameters[[i]] != 0}];
     ];
  syscond = Apply[And,syscondlist];
  Print[" "];
  Print["Starting the search for compatibility conditions."];
  Print[" "];
  Print["The program will try all possible choices for non-vanishing"];
  Print["c[i]."];
  Print[" "];
  inputlist = Complement[unknownlist,analyzelist];
  If[inputlist === {},
    Print["All of the coefficients c[i] in the density have to vanish."];
    Print["The search for compatibility conditions ends."];
    ];
  While[inputlist =!= {},
   inputpart = Part[inputlist,1];
   inputval = inputpart == 1;
   inputlist = Complement[inputlist,{inputpart}];
   inputrule = ToRules[inputval];
   Print[" "];
   Print["-------------------------------------------------------"];
   Print[" "];
   Print["* Setting ",inputrule[[1]]," :"];
   Print[" "];
   Print["Computation of the compatibility conditions."];
   comcond = Eliminate[maineqlist /. inputrule,unknownlist];
   Print[" "];
   comcondfac = MapAll[Factor,comcond];
   (* Remove possible duplicates *)
   If[Head[comcondfac] === Or,
     comcondfactab = Table[Part[comcondfac,k], {k,1,Length[comcondfac]}];
     comcondfac = Union[comcondfactab]
     ];
   Print["This is the compatibility condition:"];
   Print[" "];
   Print[comcondfac];
   sol = Reduce[comcond && syscond,parameters];
   If[comcondfac === True,
     Print[" "];
     Print["The compatibility condition is satisfied without constraints"];
     Print["on the parameters."]
     ];
   If[sol =!= False,
     If[Not[FreeQ[sol,Or]],
       lengthsol = Length[sol],
       lengthsol = 1;
       ];
     For[j = 1,j <= lengthsol,j++,
        If[lengthsol == 1,
          rules = ToRules[sol],
          rules = ToRules[Part[sol,j]];
          ];
        If[comcondfac =!= True,
          For[k = 1,k <= lengthpar,k++,
             parameter = parameters[[k]] /. rules;
             Print[" "];
             Print["For ",parameters[[k]]," = ",parameter];
             ]
          ];
        newmaineqlist = maineqlist /. rules;
        solc = Flatten[Solve[newmaineqlist /. inputrule,unknownlist]];
        solc = Simplify[Union[solc,inputrule]];
        Print[" "];
        Print["The solution of the system is:"];
        Print[" "];
        Print[solc];
        solrules = Union[rules,solc];
        evaluate[formrho,rhot,solrules];
        inputlist = Intersection[inputlist,
             Union[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc]]];
        Clear[rules,solc,newmaineqlist];
        ], (* closes for *)
     Print[" "];
     Print["The system becomes inconsistent, or the compatibility conditions"];
     Print["require that one or more of the parameters are zero."];
     Print["Not acceptable!"];
     ] ] (* closes While *),
  solc = Flatten[Solve[maineqlist,unknownlist]];
  Print[" "];
  Print["Solution of the system:"];
  Print[" "];
  Print[solc];
  evaluate[formrho,rhot,solc];
  Clear[solc];
  ]; ]; (* end main *)

(* ------------------------------------------------------------------------- *)

(* set-up for collecting the data in a log file, transcript of computations *)

logfile="";
OpenLog[filename_String]:=(logfile = OpenWrite[filename];
  If[logfile === $Failed, Return[]];
  AppendTo[$Echo,logfile]; AppendTo[$Output,logfile];);
CloseLog[]:=($Echo = Complement[$Echo,{logfile}];
  $Output = Complement[$Output,{logfile}];
  Close[logfile];);

(* End of the auxiliary functions and procedures *)
(* Start of the executable part of the program *)

menu;
OpenLog[myfile];
commentinter[];
If[formrho === {}, formrho = 0];
If[ListQ[formrho] && formrho =!= {}, formrho=Part[formrho,1]]; 
Print[" "];
Print["Working with the data file for the ",name,"."];
Print[" "];

For[i = 1,i <= noeqs,i++,
   Print[" "];
   Print["Equation ",i," of the system with ",noeqs," equation(s):"];
   Print[" "];
   Print["  ", SequenceForm[u,Subscript[i],Subscript[","],Subscript[n]]',
         " = ", subscriptform[u[i][n]'[t]]];
   ]; (* end for *)

If[formrho === 0, (* point rhogiven *)
  eqlist = Flatten[Table[u[i][n]'[t],{i,1,noeqs}]];
  scalerules = scaling[eqlist];
  counternegweight = 0;
  listfreeweights = {};
  varscalelist = {};
  allchoices = {};
  Print[" "];
  Print["Program determines the weights of the variables (and parameters)."];
  For[i = 1,i <= noeqs,i++,
     weightu[i] = weightu[i] /. scalerules;
     If[i == 1, Print[" "]; Print["For the given system:"]; Print[" "]];
     Print["* weight of u[",i,"] is ",weightu[i],"."];
     varscalelist = Union[varscalelist,{{u[i][n][t],weightu[i]}}];
     If[Not[NumberQ[weightu[i]]],
       If[Length[weightu[i]] == 1,
         listfreeweights = Union[listfreeweights,{weightu[i]}]],
       If[weightu[i] < 0,counternegweight++];
       ];
     ];
  lenweightpars = Length[weightpars];
  For[i = 1,i <= lenweightpars,i++,
     weight[Part[weightpars,i]] =
             weight[Part[weightpars,i]] /. scalerules;
     Print["* weight of ",Part[weightpars,i]," is ",
           weight[Part[weightpars,i]],"."];
     varscalelist = Union[varscalelist,
                    {{Part[weightpars,i],weight[Part[weightpars,i]]}}];
     If[Not[NumberQ[weight[Part[weightpars,i]]]],
       If[Length[weight[Part[weightpars,i]]] == 1,
         listfreeweights =
            Union[listfreeweights,{weight[Part[weightpars,i]]}]],
       If[weight[Part[weightpars,i]] < 0,counternegweight++];
       ];
     ]; (* end for *)
  counterfreeweight = Length[listfreeweights];
  If[counternegweight =!= 0, (* begin If point 1 *)
    Print[" "];
    Print["One or more of the weights are negative."];
    Print["Negative weights are not allowed."];
    Print["Aborting the computations!"];
    CloseLog[]; Abort[],
    If[counterfreeweight =!= 0,
      If[counterfreeweight > 1,
        Print[" "];
        Print["Two or more of weights have freedom."];
        Print["Enter your values for the free weights or type Abort[];"];
        Print["Your values for the free weights can be entered by typing"];
        Print["`weightu[label] = val' or `weight[variable] = val' and"];
        Print["putting `;' between your entries."];
        Input[": "], (* else of counterfreeweight > 1 *)
        Print[" "];
        Print["Program will try to determine CHOICES for ",
              listfreeweights[[1]],":"];
        Print[" "];
        tempvarscalelist = Sort[Table[Part[varscalelist,k][[2]],
                                 {k,1,Length[varscalelist]}]];
        tempvarscalelist =
               Complement[tempvarscalelist,Select[tempvarscalelist,NumberQ]];
        lentempvarscalelist = Length[tempvarscalelist];
        Do[ (* begin Do loop *)
          Print["* CHOICE ",k];
          Print[" "];
          Print["Solving the equation: ",Part[tempvarscalelist,k]," = 1."];
          Print[" "];
          attempt = Flatten[Solve[Part[tempvarscalelist,k] == 1]];
          If[attempt =!= {},
            solfreeweight = Part[listfreeweights,1] /. attempt;
            While[solfreeweight <= 0,
                 solfreeweight = solfreeweight + 1;
                 Print["Since the weight was zero or negative,"];
                 Print["it was incremented with 1."];
                 Print[" "]
                 ];
            Print[Part[listfreeweights,1]," = ",solfreeweight];
            weightrule = Part[listfreeweights,1] -> solfreeweight;
            tempvarscalelistval = tempvarscalelist /. weightrule;
            If[MemberQ[NonNegative[tempvarscalelistval],False],
              Print["Choice is rejected!"],
              varscalelistval = varscalelist /. weightrule;
              Print[" "];
              Print["List of all the variables (and parameters) with their"];
              Print["weights:"];
              Print[" "];
              Print[varscalelistval];
              Print[" "];
              testvarscalelist = Table[Part[varscalelistval,k][[2]],
                        {k,1,Length[varscalelist]}];
              If[Union[Positive[testvarscalelist]] === {False},
                Print["Since all the weights of u[i] are zero"];
                Print["this choice is rejected!"];
                Print[" "],
                allchoices = Union[allchoices,{solfreeweight}];
                ];
              ];
            ],{k,1,lentempvarscalelist} ];  (* end do *)
        Print[" "];
        Print["Simple POSITIVE choices are considered."];
        Print[" "];
        intchoices = Select[allchoices,IntegerQ];
        fracchoices = Complement[allchoices,intchoices];
        If[intchoices =!= {},
          (* Picking the minimum integer choice *)
          intchoices = Min[intchoices];
          choicerule = Part[listfreeweights,1] -> intchoices, (* else *)
          If[fracchoices =!= {},
           (* Picking the minimum fractional choice *)
           fracchoices = Min[fracchoices];
           choicerule = Part[listfreeweights,1] -> fracchoices, (* else *)
           Print["Enter your choice for the value of ",Part[listfreeweights,1],
                 "by typing its value. NO semi-colon at the end!"];
           choicerule = Part[listfreeweights,1]-> Input[": "];
           ];
          ];
        varscalelistchoice = varscalelist /. choicerule;
        If[Length[allchoices] > 1,
          Print["List of all positive choices considered: "];
          Print[" "];
          Print[Union[Table[Part[listfreeweights,1] -> Part[allchoices,k],
                             {k,1,Length[allchoices]}]]];
          Print[" "];
          Print["Some of the weights of u[i], but not all, could be zero."];
          Print["In the data file you could enter your choices in the format"];
          Print["`weightu[label] = value;' or `weight[variable] = value;'."];
          Print[" "];
          Print["Program continues with the choice: ", choicerule];
          Print[" "];
          Print["corresponding to the weights:"];
          Print[" "];
          Print[varscalelistchoice];
          Print[" "];
          ];
        varscalelist = varscalelistchoice;
      ] ] ]; (* end if point 1 *)

  Print[" "];
  Print["The rank of rho should be an integer multiple of the lowest weight"];
  Print["of the DEPENDENT variable(s). Fractional weights are allowed."];
  Print[" "];
  If[ListQ[formrho] && formrho =!= {}, formrho=Part[formrho,1]]; 
  If[formrho === 0,
    If[Not[NumberQ[rhorank]], rhorank = Input["Enter the rank of rho: "]] 
    ];
  If[formrho === 0,
  Print["Computation of the density (and flux) of RANK = ",rhorank,"."];
  Print["  "];
  ];

  nodims = noeqs + Length[weightpars];
  varscalelist = Reverse[Sort[varscalelist,OrderedQ[{#1[[2]],#2[[2]]}]&]];
  constructformrho[nodims,varscalelist], 
  (* else point rhogiven *)
  formrho = Expand[formrho];
  If[Head[formrho] === Times,
    lenformrho = 1,lenformrho = Length[formrho];
    ];
  Print[" "];
  Print["The program will only test the given form of the density rho."];
  Print["No determination of scaling properties."];
  Print["You have given this form for the density rho :"];
  Print[" "];
  Print[subscriptform[formrho]];
  Print[" "];
  ]; (* end point rhogiven *)

unknownlist = Table[c[i],{i,1,lenformrho}];
(* setting weighted parameters to 1 *)
settoone[x_]:= x = 1;
weightpars = Map[settoone,weightpars];
weightpars = {};

Clear[i,j,k,lenformrho,counternegweight,counterfreeweight,scalerules,eqlist,
     lenweightpars,varscalelist,tempvarscalelist,attempt,listfreeweights,
     allchoices,varscalelistval,temparscalelistval,intchoices,fracchoices,
     choicerule,weightrule,nodims,highestordereq,weightu,name,weight,
     lentempvarscalelist,solfreeweight,varscalelistchoice,testvarscalelist];

main[];

Print[" "];
Print["*************************** SUMMARY ***************************"];
Print[" "];
Print["Total CPU time used in the current session is ",Round[TimeUsed[]],
      " seconds."];
Print[" "];
Print["To see the density type: rho, or subscriptform[rho]."];
Print[" "];
Print[" "];
Print["To split the density in independent pieces type: split[rho],"];
Print["or split[subscriptform[rho]]."];

(*
  %%%%%%%%%%%%%%%%
Print["To see the flux type: flux, or subscriptform[flux]."];
Print[" "];
  %%%%%%%%%%%%%%%%
*)

CloseLog[];
(* 08/04/2007 WH print statement added *)
Print["Package diffdens2007.m of August 20, 2007 successfully loaded."];

(* ******************************* END ************************************* *)