(* ::Package:: *)

(* ::Chapter:: *)
(*prolog*)


StationaryPoints::usage = "Computes the stationary points of a function f(x)."


Begin["`Private`"]


(* ::Chapter:: *)
(*main code*)


(* ::Section::Closed:: *)
(*CriticalAndStationaryPointScanner*)


(*****************************************
*        STATIONARY POINTS
******************************************)
Options[CalculateStationaryPoints] = {"SpecificPoint" -> None, "Domain" -> Automatic(*for calling from IncreasingDecreasing and passing to ? no deletecases from getting here*)};

CriticalAndStationaryPointScanner[str_, Hold[CalculateStationaryPoints[expr_, OptionsPattern[]]], _, op:OptionsPattern[]] := Module[
  	{input, xys, res, domain, capt, sp},
	If[ContainsQ[expr, CalculateData|Quantity|Unit], Return[]];

  	input = preprocessUserInput[expr];
  	xys = chooseVariables[expr, input];
  	
  	(* fighting off bad parse: "stationary points of sinx -30<x<30" *)
	If[ContainsQ[input, $Failed], Return[]];
	
	domain = OptionValue[CalculateStationaryPoints, "Domain"];
  	domain = formatdomain[domain, xys]; dPrint["domain:  ",domain];

  	res = StationaryPointFinder[input, domain, xys];
  	dPrint["stat. pts::  ",res];
	res = Which[
			res === $Failed, GrayComment[Localize["no stationary points found",4338]],
			res === {} && input =!= $Failed, GrayComment[ RowTemplate[Localize["`1` is a constant function",4339],{process[Hold[expr]]}] ],				
			True, res
		];
	domain = Which[MatchQ[domain, Automatic], {},
				   MatchQ[domain, {_, _}] && MatchQ[xys,{}], {{If[MatchQ[xys, {}], Symbol["x"], First@xys], #1, #2}}&@@domain, 
				   MatchQ[domain, _[_, First@xys, _]], {{#2, #1, #3}} & @@domain, 
				   MatchQ[domain, And[_[_, First@xys, _]..]], {{#2, #1, #3}} & @@@domain, 
				   True,  domain];
				   
	(* Don't overwrite the Result pod if this was called from IncreasingDecreasingScanner *)			   
	capt = If[MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"], Localize["Stationary points", 51813], $CaptionResult];

  	If[MatchQ[res, _GrayComment], 
    		SowResult[
     			capt,
     			res,
     			10,
     			"Buttons" -> {"ReplacePod", "Identity", Identity},
     			"PodIdentifier" -> capt,
     			"Dependencies" -> {"UpdateFrequency" -> "Weekly"}
     		];
     		If[TrueQ[NumericQ[OptionValue["SpecificPoint"]]] && !MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"](*TODO: remove once plotDiscontinuitiesWithAutoRange can be made to work how we want*),
				plotDiscontinuitiesWithAutoRange[input, OptionValue["SpecificPoint"], {}, xys[[1]], OptionValue["SpecificPoint"], op],
     			If[False && MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"],
     				createplot["ParametricPlot", input, domain, "IncDecPlot"->True(*, ImageSize->AlphaScaled[1]*)](*SowResult["Plot",LineIntegralPlot[{xys[[1]], input}, {xys[[1]], 0, 20}, D[input, xys[[1]]](*-y+z*), {xx, yy}], 99]*),
     				PlotExtrema[input, Automatic, {}, xys, "IncDecPlot"->MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"],op]
     			]
     		],
     		SowResult[
   				capt,
   				{res, input, xys, Hold[expr]},
   				10,
   				"Buttons" -> If[ContainsQ[res, Root], {"ReplacePod", "StationaryPoints", {MachinePrecision, 10}}, {"ReplacePod", "StationaryPoints", {Identity, 10}}], 
   				"PodIdentifier" -> capt,
   				"Dependencies" -> {"UpdateFrequency" -> "Weekly"}
   			];	(*Print["OptionValue[SpecificPoint]",OptionValue["SpecificPoint"],Cases[{op},a:("SpecificPoint"->_):>a,1]];*)
     		If[Length[xys] < 3, 
     			If[TrueQ[NumericQ[OptionValue["SpecificPoint"]]] && !MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"](*TODO: remove once plotDiscontinuitiesWithAutoRange can be made to work how we want*),
					plotDiscontinuitiesWithAutoRange[input, OptionValue["SpecificPoint"], {}, xys[[1]], OptionValue["SpecificPoint"], op],
     				If[False && MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"],
     					createplot["ParametricPlot", {First@xys,input}, domain, "IncDecPlot"->True(*, ImageSize->AlphaScaled[1]*)](*SowResult["Plot",LineIntegralPlot[{xys[[1]], input}, {xys[[1]], 0, 20}, D[input, xys[[1]]](*-y+z*), {xx, yy}], 99]*),
     					sp = If[#==={},"SpecificPoint"->None,#]&[Cases[{op}, a:("SpecificPoint"->b_) /; !StringQ[b] :> a, 1]];
     					PlotExtrema[input, Automatic, N[res[[All, {1, 2}]]], xys, "IncDecPlot"->MemberQ[OptionValue["PrimaryScanners"], "IncreasingDecreasingScanner"],sp]
     				]
     			]
     		];
     		(*SowResult["Annotated graph", #, 8888]&/@CalculateScan`Packages`AnnotatedGraph`AnnotatedGraph[input,xys]*)
   	]
  
]


(* ::Section::Closed:: *)
(*StationaryPointFinder*)


(*****************************************
*
*	CODE FOR FINDING STATIONARY POINTS
*
******************************************)
StationaryPointFinder[exp_, cons_, x_List] := Module[
	{
		lpt, rpt, expr, cand, dd, const, var=x[[1]]
	}, 
    If[VarList[exp] === {}, Return[{}]];
    expr = gentleComplexExpand[exp];

    const = If[MatchQ[cons, {_?usersymbolQ, _, _}], #2 <= #1 <= #3 & @@cons, cons];
 
    cand = Quiet[TimeConstrained[
      If[(MatchQ[const, HoldPattern[(Less|LessEqual)[_?NumericQ, _, _?NumericQ]]])(*True*), 
        (* then *)
  	  lpt=const[[1]]; 
  	  rpt=const[[-1]];
  	  dPrint["using NReduce!!"(*,nNReduce[D[exp, var] == 0, {var, lpt, rpt}]*)];
  	  {exp/.{var->#[[1]]}, {var -> #[[1]]}} & /@ ("IsolatedRoots" /. NReduce[D[exp, var] == 0, {var, lpt, rpt}]), 
  	  (* else *) 
  	  dPrint["not using NReduce!!"];
  	  LocalExtremaFinder[expr, const /. {Infinity->$MaxNumber, -Infinity->-$MaxNumber}, x]
      ], 3.4 $singlesteptimelimit(*7 - needed more for some, but 3,3.4 is limit in critptfinder[] *)
    ]];
Print["cand: ", cand];
    dPrint["LLLL::LL:L:L:",{{cand,x,cons}}, mMatchQ[const//FullForm, HoldPattern[(Less|LessEqual)[_?NumericQ, _, _?NumericQ]]], sNReduce[D[AiryAi[var],var] == 0, {var, -39.1, 3}](*lLocalExtremaFinder[expr, const/.{Infinity->$MaxNumber,-Infinity->-$MaxNumber}, x]*)];
    If[MatchQ[cand, $Failed | _LocalExtremaFinder | {} | $Aborted], Return[$Failed]];
    
    dd = realD[expr, {x}];
    cand = Quiet[Select[cand, VectorQ[Chop[dd /. #[[-1]] /. C[_] -> 0, 10^-6], NumericQ[#] && PossibleZeroQ[#] &]&]];
  
    (* delete duplicates is for the input "stationary points of sin[1/x]" *)
    If[cand === {}, 
  		$Failed,
  		DeleteDuplicates[
  			Map[{TimeConstrained[
					FullSimplify[#[[1]], Element[C[_], Integers]],
					1/15 $singlesteptimelimit,
					Simplify[#[[1]], Element[C[_], Integers]]
				  ], #[[-1]]} &, 
			cand]
  		]
    ]
]

realD[expr_, x_] := D[expr, x] /. Derivative[1][Abs][f_] :> f/Sign[f]


(* ::Input:: *)
(*(*In[635]:= *)StationaryPointFinder[Abs[Sqrt[x]]^(x^2 - x^3 + 1/Abs[x]), Automatic, {x}]*)
(*(*Out[635]= {{0, {x -> 0.00252567}}, {1.05099, {x -> 1.19593}}}*)*)


(* ::Input:: *)
(*(*In[636]:= *)StationaryPointFinder[Abs[1 - x^3] Abs[x]^x, Automatic, {x}]*)
(*(*Out[636]= {{2.403985, {x -> -1.6365550}}}*)*)


(* ::Input:: *)
(*(*In[637]:= *)StationaryPointFinder[Sin[x], -4 <= x <= 4, {x}]*)
(*(*Out[637]= {{-1, {x -> -(\[Pi]/2)}}, {1, {x -> \[Pi]/2}}}*)*)


(* ::Input:: *)
(*(*In[638]:= *)StationaryPointFinder[Sin[x], {x,-4,4}, {x}]*)
(*(*Out[638]= {{-1, {x -> -(\[Pi]/2)}}, {1, {x -> \[Pi]/2}}}*)*)


(* ::Input:: *)
(*(*In[639]:= *)StationaryPointFinder[(4 - x^2)^2, Automatic, {x}]*)
(*(*Out[639]= {{0, {x -> -2}}, {16, {x -> 0}}, {0, {x -> 2}}}*)*)


(* ::Input:: *)
(*(*In[640]:= *)StationaryPointFinder[Sin[x] Cos[x] - Tan[x], -10 <= x <= 10, {x}]*)
(*(*Out[640]= {{0, {x -> 0}}, {0, {x -> -3 \[Pi]}}, {0, {x -> -2 \[Pi]}}, {0, {x -> -\[Pi]}}, {0, {x -> \[Pi]}}, {0, {x -> 2 \[Pi]}}, {0, {x -> 3 \[Pi]}}}*)
(**)*)


(* ::Section::Closed:: *)
(*gentleComplexExpand*)


gentleComplexExpand[expr_] /; ContainsQ[expr, Abs] := ComplexExpand[expr, TargetFunctions -> {Re, Im}]
gentleComplexExpand[expr_] := expr 


(* ::Section::Closed:: *)
(*LocalExtremaFinder*)


Clear[LocalExtremaFinder];
Options[LocalExtremaFinder] = {"Type" -> Automatic};

LocalExtremaFinder[___] := $Failed

LocalExtremaFinder[_List, _, _List, OptionsPattern] := $Failed

LocalExtremaFinder[expr_, x_Symbol] := LocalExtremaFinder[expr, Automatic, {x}]

LocalExtremaFinder[expr_, Automatic, x_Symbol, opts___] := LocalExtremaFinder[expr, Automatic, {x}, opts]

LocalExtremaFinder[a_. Sinc[f_]^b_. + c_., Automatic, vars_List, OptionsPattern[]] := With[{type = OptionValue["Type"]},
	With[{inst = FindInstance[f== 0, vars]},

		{a+c, #}& /@ inst /; MatchQ[inst, {{__Rule}}]

	] /; MatchQ[type, Automatic|"Maximum"]
] /; VectorQ[{a, b, c}, NumericQ] && a > 0 && b > 0

LocalExtremaFinder[f_, Automatic, vars_, OptionsPattern[]] /; 
	!PolynomialQ[f, vars] && 
	TimeConstrained[
		TrueQ @ Refine[Element[f, Integers], Element[vars, Integers]], 
		0.15, 
		False
	]  := $Failed

(*
In[45]:= LocalExtremaFinder[Fibonacci[n], Automatic, {n}] // Timing
Out[45]= {0.000394, $Failed}
*)

LocalExtremaFinder[in_, cons_, vars_List, OptionsPattern[]] := Catch @ Module[
	{expr, type, solns, xpts, nsol, sensible, eqns, tryNsol = False}, 
Print["in LocalExtremaFinder...", cons];
  	(* univariate case *)
  	type = OptionValue["Type"];
  	
  	expr = gentleComplexExpand[in];

	Which[
		cons === Automatic && MatchQ[expr,HeavisideTheta[__]],
			Return[{}],
		cons === Automatic,
	  		eqns = Thread[D[expr, {vars}] == 0] /. Automatic -> True;
	  		If[FreeQ[eqns, Derivative],
  				Quiet[
   				solns = CalculateTimeConstrained[
   							Reduce[eqns, vars, Reals, Backsubstitution -> True],
   							2 $singlesteptimelimit
   						];
   				If[MatchQ[Head[solns],Reduce] && ContainsQ[eqns,HoldPattern[((*TriangleWave|*)SquareWave)[__]]], solns = False]
   				],
   				solns = $Failed
   			];
   			(*solns = False;*)
   			print["Reduce solution is :: ", solns, rReduce[eqns, vars, Reals, Backsubstitution -> True]],
   		MatchQ[cons, Abs[_Symbol] < _?NumericQ | (Less | LessEqual)[_?NumericQ, _Symbol, _?NumericQ]], (* and possibly some more... *)
	  		eqns = Thread[D[expr, {vars}] == 0] /. Automatic -> True;
  			Quiet[
  				If[FreeQ[eqns, Derivative],
   					solns = CalculateTimeConstrained[
   							Reduce[Flatten[{eqns, cons}], vars, Reals, Backsubstitution -> True],
   							2 $singlesteptimelimit
   						];
   					If[MatchQ[Head[solns],Reduce] && ContainsQ[eqns,HoldPattern[((*TriangleWave|*)SquareWave)[__]]], solns = False]
   					,
   					solns = $Failed
  				]
   			];
   			print["Reduce solution is :: ", solns],   			
   		True,
   			solns = False; Print["solns: ", solns]
	];
    
	If[solns === False && cons === Automatic, Return[$Failed]];
	
  	If[MatchQ[solns, {} | {{}} | $Failed | $Aborted | _Reduce | False], 
  		If[cons === Automatic && elementaryQ[expr, vars], 
  			print["going back with a interval range"]; 
  			Return[LocalExtremaFinder[expr, Last[Get1DRange[expr, vars[[1]]]]/.{{x_, l_, r_} :> l<=x<=r}, vars, "Type" -> type]]
  		];
  		tryNsol = True; 
  		nsol =  Quiet @ NLocalExtremaFinder[expr, cons, vars, "Type" -> type];
  		print["Extrema from NLocalExtremaFinder are :: ", nsol,nNLocalExtremaFinder[expr, cons, vars, "Type" -> type]];
  		If[nsol === $Failed,
  			Throw[$Failed],
  			solns = nsol[[All, -1]];
  			solns = {N[expr /. #], #}& /@ Select[solns, (expr /. #) \[Element] Reals &];
  			solns = DeleteDuplicates @ Cases[solns, r_ /; Count[N[r], n_ /; Abs[n] > 10^17, Infinity] == 0];
  			If[type =!= Automatic && cons === Automatic, 
  				solns = Select[solns, ExtremaTest[expr, #[[-1]]] === type&]
  			];
  			Throw[ solns ]
  		],
  		sensible = SensibleResults[solns, vars];
  		If[sensible === {}, Throw[$Failed]];
  		solns = Cases[#, _Equal, Infinity] & /@ sensible[[All, 1]];
  		solns = solns /. Equal -> Rule
	];

  	solns = Select[solns, Element[#[[All, -1]] /. C[_] -> 1, Reals] &];
  	solns = Select[solns, Length[#] == Length[vars]&];
  	solns = SortByPosition[#, First, vars]& /@ solns;
  	
  	If[solns === {} && Length[vars] != Length[VarList[eqns]],
  		Throw[$Failed]
  	];
  	
  	If[solns === {},
  		solns = FindInstance[eqns, vars, Reals];
  		If[!MatchQ[solns, HoldPattern @ {{Rule[__]..}..}], Throw[$Failed]]
  	];

  	xpts = Which[
  		type =!= Automatic && cons === Automatic,
  			If[Length[vars] === 1, 
  				Select[solns, NthDerivativeTest[expr, #] === type &],
  				Select[solns, SecondDerivativeTest[expr, #] === type &]
  			],
    	type =!= Automatic && cons =!= Automatic,
    		If[Length[vars] === 1, 
  				Select[solns, NthDerivativeTest[expr, #] === type &],
  				Select[solns, SecondDerivativeTest[expr, #] === type &]
  			],
    	True,
    		solns
    ];
    
  	If[xpts === {} && tryNsol === False && !PolynomialQ[expr, vars], 
  		nsol =  Quiet @ NLocalExtremaFinder[expr, cons, vars, "Type" -> type];
 		print["Extrema from NLocalExtremaFinder are :: ", nsol];
  		If[nsol === $Failed,
  			Throw[$Failed],
  			solns = nsol[[All, -1]];
  			solns = Quiet @ Check[ {N[expr /. #], #}& /@ Select[solns, (expr /. #) \[Element] Reals &], $Failed, General::unfl ];
  			If[ContainsQ[solns, $Failed], Throw[$Failed]];
  			solns = DeleteDuplicates @ Cases[solns, r_ /; Count[N[r], n_ /; Abs[n] > 10^17, Infinity] == 0];
  			solns = If[cons === Automatic && type =!= Automatic,
  						If[Length[vars] === 1, 
  							Select[solns, NthDerivativeTest[expr, #] === type &],
  							Select[solns, SecondDerivativeTest[expr, Last[#]] === type &]
  							],
  						solns
    				];
  			Throw[ solns ]  
  		]
  	];
  	xpts = Select[xpts, Simplify[(expr /. #) /. C[_] -> 1] \[Element] Reals &];
  	(* no RootReduce if not necessary, ex: stationary points x+x^2/2+x^4/4 *)
  	xpts = With[{simple = Simplify[expr /. #, Element[C[_], Integers]]}, Chop[{If[ContainsQ[simple, Root], RootReduce@simple, simple], #}]]& /@ xpts;
  	(* Adding N in Count below for example: LocalExtremaFinder[(x^5 + x^9 - x - 1)^3, Automatic, {x}] *)
  	DeleteDuplicates @ Cases[xpts, r_ /; Count[N[r], n_ /; Abs[n] > 10^17, Infinity] == 0]
]
  

(*
In[380]:= LocalExtremaFinder[-x^4 + 15 x^2 + 1 + -y^4 + 15 y^2 + 1, Automatic, {x, y}, "Type" -> "Minimum"]

Out[380]= {{2, {x -> 0, y -> 0}}}

In[381]:= LocalExtremaFinder[-x^4 + 15 x^2 + 1 + -y^4 + 15 y^2 + 1, Automatic, {x, y}, "Type" -> "Maximum"]

Out[381]= {{229/2, {x -> -Sqrt[(15/2)], y -> -Sqrt[(15/2)]}}, {2292, {x -> -Sqrt[(15/2)], y -> Sqrt[15/2]}}, 
{2292, {x -> Sqrt[15/2], y -> -Sqrt[(15/2)]}}, {2292, {x -> Sqrt[15/2], y -> Sqrt[15/2]}}}

In[382]:= LocalExtremaFinder[Sin[x + Sqrt[2] Sin[x^2]], -\[Pi] < x < \[Pi], {x},"Type" -> "Minimum"] // N

Out[382]= {{-1., {x -> -2.87307}}, {-1., {x -> -2.68057}}, {-1., {x-> -1.73858}}, 
{-0.176321, {x -> -0.356426}}, {0.40943, {x -> 1.35457}}, {0.672405, {x -> 2.13209}}, {-0.884908, {x -> 2.8248}}}

In[280]:= LocalExtremaFinder[Sin[x^3 y^3], Automatic, {x, y}]
Out[280]= {{0, {x -> 0, y -> 0}}}

In[376]:= LocalExtremaFinder[Exp[x] Sin[y], x^2 + y^2 == 1, {x, y}, "Type" -> "Maximum"]
Out[376]= {{1.32116, {x -> 0.673612, y -> 0.739085}}}

In[183]:= LocalExtremaFinder[(9 (4 x^2 Sin[x^2] - 2 Cos[x^2]))/5000, 0 < x < 6, {x}, "Type" -> "Maximum"]
Out[183]= {{0.0149327, {x -> 1.47466}}, {0.147422, {x -> 4.52697}}, {0.23775, {x -> 5.74736}}}

In[140]:= LocalExtremaFinder[Re[(1 + I x^2)/(1 - 3 I x)], Automatic, {x}, "Type" -> "Maximum"]
Out[140]= {{1., {x -> -7.36548*10^-9}}}

In[173]:= LocalExtremaFinder[Im[(1 + I x^2)/(1 - 3 I x)] + Re[(1 + I x^2)/(1 - 3 I x)], Automatic, {x}]
Out[173]= {{0.0383277, {x -> -0.509137}}, {1.21676, {x -> 0.141665}}}

In[54]:= LocalExtremaFinder[x y z, x + y + z == 100, {x, y, z}, "Type" -> "Maximum"]
Out[54]= {{37037., {x -> 33.3333, y -> 33.3333, z -> 33.3333}}}

In[160]:= LocalExtremaFinder[(4 E^-t^2 Erf[t])/Sqrt[\[Pi]], Automatic, {t}, "Type" -> "Maximum"]
Out[160]= {{0.951747, {t -> 0.620063}}}

In[161]:= LocalExtremaFinder[(4 E^-t^2 Erf[t])/Sqrt[\[Pi]], Automatic, {t}, "Type" -> "Minimum"]
Out[161]= {{-0.951747, {t -> -0.620063}}}
*)

elementaryQ[expr_, vars_] := Simplify`FunctionSingularities[expr, vars, "ELEM"] =!= $Failed["NONELEM"]

NLocalExtremaFinder[f_, const_, vars_List, "Type" -> "Minimum"] := Module[
	{findMinimumResult, otherFindMinimumResultPos, otherFindMinimumResultNeg, time, nMinimizeResult, 
		otherMinimizeResult, minima, others = {}},
  	findMinimumResult = Quiet @ If[const === Automatic,
  					Check[FindMinimum[f, ##], $Failed]& @@ Partition[vars, 1],
  					Check[FindMinimum[{f, const}, ##], $Failed]& @@ Partition[vars, 1]
  			];
  	
  	otherFindMinimumResultPos = If[const === Automatic,
  		Quiet @ Check[FindMinimum[f, Thread[{vars, 1 + RandomReal[]/100}]], $Failed],
  		Quiet @ Check[FindMinimum[{f, const}, Thread[{vars, 1 + RandomReal[]/100}]], $Failed]
  	];
  	
  	otherFindMinimumResultNeg = If[const === Automatic,
  		Quiet @ Check[FindMinimum[f, Thread[{vars, -1 - RandomReal[]/100}]], $Failed],
  		Quiet @ Check[FindMinimum[{f, const}, Thread[{vars, -1 - RandomReal[]/100}]], $Failed]
  	];
  	
  	print["FindMinimum results ", {findMinimumResult, otherFindMinimumResultPos, otherFindMinimumResultNeg}];

	{time, nMinimizeResult} = Timing @ Quiet @ Check[NMinimize[{f, const} /. Automatic -> True, vars], $Failed];
	(* Try NMinimize with another setting and hope it will find a different minimum. SamB 0410 *)
	otherMinimizeResult = If[time < 0.15,
						TimeConstrained[ 
							Quiet @ Check[
										NMinimize[{f, const} /. Automatic -> True, vars, 
											Method -> {"DifferentialEvolution", "CrossProbability" -> 0.125, "ScalingFactor" -> 2}],
										$Failed
									],
						$singlesteptimelimit/3,
						{}],
						{}
					];

	print["NMinimize results ", {nMinimizeResult, otherMinimizeResult}];

	If[Length[vars] == 1,
		others = extremaFromNumericalRootFinding[f, const, vars];
		others = Select[others, NthDerivativeTest[f, #] === "Minimum" &]
	];
	print["others", others];
	minima = Union[
			Cases[
				{findMinimumResult, otherFindMinimumResultPos, otherFindMinimumResultNeg, nMinimizeResult, otherMinimizeResult, others} /. {} -> Sequence[],
				{_, {__Rule}},
				Infinity
			],
			SameTest -> (SetPrecision[#1, 3] === SetPrecision[#2, 3]&)
		];

	minima = DeleteCases[minima, FindMinimum[__] | $Failed];
  	minima = DeleteCases[minima, {_?(!NumericQ[#] || !Abs[#] < $MaxMachineNumber &), {__Rule}}&];
  	If[MatchQ[minima, {{_, {__Rule}}..}] && FreeQ[minima, Indeterminate],
  		If[const === Automatic && Length[vars] === 1 && Length[minima] === 1 && EvenFnQ[f, vars],
   			Sort[ AppendTo[minima, {First[minima], {First[vars] -> -minima[[-1, 1, -1]]}}] ],
   			minima
  		],
   		$Failed
   	]
]


(* ::Chapter:: *)
(*epilog*)


End[]
