(* File DDEHOPER.M: all pieces for discrete homotopy operator in one file! *)
(* WH 10/04/2003 *)
(* Based on code of 10/01/2003 *)

(* ************************* BEGIN of all ***************************** *)

(* ************* BEGIN loadDiscreteHomotopyFunctions.m ************* *)

(* Time-stamp: <Fri May 30 14:16:14 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* loadDiscreteHomotopyFunctions                                         *)
(* Purpose: To load all the functions needed for the Discrete Homotopy  *)
(*          case.                                                       *)
(* Authors: Frances Martin, Adam T. Ringler                             *)
(* Input: none                                                          *)
(* Output: none                                                         *)
(* Created: May 28, 2003 @ 1:04 PM @ CSM                                *)
(* Last Modified: May 30, 2003 @ 2:14 PM @ CSM                          *)
(************************************************************************)

(* load all necessary functions *)
(* 
<<shiftExpression.m; 
<<shiftExtract.m; 
<<maxDownShift.m;
<<maxUpShift.m;
<<discreteEulerOperator.m;
<<discreteInteriorProduct.m;
<<discreteLambdaReplace.m;
<<discreteHomotopyOperator.m;
*)

(* ************ END loadDiscreteHomotopyFunctions.m ****************** *)

(* *************** BEGIN discreteEulerOperator.m ********************** *)

(* Time-stamp: <Tue Jun  3 19:14:20 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* discreteEulerOperator[expression_, orderOfOperator_, dependVarHead_, *)
(*                       dependVarNum_];                                *)
(* Purpose: To calculate nth order Euler Operators for DDE's            *)
(* Authors: Ingo Kabirschke, Ryan Sayers                                *)
(* Input: An expression with no negative shifts: expression,            *)
(*        the order of the Euler Operator: orderOfOperator,             *)
(*        a dependent variable: dependVar                               *)
(* Output: The result of applying the nth order Euler Operator to a DDE *)
(* Created: May 28, 2003 at CSM                                         *)
(* Last Modified: June 3, 2003 at 7:12 PM at CSM                        *)
(************************************************************************)

(* Debug Flags for the DiscreteEulerOperator *)
debugDiscreteEulerOpInput = False;
debugEulerMaxShift = False;
debugEulerSum = False;
debugDiscreteEulerResult = False;

discreteEulerOperator[expression_, orderOfOperator_, 
                      dependVarHead_, dependVarNum_] :=
  Module[{maxShift, discreteEulerSum, discreteEulerResult},

    (* Debugging input *)
    If[debugDiscreteEulerOpInput,
       Print["Computing the ", orderOfOperator,
             " order discreteEulerOperator on "];
       Print[expression];
       Print["with respect to ", dependVarHead];
    ];

    (* Checks for negative shifts - if found, returns error message *)
    (* and aborts calculation                                       *)
    If[maxDownShift[expression, dependVarHead] < 0,
       Print["The discrete Euler operator needs an expression without", 
             " negative shifts."];
       Return[False];
    ];

    (* Finds the maximum positive shift in the expression *)
    maxShift = maxUpShift[expression, dependVarHead];

    (* Debugging maximum positive shift *)
    If[debugEulerMaxShift, 
      Print["The max shift is ", maxShift, "."]
    ];
   
    (* Calculate the Euler Operator before taking the partial     *)
    (* derivative - taken from Hereman's Notes (Newest), Page 31  *)
    discreteEulerSum = Sum[Binomial[k, orderOfOperator] * 
                           shiftExpression[expression, u, -k],
                          {k, orderOfOperator, maxShift}
                       ];

    (* Debugging result before partial derivative *)
    If[debugEulerSum, 
      Print["The Euler sum is "]; 
      Print[discreteEulerSum];
    ];

    (* Take the partial derivative with respect to the dependent   *)
    (* variable to complete the calculation of the Euler Operator. *)
    (* Also taken from Hereman's notes (Newest), page 31           *)
    discreteEulerResult = Expand[ D[discreteEulerSum, 
                                   dependVarHead[dependVarNum][n, t]]
                          ];

    (* Debugging result *)
    If[debugDiscreteEulerResult, 
       Print["The result of applying the ", orderOfOperator,
             " order discrete Euler operator with respect to ",
             dependVarHead[dependVarNum], " to"];
       Print[expression];
       Print["is"];
       Print[discreteEulerResult];
    ];

    (* Return nth order Euler Operator *)
    Return[discreteEulerResult];
];

(* Print["discreteEulerResult loaded successfully."]; *)

(* *************** END discreteEulerOperator.m ********************** *)

(* ***** ********* BEGIN discreteInteriorProduct.m ************* *)

(* Time-stamp: <Fri May 30 15:08:44 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* discreteInteriorProduct[expression_, dependVarHead_, dependVarNum_]  *)
(* Purpose: To calculate the discrete interior product of a given       *)
(*          expression with respect to the dependent variable           *)
(* Authors: Frances Martin, Adam Ringler, Ryan Sayers                   *)
(* Input: An expression: expression                                     *)
(*        A dependent variable head: dependVarHead                      *)
(*          Example: for u[1], u[2], u[3], the head is u.               *)
(*        The index of the dependent variable: dependVarNum             *)
(*          Example: for u[1], the dependVarNum is 1.                   *)
(* Output: The interior product of the expression                       *)
(* Created: May 28, 2003, @ 1:43 PM @ CSM                               *)
(* Last Modified: May 30, 2003 @ 3:08 PM @ CSM                          *)
(************************************************************************)

(* Debug Flags for the DiscreteInteriorProduct *)
debugDiscreteIntProdInput = False;
debugDiscreteIntProdTerms = False;
debugDiscreteIntProdResult = False;

discreteInteriorProduct[expression_, dependVarHead_, dependVarNum_] :=
  Module[{i, l, intProdResult, maxShift},

    (* Debugging the input *)
    If[debugDiscreteIntProdInput,
       Print["Applying the Discrete Interior Product to "];
       Print[expression];
    ];

    (* Finds the maximum positive shift in the expression. This will *)
    (* be used to limit the sum in the interior product calculation. *)
    maxShift = maxUpShift[expression, dependVarHead];

    (* Calculating the interior product using page 32 of Dr. Hereman notes *)
    (* expressing (D - I)^i as a binomial expansion                        *)
    intProdResult =
      Expand[
        Sum[
          Sum[
            (-1)^(i-l) * Binomial[i, l] * 
            shiftExpression[
              dependVarHead[dependVarNum][n, t] * 
              discreteEulerOperator[expression, i + 1, dependVarHead, 
                                    dependVarNum],
              dependVarHead, l],
            {l, 0, i}],
          {i, 0, maxShift - 1}]
      ];         
                         
    (* Debugging the terms from calculating the interior product *)
    (* Outputs the terms in tabular form                         *)
    If[debugDiscreteIntProdTerms,
      Print["The terms of the interior product are: "];
      Print[
        Table[
          Sum[
            (-1)^(i - l) * Binomial[i, l] * 
            shiftExpression[dependVarHead[dependVarNum][n, t] * 
                            discreteEulerOperator[expression, 
                                                  i + 1, dependVarHead, 
                                                  dependVarNum],
                            dependVarHead, l],
             {l, 0, i}],
          {i, 0, Max[{0, maxShift - 1}]}
        ]
      ]
    ];

    (* Debugging the result of applying the interior product *)
    If[debugDiscreteIntProdResult,
       Print["The Interior Product of "];
       Print[expression];
       Print[" with respect to ", dependVarHead, "[", dependVarNum, 
             "] is: "];
       Print[intProdResult];
    ];

    (* Return the interior product *)
    Return[intProdResult];

];

(* Print["discreteInteriorProduct loaded successfully."]; *)

(* *************** END discreteInteriorProduct.m ************* *)

(* ********************* BEGIN shiftExpression.m ******************** *)

(* Time-stamp: <Fri May 30 14:39:16 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* shiftExpression[expression_, dependVarHead_, shiftAmt_]              *)
(* Purpose: To shift a discrete expression by some finite amount        *)
(* Authors: Ingo Kabirschke, Ryan Sayers                                *)
(* Input: An expression: expression,                                    *) 
(*        a dependent variable head: dependVarHead,                     *)
(*          (For example, for u[1], u[2], and u[3] the head would be u) *)
(*        the distance to shift the variables: shiftAmt                 *) 
(* Output: The expression shifted by shiftAmt                           *)
(* Created: May 28, 2003                                                *)
(* Last Modified: May 30, 2003, at 2:35 PM at CSM                       *)
(************************************************************************)

(* DEBUG FLAGS *)
debugShiftInput = False;
debugShiftResult = False;

shiftExpression[expression_, dependVarHead_, shiftAmt_] := 
  Module[{shiftReplaceRule, shiftResult},
    
    (* Debugging input *)
    If[debugShiftInput, 
      Print["Shifting "];
      Print[expression];
      Print[" by ", shiftAmt];
    ];
      
    (* Define shifting rule: replace n with n + shift *)
    shiftReplaceRule = {dependVarHead[i_][n_, t_] -> 
                        dependVarHead[i][n + shiftAmt, t]};

    (* Apply shifting rule *)
    shiftResult = expression /. shiftReplaceRule;

    (* Debugging result *)
    If[debugShiftResult,
      Print["The result of shiftExpression is "];
      Print[expression];
    ];

    (* Return shifted expression *)
    Return[shiftResult];
  ];

(* Print["shiftExpression loaded successfully."]; *)

(* ********************* END shiftExpression.m ******************** *)

(* ****************** BEGIN discreteLambdaReplace.m **************** *)

(* Time-stamp: <Fri May 30 15:25:30 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* discreteLambdaReplace[expression_, dependVarHead_]                   *)
(* Purpose: To replace the dependent variable with lambda * dependVar   *)
(* Authors: Frances Martin, Kara Namanny and Ingo Kabirschke            *)
(* Input: An expression: expression,                                    *) 
(*        a dependent variable head: dependVarHead,                     *)
(*          (For example, for u[1], u[2], and u[3], the head would be u)*)
(* Output: The expression with lambda replacement                       *)
(* Created: May 29, 2003 10:35 AM at CSM                                *)
(* Last Modified: May 30, 2003 3:25 PM at MLRC                          *)
(************************************************************************)

(* Debug Flags for DiscreteLambdaReplace *)
debugDiscreteLambdaReplaceInput = False;
debugDiscreteLambdaReplaceResult = False;

discreteLambdaReplace[expression_, dependVarHead_] :=
  Module[{i, discreteLambdaReplaceRule, discreteLambdaReplaceResult},

    (* Debugging input *)
    If[debugDiscreteLambdaReplaceInput, 
      Print["Replacing the dependent variables, ", 
            dependVarHead, "[i], in the expression, "];
      Print[expression]; 
      Print[", with lambda * ", dependVarHead, "[i]"];
    ];

    (* Define the replacement rule - replaces each occurance of the    *)
    (* dependent variable head with lambda * dependVarHead             *)
    discreteLambdaReplaceRule =
      dependVarHead[i_][n_, t] -> lambda * dependVarHead[i][n, t];

    (* Apply replacement rule *) 
    discreteLambdaReplaceResult = Expand[expression /. 
                                         discreteLambdaReplaceRule];

    (* Debugging result *)
    If[debugDiscreteLambdaReplaceResult, 
      Print["The result of replacement is: "];
      Print[discreteLambdaReplaceResult];
    ];

    (* Return new expression with lambda *)
    Return[discreteLambdaReplaceResult];
  ];

(* Print["Discrete Lambda Replace loaded successfully."]; *)

(* ****************** END discreteLambdaReplace.m **************** *)

(* ********************* BEGIN discreteHomotopyOperator.m ************ *)

(* Time-stamp: <Wed Jun 11 15:46:02 MDT 2003>        *)
(************************************************************************)
(* discreteHomotopyOperator[expression_, dependVarHead_, numDependVar_] *)
(* Purpose: To calculate the Discrete Homotopy Operator of              *)
(*          an expression                                               *)
(* Authors: Ingo Kabirschke, Frances Martin, Kara Namanny               *)
(* Input: An expression: expression,                                    *) 
(*        a dependent variable head: dependVarHead,                     *)
(*          (For example, for u[1], u[2], and u[3] the head would be u) *)
(*        the # of elements in dependent variable vector: numDependVar  *) 
(* Output: The Discrete Homotopy Operator of the expression             *)
(* Created: May 29, 2003 @ 11:42 AM @ CSM                               *)
(* Last Modified: June 3, 2003, at 8:14 PM at CSM                       *)
(************************************************************************)

(* Debug Flags for the DiscreteHomotopyOperator *)
debugDiscreteHomotopyOpInput = False;
debugShiftedExpression = False;
debugDiscreteHomotopyOpResult = False;

discreteHomotopyOperator[expression_, dependVarHead_, numDependVar_] := 
  Module[{i, maximumDownShift, shiftedExpression, fullDiscreteIntProd, 
          lambdaExpression, integrand, shiftedHomotopyOpResult, 
          homotopyOpResult},

    (* Debugging input *)
    If[debugDiscreteHomotopyOpInput,
       Print["Applying the Discrete Homotopy Operator to "];
       Print[expression];
       Print["The dependent variable vector has head ", dependVarHead,
             " and ", numDependVar, " components."];
    ];

    (* Find the maximum downshift *)
    maximumDownShift = maxDownShift[expression, dependVarHead];
    
    (* Remove any negative shifts from expression, if needed. *)
    (* Ie. Shift up the amount of the maximum negative shift. *)
    If[maximumDownShift < 0,
       shiftedExpression = shiftExpression[expression, dependVarHead, 
                                    Abs[maximumDownShift]],
       (* Else (if no negative shifts) *)
       shiftedExpression = expression;
    ];

    (* Debugging the shifted expression *)
    If[debugShiftedExpression,
       Print["The shifted expression is: "];
       Print[shiftedExpression];
    ];

    (* Checking for integrability (Zeroth Euler Operator = 0)        *)
    (* If the expression is not integrable, outputs an error message *)
    (* and aborts further caculations. (Returns false)               *)
    For[i = 1, i <= numDependVar, i++,
      If[discreteEulerOperator[shiftedExpression, 0, dependVarHead, i] =!= 0,
         (* Then *)
         Print["The expression input (after shifting), "];
         Print[shiftedExpression];
         Print[" is not integrable."];
         Print["Aborting calculations."];
         Return[False];
      ];
    ];

    (* The steps below are taken from Hereman's notes (newest) page 26   *)
    (*********************************************************************)

    (* Find the interior product (little j) of the shifted expression *)
    fullDiscreteIntProd = Sum[     
                            discreteInteriorProduct[shiftedExpression, 
                                                    dependVarHead, i],
                            {i, 1, numDependVar}
                          ];

    (* Replace dependent variable with lambda * dependent variable *)
    lambdaExpression = 
      discreteLambdaReplace[fullDiscreteIntProd, dependVarHead];

    (* Divide new expression by lambda *)
    integrand = Expand[lambdaExpression / lambda];

    (* Integrate with respect to lambda from 0 to 1 to find flux (big J) *)
    shiftedHomotopyOpResult = Expand[Integrate[integrand, {lambda, 0, 1}]];

    (* If original expression was shifted to remove negative shifts, *)
    (* shift back down the same amount to get final result.          *)
    If[maximumDownShift < 0,
       homotopyOpResult = shiftExpression[shiftedHomotopyOpResult, 
                                          dependVarHead, maximumDownShift],
       (* Else, do not shift the final result *)
       homotopyOpResult = shiftedHomotopyOpResult;
    ]; 

    (* Debugging output *)
    If[debugDiscreteHomotopyOpResult,
       Print["The result of applying the Discrete Homotopy Operator is "];
       Print[homotopyOpResult];
    ];

    (* Return the result of applying the Homotopy Operator to expression *)
    Return[homotopyOpResult];
];

(* Print["Discrete Homotopy Operator loaded successfully."]; *)

(* ********************* END discreteHomotopyOperator.m ************ *)

(* ********************* BEGIN shiftExtract.m ******************** *)

(* Time-stamp: <Fri May 30 14:20:49 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* shiftExtract[expression_, dependVarHead_]                            *)
(* Purpose: To extract a list of the shifts in expression for           *)
(*          dependVarHead                                               *)
(* Authors: Ingo Kabirschke and Ryan Sayers                             *)
(* Input: An expression: expression,                                    *)
(*        the dependent variable head: dependVarHead                    *)
(*          Example: for u[1], u[2], u[3] this would be u               *)
(* Output: The list of shifts of dependVarHead in expression            *)
(* Created: May 28, 2003 at CSM                                         *)
(* Last Modified: May 30, 2003 at 2:16 PM at CSM                        *)
(************************************************************************)

(* DEBUG FLAGS for shiftExtract *)
debugShiftExtractInput = False;
debugShiftExtractOutput = False;

shiftExtract[expression_, dependVarHead_] := Module[
  {shiftExtractRule, shiftExtractCases, shiftList},

  (* Debugging input *)  
  If[debugShiftExtractInput,
     Print["The input to shiftExtract is "];
     Print[expression];
  ];

  (* Defining extraction rules that change an occurance of the dependent*)
  (* variable to a number, representing the dependent variable's shift  *)
  (* Example: u[1][n + 3, t] would be replaced by 3                     *)
  shiftExtractRule = 
    {dependVarHead[i_][n + shift_, t_] -> shift,
     dependVarHead[i_][n, t_] -> 0};

  (* Makes a list of every occurance of the dependent variables *)
  (* in the expression                                          *)
  shiftExtractCases = 
    Cases[expression, dependVarHead[i_][n_, t_], {0, Infinity}];

  (* Applying extraction rule - this results in a list of numbers that *)
  (* represent the shifts of each occurance of the dependent variable. *)
  shiftList = shiftExtractCases /. shiftExtractRule;

  (* Debugging result *)
  If[debugShiftExtractOutput,
     Print["The output of shiftExtract is "];
     Print[shiftList];
  ];

  (* Return list of shifts *)
  Return[shiftList];
];

(* Print["shiftExtract loaded successfully."]; *)

(* ********************* END shiftExtract.m ******************** *)

(* ********************* BEGIN maxDownShift.m ******************** *)

(* Time-stamp: <Fri May 30 14:32:11 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* maxDownShift[expression_, dependVarHead_]                            *)
(* Purpose: To find the maximum negative shift in a given expression    *)
(* Authors: Frances Martin, Kara Namanny, Adam Ringler                  *)
(* Input: a discrete expression: expression                             *)
(*        a dependent variable head: dependVarHead                      *)
(*          Example: for the variables u[1], u[2], u[3] the head is u   *)
(* Output: maximum negative shift in expression                         *)
(* Created: May 28, 2003 @ 10:56 AM @ CSM                               *)
(* Last Modified: May 30, 2003 @ 2:30 PM @ CSM                          *)
(************************************************************************)

(* Debug Flags for MaxDownShift *)
debugMaxDownShiftInput = False;
debugMaxDownShiftResult = False;

maxDownShift[expression_, dependVarHead_] := 
  Module[{maxNegativeShift},

  (* Debugging input *)
  If[debugMaxDownShiftInput,
    Print["Finding the maximum negative shift in "];
    Print[expression];
  ];

  (* Find the maximum negative shift in a given expression, using  *)
  (* Mathematica's Min function                                    *)
  maxNegativeShift = Min[shiftExtract[expression, dependVarHead]];

  (* Debugging result *)
  If[debugMaxDownShiftResult,
    Print["The maximum negative shift in "];
    Print[expression];
    Print[" is: ", maxNegativeShift];
  ];

  (* Return maximum negative shift in expression *)
  Return[maxNegativeShift];
];

(* Print["maxDownShift loaded successfully."]; *)

(* ********************* END maxDownShift.m ******************** *)

(* ********************* BEGIN maxUPShift.m ******************** *)

(* Time-stamp: <Fri May 30 14:29:16 Mountain Daylight Time 2003>        *)
(************************************************************************)
(* maxUpShift[expression, dependVarHead]                                *)
(* Purpose: To find the maximum positive shift in a given expression    *)
(* Authors: Frances Martin, Kara Namanny, Adam Ringler                  *)
(* Input: a discrete expression: expression                             *)
(*        a dependent variable head: dependVarHead                      *)
(*          Example: for u[1], u[2], u[3] the head is u                 *)
(* Output: The maximum positive shift in expression                     *)
(* Created: May 28, 2003 @ 10:56 AM @ CSM                               *)
(* Last Modified: May 30, 2003 @ 2:29 PM @ CSM                          *)
(************************************************************************)

(* Debug Flags for MaxUpShift *)
debugMaxUpShiftInput = False;
debugMaxUpShiftResult = False;

maxUpShift[expression_, dependVarHead_] := Module[{maxPositiveShift},

  (* Debugging input *)
  If[debugMaxUpShiftInput,
    Print["Finding the maximum positive shift in "];
    Print[expression];
  ];

  (* Finding the maximum positive shift in the given expression, using *)
  (* Mathematica's Max function                                        *)
  maxPositiveShift = Max[shiftExtract[expression, dependVarHead]];

  (* Debugging result *)
  If[debugMaxUpShiftResult,
    Print["The maximum positive shift in "];
    Print[expression];
    Print[" is: ", maxPositiveShift];
  ];

  (* Returning the maximum positive shift in the given expression *)
  Return[maxPositiveShift];
];

(* Print["maxUpShift loaded successfully."]; *)

(* ********************* END maxUPShift.m ******************** *)

Print["Discrete Homotopy Operator code loaded successfully."];

(* ************************* END of all ***************************** *)