BeginPackage["DrawLowerEnvelope`"]


(*Function BinarySearchSolveIncreasing solves all increasing functions*)

  BinarySearchSolveIncreasing[f_, n_, maxError_, variableName_] := 
  Module[{},
  (*Set the range of real numbers in which we will find our answer*)
      rbegin = 0;
    rend = 100;
    (*Choose our initial guess*)index = (rend + rbegin)/2;
    val = ReplaceAll[f, variableName -> index];
    (*Search until we find an answer*)
      While[(rbegin != rend),(*If the search range was of length 1, 
        make it length 0*)If[(rend - rbegin) <= maxError, rbegin = rend, 
        If[val > n, rend = index, rbegin = index]];
      index = (rend + rbegin)/2;
      val = ReplaceAll[f, variableName -> index];];
    index]

(*Function BinarySearchSolveDecreasing solves all decreasing functions************************************************************************)


BinarySearchSolveDecreasing[f_, n_, maxError_, variableName_] := 
  Module[{},
  
  (*Set the range of real numbers in which we will find our answer*)
      rbegin = 0;
    rend = 100;
    (*Choose our initial guess*)
    index = (rend + rbegin)/2;
    val = ReplaceAll[f, variableName -> index];
    (*Search until we find an answer*)
      While[(rbegin != rend),(*If the search range was of length 1, 
        make it length 0*)If[(rend - rbegin) <= maxError, rbegin = rend, 
        If[val > n, rbegin = index, rend = index]];
      index = (rend + rbegin)/2;
      val = ReplaceAll[f, variableName -> index];];
    index]

(*Function BinarySearchSolve takes as parameters a function f and an integer n, 
  and returns either a real number x such that f(x) = 
    n or - 1 if no such non - 
      negative integer exists.************************************************************************)




BinarySearchSolve[f_, n_, maxError_, variableName_] := 
  Module[{}, result = -1;
    If[ReplaceAll[f, variableName -> 99] > ReplaceAll[f, variableName -> 1], 
      result = BinarySearchSolveIncreasing[f, n, maxError, variableName], 
      result = BinarySearchSolveDecreasing[f, n, maxError, variableName]];
    result]

(*Function NextEqualsInConfig finds the beginning of a closed group, 
  starting from a given index and going backwards.************************************************************************)


(*NextEqualsInConfig[index_, config_] := Module[{},
    
    If[index == 1, Return[1]];
    
    If[config[[index - 1]] == "=", Return[index],
      If[index == 1, Return[1], 
        For[i = index - 1, i > 0, i--, 
          If[config[[i]] == "=" || config[[i]] == "<", Return[i], 
            Return[1]]]]
      ]
    ]
*)

(* NextEqualsInConfig[index_, config_] := Module[{i}, 
    
       
          If[index == 1, Return[1], 
                For[i= index - 1,i> 0, i--, 
                        If[config[[i]] == "=" || config[[i]] == "<", Return[i+1]]
                ]; 
                Return[1];
         ]; 
       
      
    ]
*)

NextEqualsInConfig[index_, config_] := Module[{i}, 
    
       
         For[i= index - 1,i> 0, i--, 
                If[config[[i]] == "=" || config[[i]] == "<", Return[i+1]]
         ]; 
         Return[1];     
    ] 



(*This function calculates the flow time************************************************************************)

CalculateFlowTime[InputPlast_, currentConfig_] := Module[{}, flowTime = 0;
    
    runningTimesX = {};
    completionTimes = {};
    allInTermsOfPlast = {};
    For[i = Length[currentConfig], i > 0, i--, 
      runningTimesX = Append[runningTimesX, -1]];
    For[i = Length[currentConfig], i > 0, i--, 
      completionTimes = Append[completionTimes, -1]];
    For[i = 1, i <= Length[currentConfig], i++, 
      allInTermsOfPlast = Append[allInTermsOfPlast, GlobalListOfJobs[[i]]]];
    allInTermsOfPlast[[Length[currentConfig]]] = Plast;
    Plast = InputPlast;
    (*Determine all Pi's in terms of Plast*)(*Evaluate "less than" case, 
      since it's easiest*)For[i = Length[currentConfig], i > 0, i--, 
      If[currentConfig[[i]] == "<", allInTermsOfPlast[[i]] = Plast]];
    (*Evaluate "greater than" case*)For[i = Length[currentConfig], i > 0, i--,
       If[currentConfig[[i]] == ">", 
        allInTermsOfPlast[[i]] = allInTermsOfPlast[[i + 1]] + Plast]];
    allInTermsOfPlast;
    (*Evaluate "equals" case*)For[loopCounter = Length[currentConfig], 
      loopCounter > 0, loopCounter--, 
      If[currentConfig[[loopCounter]] == "=", 
        Module[{}, currentIndex = loopCounter;
          nextEquals = NextEqualsInConfig[currentIndex, currentConfig];
          
          rDifference = 
            JobReleases[[loopCounter + 1]] - JobReleases[[nextEquals]];
          
          summation = 
            Sum[allInTermsOfPlast[[j]]^(-1/alpha), {j, nextEquals, 
                loopCounter}];
          temp = allInTermsOfPlast[[loopCounter]];
          
          bs = BinarySearchSolve[summation, rDifference, errorForBinarySearch,
               allInTermsOfPlast[[loopCounter]]];
          allInTermsOfPlast = ReplaceAll[allInTermsOfPlast, temp -> bs];]]];
    (*REPORT NEW allInTermsOfPlast*)allInTermsOfPlast;
    (*CALCULATE RUNNING TIMES Xi's*)For[i = 1, i <= Length[allInTermsOfPlast],
       i++, runningTimesX[[i]] = (allInTermsOfPlast[[i]]^(-1/alpha))];
    (*REPORT CALCULATES runningTimesX*)runningTimesX;
    (*CALCULATE COMPLETION TIMES*)completionTimes[[1]] = runningTimesX[[1]];
    For[i = 2, i <= Length[runningTimesX], i++, 
      completionTimes[[i]] = 
        runningTimesX[[i]] + 
          Max[completionTimes[[i - 1]], JobReleases[[i]]]];
    completionTimes;
    (*CALCULATE TOTAL FLOW TIME*)For[i = 1, i <= Length[completionTimes], i++,
       flowTime = completionTimes[[i]] + flowTime];
    flowTime]

(*This function calculates completion times and powers************************************************************************)



FindAllPowers[InputPlast_, currentConfig_] :=
  
  
  Module[{}, allInTermsOfPlast = {};
    For[i = 1, i <= Length[currentConfig], i++, 
      allInTermsOfPlast = Append[allInTermsOfPlast, GlobalListOfJobs[[i]]]];
    allInTermsOfPlast[[Length[currentConfig]]] = Plast*1.0;
    Plast = InputPlast;
    (*Determine all Pi's in terms of Plast*)(*Evaluate "less than" case, 
      since it's easiest*)
    For[i = Length[currentConfig], i > 0, i--, 
      If[currentConfig[[i]] == "<", allInTermsOfPlast[[i]] = Plast*1.0]];
    
    (*Evaluate "greater than" case*)For[i = Length[currentConfig], i > 0, i--,
       If[currentConfig[[i]] == ">", 
        allInTermsOfPlast[[i]] = allInTermsOfPlast[[i + 1]]*1.0 + Plast*1.0]];
    
    (*Evaluate "equals" case*)
    
    For[loopCounter = Length[currentConfig], loopCounter > 0, loopCounter--, 
      If[currentConfig[[loopCounter]] == "=", 
        Module[{}, currentIndex = loopCounter;
          nextEquals = NextEqualsInConfig[currentIndex, currentConfig];
          
          
          rDifference = 
            JobReleases[[loopCounter + 1]] - JobReleases[[nextEquals]];
          
          
          summation = 
            Sum[allInTermsOfPlast[[j]]^(-1/alpha), {j, nextEquals, 
                loopCounter}];
          
          temp = allInTermsOfPlast[[loopCounter]];
          
          bs = BinarySearchSolve[summation, rDifference, errorForBinarySearch,
               allInTermsOfPlast[[loopCounter]]];
          
          allInTermsOfPlast = ReplaceAll[allInTermsOfPlast, temp -> bs*1.0];]]];
    
    (*Return allInTermsOfPlast*)
    Return[allInTermsOfPlast];
    ]

FindCompletionTimes[allInTermsOfPlast_] := Module[{}, runningTimesX = {};
    
    completionTimes = {};
    For[i = Length[currentConfig], i > 0, i--, 
      runningTimesX = Append[runningTimesX, -1]];
    For[i = Length[currentConfig], i > 0, i--, 
      completionTimes = Append[completionTimes, -1]];
    (*CALCULATE RUNNING TIMES Xi's*)For[i = 1, i <= Length[allInTermsOfPlast],
       i++, runningTimesX[[i]] = (allInTermsOfPlast[[i]]^(-1/alpha))];
    (*REPORT CALCULATES runningTimesX*)runningTimesX;
    (*CALCULATE COMPLETION TIMES*)
    completionTimes[[1]] = runningTimesX[[1]]/1.0;
    
    For[i = 2, i <= Length[runningTimesX], i++, Module[{},
          (*
            If[runningTimesX[[i]] \[Element] Reals, Print["real"]
                  Module[{},
                    Print["runningTimesX[[i]]=", runningTimesX[[i]]];
                    
                    Print["completionTimes[[i-1]]=", 
                      completionTimes[[i - 1]]];
                    Print["JobReleases[[i]]=", JobReleases[[i]]];
                    ]
                ]
              
              If[runningTimesX[[i]] \[Element] Reals, 1 + 1,
                  Print["runningTimesX[[", i, "]]=", runningTimesX[[i]]/1.0]
                  ];*)
          
          
          
          completionTimes[[i]] = 
            runningTimesX[[i]] + 
              Max[completionTimes[[i - 1]]/1.0, JobReleases[[i]]]];
      ];
    Return[completionTimes];
    
    ]


(*Function calculates the energy level for a given set of Pi's*******************************************************************************)



CalculateEnergyUsed[InputPlast_, currentConfig_] := 
  Module[{}, runningTimesX = {};
    
    completionTimes2 = {};
    allInTermsOfPlast2 = {};
    For[i = Length[currentConfig], i > 0, i--, 
      runningTimesX = Append[runningTimesX, -1]];
    For[i = Length[currentConfig], i > 0, i--, 
      completionTimes2 = Append[completionTimes2, -1]];
    For[i = 1, i <= Length[currentConfig], i++, 
      allInTermsOfPlast2 = Append[allInTermsOfPlast2, GlobalListOfJobs[[i]]]];
    allInTermsOfPlast2[[Length[currentConfig]]] = Plast;
    Plast = InputPlast;
    (*Determine all Pi's in terms of Plast*)(*Evaluate "less than" case, 
      since it's easiest*)For[i = Length[currentConfig], i > 0, i--, 
      If[currentConfig[[i]] == "<", allInTermsOfPlast2[[i]] = Plast]];
    (*Evaluate "greater than" case*)For[i = Length[currentConfig], i > 0, i--,
       If[currentConfig[[i]] == ">", 
        allInTermsOfPlast2[[i]] = allInTermsOfPlast2[[i + 1]] + Plast]];
    allInTermsOfPlast2;
    (*Evaluate "equals" case*)For[loopCounter = Length[currentConfig], 
      loopCounter > 0, loopCounter--, 
      If[currentConfig[[loopCounter]] == "=", 
        Module[{}, currentIndex = loopCounter;
          nextEquals = NextEqualsInConfig[currentIndex, currentConfig];
          
          rDifference = 
            JobReleases[[loopCounter + 1]] - JobReleases[[nextEquals]];
          
          summation = 
            Sum[allInTermsOfPlast2[[j]]^(-1/alpha), {j, nextEquals, 
                loopCounter}];
          temp = allInTermsOfPlast2[[loopCounter]];
          
          bs = BinarySearchSolve[summation, rDifference, errorForBinarySearch,
               allInTermsOfPlast2[[loopCounter]]];
          allInTermsOfPlast2 = ReplaceAll[allInTermsOfPlast2, temp -> bs];]]];
    energyUsed = 
      Sum[allInTermsOfPlast2[[i]]^((alpha - 1)/alpha), {i, 
          Length[allInTermsOfPlast2]}];
    energyUsed]


FindNewConfiguration[completionTimes_, releaseTimes_, allInTermsOfPlast_, 
    currentConfiguration_] := Module[{},
 
    newConfig = currentConfiguration;
    For[i = 1, i <= Length[currentConfiguration] - 1, i++, Module[{},
          
          If[newConfig[[i]] == "=" && 
              allInTermsOfPlast[[i]] > (allInTermsOfPlast[[i + 1]] + 
                    allInTermsOfPlast[[Length[allInTermsOfPlast]]]), 
            newConfig[[i]] = ">"];
          
          
          If[newConfig[[i]] == "<" && 
              completionTimes[[i]] > releaseTimes[[i + 1]], 
            newConfig[[i]] = "="];
          
          ];
      ];
    
    
    Return[newConfig];
    ]



FindOptimalChangePoint[endPoint_, currentConfig_] := Module[{},
    ending = endPoint;
    For[i = 1, i<=Length[currentConfig], i++, newConfig[[i]]=currentConfig[[i]]];

      
      While[newConfig == currentConfig, Module[{},
      
      inTermsOfPlast = FindAllPowers[ending, currentConfig];
      completionTimes = FindCompletionTimes[inTermsOfPlast];


          runningX = {};
          jobSpeeds = {};
          jobEnergies = {};


          (*CALCULATE RUNNING TIMES Xi's*)
          For[mm = Length[currentConfig], mm > 0, mm--, runningX =Append[runningX, -1]];
          For[mm = 1, mm <= Length[inTermsOfPlast],mm++, runningX[[mm]] =(inTermsOfPlast[[mm]]^(-1.0/alpha))]; 
          (*CALCULATE speeds Si's*)
          For[i = Length[currentConfig], i > 0, i--, jobSpeeds = Append[jobSpeeds,-1]];
          For[i = 1, i <= Length[runningX],i++, jobSpeeds[[i]] =(1.0/runningX[[i]])];  

          (*CALCULATE job energies Ei's*)
          For[i = Length[currentConfig], i > 0, i--, jobEnergies = Append[jobEnergies,-1]];
          For[i = 1, i <= Length[jobSpeeds],i++, jobEnergies[[i]] =(jobSpeeds[[i]]^(alpha-1.0))];  

        enAll=CalculateEnergyUsed[ending, currentConfig]*1.0;

        Round[enAll*1000.00]/1000.00>>>C:\output.txt;
        currentConfig>>>C:\output.txt;
        Round[inTermsOfPlast*1000.00]/1000.00>>>C:\output.txt;
        Round[completionTimes*1000.00]/1000.00>>>C:\output.txt;
        Round[runningX*1000.00]/1000.00>>>C:\output.txt;
        Round[jobSpeeds*1000.00]/1000.00>>>C:\output.txt;
        Round[jobEnergies*1000.00]/1000.00>>>C:\output.txt;
        
        
          
          newConfig =  FindNewConfiguration[completionTimes, JobReleases,inTermsOfPlast, currentConfig];
          ending = ending - howClose;
          If[ending <= 0, Return[0];];
          ]
         ];
    
    Return[ending];
   ] 






LowerEnvelope[] := Module[{},
    
    (*Always add the end point to start with, so it will have n + 1 points*)
    intervalsOfOptimalConfigurations = 
      Append[intervalsOfOptimalConfigurations, optPoint];
    allOptimalConfigurations = 
      Append[allOptimalConfigurations, currentConfig ];
    For[i = Length[currentConfig], i > 0, i--, 
      newConfig = Append[newConfig, -1]];
    
    If[lowerEnergy<0.0001,lowerEnergy=0.0001];
    While[optPoint > lowerEnergy, Module[{},
    
    allEnergiesForOptimalConfigurations = Append[allEnergiesForOptimalConfigurations,FindAllPowers[optPoint, currentConfig]]; 
    allEnergiesUsed = Append[allEnergiesUsed,CalculateEnergyUsed[optPoint*1.0, currentConfig]]; 
    
    optPoint = FindOptimalChangePoint[optPoint, currentConfig]*1.0;
            
            
            
            
            currentConfig  = newConfig;
            
            allOptimalConfigurations = 
              Append[allOptimalConfigurations, currentConfig ];
            
            intervalsOfOptimalConfigurations = 
              Append[intervalsOfOptimalConfigurations, optPoint];
            ];
        ];
      
      intervalsOfOptimalConfigurations[[Length[intervalsOfOptimalConfigurations]]] = lowerEnergy*1.0;
    
    
    For[k = 1, k < Length[intervalsOfOptimalConfigurations], k++,
      globalListOfPlots = 
          Append[globalListOfPlots, 
            ParametricPlot[{CalculateEnergyUsed[t, 
                  allOptimalConfigurations[[k]]], 
                CalculateFlowTime[t, allOptimalConfigurations[[k]]]}, {t, 
                intervalsOfOptimalConfigurations[[k]], 
                intervalsOfOptimalConfigurations[[k + 1]]}, 
              PlotStyle -> {colorList[[Mod[k, 2] + 1] ]}, 
              DisplayFunction -> Identity]
            ];
      ]
    
    ]



ClearAllGlobalVariables[] := Module[{},
        newConfig = {};
             allOptimalConfigurations = {};
        intervalsOfOptimalConfigurations = {};
     k = 1;
     globalListOfPlots = {};
    
    ]

(*Up to 20 jobs is supported*)
GlobalListOfJobs = {P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12, P13, P14, P15, P16, P17, P18, P19, P20};
colorList = {RGBColor[1, 0, 0], RGBColor[0, 0, 1]};
newConfig = {};
allOptimalConfigurations = {};
intervalsOfOptimalConfigurations = {};
allEnergiesForOptimalConfigurations = {};
allEnergiesUsed = {};
k = 1;
globalListOfPlots = {};
"start">>C:\output.txt;
EndPackage[ ]


(*Save in C:\Program Files\Wolfram Research\Mathematica\4.1\DrawLowerEnvelope.m*)