(* Function to find attractors and ALL basins. *)
(* Usage: BartCA[105,3] *)

BartCA[rule_,bitLength_] := 
	Module[{n,l,DynamicsList,p,initPattern,x}, 
		    (* Declare some utility functions first *)
		    (* The function below converts the input in decimal to a *)
	           (* list of its binary equivalent suitable for use in CellularAutomaton *)
		    DecimalToBinaryList[n_] := IntegerDigits[n, 2, bitLength];
		    		   
		    (* Convert a binary number (represented as a list) to a decimal integer *)
		    BinaryToDecimal[l_] := FromDigits[l, 2];
		    
		    (* The function below can be used to convert the output from CellularAutomaton to integers *)
		    BinaryToDecimalEvolution[DynamicsList_] := Map[BinaryToDecimal, DynamicsList];

		    (* Functions below puts the output in the format we want: {basin,{attractor}} OR {attractor} *)
		    (* getEvolutionListAndAttractorPositions works because the last element in each evolution list is *)
		    (* repeated again somewhere in the evolution list.  Hence, if we will get an array of two elements if we *)
		    (* find the positions of the last element in each evolution list.  *)
		    getEvolutionListAndAttractorPositions[l1_] := Join[{l1}, {Flatten[Position[l1, Last[l1]]]}];

		    (* Now, lets say we have the following evolution list and attractor positions in the list: *)
		    (* LP = {{{0,0},{1,2}},{{1,3,7,13,7},{3,5}}} *)
		    (* We want: {{{},{0}},{{1,3},{7,13,7}}} *)
	          evolutionFormat[LP_] := {Flatten[Take[LP[[1]],LP[[2]][[1]]-1]],Take[LP[[1]],LP[[2]]]};

	           (* Function below is obvious: it computes one iteration of an evolution rule *)
		   (* by running CellularAutomaton once *)
		    OneIteration[p_] := Last[CellularAutomaton[rule,p,1]];
	           (*  Run until you hit an attractor using the NestWhileList function *)
		   CARunUntilAttractor[initPattern_] := NestWhileList[OneIteration,initPattern, UnsameQ, All];

		    x = Table[i,{i,0,(2^bitLength)-1}];

		    Map[evolutionFormat,Map[getEvolutionListAndAttractorPositions,Map[BinaryToDecimalEvolution,Map[CARunUntilAttractor,(Map[DecimalToBinaryList,x])]]]]];


(* Function to find unique attractors *)
(* Format: {{attractor 1},{attractor 2},...} *)
(* Example: {{0},{3,12},{6,9}} *)
(* Usage:  GetUniqueAttractors[BartCA[105,3]]*)

GetUniqueAttractors[evolutionList_] :=
	Module[{tempList},
		(* Utility Function to check for empty list *)
		CheckIfEmptyQ[l_] := Length[l[[1]]] == 0;

		(* Mathematica example function to take union of elements WITHOUT sorting *)
		UnsortedUnion[x_]:=Module[{f},f[y_]:=(f[y]=Sequence[];y);f/@x];

		(* This utility function finds the position of the smallest number in each attractor *)
		(* We use this as a predicate in NestWhile to rotate the attractor until the smallest *)
		(* number is the head of the attractor *)
		FirstPosSmallestQ[l_] := Position[l, Min[l]] != {{1}};

		RotateUntil[l_] := NestWhile[RotateLeft, l, FirstPosSmallestQ];

		(* If the first element is empty, then we know that we have an attractor *)
		(* Then, we just extract the second element from each list, that is the attractor. *)
		(* Next, we flatten each element in the list for nicer formatting *)
		(* By taking an UnsortedUnion, we remove the last duplicate entry in each attractor *)
		(* We then sort the list so that attractors of the same period are grouped *)
		(* Then, we just use RotateUntil so that we can finally take the Union to get the unique attractors *)		
		
		Union[Map[RotateUntil,Sort[Map[UnsortedUnion,Map[Flatten,Map[Rest,Select[evolutionList,CheckIfEmptyQ]]]]]]]];

(* Function to find attractor periods *)
(* Output format:  {{no-of-attractors of period p1, {Period, p1}},{no-of-attractors of period p2,{Period, p2}},...*)
(* Example:  {{1, {Period, 1}},{2,{Period, 2}}} *)
(* Usage: FindAttractorPeriods[GetUniqueAttractors[BartCA[105,3]]] *)

FindAttractorPeriods[UniqueAttractorList_] := 
	Module[{},
	OutputFormat[l_] := {Length[l],{Period,First[l]}};
	Map[OutputFormat,Split[Map[Length,UniqueAttractorList]]]];

(* Function to find if we only have Isles of Eden in the evolution list *)
(* Output:  True or False *)
(* Usage: AllIslesOfEdenQ[BartCA[105,3]] *)

AllIslesOfEdenQ[evolutionList_] :=
	Module[{},
	(* Utility Function to check if transient is empty *)
	CheckIfEmptyQ[l_] := Length[l[[1]]] == 0;
	
	(* First apply function above to every initial condition *)
	(* Next take Union of the result *)
	(* If we only have Isles of Eden, then the result would be the singleton set {True} *)
	(* else it will be {True,False} (or {False,True}) *)
	Union[Map[CheckIfEmptyQ,evolutionList]] == {True}];
	
BartCAInitialCondition[rule_,bitLength_,init_,t_]:=Module[{x},
      (*Declare some utility functions first*)
      (*The function below converts the input in decimal to a*)
      (*list of its binary equivalent suitable for use in CellularAutomaton*)
      DecimalToBinaryList[n_]:=IntegerDigits[n,2,bitLength];
      (*Convert a binary number (represented as a list) to a decimal integer*)
      BinaryToDecimal[l_]:=FromDigits[l,2];
      Map[BinaryToDecimal,CellularAutomaton[rule,DecimalToBinaryList[init],t]]];

(* FUNCTION BELOW DOES NOT WORK YET *)

BartCARandomInitialCondition[rule_,bitLength_] := 
	Module[{x}, 
 	        (* Declare some utility functions first *)
		(* The function below converts the input in decimal to a *)
	        (* list of its binary equivalent suitable for use in CellularAutomaton *)
		DecimalToBinaryList[n_] := IntegerDigits[n, 2, bitLength];
		    		   
		(* Convert a binary number (represented as a list) to a decimal integer *)
		BinaryToDecimal[l_] := FromDigits[l, 2];
		    
		(* Function below is obvious: it computes one iteration of an evolution rule *)
		(* by running CellularAutomaton once *)
		OneIteration[p_] := Last[CellularAutomaton[rule,p,1]];
	        (*  Run until you hit an attractor using the NestWhileList function *)
		CARunUntilAttractor[initPattern_] := NestWhileList[OneIteration,initPattern, UnsameQ, All];

		(* Use Random[Integer,{0,2^bitLength-1}] to generate psuedo-random integer in the given range*)
		x = Random[Integer,{0,2^bitLength-1}];

		Map[BinaryToDecimal,CARunUntilAttractor[DecimalToBinaryList[x]]]];