(* ::Package:: *)

(* :Name: IteratedFunctionSystems` *)
(* :Author: Mark McClure, Feb 2002 *)
(* :Copyright: Copyright Mark McClure, 2002-7 *)
(* :Package Version: 2.0 *)
(* :Mathematica Version: 6.0 *)
(* :Summary:
  Generates images of self affine sets. 
*)

BeginPackage["FractalGeometry`IteratedFunctionSystems`",
  (* "Utilities`FilterOptions`", *)
  "FractalGeometry`Common`"]

IteratedFunctionSystems::usage = "IteratedFunctionSystems is a package \
defining several functions used to generate images of self-affine sets."

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

ShowIFSStochastic::usage = "ShowIFSStochastic[IFS_,numPoints_] generates \
numPoints points approximating the self similar sets defined by IFS using a \
stochastic algorithm.  The IFS is represented by a list of affine functions \
written in the form {{{a1,a2},{b1,b2}},{x0,y0}}."

FindProbabilities::usage = "FindProbabilities[IFS_] estimates a \
probability list for IFS which generates a uniform distribution \
of points over the self-affine set determined by IFS."

FindIFSDimension::usage = "FindDimension[IFS_] computes a numerical \
estimate of the dimension of the corresponding self-affine set."

(*
Initiator::usage = "Initiator is an option for ShowIFS \
indicating an initial list of Graphics primitives for the IFS \
to operate on."
*)

Probabilities::usage = "Probabilities is an option for ShowIFSStochastic \
indicating a list of probabilities for the functions in the IFS. In \
ShowIFSStochastic[IFS_, numPoints_, \
Probabilities->probabilities], probabilities should be a list of positive \
numbers whose sum is one."

SierpinskiPedalTriangleIFS::usage = "SierpinskiPedalTriangleIFS[{A_,B_,C_}] 
returns an IFS to generate the Sierpinski pedal triangle with vertices \
A, B, and C."

Begin["`Private`"]

(* Set the Options *)
(* For ShowIFS *)
ShowIFSOptions = Join[{Initiator -> Point[{0,0}], 
  Colors -> False, PlotStyle -> {}},
  Options[Graphics]];
keywords = Complement[Union[First /@ ShowIFSOptions], 
  {DisplayFunction, DefaultFont}];
vals = keywords /.  ShowIFSOptions;
special = {
  DisplayFunction :> $DisplayFunction, 
  DefaultFont :> $DefaultFont};
Options[ShowIFS] = 
  Union[Apply[Rule,Transpose[{keywords,vals}],{1}], special];

Options[ShowIFS] = Union[{AspectRatio -> Automatic,
  Initiator -> Point[{0,0}], Colors -> False, 
  PlotStyle -> {}}, Options[Graphics],
  SameTest -> (First[#1]===First[#2]&)];

(* For ShowIFSStochastic *)
ShowIFSStochasticOptions = Join[{AspectRatio -> Automatic, Axes -> False,
  Colors -> False,
  Probabilities -> Automatic, PlotStyle -> AbsolutePointSize[0.005]}, 
  Options[Graphics]];
keywords = Complement[Union[First /@ ShowIFSStochasticOptions], 
  {DisplayFunction, DefaultFont}];
vals = keywords /.  ShowIFSStochasticOptions;
special = {DisplayFunction :> $DisplayFunction, DefaultFont :> $DefaultFont};
Options[ShowIFSStochastic] = 
  Union[Apply[Rule,Transpose[{keywords,vals}],{1}], special];

Options[ShowIFSStochastic] = Join[{Colors -> False,
  Probabilities -> Automatic, PlotStyle -> AbsolutePointSize[0.4]},
  Options[Graphics]];

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


(* The Functions *)
ShowIFS[IFS_, depth_Integer?(# >= 0 &), opts___] := Module[
   {initiator, plotStyle, colors, toFunc, funcs, F,
    attractor, at, x},
   initiator = Initiator /. {opts} /. Options[ShowIFS];
   plotStyle = PlotStyle /. {opts} /. Options[ShowIFS];
   colors = Colors /. {opts} /. Options[ShowIFS];
   
   Which[colors === Automatic || colors === True,
    colors = 
     ColorData["Rainbow"] /@ 
      Range[0., 1 - 1./Length[IFS], 1./Length[IFS]],
    Head[colors] === String,
    colors = 
     ColorData[colors] /@ 
      Range[0., 1 - 1./Length[IFS], 1./Length[IFS]],
    colors =!= None && colors =!= False && Head[colors] =!= List,
    colors = 
     colors /@ Range[0., 1 - 1./Length[IFS], 1./Length[IFS]]];
   
   toFunc[{A_, b_}] := Module[{cfOut, fOut},
     cfOut = Compile[{{v, _Real, 1}}, A.v + b];
     fOut[{x_?NumericQ, y_?NumericQ}] := cfOut[{x, y}];
     fOut[x_List] := fOut /@ x;
     fOut[Point[pts_]] := Point[fOut[pts]];
     fOut[Line[x_]] := Line[fOut[x]];
     fOut[Arrow[x_]] := Arrow[fOut[x]];
     fOut[Polygon[x_, pOpts___]] := Polygon[fOut[x], pOpts];
     fOut[x_] := x;
     fOut];
   
   funcs = toFunc /@ IFS;
   F[Point[pt : {_?NumericQ, _?NumericQ}]] := 
    Point[Table[f[pt], {f, funcs}]];
   F[Point[pts : {{_?NumericQ, _?NumericQ} ..}]] :=
    Point[Flatten[Table[f /@ pts, {f, funcs}], 1]];
   F[Line[pts : {{_?NumericQ, _?NumericQ} ..}]] :=
    Line[Table[f /@ pts, {f, funcs}]];
   F[Line[pts : {{{_?NumericQ, _?NumericQ} ..} ..}]] :=
    Line[Flatten[Table[Map[f, pts, {2}], {f, funcs}], 1]];
   F[Arrow[pts : {{_?NumericQ, _?NumericQ} ..}, s___]] :=
    Table[Arrow[f /@ pts, s], {f, funcs}];
   F[Polygon[pts : {{_?NumericQ, _?NumericQ} ..}]] :=
    Polygon[Table[f /@ pts, {f, funcs}]];
   F[Polygon[pts : {{{_?NumericQ, _?NumericQ} ..} ..}]] :=
    Polygon[Flatten[Table[Map[f, pts, {2}], {f, funcs}], 1]];
   F[Polygon[pts : {{_?NumericQ, _?NumericQ} ..}, 
      VertexColors -> vc_]] :=
    Polygon[Table[f /@ pts, {f, funcs}], 
     VertexColors -> Table[vc, {Length[funcs]}]];
   F[Polygon[pts : {{{_?NumericQ, _?NumericQ} ..} ..}, 
      VertexColors -> vc_]] :=
    Polygon[Flatten[Table[Map[f, pts, {2}], {f, funcs}], 1],
     VertexColors -> Flatten[Table[vc, {Length[funcs]}], 1]];
   F[ll_List] := F /@ ll;
   F[x_] := x;
   If[colors =!= False && colors =!= None,
    attractor = Nest[F, initiator, depth - 1];
    attractor = at /@ Table[f[attractor], {f, funcs}];
    attractor = Inner[List, colors, attractor, List] /. 
      at[x__] -> x,
    attractor = Nest[F, initiator, depth]];
   Graphics[attractor,
    FilterRules[{opts}, Options[Graphics]],
    FilterRules[Options[ShowIFS], Options[Graphics]]]
   ];


ShowIFS[IFS_,tol_Real,opts___] := Module[
  {toFunc, funcs, attractor, color, colors,
   init, initiator, coloredAttractor, x, plotStyle,
   twoNorms, extend, f, ff, extendGPrimitives, extendedFuncs},
  initiator = Initiator /. {opts} /. Options[ShowIFS];
  plotStyle = PlotStyle /. {opts} /. Options[ShowIFS];
  colors = Colors  /. {opts} /. Options[ShowIFS];

  (* Generate the list of colors *)
  Which[
    colors === Automatic || colors === True,
      colors =  ColorData["Rainbow"] /@ 
        Range[0.,1 - 1./Length[IFS], 1./Length[IFS]],
    Head[colors] === String,
      colors = ColorData[colors] /@  
        Range[0.,1 - 1./Length[IFS], 1./Length[IFS]],
    colors =!= None && colors =!= False && Head[colors] =!= List,
      colors = colors /@
        Range[0.,1 - 1./Length[IFS], 1./Length[IFS]]];


  toFunc[{A_, b_}] := Module[{cfOut, fOut},
     cfOut = Compile[{{v, _Real, 1}}, A.v + b];
     fOut[{x_?NumericQ, y_?NumericQ}] := cfOut[{x, y}];
     fOut[x_List] := fOut /@ x;
     fOut[Point[pts_]] := Point[fOut[pts]];
     fOut[Line[x_]] := Line[fOut[x]];
     fOut[Arrow[x_]] := Arrow[fOut[x]];
     fOut[Polygon[x_, pOpts___]] := Polygon[fOut[x], pOpts];
     fOut[x_] := x;
     fOut];
  funcs = toFunc /@ IFS;

  extendGPrimitives[ff_List,Polygon[pp:{{_,_}..}]] :=
    Polygon[Table[f[pp],{f,ff}]];
  extendGPrimitives[ff_List,Polygon[pp:{{{_,_}..}..}]] :=
    Polygon[Flatten[Table[f[pp],{f,ff}],1]];
  extendGPrimitives[ff_List,Line[pp:{{_,_}..}]] :=
    Line[Table[f[pp],{f,ff}]];
  extendGPrimitives[ff_List,Line[pp:{{{_,_}..}..}]] :=
    Line[Flatten[Table[f[pp],{f,ff}],1]];
  extendGPrimitives[ff_List,Point[pp:{_,_}]] :=
   Point[Table[f[pp],{f,ff}]];
  extendGPrimitives[ff_List,Point[pp:{{_,_}..}]] :=
    Point[Flatten[Table[f[pp],{f,ff}],1]];
  extendGPrimitives[ff_List,ll_List] := 
    extendGPrimitives[ff,#]& /@ ll;
  extendGPrimitives[ff_List,x_] := x;

  twoNorms = Max[SingularValueList[#]]& /@
    N[First /@ IFS];
  extend[ff[f_,r_Real]] := If[r<tol, ff[f,r],
    Table[ff[Composition[f,funcs[[i]]],r*twoNorms[[i]]],{i,Length[IFS]}]];
  extend[ffs_List] := extend /@ ffs;

(*  extendedFuncs = Map[First, Flatten /@ 
    FixedPoint[extend,ff[#&,1.0]], {2}];
*)  extendedFuncs = FixedPoint[extend,ff[#&,1.0]];
  extendedFuncs = If[Head[#]===ff,{#},#]& /@ extendedFuncs;
  extendedFuncs = Map[First, Flatten /@ extendedFuncs, {2}];
  attractor = Table[extendGPrimitives[ff,initiator],
    {ff,extendedFuncs}];
  If[colors =!= False && colors =!= None, 
    coloredAttractor = Table[{colors[[i]],attractor[[i]]},
      {i,Length[IFS]}],
    coloredAttractor = attractor 
  ];

  Graphics[{Sequence@@Flatten[{plotStyle}], coloredAttractor},
    FilterRules[{opts}, Options[Graphics]],
    FilterRules[Options[ShowIFS], Options[Graphics]]]
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, IFS]] ||
    Message[IteratedFunctionSystems::badIFS, ShowIFS] ) &&
  ( And @@ Map[OptionQ, {opts}] ||
      Message[IteratedFunctionSystems::badOpt,ShowIFS] )



ShowIFSStochastic[IFS_, n_, opts___] := Module[
  {funcs, toFunc, coords, chooser, randomSequence,
   color, colors, coloredPoints, d, pointColors,
   pList, pSums, twoNorms, approximateDim, plotStyle},
        
  plotStyle = PlotStyle /. {opts} /. Options[ShowIFSStochastic];
  colors = Colors /. {opts} /. Options[ShowIFSStochastic];
  pList = Probabilities /. {opts} /. Options[ShowIFSStochastic];

  (* Generate the list of colors *)
  Which[
    colors === Automatic || colors === True,
      colors =  Hue[#,.7,.7]& /@ 
        Range[0.,1 - 1./Length[IFS], 1./Length[IFS]],
    Head[colors] === String,
      colors = ColorData[colors] /@ 
        Range[0.,1 - 1./Length[IFS], 1./Length[IFS]],
    colors =!= None && colors =!= False && Head[colors] =!= List,
      colors = colors /@
        Range[0.,1 - 1./Length[IFS], 1./Length[IFS]]];

  toFunc[{A_,b_}] := Compile[{{v,_Real,1}},
    A.v + b];
  funcs = toFunc /@ IFS;

  If[pList == Automatic,
    twoNorms = Max[SingularValueList[#]]& /@
      N[First /@ IFS];
    approximateDim = d /. 
      FindRoot[Plus @@ (#^d & /@ twoNorms) == 1, {d,1}];
    pList = #^approximateDim & /@ twoNorms
  ];

  randomSequence = RandomChoice[pList -> Range[Length[IFS]], n];
  coords = Developer`ToPackedArray[Drop[ComposeList[
    Table[funcs[[ randomSequence[[i]] ]], 
      {i,1,n}], {0,0}], Min[10, n]]];

  If[colors =!= False && colors =!= None, 
    coloredPoints = {colors[[#]],Point[coords[[Flatten[
      Position[Drop[randomSequence, Min[9, n]],#]]]]]}& /@
        Range[Length[IFS]];,
    coloredPoints = Point[coords]];

  Show[Graphics[{Sequence@@Flatten[{plotStyle}], coloredPoints},
	FilterRules[{opts}, Options[Graphics]],
    FilterRules[Options[ShowIFSStochastic], Options[Graphics]]]]
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, IFS]] ||
    Message[IteratedFunctionSystems::badIFS, ShowIFSStochastic] )  &&
  ( (IntegerQ[n] && n >= 0 ) ||
    Message[IteratedFunctionSystems::badInt, ShowIFSStochastic, n] ) &&
  ( And @@ Map[OptionQ, {opts}] ||
      Message[IteratedFunctionSystems::badOpt,ShowIFSStochastic] )


FindIFSDimension[IFS_] := Module[{twoNorms, d},
  twoNorms = Max[SingularValueList[#]]& /@
    N[First /@ IFS];
  d /. FindRoot[Plus @@ (#^d & /@ twoNorms) == 1, {d,1}]
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, IFS]] ||
    Message[IteratedFunctionSystems::badIFS, 
      FindIFSDimension] )


FindProbabilities[IFS_] := Module[{d, twoNorms, approximateDim},
  twoNorms = Max[SingularValueList[#]]& /@
    N[First /@ IFS];
  approximateDim = d /. 
    FindRoot[Plus @@ (#^d & /@ twoNorms) == 1, {d,1}];
  #^approximateDim & /@ twoNorms
] /;
  (And @@ Flatten[Map[MatchQ[#,
    {{{a_?NumericQ, b_?NumericQ}, {c_?NumericQ, d_?NumericQ}},
    {e_?NumericQ, f_?NumericQ}}] &, IFS]] ||
    Message[IteratedFunctionSystems::badIFS, 
      FindProbabilities] )
      
SierpinskiPedalTriangleIFS[{A_,B_,C_}] := Module[
  {A1, B1, C1, M, v, a, b, c, d, x, y, eqs,
     fA, fB, fC}, 

  A1 = C + ((B-C).(A-C))/((B-C).(B-C)) (B-C);
  B1 = C + ((A-C).(B-C))/((A-C).(A-C)) (A-C);
  C1 =A+ ((B-A).(C-A))/((B-A).(B-A)) (B-A);
  M = {{a,b},{c,d}}; 
    v = {x,y};
  eqs = {M.A + v == A, M.B + v == B1, M.C + v == C1};
  {fA} = {{{a,b},{c,d}},{x,y}} /. 
    Solve[eqs,{a,b,c,d,x,y}];
  eqs = {M.A + v == A1, M.B + v == B, M.C + v == C1};
  {fB} = {{{a,b},{c,d}},{x,y}} /. 
    Solve[eqs,{a,b,c,d,x,y}];
  eqs = {M.A + v == A1, M.B + v == B1, M.C + v == C};
  {fC} = {{{a,b},{c,d}},{x,y}} /. 
    Solve[eqs,{a,b,c,d,x,y}];
  {fA,fB,fC}
];
      
End[]  (* End Private Context *)

Protect[IteratedFunctionSystems, ShowIFS, ShowIFSStochastic, 
  FindProbabilities, FindIFSDimension,
  SierpinskiPedalTriangleIFS, Colors,
  Initiator, Probabilities]

EndPackage[]
