(* :Title: Bijective Rules *)

(* :Context: BijectiveRules` *)

(* :Author: Eric Rowland *)
(*          http://math.rutgers.edu/~erowland *)

(* :Date: {2008, 7, 4} *)

(* :Package Version: 1.10 *)

(* :Mathematica Version: 6.0 *)

(* :Discussion:
	This package contains functions for studying right bijective and left bijective cellular automata.

	As in CellularAutomaton, rules are written in the form {n, k, r}, where n is the rule number,
	k is the number of colors, and r is the radius (so that 2 r + 1 is the number of cells on which the rule depends).

	Leftful initial conditions are represented as {listofcells, backgroundcolor}, as in CellularAutomaton.
	Rightful initial conditions are represented as {backgroundcolor, listofcells}.
*)


BeginPackage["BijectiveRules`"]

(* \[Lambda] *)
AgreementLength::usage =
"AgreementLength[list1, list2] gives the number of positions at the beginning of two lists for which they contain identical elements."

ApplyRule::usage =
"ApplyRule[{n, k, r}, cells] gives the color of the cell yielded by the sequence cells under the given rule."

ApplyRightBijectiveRuleInverse::usage =
"ApplyRightBijectiveRuleInverse[{n, k, r}, cells, rightcell] gives Append[Rest[cells], c], where c is such that Append[cells, c] yields rightcell under the given rule."

BijectiveQ::usage =
"BijectiveQ[rule, p] returns True if rule is bijective in position p and False otherwise."
LeftBijectiveQ::usage =
"LeftBijectiveQ[rule] returns True if rule is left bijective and False otherwise."
RightBijectiveQ::usage =
"RightBijectiveQ[rule] returns True if rule is right bijective and False otherwise."

(* \[Iota] *)
BorderBlockLength::usage =
"BorderBlockLength[list, backgroundcolor] gives the length of the first single-color block following the initial block of color backgroundcolor in list.
BorderBlockLength[list] uses background color 0."

ColorCycle::usage =
"ColorCycle[rule, c] gives the eventual cycle of tail colors obtained from tail color c under rule.
If c is present in this cycle, c appears as the first entry in the cycle; otherwise the cycle is given in its canonical rotation."

ColorCycles::usage =
"ColorCycles[rule] lists the possible cycles of tail colors under rule.
Each cycle is given in its canonical rotation."

(* a *)
ConvergenceSequence::usage =
"ConvergenceSequence[rule, init, power, l] gives the agreement lengths for rows offset by successive powers of l in the evolution of the cellular automaton with the specified rule from initial condition init.
ConvergenceSequence[rule, init, power] uses l = LCM[1, 2, ..., k]."

DependenceStrengths::usage =
"DependenceStrengths[rule] gives, for each position in the rule, the probability that changing the input in that position changes the output."

EquivalentRules::usage =
"EquivalentRules[{n, k, r}] gives the rule numbers of all k\[Hyphen]color, radius r rules that are equivalent to the given rule under a permutation of colors."

LeftBijectiveInverse::usage =
"LeftBijectiveInverse[rule] gives the radius 1/2 rule that computes the inverse of the given left bijective rule."
RightBijectiveInverse::usage =
"RightBijectiveInverse[rule] gives the radius 1/2 rule that computes the inverse of the given right bijective rule."

LeftBijectiveRules::usage =
"LeftBijectiveRules[k, r] lists the rule numbers of all left bijective k\[Hyphen]color, radius r rules."
RightBijectiveRules::usage =
"RightBijectiveRules[k, r] lists the rule numbers of all right bijective k\[Hyphen]color, radius r rules."
LeftOrRightBijectiveRules::usage =
"LeftOrRightBijectiveRules[k, r] lists the rule numbers of all k\[Hyphen]color, radius r rules that are either left bijective or right bijective."

LeftfulPredecessor::usage =
"LeftfulPredecessor[rule, {row, tailcolor}] gives the predecessor of a row with an infinite leftful history under the given left bijective rule."
RightfulPredecessor::usage =
"RightfulPredecessor[rule, {tailcolor, row}] gives the predecessor of a row with an infinite rightful history under the given right bijective rule."

LeftfulSuccessor::usage =
"LeftfulSuccessor[rule, {row, tailcolor}] gives the successor of a row under the given rule."
RightfulSuccessor::usage =
"RightfulSuccessor[rule, {tailcolor, row}] gives the successor of a row under the given rule."

LeftRightReflection::usage =
"LeftRightReflection[rule] gives the left\[Dash]right reflection of the given rule."

RandomLeftBijectiveRule::usage =
"RandomLeftBijectiveRule[k, r] produces a random left bijective k\[Hyphen]color, radius r rule."
RandomRightBijectiveRule::usage =
"RandomRightBijectiveRule[k, r] produces a random right bijective k\[Hyphen]color, radius r rule."

ReverseRow::usage =
"ReverseRow[list] gives the left\[Dash]right reflection of a row of the form {lefttailcolor, row} or {row, righttailcolor}."

ReversibleCellularAutomaton::usage =
"ReversibleCellularAutomaton[rule, init, t] generates a list representing the evolution of the reversible cellular automaton with the specified rule from initial condition init for t steps.
ReversibleCellularAutomaton[rule, init, {tspec}] gives only those parts of the evolution specified by tspec."
Padding::usage =
"Padding is an option for ReversibleCellularAutomaton that determines the number of columns of padding on the left or right."

Rule30ConvergenceData::usage =
"Rule30ConvergenceData[n] gives the number of rightmost consecutive black cells on row 2^n-1 of rule 30 begun from a single black cell."

RuleTable::usage =
"RuleTable[rule] gives the local rules specified by rule."


Unprotect["BijectiveRules`*"]

Begin["`Private`"]

CanonicalRotation[list_] := First[Sort[NestList[RotateLeft, list, Length[list] - 1]]]

Classify[list_, f_, g_ : Identity] :=
	Module[{i, vals = {}, val, classes = {}},
		For[i = 1, i <= Length[list], i++,
			If[# == {},
				AppendTo[vals, val]; AppendTo[classes, {g[list[[i]]]}], 
				AppendTo[classes[[#[[1, 1]]]], g[list[[i]]]]
			] & @ Position[vals, val = f[list[[i]]], {1}, 1]
		];
		Head[list] @@@ classes
	]

FixedPointPeriod[f_, expr_] :=
Module[{list = {expr}, val = expr, i = 0},
	While[!MemberQ[list, val = f[val], {1}],
		AppendTo[list, val];
		i++
	];
	Take[list, {Position[list, val, {1}, 1][[1,1]], -1}]
]

Shear[array_List, n_Integer : 1] := MapIndexed[RotateRight[#1, (#2[[1]] - 1) n] &, array, {1}]

ToArray[array_List, padding_] /; VectorQ[array[[1,1]]] :=
	Join[#[[1]], ConstantArray[#[[2]], padding]] & /@ array
ToArray[array_List, padding_] /; VectorQ[array[[1,2]]] :=
	Join[ConstantArray[#[[1]], padding], #[[2]]] & /@ array

ToRule[n_Integer?NonNegative /; n < 256] := {n, 2, 1}
ToRule[{n_Integer?NonNegative, k : (_Integer?Positive) : 2, r : (_?NonNegative) : 1} /; IntegerQ[2 r] && n < k^k^(2 r + 1)] := {n, k, r}

UnpadLeft[list_, p_] :=
	Drop[list, Replace[Position[list, _?(! MatchQ[#, p]&), {1}, 1, Heads -> False], {{{position_}} :> position, _ -> All + 1}] - 1]

ValidRuleQ[rule_] := MatchQ[ToRule[rule], {_, _, _}]


SetAttributes[Rule30ConvergenceData, Listable]
Rule30ConvergenceData[n_Integer /; 0 <= n <= 42] :=
{1, 3, 4, 6, 7, 9, 15, 16, 24, 25, 27, 29, 34, 36, 37, 39, 41, 43, 48, 49, 51,
	54, 55, 58, 60, 63, 64, 66, 69, 70, 72, 74, 77, 79, 80, 82, 84, 86, 90, 91,
	93, 100, 103}[[n + 1]]
Rule30ConvergenceData[n_Integer /; n > 42] := Missing["Unknown"]


AgreementLength[list1_List, list2_List] :=
Module[{p},
	LengthWhile[Transpose[PadRight[{list1, list2}, Automatic, p]], SameQ @@ # &]
]

ApplyRule[{n_, k_, r_}, cells_List] := IntegerDigits[n, k, k^(2 r + 1)][[k^(2 r + 1) - FromDigits[cells, k]]]

ApplyRightBijectiveRuleInverse[{n_, k_, r_}, cells_List, rcell_] := ApplyRightBijectiveRuleInverse[{n, k, r}, cells, rcell] =
Append[
	Rest[cells],
	First[Select[Range[0, k - 1], IntegerDigits[n, k, k^(2 r + 1)][[k^(2 r + 1) - FromDigits[Append[cells, #], k]]] == rcell &, 1]]
]

BijectiveQ[rule_?ValidRuleQ, position_Integer] :=
Module[{n, k, r},
	{n, k, r} = ToRule[rule];
	Length[Union[ReplacePart[#1, position -> #2] & @@@ RuleTable[rule]]] == k^(2 r + 1)
]
LeftBijectiveQ[rule_] := BijectiveQ[rule, 1]
RightBijectiveQ[rule_] := BijectiveQ[rule, -1]

BorderBlockLength[row_List, bgcolor_ : 0] :=
With[{unpaddedrow = UnpadLeft[row, bgcolor]},
	AgreementLength[unpaddedrow, ConstantArray[First[unpaddedrow], Length[unpaddedrow]]]
]

ColorCycle[rule_?ValidRuleQ, c_] :=
Module[{n, k, r, graph},
	{n, k, r} = ToRule[rule];
	graph = IntegerDigits[n, k, k^(2 r + 1)][[k^(2 r + 1) - ((k^(2 r + 1) - 1) Range[0, k - 1])/(k - 1)]];
	(If[MemberQ[#, c], #, CanonicalRotation[#]] &)[FixedPointPeriod[graph[[# + 1]] &, c]]
]

ColorCycles[rule_?ValidRuleQ] :=
Module[{n, k, r, graph},
	{n, k, r} = ToRule[rule];
	graph = IntegerDigits[n, k, k^(2 r + 1)][[k^(2 r + 1) - ((k^(2 r + 1) - 1) Range[0, k - 1])/(k - 1)]];
	Union[Table[CanonicalRotation[FixedPointPeriod[graph[[# + 1]] &, i]], {i, 0, k - 1}]]
]

ConvergenceSequence[rule_?ValidRuleQ, init : {_, _}, power_Integer?NonNegative, l_Integer?NonNegative] :=
Module[{padding = 1, row},
	row = ReversibleCellularAutomaton[rule, init, {{{0}}}, Padding -> padding];
	AgreementLength[row, #] - padding & /@
		FoldList[
			ReversibleCellularAutomaton[rule, {First[#1], Drop[#1, padding]}, {{{#2}}}, Padding -> padding] &,
			ReversibleCellularAutomaton[rule, init, {{{1}}}, Padding -> padding],
			l^(Range[power] - 1) (l - 1)
		]
]
ConvergenceSequence[rule_?ValidRuleQ, init_List, power_Integer?NonNegative] :=
	ConvergenceSequence[rule, init, power, LCM @@ Range[ToRule[rule][[2]]]]

DependenceStrengths[rule_?ValidRuleQ] :=
Module[{n, k, r, ruletable},
	{n, k, r} = ToRule[rule];
	ruletable = RuleTable[rule];
	Function[
		position, 
		Mean[(Length /@ Union /@ Classify[ruletable, Drop[First[#], {position}] &, Last] - 1) / (k - 1)]
	] /@ Range[2 r + 1]
]
        
EquivalentRules[rule_?ValidRuleQ] :=
Module[{n, k, r, ruletable},
	{n, k, r} = ToRule[rule];
	ruletable = RuleTable[rule];
	Union[Function[
		perm,
		FromDigits[Last /@ Reverse[Sort[ruletable /. Thread[Range[0, k - 1] -> perm]]], k]
	] /@ Permutations[Range[0, k - 1]]]
]

LeftBijectiveInverse[rule_?ValidRuleQ] /; LeftBijectiveQ[rule] :=
	LeftRightReflection[RightBijectiveInverse[LeftRightReflection[rule]]]
RightBijectiveInverse[rule_?ValidRuleQ] /; RightBijectiveQ[rule] := RightBijectiveInverse[rule] =
Module[{n, k, r},
	{n, k, r} = ToRule[rule];
	{FromDigits[
		FromDigits[#, k] & /@ Flatten[Table[
			Last /@ Rest[FoldList[
				ApplyRightBijectiveRuleInverse[{n, k, r}, ##] &,
				IntegerDigits[b1, k, 2 r],
				IntegerDigits[b2, k, 2 r]
			]],
			{b1, k^(2 r) - 1, 0, -1},
			{b2, k^(2 r) - 1, 0, -1}
		], 1],
		k^(2 r)
	], k^(2 r), 1/2}
]

LeftBijectiveRules[k_, r_] /; ValidRuleQ[{0, k, r}] :=
	Sort[
		First[LeftRightReflection[{#, k, r}]] & /@ RightBijectiveRules[k, r]
	]
RightBijectiveRules[k_, r_] /; ValidRuleQ[{0, k, r}] :=
	FromDigits[Join @@ #, k] & /@
		Tuples[Permutations[Range[0, k - 1]], k^(2 r)]
LeftOrRightBijectiveRules[k_, r_] /; ValidRuleQ[{0, k, r}] :=
	Union[Flatten[
		{#, First[LeftRightReflection[{#, k, r}]]} & /@ RightBijectiveRules[k, r]
	]]

LeftfulPredecessor[rule_?ValidRuleQ, {row_List, tail_Integer}] /;
		LeftBijectiveQ[rule] && MemberQ[ColorCycle[rule, tail], tail] :=
	ReverseRow[RightfulPredecessor[LeftRightReflection[rule], ReverseRow[{row, tail}]]]
(* "Last /@" here makes it inefficient for long lists; I should just Sow each cell *)
RightfulPredecessor[rule_?ValidRuleQ, {tail_Integer, row_List}] /;
		RightBijectiveQ[rule] && MemberQ[ColorCycle[rule, tail], tail] :=
Module[{n, k, r, newtail},
	{n, k, r} = ToRule[rule];
	newtail = Last[ColorCycle[rule, tail]];
	{newtail, Last /@ Rest[FoldList[
		ApplyRightBijectiveRuleInverse[{n, k, r}, ##] &,
		ConstantArray[newtail, 2 r],
		row
	]]}
]

LeftfulSuccessor[rule_?ValidRuleQ, {row_List, tail_Integer}] :=
	ReverseRow[RightfulSuccessor[LeftRightReflection[rule], ReverseRow[{row, tail}]]]
RightfulSuccessor[rule_?ValidRuleQ, {tail_Integer, row_List}] :=
Module[{n, k, r},
	{n, k, r} = ToRule[rule];
	{
		ApplyRule[{n, k, r}, ConstantArray[tail, 2 r + 1]], 
		ApplyRule[{n, k, r}, #] & /@ Partition[Join[ConstantArray[tail, 2 r], row], 2 r + 1, 1]
	}
]

LeftRightReflection[rule_?ValidRuleQ] :=
Module[{n, k, r},
	{n, k, r} = ToRule[rule];
	{FromDigits[Reverse[Last /@ Sort[{Reverse[#1], #2} & @@@ RuleTable[rule]]], k], k, r}
]

RandomLeftBijectiveRule[k_, r_] := LeftRightReflection[RandomRightBijectiveRule[k, r]]
RandomRightBijectiveRule[k_, r_] /; ValidRuleQ[{0, k, r}] :=
	{FromDigits[Join @@ Table[RandomSample[Range[0, k - 1]], {k^(2 r)}], k], k, r}

ReverseRow[{row_List, tail_Integer}] := {tail, Reverse[row]}
ReverseRow[{tail_Integer, row_List}] := {Reverse[row], tail}

evolveForward[{n_, k_, r_}, init : {row_List, _Integer}, tspec_] :=
	{Most[#], Last[#]} & /@
		CellularAutomaton[
			{n, k, List /@ Range[0, 2 r]},
			init,
			{tspec, {0, Length[row]}}
		]
evolveForward[{n_, k_, r_}, init : {_Integer, row_List}, tspec_] :=
	{First[#], Rest[#]} & /@
		CellularAutomaton[
			{n, k, List /@ Range[-2 r, 0]},
			Reverse[init],
			{tspec, {-1, Length[row] - 1}}
		]
evolveBackward[rule_?LeftBijectiveQ, init : {_List, _Integer}, tspec_] :=
	ReverseRow /@ evolveBackward[LeftRightReflection[rule], ReverseRow[init], tspec]
evolveBackward[rule : {n_, k_, r_}?RightBijectiveQ, init : {bg_Integer, row_List}, {t1_?Negative, t2_?NonPositive, dt_}] /;
		MemberQ[ColorCycle[rule, bg], bg] :=
With[{inversewidth = Ceiling[(Length[row] + 1)/(2 r)]},
	{First[#], Take[#, -Length[row]]} & /@ Flatten /@
		IntegerDigits[Take[Transpose[Shear[Transpose[Reverse[CellularAutomaton[
			Append[Most[RightBijectiveInverse[rule]], {{-1}, {0}}],
			{
				First[(Transpose[Shear[Transpose[((FromDigits[#, k] &) /@ Partition[#, 2 r] &) /@ #], -1]] &)[
					CellularAutomaton[
						{n, k, List /@ Range[-2 r, 0]},
						Reverse[init],
						inversewidth - 1,
						{All, {Length[row] - 2 r inversewidth, Length[row] - 1}}
					]
				]],
				((FromDigits[#, k] &) /@ Partition[Flatten[Table[#, {LCM[Length[#], 2 r]/Length[#]}]], 2 r] &)[
					Flatten[ConstantArray[#, 2 r] & /@ ColorCycle[rule, bg]]
				]
			},
			{inversewidth - 1 - t1, {0, inversewidth - 1}}
		]]]]], {t1 - 1, t2 - 1, dt}], k, 2 r]
]
(* all steps 0 through t *)
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, t_Integer?Negative | {t_Integer?Negative}, OptionsPattern[]] :=
	ToArray[evolveBackward[ToRule[rule], init, {t, 0, 1}], OptionValue[Padding]]
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, t_Integer?NonNegative | {t_Integer?NonNegative}, OptionsPattern[]] :=
	ToArray[evolveForward[ToRule[rule], init, t], OptionValue[Padding]]
(* a list containing only step t *)
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{t_Integer?NonPositive}}, OptionsPattern[]] :=
	ToArray[evolveBackward[ToRule[rule], init, {t, 0, -t}], OptionValue[Padding]]
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{t_Integer?Positive}}, OptionsPattern[]] :=
	ToArray[Rest[evolveForward[ToRule[rule], init, {0, t, t}]], OptionValue[Padding]]
(* step t alone *)
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{{t_Integer?Negative}}}, OptionsPattern[]] :=
	First[ToArray[evolveBackward[ToRule[rule], init, {t, 0, -t}], OptionValue[Padding]]]
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{{0}}}, OptionsPattern[]] :=
	First[ToArray[{init}, OptionValue[Padding]]]
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{{t_Integer?Positive}}}, OptionsPattern[]] :=
	Last[ToArray[evolveForward[ToRule[rule], init, {0, t, t}], OptionValue[Padding]]]
(* steps t1, t1 + dt, ... *)
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{t1_Integer?Negative, t2_Integer?NonPositive, dt : (_Integer?Positive) : 1}} /; t1 <= t2, OptionsPattern[]] :=
	ToArray[evolveBackward[ToRule[rule], init, {t1, t2, dt}], OptionValue[Padding]]
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{t1_Integer?Negative, t2_Integer?Positive, dt : (_Integer?Positive) : 1}} /; t1 <= t2, OptionsPattern[]] :=
	ToArray[Join[
		evolveBackward[ToRule[rule], init, {t1, -1, dt}],
		If[Mod[t1, dt] <= t2,
			evolveForward[ToRule[rule], init, {Mod[t1, dt], t2, dt}],
			{}
		]
	], OptionValue[Padding]]
ReversibleCellularAutomaton[rule_?ValidRuleQ, init_List, {{t1_Integer?NonNegative, t2_Integer?Positive, dt : (_Integer?Positive) : 1}} /; t1 <= t2, OptionsPattern[]] :=
	ToArray[evolveForward[ToRule[rule], init, {t1, t2, dt}], OptionValue[Padding]]
Options[ReversibleCellularAutomaton] = {Padding -> 0}

RuleTable[rule_?ValidRuleQ] :=
Module[{n, k, r},
	{n, k, r} = ToRule[rule];
	Transpose[{
		Tuples[Range[0, k - 1], 2 r + 1], 
		Reverse[IntegerDigits[n, k, k^(2 r + 1)]]
	}]
]


End[]

Protect["BijectiveRules`*"]
Unprotect[ApplyRightBijectiveRuleInverse, RightBijectiveInverse]

EndPackage[]