implemented Makie plot, extended N for #10

chakravala · Mar 9, 2021 · e3329c6 · e3329c6
1 parent e2fdfd0
commit e3329c6
Show file tree

Hide file tree

Showing 12 changed files with 242 additions and 150 deletions.
diff --git a/appveyor.yml → .appveyor.yml b/appveyor.yml → .appveyor.yml
@@ -1,7 +1,11 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1
+  - julia_version: 1.1
+  - julia_version: 1.2
+  - julia_version: 1.3
+  - julia_version: 1.4
+  - julia_version: 1.5
   - julia_version: nightly
 
 platform:
@@ -12,7 +16,6 @@ platform:
 ## (tests will run but not make your overall status red)
 matrix:
   allow_failures:
-  - julia_version: 0.7
   - julia_version: nightly
 
 branches:

diff --git a/.travis.yml b/.travis.yml
@@ -4,21 +4,17 @@ os:
   - linux
   - osx
 julia:
-  - 0.7
   - 1.0
+  - 1.1
+  - 1.2
+  - 1.3
+  - 1.4
+  - 1.5
   - nightly
 matrix:
   allow_failures:
-    - julia: 0.7
     - julia: nightly
 notifications:
   email: false
-# uncomment the following lines to override the default test script
-#script:
-#  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
-#  - julia -e 'Pkg.clone(pwd()); Pkg.build("Reduce"); Pkg.test("Reduce"; coverage=true)'
 after_success:
-  # push coverage results to Coveralls
-  - julia -e 'cd(Pkg.dir("Reduce")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
-  # push coverage results to Codecov
-  - julia -e 'cd(Pkg.dir("Reduce")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
+  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder()); Coveralls.submit(process_folder())'
diff --git a/LICENSE.md b/LICENSE.md
diff --git a/Project.toml b/Project.toml
@@ -0,0 +1,25 @@
+name = "Fatou"
+uuid = "5f923234-c850-556d-bb4f-28324fa1959a"
+authors = ["Michael Reed"]
+version = "1.1.1"
+
+[deps]
+ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+SyntaxTree = "a4af3ec5-f8ac-5fed-a759-c2e80b4d74cb"
+Reduce = "93e0c654-6965-5f22-aba9-9c1ae6b3c259"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+
+[compat]
+julia = "1"
+Reduce = "1.2"
+Requires = "1"
+SyntaxTree = "1"
+ColorSchemes = "3"
+LaTeXStrings = "1"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
diff --git a/README.md b/README.md
@@ -2,7 +2,7 @@
 
 [![Build Status](https://travis-ci.org/chakravala/Fatou.jl.svg?branch=master)](https://travis-ci.org/chakravala/Fatou.jl) [![Build status](https://ci.appveyor.com/api/projects/status/mdathjmu7jg57u77?svg=true)](https://ci.appveyor.com/project/chakravala/fatou-jl) [![Coverage Status](https://coveralls.io/repos/github/chakravala/Fatou.jl/badge.svg?branch=master)](https://coveralls.io/github/chakravala/Fatou.jl?branch=master) [![codecov.io](http://codecov.io/github/chakravala/Fatou.jl/coverage.svg?branch=master)](http://codecov.io/github/chakravala/Fatou.jl?branch=master)
 
-Julia package for Fatou sets. Install using `Pkg.add("Fatou")` in Julia. See [Explore Fatou sets & Fractals](https://github.com/chakravala/Fatou.jl/wiki/Explore-Fatou-sets-&-fractals) in Wiki for detailed *examples*. This package provides: `fatou`, `juliafill`, `mandelbrot`, `newton`, `basin`, `plot`, and `orbit`; along with various internal functionality using `Reduce` and Julia expressions to help compute `Fatou.FilledSet` efficiently. Full documentation is included. The `fatou` function can be applied to a `Fatou.Define` object to produce a `Fatou.FilledSet`, which can then be passed as an argument to the `plot` function of `PyPlot`. Creation of `Fatou.Define` objects is done via passing a `parse`-able function expression string (in variables `z`, `c`) and optional keyword arguments to `juliafill`, `mandelbrot`, and `newton`.
+Julia package for Fatou sets. Install using `Pkg.add("Fatou")` in Julia. See [Explore Fatou sets & Fractals](https://crucialflow.com/Fatou.jl) in Wiki for detailed *examples*. This package provides: `fatou`, `juliafill`, `mandelbrot`, `newton`, `basin`, `plot`, and `orbit`; along with various internal functionality using `Reduce` and Julia expressions to help compute `Fatou.FilledSet` efficiently. Full documentation is included. The `fatou` function can be applied to a `Fatou.Define` object to produce a `Fatou.FilledSet`, which can then be passed as an argument to `plot` functions of `Makie`, `PyPlot`, `ImageInTerminal`. Creation of `Fatou.Define` objects is done via passing a `parse`-able function expression string (in variables `z`, `c`) and optional keyword arguments to `juliafill`, `mandelbrot`, and `newton`.
 
 ## Background
 
@@ -22,13 +22,40 @@ The number of Julia threads available is detected at the startup and is reported
 Since each pixel is independent of any other pixel, it doesn’t matter in what order or on how many threads it is computed, the more you use the faster it is.
 The environment variable `JULIA_NUM_THREADS` can be used to enable the multi-threading for more than 1 thread.
 
-Please share your favorite fractals as `Fatou` snippet in the [discussion thread](https://discourse.julialang.org/t/ann-fatou-jl-easily-share-julia-fractals)!
+Please share with us your favorite fractals as `Fatou` code snippets!
 
-## PyPlot.jl compatability features
+## Examples
 
-The program can be initialized with `using Fatou, PyPlot` or `ImageInTerminal`.
+`Fatou.Define` provides the following optional keyword arguments:
 
-A Fatou set is a collection of complex valued orbits of an iterated function. To help illustrate this, an additional feature is a plot function designed to visualize real-valued-orbits. The following is a cobweb orbit plot of a function:
+```Julia
+Q::Expr 	= :(abs2(z)),           # escape criterion, (z, c) -> Q
+C::Expr 	= :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
+∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
+n::Integer  = 176,                  # vertical grid points
+N::Integer  = 35,                   # max. iterations
+ϵ::Number   = 4,                    # basin ϵ-Limit criterion
+iter::Bool  = false,                # toggle iteration mode
+p::Number   = 0,                    # iteration color exponent
+newt::Bool  = false,                # toggle Newton mode
+m::Number   = 0,                    # Newton multiplicity factor
+mandel::Bool= false,                # toggle Mandelbrot mode
+seed::Number= 0.0+0.0im,            # Mandelbrot seed value
+x0          = nothing,              # orbit starting point
+orbit::Int  = 0,                    # orbit cobweb depth
+depth::Int  = 1,                    # depth of function composition
+cmap::String= ""                    # imshow color map
+```
+
+A Fatou set is a collection of complex valued orbits of an iterated function. To help illustrate this, an additional feature is a plot function designed to visualize real-valued-orbits.
+The program can be initialized with `using Fatou, PyPlot` or `Makie` or `ImageInTerminal`.
+For `PyPlot` the `imshow` and `plot` methods can be used, while for `Makie` the `heatmap`, `contour`, `surface`, and `arrows` methods can be used.
+
+When `using ImageInTerminal`, the display of a `Fatou.FilledSet` will be plotted automatically in the terminal.
+The `orbit` method also has optional `UnicodePlots` compatibility.
+Additional plotting support can be added via Pull-Request by adding another `Requires` script to the `__init__()` function definition.
+
+The following is a cobweb orbit plot of a function:
 
 ```Julia
 juliafill(:(z^2-0.67),∂=[-1.25,1.5],x0=1.25,orbit=17,depth=3,n=147) |> orbit
@@ -88,15 +115,7 @@ basin(nf,2)
 
 ![D2(ϵ)](http://latex.codecogs.com/svg.latex?D_2(\epsilon)%20=%20\left\\{z\in\mathbb{C}:\left|\\,z%20-%20\frac{1}{\cos{\left%20(z%20\right%20)}}%20\left(1%20+%20i\right)%20\left(\sin{\left%20(z%20\right%20)}%20-%201\right)%20-%20\frac{\left(1%20+%20i\right)%20\left(\sin{\left%20(z%20-%20\frac{1}{\cos{\left%20(z%20\right%20)}}%20\left(1%20+%20i\right)%20\left(\sin{\left%20(z%20\right%20)}%20-%201\right)%20\right%20)}%20-%201\right)}{\cos{\left%20(z%20-%20\frac{1}{\cos{\left%20(z%20\right%20)}}%20\left(1%20+%20i\right)%20\left(\sin{\left%20(z%20\right%20)}%20-%201\right)%20\right%20)}}%20-%20r_i\\,\right|%3C\epsilon,\\,\forall%20r_i(\\,f(r_i)=0%20)\right\\})
 
-## ImageInTerminal.jl compatibility
-
-When `using ImageInTerminal`, the display of a `Fatou.FilledSet` will be plotted automatically in the terminal.
-The `orbit` method also has optional `UnicodePlots` compatibility.
-Additional plotting support can be added via Pull-Request by adding another `Requires` script to the `__init__()` function definition.
-
-## Detailed Explanation
-
-View [Explore Fatou sets & Fractals](https://github.com/chakravala/Fatou.jl/wiki/Explore-Fatou-sets-&-fractals) in Wiki for detailed *examples*.
+View [Explore Fatou sets & Fractals](https://crucialflow.com/Fatou.jl) in Wiki for detailed *examples*.
 
 ### Troubleshooting on Julia 1.0.1+
 

diff --git a/REQUIRE b/REQUIRE
diff --git a/src/Fatou.jl b/src/Fatou.jl
@@ -37,7 +37,7 @@ struct Define{FT<:Function,QT<:Function,CT<:Function,M,N,P,D}
     C::CT # complex fixed point coloring
     ∂::Array{Float64,1} # bounds
     n::UInt16 # number of grid points
-    N::UInt8 # number of iterations
+    N::UInt16 # number of iterations
     ϵ::Float64 # epsilon Limit criterion
     iter::Bool # toggle iteration mode
     p::Float64 # iteration color exponent
@@ -76,7 +76,7 @@ struct Define{FT<:Function,QT<:Function,CT<:Function,M,N,P,D}
         (f = genfun(newton_raphson(E,m),[:z,:c]); q = genfun(Expr(:call,:abs,E),[:z,:c]))
         c = genfun(C,[:z,:n,:p])
         e = typeof(E) == String ? parse(E) : E
-        return new{typeof(f),typeof(q),typeof(c),mandel,newt,plane,disk}(e,f,q,c,convert(Array{Float64,1},∂),UInt16(n),UInt8(N),float(ϵ),iter,float(p),newt,m,mandel,seed,x0,orbit,depth,cmap,plane,disk)
+        return new{typeof(f),typeof(q),typeof(c),mandel,newt,plane,disk}(e,f,q,c,convert(Array{Float64,1},∂),UInt16(n),UInt16(N),float(ϵ),iter,float(p),newt,m,mandel,seed,x0,orbit,depth,cmap,plane,disk)
     end
 end
 
@@ -88,7 +88,7 @@ Compute the `Fatou.FilledSet` set using `Fatou.Define`.
 struct FilledSet{FT,QT,CT,M,N,P,D}
     meta::Define{FT,QT,CT,M,N,P,D}
     set::Matrix{Complex{Float64}}
-    iter::Matrix{UInt8}
+    iter::Matrix{UInt16}
     mix::Matrix{Float64}
     function FilledSet{FT,QT,CT,M,N,P,D}(K::Define{FT,QT,CT,M,N,P,D}) where {FT,QT,CT,M,N,P,D}
         (i,s) = Compute(K)
@@ -269,37 +269,80 @@ disk(z::Complex) = (2z.re/(z.re^2+(1+z.im)^2))+im*(z.re^2+z.im^2-1)/(z.re^2+(1+z
 # define function for computing orbit of a z0 input
 function orbit(K::Define{FT,QT,CT,M,N,P,D},z0::Complex{Float64}) where {FT,QT,CT,M,N,P,D}
     M ? (z = K.seed) : (z = P ? plane(z0) : z0)
-    zn = 0x00
+    zn = 0x0000
     while (N ? (K.Q(z,z0)::Float64>K.ϵ)::Bool : (K.Q(z,z0)::Float64<K.ϵ))::Bool && K.N>zn
         z = K.F(z,z0)::Complex{Float64}
-        zn+=0x01
+        zn+=0x0001
     end; #end
     # return the normalized argument of z or iteration count
-    return (zn::UInt8,(D ? disk(z) : z)::Complex{Float64})
+    return (zn::UInt16,(D ? disk(z) : z)::Complex{Float64})
 end
 
+gridy(K::Define) = round(UInt16,(K.∂[4]-K.∂[3])/(K.∂[2]-K.∂[1])*K.n)
+
+ranges(K::FilledSet) = ranges(K.meta)
+function ranges(K::Define)
+    x = range(K.∂[1]+0.0001,stop=K.∂[2],length=K.n)
+    y = range(K.∂[4],stop=K.∂[3],length=gridy(K))
+    return x,y
+end
+
+function grid(K::Define) # generate coordinate grid
+    x, y = ranges(K); x' .+ im*y
+end # apply Newton-Orbit function element-wise to coordinate grid
+
 """
     Compute(::Fatou.Define)::Union{Matrix{UInt8},Matrix{Float64}}
 
 `Compute` the `Array` for `Fatou.FilledSet` as specefied by `Fatou.Define`.
 """
-function Compute(K::Define{FT,QT,CT,M,N,D})::Tuple{Matrix{UInt8},Matrix{Complex{Float64}}} where {FT,QT,CT,M,N,D}
-    # generate coordinate grid
-    Kyn = round(UInt16,(K.∂[4]-K.∂[3])/(K.∂[2]-K.∂[1])*K.n)
-    x = range(K.∂[1]+0.0001,stop=K.∂[2],length=K.n)
-    y = range(K.∂[4],stop=K.∂[3],length=Kyn)
-    Z = x' .+ im*y # apply Newton-Orbit function element-wise to coordinate grid
-    (matU,matF) = (Array{UInt8,2}(undef,Kyn,K.n),Array{Complex{Float64},2}(undef,Kyn,K.n))
-    @time @threads for j = 1:length(y); for k = 1:length(x);
-        (matU[j,k],matF[j,k]) = orbit(K,Z[j,k])::Tuple{UInt8,Complex{Float64}}
+function Compute(K::Define{FT,QT,CT,M,N,D})::Tuple{Matrix{UInt16},Matrix{Complex{Float64}}} where {FT,QT,CT,M,N,D}
+    Z,Kyn = grid(K),gridy(K)
+    (matU,matF) = (Array{UInt16,2}(undef,Kyn,K.n),Array{Complex{Float64},2}(undef,Kyn,K.n))
+    @time @threads for j = 1:size(Z)[1]; for k = 1:size(Z)[2];
+        (matU[j,k],matF[j,k]) = orbit(K,Z[j,k])::Tuple{UInt16,Complex{Float64}}
     end; end
     return (matU,matF)
 end
 
+# determine if plot is Iteration, Roots, or Limit
+typeplot(K::FilledSet) = typeof(K.meta.iter ? K.iter : K.mix) == Matrix{UInt16} ? "iter." : K.meta.m==1 ? "roots" : "limit"
+
+function Base.String(K::FilledSet) # annotate title using LaTeX
+    text,t = "f : z ↦ $(K.meta.E),",typeplot(K)
+    K.meta.newt ? "$text m = $(K.meta.m), $t" : "$text $t"
+end
+
 function __init__()
     println("Fatou detected $(Threads.nthreads()) julia threads.")
+    @require ColorSchemes="35d6a980-a343-548e-a6ea-1d62b119f2f4" begin
+        nonan(x) = isnan(x) ? 0.0 : x
+        function (C::ColorSchemes.ColorScheme)(K::FilledSet)
+            S = size(K.iter)
+            H = zeros(ColorSchemes.RGB{Float64},S...)
+            if K.meta.iter
+                M = length(C)/(maximum(K.iter)+1)
+                for x ∈ 1:S[1], y ∈ 1:S[2]
+                    H[x,y] = C[round(Int,M*(K.iter[x,y]+1),RoundUp)]
+                end
+            else
+                for x ∈ 1:S[1], y ∈ 1:S[2]
+                    H[x,y] = get(C,nonan(K.mix[x,y]))
+                end
+            end
+            return H
+        end
+    end
+    @require ImageInTerminal="d8c32880-2388-543b-8c61-d9f865259254" begin
+        import ColorSchemes
+        function Base.show(io::IO,K::FilledSet;c::String="",bare::Bool=false)
+            isempty(c) && (c = K.meta.cmap)
+            display(getproperty(ColorSchemes, isempty(c) ? :balance : Symbol(c))(K))
+            !bare && print(io,String(K))
+        end
+    end
+    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("makie.jl")
     @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("pyplot.jl")
-    @require ImageInTerminal="d8c32880-2388-543b-8c61-d9f865259254" include("term.jl")
     @require UnicodePlots="b8865327-cd53-5732-bb35-84acbb429228" include("uniplots.jl")
 end
 

diff --git a/src/makie.jl b/src/makie.jl
@@ -0,0 +1,79 @@
+#   This file is part of Fatou.jl. It is licensed under the MIT license
+#   Copyright (C) 2019 Michael Reed
+
+import ColorSchemes
+
+Makie.plot(K::FilledSet;bare=false,args...) = Makie.heatmap(K;bare=bare,args...)
+
+for plt ∈ (:contour,:contourf,:heatmap)
+    plt! = Symbol(plt,:!)
+    @eval function Makie.$plt(K::FilledSet;bare::Bool=false,args...)
+        scene, layout = Makie.layoutscene()
+        ax = layout[1,1] = if K.meta.newt 
+            Makie.Axis(scene,title=String(K),ylabel="Fatou set: z ↦ z-m×f(z)/f'(z)")
+        else
+            Makie.Axis(scene,title=String(K))
+        end
+        plt = Makie.$plt!(ax,K;args...)
+        layout[1,2] = Makie.Colorbar(scene,plt;width=30)
+        return scene
+    end
+    @eval function Makie.$plt!(ax,K::FilledSet;bare::Bool=false,args...)
+        !haskey(args,:colormap) && !isempty(K.meta.cmap) && (return Makie.$plt!(ax,K;bare=bare,colormap=Symbol(K.meta.cmap),args...))
+        r1,r2 = ranges(K)
+        Z = $(plt≠:heatmap ? :(transpose(K.meta.iter ? K.iter : K.mix)) : :(reverse(transpose(K.meta.iter ? K.iter : K.mix),dims=2)))
+        if bare
+            Makie.$plt!(ax,r1,r2,Z;args...)
+        else
+            Makie.$plt!(ax,r1,r2,Z;axis=(title=String(K),),args...)
+        end
+    end
+end # center!(scene)
+
+function Makie.surface(K::FilledSet;bare::Bool=false,args...)
+    !haskey(args,:colormap) && !isempty(K.meta.cmap) && (return Makie.surface(K;bare=bare,colormap=Symbol(K.meta.cmap),args...))
+    r1,r2 = ranges(K)
+    Z = reverse(reverse(transpose(K.meta.iter ? K.iter : K.mix),dims=1),dims=2)
+    if bare
+        Makie.surface(reverse(r1),r2,Z;(xreversed=true,)args...)
+    else
+        Makie.surface(reverse(r1),r2,Z;axis=(title=String(K),xreversed=true),args...)
+    end
+end
+
+function Makie.arrows(K::FilledSet;bare::Bool=false,args...)
+    r1,r2 = ranges(K)
+    Z = transpose(K.set)
+    Makie.arrows(r1,r2,real.(Z),imag.(Z);args...)
+end
+
+function orbit(E,f::Function,bi::Matrix{Float64},orb::Int=0,depth::Int=1,incr::Int=384)
+    x,N,N2,orbit,bis = real_orb(E,f,bi,orb,depth,incr)
+    # setup plot
+    funs,funt = ~(orb == 0) ? (["ϕ(x₀₋ₙ)"],", IC: x₀ = $(bis[3]), n∈0:$orb") : ([],"")
+    scene, layout = Makie.layoutscene()
+    ax = layout[1,1] = Makie.Axis(scene,title="x ↦ $E$funt")
+    # plot background lines
+    plt = [Makie.lines!(ax,x[:],N[:,1];linestyle=:dash),
+    Makie.lines!(ax,x[:],N[:,2]),
+    # plot orbit cobweb path
+    Makie.lines!(ax,orbit[:,1],orbit[:,2];color=:red)]
+    # plot f^2,f^3,f^4,...
+    for h ∈ 3:depth+1
+        push!(plt,Makie.lines!(ax,x,N[:,h];linewidth=1,color=ColorSchemes.tab10[h>5 ? h-1 : h-2]))
+    end
+    if orb≠0
+        ran = range(bi[1],stop=bi[2],length=length(N2))
+        plt = [plt...,[
+            Makie.lines!(ax,ran,N2[:];color=:gray,linestyle=:dot,linewidth=1),
+            Makie.scatter!(ax,ran,N2[:];color=:gray,marker=:x)]]
+    end
+    # trim graph
+    d=1.07
+    Makie.xlims!(ax,bi[1],bi[2])
+    Makie.ylims!(ax,minimum([d*minimum(N[:,2]),0]),maximum([d*maximum(N[:,2]),0]))
+    # set legend
+    layout[1,2] = Legend(scene,plt,vcat(["y=x","ϕ(x)","(xₙ,ϕ(xₙ))"],["ϕ^$x(x)" for x ∈ 2:depth],funs))
+    return scene
+end
+
diff --git a/src/orbitplot.jl b/src/orbitplot.jl
@@ -15,7 +15,7 @@ julia> juliafill("z^2-0.67",∂=[-1.25,1.5],x0=1.25,orbit=17,depth=3) |> orbit
 ```
 """
 function orbit(K::Define{FT,QT,CT}) where {FT,QT,CT}
-    !isdefined(Fatou,:UnicodePlots) && !isdefined(Fatou,:PyPlot) && throw(error("Requires `using PyPlot` or `using UnicodePlots`"))
+    !isdefined(Fatou,:UnicodePlots) && !isdefined(Fatou,:PyPlot) && !isdefined(Fatou,:Makie) && throw(error("Requires `using PyPlot` or `using UnicodePlots` or `using Makie`"))
     K.x0 == nothing ? (bi = K.∂[1:2]') : (bi = [K.∂[1:2]...,K.x0]')
     orbit(K.E,z->K.F(z,0),convert(Array{Float64},bi),K.orbit,K.depth,Int(K.n))
 end