Juliaで「グモウスキー・ミラの写像」を描いた

[グモウスキー・ミラの写像](https://ja.wikipedia.org/wiki/%E3%82%B0%E3%83%A2%E3%82%A6%E3%82%B9%E3%82%AD%E3%83%BC%E3%83%BB%E3%83%9F%E3%83%A9%E3%81%AE%E5%86%99%E5%83%8F)をJulia(version 1.1.x)で描きます。

# 初版
イテレータについては[お気楽 Julia プログラミング超入門 Julia の基礎知識]
(http://www.nct9.ne.jp/m_hiroi/light/juliaa02.html)を参考にしました。

```gm.jl
struct DynamicalSystem
    initial::Array{Float64,1}
    map::Function
end
Base.iterate(ds::DynamicalSystem, pt=ds.initial ) = (pt, ds.map(pt) )

function push_points(initial::Array{Float64,1} ,
                    map::Function ,
                    drop::Int64  ,
                    num::Int64 ,
                    plot_x::Array{Float16,1} , plot_y::Array{Float16,1})
    for pt in Iterators.take(
                Iterators.drop(DynamicalSystem(initial,map ),drop ) , num)
        push!(plot_x,pt[1])
        push!(plot_y,pt[2])
    end
end

@enum GMType f1g1 f1g2 f2g1 f2g2

struct GumowskiMira
    gmtype::GMType
    α::Float64
    σ::Float64
    μ::Float64

    # g1,g2 : 有界非線形性 (bounded nonlinerity)
    g1::Function
    g2::Function

    f1::Function
    f2::Function

    function GumowskiMira(gmtype::GMType ;α::Float64=0. ,σ::Float64=0. ,μ::Float64)
        g1=x->μ*x + (2.0(1.0-μ)x^2 )/(1.0 + x^2)
        g2=x->μ*x + (1.0-μ)x^2*exp((1.0 - x^2)/4.0)
        f1=g->function(pt::Array{Float64,1})
            xx =  pt[2] + g(pt[1])
            yy = -pt[1] + g(xx)
            [xx,yy]
        end
        f2=g->function(pt::Array{Float64,1})
            xx =  pt[2] + α*pt[2]*( 1.0 - σ*pt[2]^2 ) + g(pt[1])
            yy = -pt[1] + g(xx)
            [xx,yy]
        end
        new(gmtype,α,σ,μ,g1,g2,f1,f2)
    end

end

gmmap(gm::GumowskiMira)=
                    if      (gm.gmtype==f1g1)
                        gm.f1(gm.g1)
                    elseif (gm.gmtype==f1g2)
                        gm.f1(gm.g2)
                    elseif (gm.gmtype==f2g1)
                        gm.f2(gm.g1)
                    elseif (gm.gmtype==f2g2)
                        gm.f2(gm.g2)
                    else
                        error("gmtypeが不正")
                    end


#gm = GumowskiMira(f2g1,α=0.008 , σ=0.05, μ=-0.496) # 「神話の鳥 (mythic bird)」3枚羽の翼
gm = GumowskiMira(f2g1,α=0.009 , σ=0.05, μ=-0.801) # 「神話の鳥 (mythic bird)」5枚羽の翼
#gm = GumowskiMira(f2g2,α=0.0083, σ=0.1 , μ=-0.38 ) #
#gm = GumowskiMira(f2g2,α=0.01  , σ=0.1 , μ=0.8 ) #
#gm = GumowskiMira(f1g1,μ=0.39 ) #
#gm = GumowskiMira(f1g1,μ=0.365 ) #
#gm = GumowskiMira(f1g1,μ=0.34 ) #

initial = [1.1,1.1]
plot_x=Float16[]
plot_y=Float16[]
push_points(initial,gmmap(gm),20000,200000,plot_x,plot_y)

using Plots
pyplot()
scatter(plot_x,plot_y,markersize=1,markercolor=:black)
```
![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/402919/10785773-d10d-77b7-51e9-83671e59cf1d.png)

# 第二版(色をつける、など)
- moduleに分けました
- 型パラメータを使用しました
- パラメータの指定に名前付きタプルを使用しました
- 色つきとモノクロの二種類に対応しました(with_rgbパラメータ)

```DynamicalSystem.jl
module DynamicalSystem

using ColorTypes

export make_seq , make_initial_seq , push_points!
struct Seq{P}
    initial::P
    map::Function
end
Base.iterate(ds::Seq, pt=ds.initial ) = (pt, ds.map(pt) )

function make_seq(initial::P ,
                    map::Function ;
                    ndrop::Int64 = 0,
                    ntake::Int64 = 0) where{P} 
    Iterators.take(
        Iterators.drop(Seq{P}(initial,map ),ndrop ) , ntake)
end

function make_initial_seq(;
            density::Float64,radius::Float64=1.0) ::Base.Generator
    length =2*Int(round(density*radius))
    range_x=range(-radius,radius,length=length)
    range_y=range(-radius,radius,length=length)
    ((x,y) for x in range_x,y in range_y)
end

# push 2d points
function push_points!(plot_x::Array{T,1} ,
                       plot_y::Array{T,1} ;
                       initial_seq::Base.Generator ,
                       map::Function ,
                       ndrop::Int64  ,
                       ntake::Int64  ) where {T,U}
    for initial in initial_seq
        for pt in make_seq(initial,map,ndrop=ndrop,ntake=ntake)
            push!(plot_x,pt[1])
            push!(plot_y,pt[2])
        end
    end
end

# push 2d points with color
function push_points!(plot_x::Array{T,1} ,
                       plot_y::Array{T,1} ,
                       rgbs  ::Array{RGB{U},1} ;
                       make_color::Function ,
                       initial_seq::Base.Generator ,
                       map::Function ,
                       ndrop::Int64  ,
                       ntake::Int64  ) where {T,U}
    for initial in initial_seq
        rgb::RGB{U}=make_color(initial)
        for pt in make_seq(initial,map,ndrop=ndrop,ntake=ntake)
            push!(plot_x,pt[1])
            push!(plot_y,pt[2])
            push!(rgbs,rgb)
        end
    end
end


end # module DynamicalSystem

```

```GumowskiMira.jl
module GumowskiMira
export gmmap,GMType

@enum GMType f1g1 f1g2 f2g1 f2g2

function g1(x::Float64 ; μ::Float64) ::Float64
    μ*x + (2.0(1.0-μ)x^2 )/(1.0 + x^2)
end

function g2(x::Float64 ; μ::Float64) ::Float64
    μ*x + (1.0-μ)x^2*exp((1.0 - x^2)/4.0)
end

function f1(g::Function,pt::Tuple{Float64,Float64};
            μ::Float64) ::Tuple{Float64,Float64}
    xx =  pt[2] + g(pt[1],μ=μ)
    yy = -pt[1] + g(xx   ,μ=μ)
    (xx,yy)
end

function f2(g::Function,pt::Tuple{Float64,Float64};
            α::Float64,σ::Float64,μ::Float64)::Tuple{Float64,Float64}
    xx =  pt[2] + α*pt[2]*( 1.0 - σ*pt[2]^2 ) + g(pt[1],μ=μ)
    yy = -pt[1] + g(xx,μ=μ)
    (xx,yy)
end

function gmmap(gmtype::GMType;
                α::Float64=0.,σ::Float64=0.,μ::Float64) ::Function
    if      (gmtype==f1g1) pt->f1(g1,pt,μ=μ)
    elseif (gmtype==f1g2)  pt->f1(g2,pt,μ=μ)
    elseif (gmtype==f2g1)  pt->f2(g1,pt,α=α,σ=σ,μ=μ)
    elseif (gmtype==f2g2)  pt->f2(g2,pt,α=α,σ=σ,μ=μ)
    else
        error("gmtypeが不正")
    end
end

end # module GumowskiMira
```

```gm.jl
using Colors
using FixedPointNumbers
using Plots

# my modules
import DynamicalSystem
import GumowskiMira
# 

function calc_points(param::NamedTuple)
    
    function make_hsv(initial::Tuple{T,T}) where{T}
        hue= ( hypot(initial...) * 360 ) % 360
        HSV(hue,1,1)
    end
    
    initial_seq=DynamicalSystem.make_initial_seq(
                    density=param.initial_density,
                    radius =param.initial_radius)

    plot_x=Float16[]
    plot_y=Float16[]
    rgbs=RGB{N0f8}[]
    if param.with_rgb
        DynamicalSystem.push_points!(plot_x,plot_y,rgbs,
                        make_color  =make_hsv     ,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,rgbs
    else
        DynamicalSystem.push_points!(plot_x,plot_y,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,nothing
    end
end

function plot_points(param::NamedTuple ,
                plot_x::Array{T,1}  ,
                plot_y::Array{T,1}  ,
                rgbs::Union{Array{RGB{U},1},Nothing} 
                ) where{ T<:AbstractFloat , U }
    param.backend()
    if !isnothing(rgbs)
        scatter(plot_x,plot_y
                ,background_color=:black
                ,aspect_ratio=:equal
                ,size=param.size
                ,markerstrokewidth=0
                ,markeralpha=param.markeralpha
                ,markersize=param.markersize
                ,markercolor=rgbs)

    else
        scatter(plot_x,plot_y
                ,background_color=:black
                ,aspect_ratio=:equal
                ,size=param.size
                ,markerstrokewidth=0
                ,markeralpha=param.markeralpha
                ,markersize =param.markersize
                ,markercolor=:white)
    end
end

gm_map_1=GumowskiMira.gmmap(GumowskiMira.f2g1,α=0.008 , σ=0.05, μ=-0.496) # 「神話の鳥 (mythic bird)」3枚羽の翼
gm_map_2= GumowskiMira.gmmap(GumowskiMira.f2g1,α=0.009 , σ=0.05, μ=-0.80) #「神話の鳥 (mythic bird)」5枚羽の翼
gm_map_3= GumowskiMira.gmmap(GumowskiMira.f2g2,α=0.0083, σ=0.1 , μ=-0.38 ) #
gm_map_4= GumowskiMira.gmmap(GumowskiMira.f2g2,α=0.01  , σ=0.1 , μ=0.8 ) #
gm_map_5= GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.39 ) #
gm_map_6= GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.365) #
gm_map_7= GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.34 ) #

gm_map=gm_map_2

param_0_0=(
    backend        = pyplot ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = true  ,markeralpha    =1 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=0 , ntake=1 )

param_0_1=(
    backend        = unicodeplots ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = false ,markeralpha    =1 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=1000 , ntake=1 )

param_1=(
    backend        = pyplot ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = false  ,markeralpha    =0.5 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=1000 , ntake=100 )

param_2=(
    backend        = pyplot ,
    map            = gm_map ,
    size           = (700,700),
    with_rgb       = true  ,markeralpha    =0.1 , markersize =1 ,
    initial_density= 5.0   ,initial_radius =3.0 , 
    ndrop=5000 , ntake=1000 )

plot_points(calc_points( param_2 )...)

```
![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/402919/2380e8b3-7722-91cb-98e1-3a2f0ce6fb27.png)

with_rgbをfalseにするとモノクロで描画します
![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/402919/5f5b84f0-e118-2fe4-e41f-1e188414d1c5.png)

反復回数の最低値ndropを0に、取得数ntakeを1にすると初期値を描画します。
![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/402919/cc86c563-47f2-9bad-d503-779f81875fd0.png)

#第三版（複数描画
``` gm.jl
using Colors
using FixedPointNumbers
using Plots

# my modules
import DynamicalSystem
import GumowskiMira
#

function calc_points(param::NamedTuple)

    function make_hsv(initial::Tuple{T,T}) where{T}
        hue= ( hypot(initial...) * 360 ) % 360
        HSV(hue,1,1)
    end

    initial_seq=DynamicalSystem.make_initial_seq(
                    density=param.initial_density,
                    radius =param.initial_radius)

    plot_x=Float16[]
    plot_y=Float16[]
    rgbs=RGB{N0f8}[]
    if param.with_rgb
        DynamicalSystem.push_points!(plot_x,plot_y,rgbs,
                        make_color  =make_hsv     ,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,rgbs
    else
        DynamicalSystem.push_points!(plot_x,plot_y,
                        initial_seq =initial_seq  ,
                        map         =param.map    ,
                        ndrop       =param.ndrop  ,
                        ntake       =param.ntake        )
        param,plot_x,plot_y,nothing
    end
end

function plot_points(param::NamedTuple ,
                plot_x::Array{T,1}  ,
                plot_y::Array{T,1}  ,
                rgbs::Union{Array{RGB{U},1},Nothing}
                )::Plots.Plot where{ T<:AbstractFloat , U }
    if !isnothing(rgbs)
        scatter(plot_x,plot_y , markercolor=rgbs ,
                markeralpha=param.markeralpha    ,
                title =param.title)

    else
        scatter(plot_x,plot_y ,
                markeralpha=param.markeralpha ,
                title =param.title)
    end
end

param_f2g1_1=(
    title="f2g1:α=0.008,σ=0.05,μ=-0.496" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g1,
                        α=0.008 , σ=0.05, μ=-0.496) ,
                        # 「神話の鳥 (mythic bird)」3枚羽の翼
    with_rgb       = false  ,markeralpha   =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f2g1_2=(
    title="f2g1:α=0.009,σ=0.05,μ=-0.80" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g1 ,
                        α=0.009 , σ=0.05, μ=-0.80) ,
                        #「神話の鳥 (mythic bird)」5枚羽の翼
    with_rgb       = false ,markeralpha    =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f2g2_1=(
    title="f2g2:α=0.0083,σ=0.1,μ=-0.38" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g2,
                        α=0.0083, σ=0.1 , μ=-0.38 )  ,
    with_rgb       = false ,markeralpha   =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f2g2_2=(
    title="f2g2:α=0.01,σ=0.1,μ=0.8" ,
    map=GumowskiMira.gmmap(GumowskiMira.f2g2,
                        α=0.01  , σ=0.1 , μ=0.8 )  ,
    with_rgb       = false ,markeralpha    =0.3 ,
    initial_density= 10.   ,initial_radius =1. ,
    ndrop=1000 , ntake=50 )
param_f1g1_1=(
    title="f1g1:μ=0.39" ,
    map=GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.39 )  ,
    with_rgb       = false   ,markeralpha   =0.1 ,
    initial_density= 10.0   ,initial_radius =3.0 ,
    ndrop=1000 , ntake=10 )
param_f1g1_2=(
    title="f1g1:μ=0.365" ,
    map=GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.365)  ,
    with_rgb       = false  ,markeralpha    =0.1 ,
    initial_density= 10.0   ,initial_radius =3.0 ,
    ndrop=1000 , ntake=10 )
param_f1g1_3=(
    title="f1g1:μ=0.34" ,
    map=GumowskiMira.gmmap(GumowskiMira.f1g1,μ=0.34 )  ,
    with_rgb       = false  ,markeralpha    =0.1 ,
    initial_density= 10.0   ,initial_radius =3.0 ,
    ndrop=1000 , ntake=10 )


pyplot()
plot(
    plot(
        plot_points(calc_points( param_f2g1_1 )...) ,
        plot_points(calc_points( param_f2g1_2 )...) ,
        layout=(1,2) 
    ),
    plot(
        plot_points(calc_points( param_f2g2_1 )...) ,
        plot_points(calc_points( param_f2g2_2 )...) ,
        layout=(1,2) 
    ),
    plot(
        plot_points(calc_points( param_f1g1_1 )...) ,
        plot_points(calc_points( param_f1g1_2 )...) ,
        plot_points(calc_points( param_f1g1_3 )...) ,
        layout=(1,3) 
    ),
    layout=(3,1),
    aspect_ratio=:equal,
    markercolor=:black ,
    markerstrokewidth=0 ,
    markersize =1 ,
    size=(700,800)
)
```

![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/402919/255075b8-74fb-1d93-a232-e4aeb69173de.png)


# 参考
- [2次元カオス図形を芸術的に描画する(石立 喬)] (https://codezine.jp/article/detail/327)
- [Julia入門のための参考資料まとめ](https://qiita.com/yamasushi/items/e3eded4bfa42dcc4ec42)