(* ::Package:: *)

(* :Title: Digraph Fractals *)
(* :Context: FractalGeometry`DigraphFractals` *)
(* :Author: Mark McClure *)
(* :Copyright: Copyright 1999, 2003, 2008 Mark McClure *)
(* :Package Version: 3.0 *)
(* :Mathematica Version: 6.0 *)
(* :Summary:
    Generates images of digraph self-similar and
    digraph self-affine sets.
*)
(* :Keywords: fractal, self-affine, self-similar, 
    digraph self-affine, digraph self-similar,
    mixed self-affine, mixed self-similar,
    digraph IFS
*)


(* :History: Original version written 1999. 

  Version 3.0 January 2008
  Now V6 compatible.  In particular, mult-point
  objects are used for efficiency.

  Version 2.0 October 2003. 
  Major enhancements include:
    * Moved to the FractalGeometry` context.
      Primarily to allow interaction with the
      SelfAffineTiles package.
    * Added the Colors option to 
      ShowDigraphFractalsStochastic.
    * Added the StronglyConnectedDigraphQ function.
*)


(* :Sources:
  Bandt, C. 1989. "Self-similar sets III. Constructions with
  sofic systems."  Monatsh. Math. 108:89-102.
  
  Mauldin, R.D. and Williams, S.C., 1988.  "Hausdorff dimension 
  in graph directed constructions." Trans. Amer. Math. Soc. 
  309:811-829.
  
  McClure, M. 2000. "Directed-graph iterated function systems."
  Mathematica in Education and Research. 9(2).  

  McClure, M. 2008. "Fractals 6.0."  MiER 13(1).
*)


BeginPackage["FractalGeometry`DigraphFractals`",
  "FractalGeometry`Common`"]


DigraphFractals::usage = "DigraphFractals is a package defining \
several functions used to generate images of digraph Self \
Similar sets."

ShowDigraphFractals::usage = "ShowDigraphFractals[digraph_,depth_] \
generates the approximation to the digraph self similar sets \
defined by digraph to order depth using a deterministic algorithm. \
The digraph is represented by a matrix of lists of affine functions \
written in the form {{{a1,a2},{b1,b2}},{x0,y0}}."

ShowDigraphFractalsStochastic::usage = \
"ShowDigraphFractalsStochastic[digraph_,numPoints_] \
generates numPoints points approximating the digraph Self Similar \
sets defined by digraph using a stochastic algorithm. \
The digraph is represented by a matrix of lists of affine functions \
written in the form {{{a1,a2},{b1,b2}},{x0,y0}}."

ShowDigraph::usage = \
"ShowDigraph[digraphIFS] uses GraphPlot show an image of the digraph \
corresponding to digraphIFS."

ComputePMatrix::usage = "ComputePMatrix[digraph_] estimates a \
probability matrix for digraph which generates a uniform distribution \
of points over each digraph Fractal."

FindDigraphDimension::usage = "FindDigraphDimension[digraph_] computes a \
numerical estimate of the dimension of the corresponding digraph \
fractals."

StronglyConnectedDigraphQ::usage = "StronglyConnectedDigraphQ[digraph_] \
test the digraph to see if it is strongly connected."

Initiators::usage = "Initiators is an option for ShowDigraphFractals \
indicating an initial list of Graphics primitives for the digraph \
to operate on.  Initiators should be a list of lists of Graphics \
primitives with length equal to the dimension of the digraph."

PMatrix::usage = "PMatrix is an option for \
ShowDigraphFractalsStochastic indicating a matrix of probabilities \
for the functions in the digraph. in \
ShowDigraphFractalsStochastic[digraph_, numPoints_], PMatrix->pMatrix], \
pMatrix should be a matrix of lists of positive numbers with the same \
shape as digraph.  The sum of the sums of the lists in any column of \
PMatrix should be one."

Begin["`Private`"]

(* Set the Options *)

(* For ShowDigraphFractals *)
Options[ShowDigraphFractals] = Join[{Initiators -> Automatic, 
  PlotStyle -> {}}, Options[Graphics]];

(* For ShowDigraphFractalsStochastic *)
Options[ShowDigraphFractalsStochastic] = Join[{PMatrix -> Automatic, 
  PlotStyle -> {AbsolutePointSize[0.4]}, Colors -> False}, Options[Graphics]];

(* The Functions *)

(* Error Messages *)
DigraphFractals::optx = "Unknown Option `1`."
DigraphFractals::badDigraph = "The first argument in `1` must be a
matrix of lists of affine lists."
DigraphFractals::badInt = "The second argument in `1` must be a
non-negative integer."
DigraphFractals::badOpt = "Options expected as optional arguments in `1`."

Options[ShowDigraph] = Options[GraphPlot];
ShowDigraph[digraphIFS_, opts___] := Module[
  {adj},
  adj = Map[Length, digraphIFS, {2}];
  GraphPlot[adj, opts,
    SelfLoopStyle -> 0.2, MultiedgeStyle -> 0.3,
    DirectedEdges -> True, VertexLabeling -> True]];

ShowDigraphFractals[digraphIFS_, depth_, opts___] := Module[
  {initiators, plotStyle, toFunc, toGraphicFunction,
  step, funcMatrix1, funcMatrix2, nestedFuncMatrix, 
  gfMatrix, attractors,h,x},

  initiators=Initiators/.{opts}/.Options[ShowDigraphFractals];
  plotStyle=PlotStyle/.{opts}/.Options[ShowDigraphFractals];
  If[initiators===Automatic,
    initiators=Table[{Point[{0,0}]},{Length[digraphIFS]}]];

  toFunc[{A_,b_}]:=Module[{cfOut,fOut},
    cfOut=Compile[{{v,_Real,1}},A.v+b];
    fOut[{x_?NumericQ,y_?NumericQ}]:=cfOut[{x,y}];
    fOut[Point[{x_?NumericQ, y_?NumericQ}]]:=Point[{fOut[{x,y}]}];
    fOut[Point[pts:{{_?NumericQ,_?NumericQ}..}]] := Point[fOut/@pts];
    fOut[Line[pts:{{_?NumericQ,_?NumericQ}..}]]:=Line[{fOut/@pts}];
    fOut[Line[pts:{{{_?NumericQ,_?NumericQ}..}..}]]:=Line[Map[fOut,pts,{2}]];
    fOut[Arrow[pts:{{_?NumericQ,_?NumericQ}..},s___]]:=Arrow[fOut/@pts,s];
    fOut[Polygon[pts:{{_?NumericQ,_?NumericQ}..}, vc___]]:=
      Polygon[{fOut/@pts}, vc];
    fOut[Polygon[pts:{{{_?NumericQ,_?NumericQ}..}..}, vc___]]:=
      Polygon[Map[fOut,pts,{2}], vc];
    fOut[x_] := x;
    fOut];

  toGraphicFunction[fs_List] := Module[
    {gfOut},
    gfOut[Point[pt:{_?NumericQ,_?NumericQ}]] := 
      Point[Through[fs[pt]]];
    gfOut[Point[pts:{{_?NumericQ,_?NumericQ}..}]] := 
      Point[Flatten[Through[fs[#]]& /@ pts,1]];
    gfOut[Line[pts:{{_?NumericQ,_?NumericQ}..}]] := 
      Line[Table[f/@pts,{f,fs}]];
    gfOut[Arrow[pts:{{_?NumericQ,_?NumericQ}..},s___]] := 
      Table[Arrow[f/@pts,s],{f,fs}];
    gfOut[Line[pts:{{{_?NumericQ,_?NumericQ}..}..}]] := 
      Line[Flatten[Table[Map[f,pts,{2}],{f,fs}],1]];
    gfOut[Polygon[pts:{{_?NumericQ,_?NumericQ}..}]] := 
      Polygon[Table[f/@pts,{f,fs}]];
    gfOut[Polygon[pts:{{{_?NumericQ,_?NumericQ}..}..}]] := 
      Polygon[Flatten[Table[Map[f,pts,{2}],{f,fs}],1]];
    gfOut[ll_List] := gfOut /@ ll;
    gfOut[h[x_]] := h[gfOut[x]];
    gfOut[x_] := x;
    gfOut];

  funcMatrix1 = Map[toFunc, digraphIFS, {3}];
  funcMatrix2 = Apply[h, funcMatrix1, {2}];
  step[matrix_] := Map[Flatten,Inner[Tuples[Composition[##]]&, 
    funcMatrix2, Apply[h, matrix, {2}], List],{2}];
  If[depth > 0, 
    nestedFuncMatrix = Nest[step,funcMatrix1,depth-1];
    gfMatrix = Map[toGraphicFunction,nestedFuncMatrix,{2}];
    attractors = Inner[#1[#2]&, gfMatrix, h/@initiators, List];
    attractors = attractors /. h[x_] -> x,
    attractors = initiators];
  Graphics[{plotStyle,#}, Sequence@@Flatten[{
    FilterRules[{opts}, Options[Graphics]],
    FilterRules[Options[ShowDigraphFractals],
    Options[Graphics]]}]]&/@attractors] /; 
(And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, digraphIFS, {3}]] ||
    Message[DigraphFractals::badDigraph, ShowDigraphFractals] )  &&
  ( (IntegerQ[depth] && depth >= 0 ) ||
    Message[DigraphFractals::badInt, ShowDigraphFractals, depth] ) &&
  ( And @@ Map[OptionQ, {opts}] ||
      Message[DigraphFractals::badOpt,ShowDigraphFractals] );


ShowDigraphFractalsStochastic[digraph_, numPoints_, opts___] :=
  Module[{matrices, eigenvalues, s, sMatrix, cnt, point,
    spectralRadius, spectralRadius0, pMatrix, pointFunction,
    approximateDim, dimMatrix, perronNumbers,
    pMatrix1, pMatrix2, pMatrixNormalizer, i, j,
    toFuncs, digraphFuncs, v, v1, v2, pointColors,
    currentPoint, points, choose, colors, plotStyle,
    valid = First /@ Options[ShowDigraphFractalsStochastic]},

  Scan[If[!MemberQ[valid, First[#]],
    Message[DigraphFractals::optx, ToString[First[#]]]]&,
    Flatten[{opts}]
  ];
  
  colors = Colors /. {opts} /. Options[ShowDigraphFractalsStochastic];
  plotStyle = PlotStyle /. {opts} /. Options[ShowDigraphFractalsStochastic];
  pMatrix = PMatrix /. {opts} /. Options[ShowDigraphFractalsStochastic];

  Which[
    colors === Automatic || colors === True,
      colors =  Hue[#,.7,.7]& /@ 
        Range[0.,1 - 1./Length[digraph], 1./Length[digraph]],
    Head[colors] === String,
      colors = ColorData[colors] /@  
        Range[0.,1 - 1./Length[digraph], 1./Length[digraph]],
    colors =!= None && colors =!= False && Head[colors] =!= List,
      colors = colors /@
        Range[0.,1 - 1./Length[digraph], 1./Length[digraph]]];


  (* Computing the Probabilities *)  
  pMatrixNormalizer[subPList_List] := If[
    Length[subPList] > 0,
      FoldList[Plus, 0, subPList]/Plus @@ subPList,
      subPList];

  If[pMatrix === Automatic,
      pMatrix = ComputePMatrix[digraph]];
    
  pMatrix = N[pMatrix];
  pMatrix2 = Transpose[Map[pMatrixNormalizer, pMatrix, {2}]];
  pMatrix1 = Map[Plus @@ # & , Transpose[pMatrix], {2}];
  pMatrix1 = FoldList[Plus, 0, #] & /@ pMatrix1;
  
  toFuncs[{A_,b_}] := Compile[{{v,_Real,1}},
    A.v + b];
  digraphFuncs = Map[toFuncs, N[digraph], {3}];
  v = 1; v1 = pMatrix1[[v]]; v2 = pMatrix2[[v]];
  currentPoint = {0., 0.};
  points = Table[{}, {Length[digraph]}];
  choose[v1_, v2_] := Module[
    {choice1, choice2, chooser},
    chooser = RandomReal[];
    choice1 = Length[Select[v1, (# < chooser) &]];
    chooser = RandomReal[];
    choice2 = Length[Select[v2[[choice1]], (# < chooser) &]];
    {choice1, choice2}];

  Do[{
    {i, j} = choose[v1, v2];
    currentPoint = 
      digraphFuncs[[i, v, j]][currentPoint];
    points[[i]] = {points[[i]], {v, currentPoint}};
    v = i;  v1 = pMatrix1[[v]];  v2 = pMatrix2[[v]];
    }, {numPoints}];
  If[colors === False || colors === None,
    points = points /. {v_Integer, p_} :>  p;
    points = Point/@(Partition[Flatten[#],2]& /@ points),
    points = Transpose[Table[{colors[[i]],Point[Last/@
      Cases[#,{i,{_,_}},Infinity]]}& /@ points,
      {i,1,Length[colors]}]]];

  Show[Graphics[{Sequence@@Flatten[{plotStyle}], #}], 
    FilterRules[{opts}, Options[Graphics]],
    FilterRules[Options[ShowDigraphFractalsStochastic], Options[Graphics]]] & /@ points
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, digraph, {3}]] ||
    Message[DigraphFractals::badDigraph, 
      ShowDigraphFractalsStochastic] )  &&
  ( (IntegerQ[numPoints] && numPoints >= 0 ) ||
    Message[DigraphFractals::badInt, 
      ShowDigraphFractalsStochastic, numPoints] ) &&
  ( And @@ Map[OptionQ, {opts}] ||
      Message[DigraphFractals::badOpt,ShowDigraphFractalsStochastic] )

ComputePMatrix[digraph_] := Module[
  {matrices, eigenvalues, s, sMatrix,
   spectralRadius, spectralRadius0, eigenMatrix,
   approximateDim, dimMatrix, perronNumbers,
   pMatrix1, pMatrix2, pMatrixNormalizer,
   a,b,c,d},
    
    pMatrixNormalizer[subPList_List] := If[
      Length[subPList] > 0,
        subPList/Plus @@ subPList,subPList];
    matrices = Map[First, Transpose[digraph], {3}] // N;
    eigenMatrix = matrices /. 
      {{a_Real, b_Real}, {c_Real, d_Real}} -> 
      Max[Abs /@ Eigenvalues[{{a, b}, {c, d}}]]^s;
    sMatrix = Map[Plus @@ # &, eigenMatrix, {2}] // N;
    eigenvalues = Eigenvalues[sMatrix];
    spectralRadius0 = Max[Chop[Abs@N[eigenvalues /. s -> 0]]];
    spectralRadius = Select[eigenvalues, 
      ((# /. s -> 0) == spectralRadius0) &][[1]];
    approximateDim = Chop[s /.
      FindRoot[spectralRadius == 1, {s, 1, 2}]];
    dimMatrix = sMatrix /. s -> approximateDim;
    perronNumbers = Eigensystem[dimMatrix][[2]][[1]];
    pMatrix1 = Inner[Times, dimMatrix,
      perronNumbers, List]/perronNumbers;
    pMatrix2 = eigenMatrix /. s -> approximateDim;
    pMatrix2 = Map[pMatrixNormalizer, pMatrix2, {2}];
    Transpose[pMatrix1 pMatrix2]
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, digraph, {3}]] ||
    Message[DigraphFractals::badDigraph, 
      ComputePMatrix] )


FindDigraphDimension[digraph_] := Module[
  {matrices, eigenvalues, s, sMatrix,
    spectralRadius, spectralRadius0,
    approximateDim, eigenMatrix, a, b, c, d},

  matrices = Map[First, Transpose[digraph], {3}] // N;
  eigenMatrix = matrices /. 
    {{a_Real, b_Real}, {c_Real, d_Real}} -> 
    Max[Abs /@ Eigenvalues[{{a, b}, {c, d}}]]^s;
  sMatrix = Map[Plus @@ # &, eigenMatrix, {2}] // N;
  eigenvalues = Eigenvalues[sMatrix];
  spectralRadius0 = Max[Chop[Abs@N[eigenvalues /. s -> 0]]];
  spectralRadius = Select[eigenvalues, 
    ((# /. s -> 0) == spectralRadius0) &][[1]];
  approximateDim = Chop[s /.
    FindRoot[spectralRadius == 1, {s, 1, 2}]]
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, digraph, {3}]] ||
    Message[DigraphFractals::badDigraph, 
      FindDigraphDimension] )

StronglyConnectedDigraphQ[digraph_] := Module[
  {numericDigraph, k},
  
  numericDigraph = Map[Length, digraph, {2}];
  Length[Cases[
    Sum[MatrixPower[numericDigraph, k],
          {k, 1, Length[numericDigraph]}], 0, {2}]] == 0
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, digraph, {3}]] ||
    Message[DigraphFractals::badDigraph, 
      StronglyConnectedDigraphQ] )

End[]  (* End Private Context *)

Protect[DigraphFractals, ShowDigraphFractals, 
  ShowDigraphFractalsStochastic, ShowDigraph,
  ComputePMatrix, FindDigraphDimension,
  Initiators, PMatrix]

EndPackage[]
