(* ::Package:: *)

BeginPackage["KnotTheory`OptimizedKnotPresentation`", {"KnotTheory`"}]


OptimizedKnotPresentation::about = "
	OptimizedKnotPresentation was written by Yun-chi Tang at the University of Toronto in summer 2021 under the guidance of Dr. Dror Bar-Natan.
"


OptimizedKnotPresentation::usage = "
	OptimizedKnotPresentation[pd] takes a knot in PD notation and through the greedy algorithm and local optimization techniques 
	determine a permutation of the knot that will be efficient in calculating knot invariants. The function will, from left to right, 
	at each time selecting a group of consecutive crossings to find an optimal order of depending on the knot invariant users intend to calculate. 
	OptimizedKnotPresentation[pd, options], where options is a list of rules, alows the user to control some of the parameters. 
	GreedyCriteria -> Function that help evaluate samples generated by the greedy algorithm to obtain the best greedy algorithm presentation. 
	InvariantCriteria -> Function sets the criteria, pertaining to the specific knot invariant the user intend to calculate, used 
	to conduct local optimization by changing the order of specific consecutive crossings, and Size -> n allows the user to change the 
	number of crossings that we intend to find an optimal order for every time we conduct local optimization.  
"


Begin["`Private`"]


greedyRep[pd_PD] := Module[{todo, front, x,v, len, rep},
todo = List @@ pd;
front = {};
len = 0;
v[x_X] :=Length[front  \[Intersection] (List @@ x)];
rep = {};
While[Length[todo] > 0,
x = RandomChoice[MaximalBy[todo, v]];
rep = Join[rep, {x}];
todo = DeleteCases[todo,x];
front = Complement[front \[Union] (List @@ x), front \[Intersection] (List @@ x)];
len = Max[len, Length[front]]
];
rep]


Protect[factor, tries];


Options[bestGreedy] = {factor -> maxWidthPresentationList, tries -> 1000};


bestGreedy[pd_PD, opts: OptionsPattern[bestGreedy]] := SortBy[Parallelize[Table[greedyRep[pd], {i, OptionValue[tries]}]], OptionValue[factor][#]&][[1]]


JonesResult[pdList_] := Module[{widths, JonesQuantity}, 
widths = Delete[Delete[widthPresentationList[pdList],1],-1];
JonesQuantity = Total[Range[Length[widths]] * widths * (CatalanNumber[#/2]&/@ widths)];
JonesQuantity]


maxWidthPresentationList[pdList_] := Max[Length/@FoldList[Complement [#1 \[Union] #2, #1 \[Intersection] #2] &, {}, List @@@ pdList]]


widthPresentationList[crossings_] := widthPresentationList[crossings] = Length/@FoldList[Complement [#1 \[Union] #2, #1 \[Intersection] #2] &, {}, List @@@ crossings]


maxWidth[startWidth_, crossings_, index_] := SortBy[Permutations[crossings], Max[Length /@ FoldList[Complement[#1\[Union] #2, #1 \[Intersection] #2] &, startWidth, List @@@ #]]&][[1]];


PartialJonesOptimize[startWidth_, crossings_, index_, portion_]:=
Module[{permutations,widths, catalanWidths, optimizeResult, minimumPermutation},
permutations = Join[RandomSample[Permutations[crossings], Ceiling[portion* Factorial[Length[crossings]]]], {crossings}];
widths =Delete[(Length /@ FoldList[Complement[#1\[Union] #2, #1 \[Intersection] #2] &, startWidth, List @@@ #]),1]&/@ permutations;
catalanWidths = (CatalanNumber[#/2]&/@ #)&/@ widths;
optimizeResult = Total[#] &/@( ConstantArray[Range[index,index + Length[crossings] - 1] , Length[permutations]] * widths * catalanWidths);
minimumPermutation = permutations[[First[First[Position[optimizeResult, Sort[optimizeResult][[1]]]]]]];
minimumPermutation]


Protect[startPt, invariantCriteria, optSize, partialPortion, rounds];


Options[backwardsPLO] := {startPt -> convertToList, invariantCriteria -> PartialJonesOptimize, optSize -> 3, partialPortion -> 0.1, rounds -> 3};


backwardsPLO[pd_PD, opts: OptionsPattern[backwardsPLO]] := 
Module[{s, k, sBest, optimizeList, optimizeResult, widths, width, round},
s = OptionValue[startPt][pd];
k = Length[s] - OptionValue[optSize] + 1;
round = 1;
sBest = s;
While[round <= OptionValue[rounds],
widths = FoldList[Complement [#1 \[Union] #2, #1 \[Intersection] #2] &, {}, List @@@ sBest];
While[k >= 1,
optimizeList = sBest[[k;;k+OptionValue[optSize] - 1]];
width =  widths[[k]];
optimizeResult = OptionValue[invariantCriteria][width, optimizeList,k, OptionValue[partialPortion]];
sBest = Join[If[k == 1, {}, sBest[[1;;k-1]]], optimizeResult, If[k == Length[s] - OptionValue[optSize] + 1, {}, sBest[[k + OptionValue[optSize];;-1]]]];
k = k - 1;
];
round = round + 1;
k  = Length[s] - OptionValue[optSize] + 1;
];
sBest]


convertToList[pd_PD] := List @@ pd;


Protect[GreedyCriteria, InvariantCriteria, OptSize, ComparedSamples];


Options[OptimizedKnotPresentation] := {GreedyCriteria -> JonesResult, InvariantCriteria -> PartialJonesOptimize, OptSize -> Min[Floor[Length[List @@ pd]/150] + 5, 10], ComparedSamples -> 1};


OptimizedKnotPresentation[pd_PD, opts: OptionsPattern[OptimizedKnotPresentation]] := 
(SortBy[Table[backwardsPLO[PD @@ bestGreedy[pd, {factor -> OptionValue[GreedyCriteria], tries -> 1000}], 
{startPt -> convertToList, invariantCriteria -> OptionValue[InvariantCriteria], optSize -> OptionValue[OptSize], partialPortion -> 0.1, rounds -> 10}], {i, OptionValue[ComparedSamples]}], JonesResult[#]&][[1]]) /. (Global`X -> X)


End[]


EndPackage[]
