fixed treeplot() after thanksgiving refactor

Author: jeffwack

Project.toml

@@ -1,7 +1,7 @@
 name = "Mandelbrot"
 uuid = "52854307-bbd1-4ac7-933f-adc0936abe8a"
 authors = ["Jeffrey Wack"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"

README.md

@@ -5,6 +5,6 @@ This package implements several algorithms related to complex quadratic dynamics
 - The external angles of a hyperbolic component can be calculated from an angled internal address, describing the path to this component from zero. [Read about internal addresses.](https://arxiv.org/abs/math/9411238)
 - A combinatorial description of a Hubbard trees can be generated from a kneading sequence, and when oriented in the plane can produce external angles. [Read about Hubbard trees.](https://www.mat.univie.ac.at/~bruin/papers/bkafsch.pdf)
 
-This package is a work in progress. Currently, KneadingSequence, AngledInternalAddress, HubbardTree, and OrientedHubbardTree can be constructed. To run the spider algorithm use parameter().
+This package is a work in progress. To see what it can do, run treeplot(theta) where theta is a rational number with an odd denominator. The external ray of the Mandelbrot set with angle theta lands at a hyperbolic component, and this will plot the Hubbard tree corresponding to this hyperbolic component. Note that in Julia, rational numbers are declared with two slashes, as "3//5".
 
 [![Build Status](https://github.com/jeffwack111/Spiders.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/jeffwack111/Spiders.jl/actions/workflows/CI.yml?query=branch%3Amain)

src/HubbardTrees.jl

@@ -3,7 +3,7 @@
 #Henk Bruin, Alexandra Kafll, Dierk Schleicher
 abstract type AbstractHubbardTree <: Graph end
 
-function criticalpoint(H::AbstractHubbardTree)
+function criticalpoint(H)
     return first(filter(x->x[1]==KneadingSymbol('*'),vertices(H)))
 end
 
@@ -154,8 +154,9 @@ function findone(f,A)
     end
 end
 
-function isbetween(htree::AbstractHubbardTree,a,b,c)
-    zero = htree.criticalpoint
+#is this for some reason better than a graph search?
+function isbetween(htree,a,b,c) 
+    zero = criticalpoint(htree)
     K = shift(zero)
     (type,vertex) = iteratetriod(K,(a,b,c))
     if type == "flat" && vertex == a

src/Mandelbrot.jl

@@ -42,5 +42,4 @@ include("showrays.jl")
 include("showspider.jl")
 include("showtree.jl")
 
-
 end

src/angledoubling.jl

@@ -340,7 +340,7 @@ function numembeddings(theta::Rational)
     return numembeddings(internaladdress(theta))
 end
 
-function period(theta::Rational)
+function period(theta::RationalAngle)
     return orbit(theta).period
 end
 

src/dynamicrays.jl

@@ -46,7 +46,7 @@ function dynamicrays(c::Complex,angle::BinaryExpansion,R::Real,res::Int,depth::I
     radii = collect(LinRange(R,sqrt(R),res))
 
     #similar to standard spider but is outer
-    rays = Sequence{Vector{ComplexF64}}([exp(2im*pi*Rational(theta)).*radii for theta in orb.items],orb.preperiod)
+    rays = Sequence{Vector{ComplexF64}}([exp(2im*pi*RationalAngle(theta).value).*radii for theta in orb.items],orb.preperiod)
 
     for jj in 1:depth
         for (target,source) in goesto(rays)

src/embedtrees.jl

@@ -1,9 +1,7 @@
 struct HyperbolicComponent
-    criticalpoint::Sequence
-    adj::Dict{KneadingSequence,Vector{KneadingSequence}}
-    boundary::Vector{KneadingSequence}
+    htree::OrientedHubbardTree
     onezero::Dict{KneadingSequence, Union{Nothing,Digit{2}}}
-    angle::Rational
+    angle::RationalAngle
     rays::Dict{BinaryExpansion,Vector{ComplexF64}}
     parameter::ComplexF64
     vertices::Dict{KneadingSequence, ComplexF64} #maps each vertex in the tree to a point in the plane
@@ -17,7 +15,7 @@ function HyperbolicComponent(OHT::OrientedHubbardTree)
 
     criticalorbit = orbit(OHT.criticalpoint) #This is a sequence orbit
 
-    theta = Rational(first(criticalanglesof(OZ,OHT)))
+    theta = RationalAngle(first(criticalanglesof(OZ,OHT)))
     c = parameter(standardspider(theta),250)
 
     #critical orbit
@@ -37,8 +35,8 @@ function HyperbolicComponent(OHT::OrientedHubbardTree)
         end
     end
 
-    merge!(rays,dynamicrays(c,BinaryExpansion(1//2),100,10,40))
-    merge!(rays,dynamicrays(c,BinaryExpansion(0//1),100,10,40))
+    merge!(rays,dynamicrays(c,BinaryExpansion(RationalAngle(1//2)),100,10,40))
+    merge!(rays,dynamicrays(c,BinaryExpansion(RationalAngle(0//1)),100,10,40))
 
     paramorbit = [0.0+0.0im]
     n = period(theta)
@@ -50,16 +48,24 @@ function HyperbolicComponent(OHT::OrientedHubbardTree)
 
     zvalues = Dict{Sequence,ComplexF64}()
     for node in keys(OHT.adj)
-        if star() in node.items #then it is in the critical orbit
-            idx = findone(x->x==star(),node.items)
+        if KneadingSymbol('*') in node.items #then it is in the critical orbit
+            idx = findone(x->x==KneadingSymbol('*'),node.items)
             push!(zvalues,Pair(node,paramorbit[mod1(n-idx+2,n)]))
         else
             list = [rays[angle][end] for angle in anglelist[node]]
             push!(zvalues,Pair(node,sum(list)/length(list)))
         end
     end
     E = standardedges(OHT.adj,OHT.criticalpoint,zvalues)
-    return HyperbolicComponent(OHT.criticalpoint,OHT.adj,OHT.boundary,OZ,theta,rays,c,zvalues,E)
+    return HyperbolicComponent(OHT,OZ,theta,rays,c,zvalues,E)
+end
+
+function HyperbolicComponent(theta::Rational)
+    H0 = HyperbolicComponent(OrientedHubbardTree(theta))
+    for ii in 1:5
+        refinetree!(H0)
+    end
+    return H0
 end
 
 function Base.show(io::IO,H::HyperbolicComponent)
@@ -79,8 +85,8 @@ function standardedges(adj,zero,zvalues)
 end
 
 function longpath(EHT, start, finish)
-    nodes = nodepath(EHT.adj, start, finish)
-    edges = edgepath(EHT.adj, start, finish)
+    nodes = nodepath(EHT.htree, start, finish)
+    edges = edgepath(EHT.htree, start, finish)
 
     edgevectors = EHT.edges
 

src/showtree.jl

@@ -4,24 +4,24 @@ end
 
 function Makie.plot!(myplot::TreePlot)
     EHT = myplot.EHT[]
-    (EdgeList,Nodes) = adjlist(EHT.adj)
-    criticalorbit = orbit(EHT.criticalpoint)
+    (EdgeList,Nodes) = adjlist(EHT.htree.adj)
+    criticalorbit = orbit(EHT.htree.criticalpoint)
 
     labels = []
     nodecolors = []
     for node in Nodes
         firstchar = node.items[1]
-        if firstchar == star()
+        if firstchar == KneadingSymbol('*')
             push!(nodecolors,"black")
-        elseif firstchar == A() #we are fully in one of the 4 regions
+        elseif firstchar == KneadingSymbol('A') #we are fully in one of the 4 regions
             if EHT.onezero[node] == Digit{2}(0)
                 push!(nodecolors,"blue")
             elseif EHT.onezero[node] === nothing
                 push!(nodecolors,"turquoise")
             elseif EHT.onezero[node] == Digit{2}(1)
                 push!(nodecolors,"green")
             end
-        elseif firstchar == B()
+        elseif firstchar == KneadingSymbol('B')
             if EHT.onezero[node] == Digit{2}(0)
                 push!(nodecolors,"red")
             elseif EHT.onezero[node] === nothing
@@ -81,6 +81,9 @@ function Makie.plot!(myplot::TreePlot)
     return myplot
 end
 
+function treeplot(theta::Rational)
+    return treeplot(HyperbolicComponent(theta))
+end
 
 function plottree!(scene, H::HubbardTree)
 

src/spidermap.jl

@@ -1,6 +1,6 @@
 struct Spider
-    angle::Rational
-    orbit::Sequence{Rational}
+    angle::RationalAngle
+    orbit::Sequence{RationalAngle}
     kneading_sequence::KneadingSequence
     legs::Vector{Vector{ComplexF64}}
 end
@@ -10,7 +10,7 @@ function Base.show(io::IO, S::Spider)
     return println("A "*repr(S.angle)*"-spider with kneading sequence "*repr(S.kneading_sequence)*" and "*repr(npoints)*" points.")
 end
 
-function standardspider(theta::Rational)
+function standardspider(theta::RationalAngle)
     orb = orbit(theta)
     K = thetaitinerary(theta,orb) #the kneading sequence is calculated like so to avoid recalculating the orbit of theta
     legs = standardlegs(orb)
@@ -21,17 +21,12 @@ function standardlegs(orbit::Sequence)
     r = collect(LinRange(100,1,10))  #NOTE - it may be better to keep this as a 'linrange' but I don't understand what that means
     legs  = Vector{ComplexF64}[]
     for theta in orbit.items
-        push!(legs,(cos(theta*2*pi)+1.0im*sin(theta*2*pi)) .* r)
+        push!(legs,(cos(theta.value*2*pi)+1.0im*sin(theta.value*2*pi)) .* r)
     end
     legs[1] = legs[1] .- legs[1][end]
     return legs
 end
 
-function standardlegs(angle::Rational)
-    info = SpiderInfo(angle)
-    return standardlegs(info)
-end
-
 function grow!(legs::Vector{Vector{ComplexF64}},scale::Real,num::Int)
     #Get the arrow pointing from the second to last point to the last point
     for leg in legs
@@ -72,11 +67,11 @@ function mapspider!(S::Spider)
     #this angle lays in region A by definition. 
 
     if a/2 < theta1 && theta1 < (a+1)/2
-        cregion = A()
-        dregion = B()
+        cregion = KneadingSymbol('A')
+        dregion = KneadingSymbol('B')
     else
-        cregion = B()
-        dregion = A()
+        cregion = KneadingSymbol('B')
+        dregion = KneadingSymbol('A')
     end
 
     for ii in 2:n-1
@@ -106,7 +101,7 @@ function mapspider!(S::Spider)
         thetau = (angle(u)/(2*pi)+1)%1
 
         if abs2(thetau - a/2) < abs2(thetau - (a+1)/2)
-            if cregion == A()
+            if cregion == KneadingSymbol('A')
                 if S.kneading_sequence.items[end] == '2'
                     newLegs[end] = path_sqrt(S.legs[1]./λ)
                 else
@@ -120,7 +115,7 @@ function mapspider!(S::Spider)
                 end
             end
         else
-            if cregion == A()
+            if cregion == KneadingSymbol('A')
                 if S.kneading_sequence.items[end] == '1'
                     newLegs[end] = path_sqrt(S.legs[1]./λ)
                 else
@@ -183,7 +178,7 @@ function parameter(S0::Spider,max_iter::Int)
     for ii in 1:max_iter
         mapspider!(S)
         c = S.legs[2][end]/2
-        println(repr(c)*" delta "*repr(abs(c-c_last)))
+        #println(repr(c)*" delta "*repr(abs(c-c_last)))
         if abs(c-c_last)<(1e-15)
             return c
         end