(* Modules for Nonlinear Unconstrained Optimization *)
(* by Zafer Gurdal *)
(* Copyright 1995  *)
(* Department of Engineering Science and Mechanics *)
(* Virginia Polytechnic Institute and State University *)
(* Blacksburg, VA 24061-0219 *)
(* 703-231-5905 *)
(* e-mail: gurdal@gurdal.esm.vt.edu (NeXT Mail Accepted) *) 
 
BeginPackage["UnconstrainedMin`"]

UnconstrainedMin::usage = 
 "Following techniques are available in the UnconstrainedMin 
  package along with associated utilities.  
 
SequentialSimplex: initialSimplex, simplexCycle, simplexPlot;  
  ConjugateDirections: PowellCycle;  
  SteepestDescent: GradientVector, NormalizedDirection; 
  ConjugateGradient: 
  NewtonsMethod: HessianMatrix;

  There are also plotting utilities, such as ProgressPlot that 
  shows the progress of the optimization"

SolutionHistory::usage = "SolutionHistory is an option of 
	optimizations which specifies whether the 
	entire set of cycles found will be returned in a List. 
	The default (False) is to return only the final solution."
	
PlotHistory::usage = "PlotHistory is an option 
 of optimizations which provides the nodes of the simplex
 for every cycle for plotting purposes. 
 The default (False) no information for plotting is returned."	
	
Epsilon::usage = "Epsilon is an option of optimizations which 
 controls the convergence criterion.  The default value of Epsilon
 is Automatic which corresponds to eps = 1.0 10^-9.  
 Use Epsilon -> eps, to change it."

OneDMin::usage = "OneDMin[f, x, f0x0, dir] performs one dimensional 
 minimization of a function f(x) starting from x0 (stores in f0x0
 along with the value of the function f0 in the form f0x0 = 
 {f[x0], x0} )."
 
AnimatePlot::usage = "AnimatePlot is an option of the ProgressPlot
 and SimplexProgressPlot.  If the option is True, the progress of 
 the design variables in the design space is animated as a function 
 of the iterations. If it is False, then the entire progress plot
 is generated as a single plot."

ProgressPlot::usage = "ProgressPlot[fcont,data] adds the plot 
 of solution point data onto the plot of the objective 
 function provided in fcont. For the Sequential Simplex algorithm
 use SimplexProgressPlot, instead of the ProgressPlot to visualize 
 the progress of the iterations.  
 Option AnimatePlot->True activates the generation of plots for 
 each iteration."

Options[ProgressPlot] = {AnimatePlot -> False}

(* SEQUENTIAL SIMPLEX *)
		
SequentialSimplex::usage = "SequentialSimplex[f,x,x0,size]
 Implements Reflection, Expansion, Contraction with appropriate 
 coefficients alpha, beta, and gamma, respectively, to find a
 minimum of the function f which is a function of x starting
 from x0. Type Options[SequentialSimplex] to see various options 
 of the algorithm"

initialSimplex::usage = "initialSimplex[x0,simSize] 
	Creates the coordinates of the remaining vertices of a 
	simplex of size simSize in n dimensional space starting from 
	the vertex x0." 
 
SimplexCycle::usage = "SimplexCycle[f,fxcom,alpha,beta,gamma] 
	Performs a one cycle simplex resulting in either reflection,
	or expansion, or contraction." 
 
simplexPlot::usage = "simplexPlot[data] prepares the graphics 
	of a two dimensional simplex for plotting." 

SimplexProgressPlot::usage = "SimplexProgressPlot[fcont,data] 
 adds the plot of solution point data onto the plot of the objective 
 function provided in fcont. For the Sequential Simplex algorithm
 use SimplexProgressPlot, instead of the ProgressPlot to visualize 
 the progress of the iterations.  
 Option AnimatePlot->True activates the generation of plots for 
 each iteration."

ReflectioN::usage = "Option of Sequential Simplex 
	algorithm which specifies the coefficient of reflecion. 
        Default value is alpha = 1.0 " 

ExpansioN::usage = "Option of Sequential Simplex 
	algorithm which specifies the coefficient of expansion.
        Default value is beta = 2.0 " 

ContractioN::usage = "Option of Sequential Simplex 
	algorithm which specifies the coefficient of contraction.
        Default value is beta = 0.5 " 

Options[SequentialSimplex] = 
	{SolutionHistory ->False, 
	 PlotHistory->False, 
	 Epsilon->Automatic, 
         MaxIterations -> 50,  
         ReflectioN -> 1.0, 
         ExpansioN -> 2.0, 
         ContractioN -> 0.5 }

(* POWEL'S CONJUGATE DIRECTIONS *)

ConjugateDirections::usage = "ConjugateDirections[f,x,x0]
 Minimizes the function f(x), starting from x0 using Powell's 
 method of Conjugate Directions.
 Type Options[ConjugateDirections] to see various options 
 of the algorithm."
 
PowellCycle::usage = "PowellCycle[f, x, f0x0, dirs]
 One cycle of Powell's conjugate directions iterations including 
 the pattern search."	
	
Options[ConjugateDirections] = 
	{SolutionHistory -> False, 
	 PlotHistory -> False, 
	 Epsilon -> Automatic, 
	 MaxIterations -> 50 }

(* STEEPEST DESCENT  *)

GradientVector::usage = "GradientVector[func, vars, pt]
 evaluates the gradient of the function func with respect to the 
 variables vars at point pt. It the point pt is not provided, 
 the function will return just the gradient vector."
 
NormalizedDirection::usage = "NormalizedDirection[ direct]
 normalizes the direction vector by its Euclidean norm and 
 returns the direction and the norm."

SteepestDescent::usage = "SteepestDescent[f,x,x0]
 Minimizes the function f of variables x starting from x0 using
 steepest descent directions. Type Options[SteepestDescent]
 to see various options of the algorithm."

Options[SteepestDescent] = 
	{SolutionHistory -> False, 
	 PlotHistory -> False, 
	 Epsilon -> Automatic, 
	 MaxIterations -> 50 }

(* FLETCHER-REEVES CONJUGATE GRADIENT DIRECTIONS  *)

ConjugateGradient::usage = "ConjugateGradient[f,x,x0]
 Minimizes the function f of variables x starting from x0 using
 steepest descent directions. Type Options[ConjugateGradient]
 to see various options of the algorithm."

Options[ConjugateGradient] = 
	{SolutionHistory -> False, 
	 PlotHistory -> False, 
	 Epsilon -> Automatic, 
	 MaxIterations -> 50 }

(* NEWTONS METHOD *)
	 
HessianMatrix::usage = "HessianMatrix[func, vars, pt]
 evaluates the Hessian matrix of the function func with respect 
 to the variables vars at point pt. If the point pt is not 
 provided, it will return the Hessian matrix."

NewtonsMethod::usage = "NewtonsMethod[f,x,x0]
 Minimizes the function f of variables x starting from x0 using
 Newton's method with exact line searches. Type  
 Options[NewtonsMethod]
 to see various options of the algorithm."

Options[NewtonsMethod] = 
	{SolutionHistory -> False, 
	 PlotHistory -> False, 
	 Epsilon -> Automatic, 
	 MaxIterations -> 50 }

(* QUASI NEWTON BFGS METHOD *)

QuasiNewtonBFGS::usage = "QuasiNewtonBFGS[f,x,x0]
 Minimizes the function f of variables x (or expression f(x))
 starting from x0 using a Quasi Newton method which implements
 BFGS update for the Hessian inverse. Type Options[NewtonsMethod]
 to see various options of the algorithm."

Options[QuasiNewtonBFGS] = 
	{SolutionHistory -> False, 
	 PlotHistory -> False, 
	 Epsilon -> Automatic, 
	 MaxIterations -> 50 }


(*  +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- *)
Begin["`Private`"]

ProgressPlot[hist_,data_,opts___?OptionQ] := 
   If[ MatrixQ[ data[[1]] ], 
       ProgressSimplex[hist, data, opts], 
       ProgressOthers[hist, data, opts] 
     ]
               
ProgressOthers[hist_,data_,opts___?OptionQ] := 
   Module[{finplot,int={ }, history = {hist} }, 
(*   Module[{finplot,int={data[[1]]}, history = {hist} }, *)
   finplot = AnimatePlot/.{opts}/.Options[ProgressPlot];
   If[ finplot, 
       FoldList[Show[#1,Graphics[{PointSize[0.02],Point[#2]}], 
                        Graphics[{Thickness[0.007],RGBColor[1,0,0],
                                      Line[ int=Append[int,#2] ]} 
                                ] 
                    ]&, history, data 
               ] ,
       Show[ history, Graphics[{PointSize[0.02],Map[Point[#]&,data]
}], 
                       Graphics[{Thickness[0.007],RGBColor[1,0,0],
                                      Line[data]} ]
           ] 
      ]
         ]

ProgressSimplex[hist_,data_,opts___?OptionQ] := 
   Module[{finplot, history = {hist} }, 
   finplot = AnimatePlot/.{opts}/.Options[ProgressPlot];
   If[ finplot,
       Map[(history = Append[history, simplexPlot[#] ];
            Show[history, DisplayFunction->$DisplayFunction])&,data
          ],
       Map[(history = Append[history, simplexPlot[#] ])&,data];
       Show[history, DisplayFunction->$DisplayFunction]
     ]
         ]

initialSimplex[x0_,asize_] := 
  Module[ {n=Length[x0], p, q, idt, mult},
    p    = asize / Sqrt[2.0] / n (Sqrt[n+1.0]+n-1) ;
    q    = asize / Sqrt[2.0] / n (Sqrt[n+1.0]-1) ;
    idt = IdentityMatrix[ n ];
    mult = Fold[ Plus, x0, q idt ]; 
    Insert[ Map[mult + p # - q # &, idt], x0, 1]
    ]  
    
simplexPlot[data_] := 
   Module[{int=data}, 
   Graphics[ Flatten[{ Line[Append[int,int[[1]] ] ],
             Map[{PointSize[0.02], Point[#]}&, int] } ,1]
           ]
         ]
         

(*----------------------------------------------------------------*)
       
SequentialSimplex[f_Symbol,x_,x0_,size_,opts___?OptionQ] := 
                              SequentialSimplex[f[x],x,x0,size,opts]

SequentialSimplex[fexp_,x_,x0_,size_,opts___?OptionQ]:= 
   Module[{init,fNo,fxbar,test,avSmpx,slnH,pltH,epep,epsCheck,MaxIt,
           counter=0, 
           slnHistory={}, pltHistory={}, epsPrecision = 1.0 10^-9},
     {maxIt, slnH, pltH, epep, alpha, beta, gamma} = 
     {MaxIterations, SolutionHistory, PlotHistory, Epsilon, 
                     ReflectioN,      ExpansioN,   ContractioN} /. 
                                {opts}/.Options[SequentialSimplex];
     init = initialSimplex[x0,size]; 
     fNo  = Map[fexp/.Thread[x->#]&, init]; 
     init = MapIndexed[Insert[init[[#2]],#,1]&,fNo];
     init = Sort[init, (Part[#2,1] < Part[#1,1])&];
     avSmpx = Level[init,1,Plus][[2]] / Length[init];
     fxbar= fexp/.Thread[x->avSmpx];
     If[pltH, 
              pltHistory=Append[pltHistory,Map[#[[2]]&,init]], 
              slnHistory=Append[slnHistory, {fxbar, avSmpx}]
             ]; 
     test = Sqrt[Apply[Plus, (Map[ #[[1]]&, init] - fxbar)^2] / 
            (Length[x0]+1)]; 
     If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];  
     While[test > epsCheck, 
           init = SimplexCycle[fexp,x,init,alpha,beta,gamma]; 
           avSmpx = Level[init,1,Plus][[2]] / Length[init];
           fxbar=fexp/.Thread[x->avSmpx]; 
           test =Sqrt[Apply[Plus,(Map[ #[[1]]&,init]-fxbar)^2]/ 
                   (Length[x0]+1) ];
           If[pltH, 
              pltHistory=Append[pltHistory,Map[#[[2]]&,init]], 
              slnHistory=Append[slnHistory, {fxbar, avSmpx}]
             ];
           If[++counter >= maxIt, (test = 0.0; Print[ "Maximum 
	Number of Iterations ", maxIt, " Exceeded "]), Null ]
          ];
     If[pltH, pltHistory, If[slnH, slnHistory, {fxbar, avSmpx}] ]
 ]
      
SimplexCycle[fun_, x_, comb_, alpha_, beta_, gamma_] :=
   Module[{fxcom=comb,dropped,xbar,xref,xexp,xcon,xshnk,
                                   fref,fexp,fcon,fshnk}, 
(* drop the node with the highest objective,fxcom is sorted *)
     dropped = Drop[fxcom, {1}];
     xbar = Level[dropped, 1, Plus][[2]] / Length[ dropped ]; 
(*               Every Cycle Starts with REFLECTION         *)
     xref = xbar + alpha ( xbar - fxcom[[1,2]] ) ; 
     fref = fun /. Thread[x->xref];
     If[ fref < Last[fxcom][[1]], 
(*       IF YES   do  EXPANSION                             *) 
(* isolate the coordinates of the simplex for shrinkage      *) 
         xexp = xbar + beta ( xref - xbar );
         fexp = fun/.Thread[x->xexp];
         If[ fexp < fref, 
(*           IF YES (following expansion) replace Xh by Xe   *) 
             fxcom = Append[ dropped, {fexp,xexp}], 
(*           IF NO  (following expansion) replace Xh by Xr   *)
             fxcom = Append[ dropped, {fref,xref}] 
           ], 
(*       IF NO    check if F(r) <= F(s)                     *) 
         If[ fref <= fxcom[[2,1]], 
(*           IF YES   replace Xh by Xr                      *) 
             fxcom = Append[ dropped, {fref,xref}], 
(*           IF NO  check if F(r) < F(h)                    *)
             If[ fref < fxcom[[1,1]],
(*               IF YES   replace Xh by Xr                  *)
                 fxcom = Prepend[ dropped, {fref,xref}], Null
                 {(* IF NO do nothing *)}
               ]; 
(*        do  CONTRACTION  (no need to check anything else  *)
             xcon = xbar + gamma ( fxcom[[1,2]] - xbar);
             fcon = fun/.Thread[x->xcon];
             If[ fcon > fxcom[[1,1]], 
(*               IF YES   Apply Shrinkage to Replace all Xi *) 
(* isolate the coordinates of the simplex for shrinkage     *)
                 xshnk = Map[ Level[#,{2}]&, fxcom]; 
(* apply shrinkage operation and compute new objectives     *)
                 xshnk = Map[# + 0.5(Last[xshnk]-#)&, xshnk];
                 fshnk = Map[fun/.Thread[x->#]&, xshnk];
(* combine new objectives with xshnk coordinates and sort   *) 
     fxcom = MapIndexed[Insert[xshnk[[#2]],#,1]&,fshnk];
     fxcom = Sort[fxcom, (Part[#2,1] < Part[#1,1])&], 
(*               IF NO    Replace Xh by Xc                  *)
                 fxcom = Prepend[ dropped, {fcon,xcon}]
               ]
           ]
       ]; 
(* do a final sort before returning back to the iterations  *)
     fxcom = Sort[fxcom, (Part[#2,1] < Part[#1,1])&] 
        ]

(*----------------------------------------------------------------*)
OneDMin[f_, x_, f0x0_, dir_] := 
   Module[ {inter, alf},
     inter = FindMinimum[f/.Thread[x->f0x0[[2]]+alf dir], {alf, 0} ];  
     {inter[[1]], f0x0[[2]] + alf dir /. inter[[2]]}
         ]
(*----------------------------------------------------------------*)

ConjugateDirections[f_Symbol,x_,x0_,opts___?OptionQ] := 
                              ConjugateDirections[f[x],x,x0,opts]
   
ConjugateDirections[fexp_,x_,x0_,opts___?OptionQ]:= 
   Module[{dirctns,desPt,slnH,pltH,epep,sol,prefval,epsCheck,maxIt,
           slnHistory={}, pltHistory={}, counter = 0, 
           test =$MaxMachineNumber, epsPrecision = 1.0 10^-9},
     {maxIt, slnH, pltH, epep} = 
     {MaxIterations, SolutionHistory, PlotHistory, Epsilon} /.
                          {opts}/.Options[ConjugateDirections];
     dirctns = IdentityMatrix[ Length[x0] ];
     desPt = {fexp/.Thread[x->x0], x0};
     If[pltH, 
              pltHistory=Append[pltHistory, x0], 
              slnHistory=Append[slnHistory, desPt]
             ];
     If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];  
     While[test > epsCheck,
           prefval = desPt[[1]]; 
           sol = PowellCycle[fexp, x, desPt, dirctns]; 
           dirctns = sol[[2]];
           desPt = {fexp/.Thread[x->sol[[1]] ], sol[[1]] };
           If[pltH, 
              pltHistory=Append[pltHistory, sol[[1]] ], 
              slnHistory=Append[slnHistory, desPt]
             ];
           test =  Abs[desPt[[1]] - prefval ];
           If[++counter >= maxIt, (test = 0.0; Print[ "Maximum 
	Number of Iterations ", maxIt, " Exceeded "]), Null ]            
          ];  
     If[pltH, pltHistory, If[slnH, slnHistory, desPt] ]
     ]
      
PowellCycle[fexp_, x_, f0x0_, dirs_] :=
   Module[{unicyc,fobjs,pos,patdir,alf,alfRHS,patsol,patpt,ndirs}, 
(* do search along each of the directions starting from {f(x0),x0} *)
     unicyc = FoldList[ OneDMin[fexp,x,#1,#2]&, f0x0, dirs]; 
     fobjs = Drop[ Map[ #[[1]]&, unicyc ], {1}]; 
     fobjs = FoldList[ Abs[#2 - #1]&, 0, fobjs ]; 
     pos = Flatten[Position[ fobjs, Max[fobjs] ]][[1]]; 
     patdir = Last[unicyc][[2]] - First[unicyc][[2]]; 
     patsol=FindMinimum[fexp/.Thread[x->First[unicyc][[2]]+alf*patdir], 
                        {alf, 0}]; 
     alf = alf /. patsol[[2]] ;
     patpt = {patsol[[1]], First[unicyc][[2]]+alf patdir};   
     alfRHS = Sqrt[Abs[(First[unicyc][[1]] - patpt[[1]]) /
                 (unicyc[[pos-1]][[1]] - unicyc[[pos]][[1]])] ]; 
     If[ Abs[alf] < alfRHS, 
         {patpt[[2]],dirs}, 
         (ndirs=ReplacePart[dirs,patdir,pos-1]; {patpt[[2]],ndirs}) ]
   ]

(*----------------------------------------------------------------*)	
GradientVector[f_, X_List] := (D[f,#]&) /@ X

GradientVector[f_, X_List, Xpt_List]:= 
			               (D[f,#]&) /@ X /. Thread[X->Xpt]
	
NormalizedDirection[ direct_] := 
   Module[{int}, 
      int = Sqrt[ Apply[ Plus, Map[ #^2 &, direct ]] ];
      {direct / int, int}
         ]
(*----------------------------------------------------------------*)

SteepestDescent[f_Symbol,x_,x0_,opts___?OptionQ] := 
                              SteepestDescent[f[x],x,x0,opts]

SteepestDescent[fexp_,x_,x0_,opts___?OptionQ] := 
   Module[{dirctn,desPt,prev,prevnorm,maxIt,epsCheck,
           slnH,pltH,epep,test =$MaxMachineNumber,counter=0,
           slnHistory={},pltHistory={},epsPrecision=1.0 10^-9}, 
     {maxIt, slnH, pltH, epep} = 
     {MaxIterations, SolutionHistory, PlotHistory, Epsilon} /.
                          {opts}/.Options[SteepestDescent];
     dirctn = GradientVector[fexp,x,x0];
     prevnorm = NormalizedDirection[dirctn]//N;
     dirctn = prevnorm[[1]];
     desPt = {fexp/.Thread[x->x0], x0};
     If[pltH, 
             pltHistory=Append[pltHistory, x0], 
             slnHistory=Append[slnHistory, desPt]
       ];
     If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];  
     While[test > epsCheck,
           prev = prevnorm[[2]]; 
           desPt = OneDMin[fexp, x, desPt, dirctn]; 
           If[pltH, 
              pltHistory=Append[pltHistory, desPt[[2]] ], 
              slnHistory=Append[slnHistory, desPt]
             ];
           dirctn = GradientVector[fexp,x,desPt[[2]] ];
           prevnorm = NormalizedDirection[dirctn]//N;
           dirctn = prevnorm[[1]];
           test =  Abs[prevnorm[[2]] - prev ];  
     If[++counter >= maxIt, (test = 0.0; Print[ "Maximum 
	Number of Iterations ", maxIt, " Exceeded "]), Null ]         
          ];  
     If[pltH, pltHistory, If[slnH, slnHistory, desPt] ]
     ]

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

ConjugateGradient[f_Symbol,x_,x0_,opts___?OptionQ] := 
                              ConjugateGradient[f[x],x,x0,opts]

ConjugateGradient[fexp_,x_,x0_,opts___?OptionQ]:= 
   Module[{dirctn,dirctn1,desPt,beta,maxIt, epsCheck,
           slnH,pltH,epep,counter = 0,test =$MaxMachineNumber,  
           slnHistory={},pltHistory={},epsPrecision=1.0 10^-9}, 
     {maxIt, slnH, pltH, epep} = 
     {MaxIterations, SolutionHistory, PlotHistory, Epsilon} /.
                          {opts}/.Options[ConjugateGradient];
     dirctn =  - GradientVector[fexp,x,x0];
     desPt = {fexp/.Thread[x->x0], x0};
     If[pltH, 
              pltHistory=Append[pltHistory, x0], 
              slnHistory=Append[slnHistory, desPt]
             ];
     If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];  
     While[test > epsCheck, 
           desPt = OneDMin[fexp, x, desPt, dirctn]; 
           If[pltH, 
              pltHistory=Append[pltHistory, desPt[[2]] ], 
              slnHistory=Append[slnHistory, desPt]
             ];
           dirctn1 = - GradientVector[fexp,x,desPt[[2]] ];
           beta =  (dirctn1 . dirctn1) / (dirctn . dirctn); 
           dirctn = dirctn1 + beta dirctn;
           test = Sqrt[Apply[Plus,Map[#^2 &,dirctn] ] ];
     If[++counter >= maxIt, (test = 0.0; Print[ "Maximum 
	Number of Iterations ", maxIt, " Exceeded "]), Null ];         
          ];  
     If[pltH, pltHistory, If[slnH, slnHistory, desPt] ]
     ]

(*----------------------------------------------------------------*)
HessianMatrix[f_, X_List] :=
				    (D[(D[f,#]&) /@ X ,#]&) /@ X 
				    
HessianMatrix[f_, X_List, Xpt_List] :=
			     (D[(D[f,#]&)/@ X ,#]&)/@ X /. Thread[X->Xpt] 
(*----------------------------------------------------------------*)

NewtonsMethod[f_Symbol,x_,x0_,opts___?OptionQ] := 
                              NewtonsMethod[f[x],x,x0,opts]

NewtonsMethod[fexp_,x_,x0_,opts___?OptionQ]:= 
   Module[{hessinv,dirctn,desPt,maxIt,epsCheck, 
           slnH,pltH,epep,test =$MaxMachineNumber, counter=0,
           slnHistory={},pltHistory={},epsPrecision=1.0 10^-9}, 
     {maxIt, slnH, pltH, epep} = 
     {MaxIterations, SolutionHistory, PlotHistory, Epsilon} /.
                          {opts}/.Options[NewtonsMethod];
     hessinv = Inverse[ HessianMatrix[ fexp, x, x0 ] ]; 
     dirctn = - hessinv . GradientVector[fexp,x,x0];     
     desPt = {fexp/.Thread[x->x0], x0};
     If[pltH, 
              pltHistory=Append[pltHistory, x0], 
              slnHistory=Append[slnHistory, desPt]
             ];
     If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];  
     While[test > epsCheck, 
           desPt = OneDMin[fexp, x, desPt, dirctn]; 
           If[pltH, 
              pltHistory=Append[pltHistory, desPt[[2]] ], 
              slnHistory=Append[slnHistory, desPt]
             ];
           hessinv = Inverse[ HessianMatrix[ fexp, x, desPt[[2]] ] ];
           dirctn= -hessinv . GradientVector[fexp,x,desPt[[2]] ];
           test =  Sqrt[Apply[Plus,Map[#^2 &,dirctn] ] ];
     If[ ++counter >= maxIt, (test = 0.0; Print[ "Maximum 
	Number of Iterations ", maxIt, " Exceeded "]), Null ];         
          ];  
     If[pltH, pltHistory, If[slnH, slnHistory, desPt] ]
     ]

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

QuasiNewtonBFGS[f_Symbol,x_,x0_,opts___?OptionQ] := 
                              QuasiNewtonBFGS[f[x],x,x0,opts]
                              
QuasiNewtonBFGS[fexp_,x_,x0_,opts___?OptionQ]:= 
   Module[{dirctn,desPt,slnH,pltH,epep,epsCheck,maxIt,
           slnHistory={}, pltHistory={}, counter = 0, 
           test =$MaxMachineNumber, epsPrecision = 1.0 10^-9,
           IDNT,BFGS,gradf,gradf1,desPt1,p0,y0,int,cnt},
     {maxIt, slnH, pltH, epep} = 
     {MaxIterations, SolutionHistory, PlotHistory, Epsilon} /.
                          {opts}/.Options[QuasiNewtonBFGS];
     IDNT = IdentityMatrix[ Length[x0] ];
     BFGS = IDNT ; 
     desPt = {fexp /. Thread[x->x0], x0};
     gradf = GradientVector[ fexp, x, x0 ] ; 
     If[pltH, 
              pltHistory=Append[pltHistory, x0], 
              slnHistory=Append[slnHistory, desPt]
             ];
     If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];  
     While[test > epsCheck,
           dirctn = -BFGS . gradf ;
           desPt1 = OneDMin[fexp, x, desPt, dirctn] ;
	   If[desPt1 == desPt, (Print[" Last  1-D  Minimization 
			Yilelds  No  Improvement"]; Break[]), Null ]; 
           gradf1 = GradientVector[ fexp, x, desPt1[[2]] ];
           If[pltH, 
              pltHistory=Append[pltHistory, desPt1[[2]] ], 
              slnHistory=Append[slnHistory, desPt1]
             ];
           p0 = desPt1[[2]] - desPt[[2]] ;
           y0 = gradf1 - gradf ; 
           cnt = p0 . y0 ; 
           int = IDNT - Outer[ Times, p0, y0 ]/cnt ;
           BFGS= int .BFGS. Transpose[int] + Outer[Times,p0,p0]/cnt;
           gradf = gradf1 ;  
           desPt = desPt1;
           test =  Sqrt[Apply[Plus,Map[#^2 &,gradf] ] ];
           If[++counter >= maxIt, (test = 0.0; Print[ "Maximum 
		Number of Iterations ", maxIt, " Exceeded "]), Null ]            
          ];  
     If[pltH, pltHistory, If[slnH, slnHistory, desPt] ]
     ]
(*----------------------------------------------------------------*)  

End[]
(* Protect[ ] *)
(* SetAttributes[{ }, Listable] *)

EndPackage[]