(* ::Package:: *)

BeginPackage["PermutationGraphs`"]

(* \[Copyright] Copyright 2016 Tom Verhoeff *)

(*
    This package is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This package is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this package.  If not, see <http://www.gnu.org/licenses/>.
*)

(*
	Permutations are lists of natural numbers, possibly with duplicates.
	The list of occurence frequencies of numbers in such permutations is called its signature.
	A signature is represented by a list s of natural numbers,
	where s[[i]] is the occurrence frequency of number i-1.
*)

multiset::usage=
	"Creates an ordered list of natural numbers of given occurrence frequencies.";

extend::usage=
	"Extends the given permutation(s) with a given tail.";

shorten::usage=
	"Shortens the given permutation(s) with a given length.";

signature::usage=
	"Determines the signature (occurrence frequencies) of a given permutation.";

neighbor::usage=
	"Whether two permutations of the same length are neighbors in the neighbor-swap graph.
	That is, whether they differ in the swap of two adjacent elements.";

swapPair::usage=
	"Swaps a pair of numbers in a permutation.";

edgeIndex::usage=
	"Returns leftmost index of pair swapped for two neighbors.";

coordinates::usage=
	"Returns coordinate list for given permutation.";


halveSignature::usage=
	"Returns the signature with frequencies halved (rounded down).";

stutterCount::usage=
	"Returns the number of stutter permutations of a given signature.";

nonStutterCount::usage=
	"Returns the number of non-stutter permutations of a given signature.";

stutters::usage=
	"List of boolean values indicating for each left-index-odd pair
	whether the elements in the pair are equal.";

stutterPermutationQ::usage=
	"Determines whether a given permutation is a stutter permutation.";

stutterFreeQ::usage=
	"Tests whether a collection of permutations contains no stutter permutations.";

stutterize::usage=
	"Converts argument into stutter permutation(s), by repeating every number";

selectOdds::usage=
	"Returns list of numbers with odd occurrence frequencies in the given signature.";

permutations::usage=
	"Constructs list of all permutations with given signature.";

stutterPermutations::usage=
	"Constructs list of all stutter permutations with given signature.";

nonStutterPermutations::usage=
	"Constructs list of all non-stutter permutations with given signature.";

selectByTail::usage=
	"Selects those permutations that end with the given tail.";

HpathQ::usage=
	"Determines whether the path is a Hamiltonian path on the non-stutter pmerutations of the given signature.";

HcycleQ::usage=
	"Determines whether the path is a Hamiltonian cycle on the non-stutter pmerutations of the given signature.";

showPath::usage=
	"Shows the given path together with some properties.";

transform::usage=
	"Transforms the permutation(s) according to the given renaming.";


conditionHpath::usage=
	"Condition for existence of Hamiltonian path on all permutations of the given signature.";

conditionHcycle::usage=
	"Condition for existence of Hamiltonian cycle on all permutations of the given signature.";

conditionHpathNS::usage=
	"Condition for existence of Hamiltonian path on all non-stutter permutations of the given signature.";

conditionHcycleNS::usage=
	"Condition for existence of Hamiltonian cycle on all non-stutter permutations of the given signature.";

Hpath::usage=
	"Constructs (if possible) a Hamiltonian path/cycle on all permutations of the given signature.";

HpathNS::usage=
	"Constructs a Hamiltonian path (or cycle) on all non-stutter permutations of the given signature.";

HpathAlt::usage=
	"Constructs a Hamiltonian path on all permutations of the given signature; serves for [k, 2, 1].";

HpathAltQ::usage=
	"Determines whether the path is a Hamiltonian path on all non-stutter and one stutter permuations of the given signature";

HpathNSAlt::usage=
	"Constructs a Hamiltonian path on all non-stutter and one stutter permuations of the given signature";


Begin["Private`"];


Needs["PathOperations`","/Users/wstomv/Documents/Mathematics/Mathematica Definitions/PathOperations.wl"]


(* multiset[freq]

	pre: freq is a list of natural numbers, representing occurrence frequencies
	ret: lexicographically smallest list with given occurrence frequencies,
		that is, number i occurs with frequency freq[[i+1]], and the list is sorted
*)

multiset[freq_List] :=
	Catenate @ MapIndexed[ConstantArray[First[#2]-1, #1]&, freq];

multiset[freq___Integer] :=
	multiset[{freq}];


(* extend[perm, or list of perms, tail]

	pre: tail is a list of numbers
	ret: tail is appended to each perm 
*)

extend[perm:{___Integer}, tail_List] :=
	Join[perm, tail];

extend[perms_List, tail_List] :=
	Map[extend[#, tail]&, perms];

extend[tail_List][perm:{___Integer}] :=
	extend[perm, tail];

extend[tail_List][perms_List] :=
	extend[perms, tail];


(* shorten[perm, or list of perms, len]

	pre: len is a natural number
	ret: each perm is shortened by dropping len items
*)

shorten[perm:{___Integer}, len_Integer] :=
	Drop[perm, -len];

shorten[perms_List, len_Integer] :=
	Map[shorten[#, len]&, perms];

shorten[len_Integer][perm:{___Integer}] :=
	shorten[perm, len];

shorten[len_Integer][perms:{___Integer}] :=
	shorten[perms, len];


(* signature[perm]

	pre: perm is a list of natural numbers
	ret: list such that i-th element is occurence frequency of i-1 in perm; without trailing zeroes
*)

signature[perm_List] :=
	signature[perm, {}];

signature[{}, s:{___Integer}] :=
	s;

signature[perm:{a_Integer, t___Integer}, s:{z___Integer}] :=
	If[Length[s] < a+1,
		signature[perm, {z, 0}],
		signature[{t}, ReplacePart[s, a+1 -> s[[a+1]]+1]]
	];


(* neighbor[p, q]

	pre: Length[p] \[Equal] Length[q]
	ret: whether p and q differ by a swap of two adjacent elements
*)

neighbor[p_List, q_List] :=
	With[
		{diff = Position[MapThread[#1 != #2 &, {p, q}], True]},
		And[
			Length[diff] == 2,
			diff[[1, 1]]+1 == diff[[2, 1]],
			Extract[p, diff] == Reverse @ Extract[q, diff] (* see Note *)
		]
	];
(* Note: this third condition is superfluous when p and q have the same signature *)

neighbor[{p_List, q_List}] :=
	neighbor[p, q];


(* swapPair[perm, i]

	pre: 1 \[LessEqual] i, j \[LessEqual] Length[perm]   or   1 \[LessEqual] i < Length[perm] if only one parameter
	ret: perm with numbers at positions i and i+1 swapped
*)

swapPair[perm_List, i_Integer, j_Integer] :=
	ReplacePart[perm, {i->perm[[j]], j->perm[[i]]}] /; i != 0 && j != 0

swapPair[perm_List, i_Integer] :=
	With[{j = i+1},
		ReplacePart[perm, {i->perm[[j]], j->perm[[i]]}]
	] /; i != 0

swapPair[i_Integer, j_Integer][perm_List] :=
	swapPair[perm, i, j];

swapPair[i_Integer][perm_List] :=
	swapPair[perm, i];


(* edgeIndex[p, q]

	pre: neighbor[p, q]
	ret: i such that swapPair[p, i] \[Equal] q
*)

edgeIndex[p_List, q_List] :=
	Position[MapThread[#1 != #2 &, {p, q}], True][[1, 1]];


(* coordinates[perm]

	pre: perm is a permutation
	ret: coordinates for embedding given permuatation
*)

coordinates[perm_, straight_] :=
	With[{pos0 = Position[perm, 0][[1, 1]]},
		{Sequence @@ Flatten[Position[If[straight, DeleteCases[perm, 0], perm], 1]], pos0}
	];


(* halveSignature[s]

	pre: s is a signature
	ret: signature s with every frequency halved (rounded down)
*)

halveSignature[s___List] :=
	Map[Quotient[#, 2]&, s];


(* stutterCount[s], nonStutterCount[s]

	pre: s is a signature
	ret: number of (non-)stutter permutations with signature s
*)

stutterCount[s___Integer] :=
	With[{list = {s}},
		If[Length[Select[list, OddQ]] >= 2,
			0,
			Multinomial @@ halveSignature[list]
		]
	];

stutterCount[{s___Integer}] :=
	stutterCount[s];

nonStutterCount[s___Integer] :=
	Multinomial[s] - stutterCount[s];

nonStutterCount[{s___Integer}] :=
	nonStutterCount[s];


(* stutters[perm]

	pre: perm is a permutation (list of natural numbers)
	ret: list with for each left-index-odd pair whether the elements in the pair are equal
*)

stutters[perm_List] :=
	Map[First[#] == Last[#] &,
		Partition[perm, 2]
	];


(* stutterPermutationQ[perm]

	pre: perm is a permutation (list of natural numbers)
	ret: whether perm is a stutter permutation
*)

stutterPermutationQ[perm_List] :=
	And @@ stutters[perm];


(* stutterFreeQ[perms]

	pre: perms is a list of permutations
	ret: whether perms is free of stutter permutations
*)

stutterFreeQ[perms_List] :=
	AllTrue [perms, ! stutterPermutationQ[#] &];


(* stutterize[perm or list of perms]

	pre: argument is a (list of) permutations
	ret: stutterized version of argument
*)

stutterize[perm:{___Integer}] :=
	Catenate @ Transpose[{perm, perm}];

stutterize[perms:{__List}] :=
	Map[stutterize, perms];


(* selectOdds[s]

	pre: s is a signature
	ret: list of numbers that occur with odd frequency
*)

selectOdds[s___Integer] :=
	selectOdds[{s}];

selectOdds[s:{___Integer}] :=
	Map[First[#]-1 & , Position[s, _?OddQ]];


(* permutations[s], stutterPermutations[s], nonStutterPermutations[s]

	pre: s is a signature
	ret: list of all stutter permutations with signature s
*)

permutations[s___Integer] :=
	Permutations[multiset[s]];

permutations[{s___Integer}] :=
	permutations[s];

stutterPermutations[s___Integer] :=
	Block[{odds, result},
		odds = selectOdds[s];
		If[Length[odds] >= 2,
			{},
			result = stutterize @ permutations @ halveSignature[{s}];
			If[Length[odds] == 1,
				extend[result, odds],
				result
			]
		]
	];

stutterPermutations[{s___Integer}] :=
	stutterPermutations[s];

nonStutterPermutations[s___Integer] :=
	Select[permutations[s], ! stutterPermutationQ[#]&];

nonStutterPermutations[{s___Integer}] :=
	nonStutterPermutations[s];


(* selectByTail[perms, tail]

	pre: perms is a list of permutations
	ret: list with those permutations in perms that end with the given tail
*)

selectByTail[perms_List, tail_List] :=
	With[{len = Length[tail]},
		Select[perms, Take[#, -len] == tail &]
	];

selectByTail[tail_List][perms_List] :=
	selectByTail[perms, tail];


(* HpathQ[perms, s], HcycleQ[perms, s]

	pre: perms is a list of permutations, and s is a signature.
	ret: whether perms is a Hamiltonian path/cycle on the non-stutter permutations of the given signature
*)

HpathQ[perms_List, s:{___Integer}] :=
	And[
		pathQ[perms, neighbor],
		Length[perms] == nonStutterCount[s],
		DuplicateFreeQ[perms],
		stutterFreeQ[perms],
		AllTrue[perms, signature[#] == s &]
	];

HcycleQ[perms_List, s:{___Integer}] :=
	HpathQ[perms, s] && cycleQ[perms, neighbor];


(* showPath[perms]

	pre: perms is a list of permutations
	ret: a grid with the path and some of its properties
*)

showPath[perms:{}] :=
	Grid[{
		{"path",perms},
		{"length",Length[perms]}
		},
		Frame->All
	];

showPath[perms:{perm_List, z___List}] :=
	With[{s = signature[perm]},
		Grid[{
			{"signature", Row[s, ","]},
			{"# vertices:", Multinomial @@ s},
			{"# stutter permutations:", stutterCount @@ s},
			{"# non-stutter permutations:", nonStutterCount @@ s},
			{"length", Length[perms]},
			{"is path?", pathQ[perms, neighbor]},
			{"is cycle?", cycleQ[perms, neighbor]},
			{"stutter free?", stutterFreeQ[perms]},
			{"path",
				If[Length[perms] <= 100,
					Grid[Prepend[{"permutation","is neighbor?", "is stutter?"}] @
						Transpose[{
							perms,
							Prepend[Null] @ Map[neighbor[Sequence @@ #]&, Partition[perms, 2, 1]],
							Map[stutterPermutationQ, perms]
						}],
						Frame->All
					],
					Grid[{{First@perms}, { ". . ."}, {Last @ perms}}, Frame->All]
				]
			}
			},
			Frame->All
		]
	];


(* transform[perm or list of perm, tr]]

	pre: tr is a list of natural numbers, and
		all numbers in the permutations are < Length[tr]
	ret: perm is transformed according to tr,
		that is, each occurence of i with 0 \[LessEqual] i < Length[tr] is replaced by tr[[i+1]]
*)

transform[perm:{___Integer}, tr_List] :=
	Map[tr[[#+1]] &, perm];

transform[perms:{___List}, tr_List] :=
	Map[transform[#, tr] &, perms];

transform[tr_List][perm:{___Integer}] :=
	transform[perm, tr];

transform[tr_List][perms:{___List}] :=
	transform[perms, tr];


(* conditionHpath[s], conditionHcycle[s], conditionHpathNS[s], conditionHcycleNS[s]

	pre: s is a signature
	ret: Whether all (non-stutter) permutations with signature s admit a Hamiltonian path/cycle
*)

conditionHpath[{s___Integer}] :=
	conditionHpath[s];

conditionHpath[s___Integer, 0, t___Integer] :=
	conditionHpath[s, t];

conditionHpath[s___Integer] :=
	With[{len = Length[{s}]},
		(len <= 1) ||
		(len == 2 && MemberQ[{s}, 1]) ||
		(Length[selectOdds[s]] >= 2)
	];

conditionHcycle[{s___Integer}] :=
	conditionHcycle[s];

conditionHcycle[s___Integer, 0, t___Integer] :=
	conditionHcycle[s, t];

conditionHcycle[s___Integer] :=
	conditionHpath[s] &&
	Length[{s}] > 2 &&
	! MatchQ[Sort[{s}, Greater], {_?EvenQ, 1, 1}];

conditionHpathNS[{s___Integer}] :=
	conditionHpathNS[s];

conditionHpathNS[s___Integer, 0, t___Integer] :=
	conditionHpathNS[s, t];

conditionHpathNS[s___Integer] :=
	True;

conditionHcycleNS[{s___Integer}] :=
	conditionHcycleNS[s];

conditionHcycleNS[s___Integer, 0, t___Integer] :=
	conditionHcycleNS[s, t];

conditionHcycleNS[s___Integer] :=
	conditionHpathNS[s] &&
	! (Length[{s}] <= 2 && Count[{s}, _?EvenQ] != 2) &&
	! MatchQ[Sort[{s}, Greater], {_?EvenQ, 1, 1}]


(* Hpath[s]

	pre: s is a signature, either with length \[LessEqual] 2 and at least one 1, or with at least two odd frequencies
	ret: Hamiltonian path (or cycle) on all permutations with signature s
*)

Hpath[{s___Integer}] :=
	Hpath[s];

Hpath[s___Integer] :=
	Missing["Nonexistent"] /; ! conditionHpath[s];

Hpath[s___Integer] :=
	With[{list = {s}},
		transform[Hpath[Sequence @@ Sort[list, Greater]], Map[#-1 &, Ordering[list, All, Greater]]]
	] /; (conditionHpath[s] && ! OrderedQ[{s}, GreaterEqual]);

(* drop all trailing 0 frequencies *)
Hpath[s___Integer, 0] :=
	Hpath[s];

Hpath[] :=
	{};

Hpath[n_Integer] :=
	{multiset[n]};

Hpath[k0_Integer, 1] :=
	Block[{result},
		result = Reverse @ IdentityMatrix[k0+1]
	];


(* HpathNS[s]

	pre: s is a signature
	ret: Hamiltonian path (or cycle) on all non-stutter permutations with signature s
*)

HpathNS[{s___Integer}] :=
	HpathNS[s];

(* If the signature is not sorted, then sort it, solve it, and transform back. *)
HpathNS[s___Integer] :=
	With[{list = {s}},
		transform[HpathNS[Sequence @@ Sort[list, Greater]], Map[#-1 &, Ordering[list, All, Greater]]]
	]/;(!OrderedQ[{s}, GreaterEqual]);

(* drop all trailing 0 frequencies *)
HpathNS[s___Integer, 0] :=
	HpathNS[s];

HpathNS[] :=
	{};

HpathNS[n_Integer] :=
	{};
(* N.B. the stutter permutation 0^n could be included. *)

(* ret: path from 0^k0 1 to 1 0^k0, if k0 odd, else from 0^(k0-1) 1 0 to 1 0^k0 *)
HpathNS[k0_Integer, 1] :=
	Block[{result},
		result = Reverse @ IdentityMatrix[k0+1];
		If[OddQ[k0], result, Rest @ result]
	];
(* N.B. If k0 is even, then the stutter permutation could be included at the start. See HpathNSAlt. *)


(* HpathNS[k0_?EvenQ, k1_?EvenQ]

	ret: cycle starting in 1^(k1-1) 0^k0 1, and containing the two specific edges
		0^(k0-2) 1 0 1^(k1-1) 0 ~ 0^(k0-1) 1^k1 0
		1^(k1-2) 0 1 0^(k0-1) 1 ~ 1^(k1-1) 0^k0 1
*)

HpathNS[k0_?EvenQ, k1_?EvenQ] :=
	(* k0 \[GreaterEqual] k1 \[GreaterEqual] 2 *)
	Block[{p0, p1},
		(* split on trailing bit x *)
		p1 = extend[HpathNS[k0, k1-1], {1}]; (* p1 is a path from 0^(k0-1) 1^(k1-1) 0 1 to 1^(k1-1) 0^k0 1 *)
		p0 = extend[HpathNS[k0-1, k1], {0}];
		If[k0 == k1, (* k0-1 < k1 *) p0 = Reverse @ p0]; (* p0 is a path from 0^(k0-1) 1^k1 0 to 1^(k1-1) 0^(k0-1) 1 0 *)

		Join[
			Reverse @ p1,
			p0
		]
	];


(* HpathNS[k0_?OddQ, k1_?EvenQ]

	ret: path from 0^k0 1^k1 to 1^(k1-1) 0^k0 1
*)

HpathNS[k0_?OddQ, k1_?EvenQ] :=
	(* k0 > k1 \[GreaterEqual] 2 *)
	Block[{p0, p1},
		(* split on trailing bit x *)
		p1 = extend[HpathNS[k0, k1-1], {1}]; (* a Ham. path from 0^k0 1^k1 to 1^(k1-1) 0^k0 1 *)
		p0 = extend[HpathNS[k0-1, k1], {0}]; (* a Ham. cycle from 1^(k1-1) 0^(k0-1) 1 0 *)

		Join[
			Most @ p1,
			RotateLeft[p0],
			{Last @ p1}
		]
	];

(* HpathNS[k0_?EvenQ, k1_?OddQ]

	ret: path from 0^(k0-1) 1^k1 0 to 1^k1 0^k0
*)

HpathNS[k0_?EvenQ, k1_?OddQ] :=
	(* k0 > k1 \[GreaterEqual] 3 *)
	Block[{p0, p1, v},
		(* split on trailing bit x *)
		p1 = extend[HpathNS[k0, k1-1], {1}]; (* a Ham. cycle containing edge 0^(k0-1) 1^(k1-1) 0 1 ~ 0^(k0-2) 1 0 1^(k1-2) 0 1 *)
		p0 = extend[HpathNS[k0-1, k1], {0}]; (* a Ham. path from 0^(k0-1) 1^k1 0 to 1^k1 0^k1 *)
		v = First @ p0;
		Join[
			{v},
			cutCycle[Reverse @ p1, swapPair[v, -2]],
			Rest @ p0
		]
	];

(* Is there also always a path that includes one stutter permutation at the end?
	Yes, see HpathNSAlt. *)


(* HpathNS[k0_?OddQ, k1_?OddQ]

	ret: path from 0^k0 1^k1 to 1^k1 0^k0
*)

HpathNS[k0_?OddQ, k1_?OddQ] :=
	(* k0 \[GreaterEqual] k1 \[GreaterEqual] 3 *)
	Block[{
			p11 = HpathNS[k0, k1-2](* p11 is a path from 0^k0 1^(k1-2) *),
			p1101 = HpathNS[k0-1, k1-3] (* this is a cycle or empty; also need spurs; p1110 is the same, and fully parallel to p1101 *),
			p0101 = HpathNS[k0-2, k1-2] (* this is a path from 0^(k0-2) 1^(k1-2); p1001, p0110, p1010 are the same, and fully parallel to p0101 *),
			p0001 = HpathNS[k0-3, k1-1] (* this is a cycle or empty; also need spurs; p0010 is the same, and fully parallel to p0001 *),
			p00 = HpathNS[k0-2, k1] (* p00 is a path; from 0^(k0-2) 1^k1m if k0-2\[GreaterEqual]k1, else reversed *),
			sp11 = extend[stutterPermutations[k0-1, k1-3], {1, 1}],
			sp00 = extend[stutterPermutations[k0-3, k1-1], {0, 0}],
			p11xy, p00xy, c11xy, c00xy, tube
		},
		(*Print[{"p11,p1101,p0101,p0001,p00,sp11,sp00 (before):",p11,p1101,p0101,p0001,p00,sp11,sp00}];*)
		(* adjust for boundary case where p1101 or p0001 is empty ; align the cycles with paths,
			and incorporate the stutter permutations *)
		c11xy = If[Length[p1101] == 0,
			Map[Join[sp11[[1]], #] &, {{0, 1}, {1, 0}}],
			p11xy = RotateLeft @ createZigZagPath[extend[cutCycle[p1101, Append[0] @ Most @ p0101[[1]]], {1, 1}], {1, 0}, {0, 1}];
			(*Print["p11xy:",p11xy];*)
			incorporateSpursInZigZag[p11xy, sp11, neighbor, 2]
		];
		c00xy = If[Length[p0001] == 0,
			Map[Join[sp00[[1]], #] &, {{1, 0}, {0, 1}}],
			p00xy = RotateLeft @ createZigZagPath[extend[cutCycle[p0001, Append[1] @ Most @ p0101[[-1]]], {0, 0}], {1, 0}, {0, 1}];
			(*Print["p00xy:",p00xy];*)
			Reverse @ incorporateSpursInZigZag[p00xy, sp00, neighbor, 2]
		];
		If[k0-2 < k1, p00 = Reverse @ p00];

		(* combine p11, 2*p1101+spikes, 4*p0101, 2*p0001+spikes, p00 into one path *)

		(* p1101 and its copy p1110 will be combined into a cycle having all w01-w10 edges, and then with their stutter permutations into one cycle *)

		(* p0001 and its copy p0010 will be combined into a cycle having all w01-w10 edges, and then with their stutter permutations into one cycle *)

		(* the four copies of p0101, p0110, p1001, p1010 are combined into a path *)
		tube = createSquareTube[Reverse @ p0101, {0, 1}, {1, 0}];
		(*Print[{"c11xy, c00xy, tube:",c11xy, c00xy, tube}];*)

		Join[extend[p11, {1, 1}], tube[[1;;3]], c00xy, tube[[4;;-3]], c11xy, tube[[-2;;-1]], extend[p00, {0, 0}]]
	];


(* [k, 1, 1]

	If k is even, then return a Hamiltonian path from 1 0^k 2 to 2 0^k 1.
	If k is odd, then return a Hamiltonian cycle, starting at 1 0^k 2.
*)

HpathNS[k_Integer, 1, 1] :=
	(* k \[GreaterEqual] 1 *)
	Block[{ (* Split off the trailing number x *)
			p2 = extend[Hpath[k, 1, 0], {2}] (* a path from 0^k 1 2 to 1 0^k 2 *),
			p1 = extend[Hpath[k, 0, 1], {1}] (* a path from 2 0^k 1 to 0^k 2 1 *),
			p0 = extend[HpathNS[k-1, 1, 1], {0}] (* by IH, a path/cycle from 1 0^k 2 (path to 2 0^k 1) *)
		},
		(*Print[{"p2, p1, p0 (before):", p2, p1, p0}];*)
		If[k == 1,
			(* reverse these, because sorting of signature reversed up the order *)
			p1 = Reverse @ p1;
			p0 = Reverse @ p0;
		];
		(*Print[{"p2, p1, p0 (after):", p2, p1, p0}];*)
		If[OddQ[k],
			Join[Take[p2, -1], p0, Reverse @ p1, Most @ p2],
			(* EvenQ[k], also need to add 0^k0 1 2 and 0^k0 2 1 *)
			Join[Take[p2, -1], p0, Reverse @ Most @ p2, p1]
		]
	];


(* HpathAlt[s]  Generic part of definition

	pre: s is a signature
	ret: alternative Hamiltonian path/cycle on all permutations of given signature;
		in view of cycle for HpathNS[k, 2, 1]

*)

HpathAlt[{s___Integer}] :=
	HpathAlt[s];

(* If the signature is not sorted, then sort it, solve it, and transform back. *)
HpathAlt[s___Integer] :=
	With[{list = {s}},
		transform[HpathAlt[Sequence @@ Sort[list, Greater]], Map[#-1 &, Ordering[list, All, Greater]]]
	]/;(!OrderedQ[{s}, GreaterEqual]);

(* drop all trailing 0 frequencies *)
HpathAlt[s___Integer, 0] :=
	HpathAlt[s];


(* HpathAlt[k, 1, 1] alternatives, in view of cycle for HpathNS[k, 2, 1]

	ret: a Hamiltonian _path_ from 1 2 0^k to 0 2 1 0^(k-1),
		including segment [k, 1]~2 from 0^k 1 2 to 1 0^k 2
*)

HpathAlt[2, 1, 1] :=
	Block[{
			p2 = extend[Hpath[2, 1, 0], {2}],
			p1 = extend[Hpath[2, 0, 1], {1}]
		},
		Join[
			{{1, 2, 0, 0}, {2, 1, 0, 0}, {2, 0 ,1, 0}},
			Reverse @ p1,
			p2,
			{{1, 0, 2, 0}, {0, 1, 2, 0}, {0, 2, 1, 0}}
		]
	];

HpathAlt[3, 1, 1] :=
	Block[{
			p2 = extend[Hpath[3, 1, 0], {2}],
			p1 = extend[Hpath[3, 0, 1], {1}]
		},
		Join[
			{{1, 2, 0, 0, 0}, {2, 1, 0, 0, 0}, {2, 0 ,1, 0, 0}, {2, 0, 0, 1, 0}},
			Reverse @ p1,
			p2,
			{{1, 0, 0, 2, 0}, {1, 0, 2, 0, 0}, {0, 1, 2, 0, 0}, {0, 1, 0, 2, 0},
				{0, 0, 1, 2, 0}, {0, 0, 2, 1, 0}, {0, 2, 0, 1, 0}, {0, 2, 1, 0, 0}
			}
		]
	];

HpathAlt[k_Integer, 1, 1] :=
	Block[{
		p2 = extend[Hpath[k, 1, 0], {2}],
		p1 = extend[Hpath[k, 0, 1], {1}],
		p20 = extend[Hpath[k-1, 1, 0], {2, 0}],
		p10 = extend[Hpath[k-1, 0, 1], {1, 0}],
		p00 = extend[HpathAlt[k-2, 1, 1], {0, 0}] (* by ind. hyp. a path from 1 2 0^(k-2) to 0 2 1 0^(k-3) *)
	},
		(*Print[{"p2, p1, p20, p10, p00:", p2, p1, p20, p10, p00}];*)
		Join[
			p00[[;;k]],
			Take[p10, -1],
			Reverse @ p1,
			p2,
			Reverse @ p20,
			Most @ p10,
			p00[[k+1;;]]
		]
	] /; k >= 4;


(* HpathNS[even, 2, 1] 

	ret: cycle on non-stutter permutations for signature [even, 2, 1]
*)

HpathNS[k_?EvenQ, 2, 1] :=
	(* k \[GreaterEqual] 2 *)
	Block[{
			p2 = extend[HpathNS[k, 2, 0], {2}] (* a cycle *),
			p1 = extend[HpathAlt[k, 1, 1], {1}], (* a path from 1 2 0^k to 0 2 1 0^(k-1)*)
			p0 = extend[If[k-1 == 1, transform[{2, 1, 0}], Reverse] @ HpathAlt[k-1, 2, 1], {0}]
			(* when k-1 \[Equal] 1, swap 0 and 2 in [k-1, 2, 1] *),
			v = Join[{1}, ConstantArray[0, k], {1, 2}],
			c
		},
		(*Print[{"p2, p1, p0:", p2, p1, p0}];*)
		c = Join[p0, p1];
		(*Print[{"c:", c}];*)
		Join[
			cutCycle[p2, swapPair[v, 1]],
			cutCycle[c, swapPair[v, -2]]
		]
	];


(* HpathAlt[odd, 2, 1] alternative, in view of cycle for HpathNS[even, 2, 1]

	ret: path from 1 2 0^k 1 to 0 2 1 0^(k-1) 1
*)

HpathAlt[k_?OddQ, 2, 1] :=
	(* k \[GreaterEqual] 3 *)
	Block[{
			p12 = extend[Hpath[k, 1, 0], {1, 2}],
			p02 = HpathNS[k-1, 2, 0] (* p20 is parallel to p02 *),
			p1 = extend[HpathAlt[k, 1, 1], {1}],
			p10 = extend[HpathAlt[k-1, 1, 1], {1, 0}],
			p00 = extend[If[k-2 == 1, transform[{2, 1, 0}], Reverse] @ HpathAlt[k-2, 2, 1], {0, 0}]
			(* when k-2 \[Equal] 1, swap 0 and 2 in [k-2, 2, 1] *),
			sp = stutterPermutations[k-1, 2] (* stutter permutations for p02, p20 *),
			cns, c2,
			split1, p1p12,
			v1 = Join[ConstantArray[0, k], {1, 2, 1}],
			v2 = Join[{1}, ConstantArray[0, k], {2, 1}],
			c0,
			p, split2
		},
		(*Print[{"p12, p02, p1, p10, p00:", p12, p02, p1, p10, p00}];*)

		cns = createZigZagPath[p02, {2, 0}, {0, 2}];
		c2 = incorporateSpursInZigZag[cns, sp, neighbor, 2];
		(*Print[{"sp, cns, c2:", sp, cns, c2}];*)

		split1 = splitPath[p1, {v1, swapPair[v2, -3]}]; (* 3 segments *)
		p1p12 = Catenate @ MapAt[createZigZagPath[shorten[#, 2], {2, 1}, {1, 2}] &, split1, 2];
		(*Print[{"split1, p1p12:", split1, p1p12}];*)

		c0 = Join[p10, p00];
		(*Print["c0:", c0];*)

		p = Join[p1p12[[1;;2]], Reverse @ RotateLeft[c0, 2], p1p12[[3;;]]];
		(*Print[{"v2 neighbor, p:", swapPair[v2, -2], p}];*)
		split2 = splitPath[p, {swapPair[v2, -2]}];
		(*Print[{"split2:", split2}];*)

		Join[
			split2[[1]],
			RotateLeft[c2, 2],
			split2[[2]]
		]
	];


(* HpathAltQ[perms, s]

	pre: perms is a list of permutations, and s is a signature.
	ret: whether perms is a Hamiltonian path on the non-stutter and one stutter permutations of the given signature
*)

HpathAltQ[perms_List, s:{___Integer}] :=
	And[
		pathQ[perms, neighbor],
		Length[perms] == nonStutterCount[s] + 1,
		DuplicateFreeQ[perms],
		stutterFreeQ[Most[perms]],
		stutterPermutationQ[Last[perms]],
		AllTrue[perms, signature[#] == s &]
	];



(* HpathNSAlt[s] alternative Hamiltonian paths that includes one stutter permutation at the end;
	generic part of definition

	pre: s is a signature
	ret: Hamiltonian path on all non-stutter permutations, plus one stutter permutation at the end
*)

HpathNSAlt[{s___Integer}] :=
	HpathNSAlt[s];

(* If the signature is not sorted, then sort it, solve it, and transform back. *)
HpathNSAlt[s___Integer] :=
	With[{list = {s}},
		transform[HpathNSAlt[Sequence @@ Sort[list, Greater]], Map[#-1 &, Ordering[list, All, Greater]]]
	]/;(!OrderedQ[{s}, GreaterEqual]);

(* drop all trailing 0 frequencies *)
HpathNSAlt[s___Integer, 0] :=
	HpathNSAlt[s];

HpathNSAlt[] :=
	Hpath[];

HpathNSAlt[n_Integer] :=
	Hpath[n];

HpathNSAlt[k0_EvenQ, 1] :=
	Reverse @ Hpath[k0, 1];


(* HpathNSAlt[k0_?EvenQ, k1_?EvenQ]

	ret: a pair of paths,
		the first from  1^(k1-1) 0^k0 1 to stutter perm 1^k1 0^k0, and
		the second from 0^(k0-1) 1^k1 0 to stutter perm 0^k0 1^k1;
		these are needed for the odd+even and even+odd cases
*)

HpathNSAlt[k0_?EvenQ, k1_?EvenQ] :=
	(* k0 \[GreaterEqual] k1 \[GreaterEqual] 2 *)
	Block[{p0, p1, q0, q1},
		(* split on trailing bit x *)
		p1 = extend[HpathNS[k0, k1-1], {1}]; (* p1 is a path from 0^(k0-1) 1^k1 0 to 1^k1 0^k0 *)
		p0 = extend[HpathNSAlt[k0-1, k1], {0}]; (* p0 is a path from 0^(k0-1) 1^k1 0 to 1^k1 0^k0 *)

		q1 = extend[HpathNSAlt[k0, k1-1], {1}]; (* q1 is a path from 1^(k1-1) 0^k0 1 to 0^k0 1^k1 *)
		q0 = extend[HpathNS[k0-1, k1], {0}];
		If[k0 == k1, (* k0-1 < k1 *) q0 = Reverse @ q0]; (* q0 is a path from 0^(k0-1) 1^k1 0 to 1^k1 0^k0 *)

		{Join[Reverse @ p1, p0], Join[q0, q1]}
	];


(* HpathNSAlt[k0_?OddQ, k1_?EvenQ]

	ret: a path from 0^k0 1^k1 to (stutter) 1^k1 0^k0
*)

HpathNSAlt[k0_?OddQ, k1_?EvenQ] :=
	(* k0 > k1 \[GreaterEqual] 2 *)
	Block[{p0, p1},
		(* split on trailing bit x *)
		p1 = extend[HpathNS[k0, k1-1], {1}]; (* a Ham. path from 0^k0 1^k1 to 1^(k1-1) 0^k0 1 *)
		p0 = extend[HpathNSAlt[k0-1, k1][[1]], {0}]; (* a path from 1^(k1-1) 0^(k0-1) 1 0 ending in stutter perm *)

		Join[p1, p0]
	];

(* HpathNSAlt[k0_?EvenQ, k1_?OddQ]

	ret: a path from 1^k1 0^k0 to (stutter) 0^k0 1^k1
*)

HpathNSAlt[k0_?EvenQ, k1_?OddQ] :=
	(* k0 > k1 \[GreaterEqual] 3 *)
	Block[{p0, p1},
		(* split on trailing bit x *)
		p1 = extend[HpathNSAlt[k0, k1-1][[2]], {1}]; (* a path from 0^k0 1^k1 to 1^(k1-1) 0^k0 1 *)
		p0 = extend[HpathNS[k0-1, k1], {0}]; (* a Ham. path from 0^(k0-1) 1^k1 0 to 1^k1 0^k1 *)
		Join[Reverse @ p0, p1]
	];


End[ ];

EndPackage[ ]
