text:
using Plots, LaTeXStrings, FourierGPE
gr(titlefontsize=12,size=(500,300),colorbar=false)

# Units:
# Units of ξ for length, 1/μ for time
# related by μ ≡ ħ²/mξ² = mc^2 where
# c = ħ/mξ is the speed of sound of
# the uniform system.

## set simulation parameters
L=(16.,16.,16.)
N=(64,64,64)
sim = Sim(L,N)
@unpack_Sim sim

γ = 0.05
tf = 4/γ
Nt = 200
t = LinRange(0.,tf,Nt)

## useful state functions
ψ0(x,y,z,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,z,0.0)/μ,0.0)+im*0.0)
healing(x,y,z,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,z,μ,g)))

## make initial state
x,y,z = X
ψi = ψ0.(x,y',reshape(z,1,1,length(z)),μ,g)
ψi = randn(N)+im*randn(N)
ϕi = kspace(ψi,sim)

@pack_Sim! sim

## Evolve in k space
sol = runsim(sim)

## plot
# function showpsi(x,y,ψ)
#     p1 = heatmap(x,y,abs2.(ψ),aspectratio=1)
#     xlabel!(L"x/\xi");ylabel!(L"y/\xi")
#     title!(L"|\psi|^2")
#     p2 = heatmap(x,y,angle.(ψ),aspectratio=1)
#     xlabel!(L"x/\xi");ylabel!(L"y/\xi")
#     title!(L"\textrm{phase}(\psi)")
#     p = plot(p1,p2,size=(600,300))
#     return p
# end

# pull out the ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg[:,:,1])

# animate a slice
# dipole decay
anim = @animate for i=1:Nt
    showpsi(x,y,xspace(sol[i],sim)[:,:,1])
end

gif(anim,"./examples/3dquench.gif",fps=30)

## animate isosurface in Makie
using Makie, AbstractPlotting

function dense(i)
    ψm = xspace(sol[i],sim)
    density = abs2.(ψm)
    pmax = maximum(density)
    return density/pmax
end

function densityfilm(Nt,saveto="examples/3dquenchisotest.gif")
scene = Scene()
tindex = Node(1)

scene = volume(lift(i -> dense(i), tindex), algorithm = :iso,show_axis=false)

record(scene, saveto, 1:Nt-10) do i
    tindex[] = i
end
end

densityfilm(Nt)


text:
### A Pluto.jl notebook ###
# v0.12.8

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using CSV, DataFrames
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.03"

# ╔═╡ 5da2b2aa-243c-11eb-395e-a3edd16d9360
data = CSV.read(datadir("purchaseData.csv"), copycols = true)

# ╔═╡ c9268920-243c-11eb-1ddf-a75015c0bb1c
typeof(data)

# ╔═╡ 802b5c32-243c-11eb-2603-778c180b1811
first(data, 6)

# ╔═╡ 802b8b4e-243c-11eb-0fb0-e31c35cc6262
describe(data)

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
df = CSV.read(datadir("purchaseDataImputed.csv"), DataFrame; delim=",")

# ╔═╡ aae3e650-23b1-11eb-330b-af6ee2f72268
size(df)

# ╔═╡ ab0e5336-23b1-11eb-2522-8fc24ac2b5b3
names(df)

# ╔═╡ ab0ecb18-23b1-11eb-0e1a-0dc82d5fe127
first(df, 6)

# ╔═╡ ab18a048-23b1-11eb-19e7-cfa7649ae028
describe(df)

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.03"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═5da2b2aa-243c-11eb-395e-a3edd16d9360
# ╠═c9268920-243c-11eb-1ddf-a75015c0bb1c
# ╠═802b5c32-243c-11eb-2603-778c180b1811
# ╠═802b8b4e-243c-11eb-0fb0-e31c35cc6262
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═aae3e650-23b1-11eb-330b-af6ee2f72268
# ╠═ab0e5336-23b1-11eb-2522-8fc24ac2b5b3
# ╠═ab0ecb18-23b1-11eb-0e1a-0dc82d5fe127
# ╠═ab18a048-23b1-11eb-19e7-cfa7649ae028
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12

ext: jl


text:
struct Interval
    start_iv
    end_iv
    interval
    Interval() = new(nothing, nothing, [])
    Interval(start_iv::Real, end_iv::Real) = new(start_iv, end_iv, start_iv:end_iv)

end

Base.minimum(iv::Interval) = iv.start_iv
Base.maximum(iv::Interval) = iv.end_iv
Base.in(x::Real, iv::Interval) = iv.start_iv isa Number && iv.end_iv isa Number ? (x ∈ iv.interval) : false
Base.isempty(iv::Interval) = iv.start_iv isa Number && iv.end_iv isa Number ? false : true
Base.issubset(x, iv::Interval) = x ⊆ iv.interval
Base.intersect(x, iv::Interval) = x ∩ iv.interval
Base.show(io::IO, iv::Interval) = iv.start_iv isa Number && iv.end_iv isa Number ? print(io, "\u301a$(iv.interval)\u301b") : print(io, "∅")


text:
function set_BLAS_threads()
    # See https://github.com/JuliaLang/julia/issues/33409, no fix even in 2021
    if (BLAS.vendor() == :openblas64)
        blas_num_threads = min(4, Int64(round(Sys.CPU_THREADS / 2)))  # even lower?
        BLAS.set_num_threads(blas_num_threads)
    end
end

function ArgParse.parse_item(::Type{Vector{Float64}}, x::AbstractString)
    return parse.(Float64, split(strip(x, [' ','\"',']','[']), ','))
end


text:
using Comonicon
using Test
using Pkg

Comonicon.disable_cache()

@testset "cache flags" begin
    @test Comonicon.Parse.CACHE_FLAG[] == false
    Comonicon.enable_cache()
    @test Comonicon.Parse.CACHE_FLAG[] == true
end

Comonicon.disable_cache()

@testset "configurations" begin
    include("configurations.jl")
end

@testset "codegen" begin
    include("codegen.jl")
end

@testset "parse" begin
    include("parse.jl")
end

@testset "build" begin
    try
        include("build.jl")
    finally
        Pkg.rm(PackageSpec(name = "Foo"))
    end
end


text:
using Test
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
splitln(str) = split(strip(str), "\n")
import ExtensibleEffects: autorun

# Task / Future
# -------------

efftask = @async @syntax_eff begin
  a = @async 3
  b = @async begin sleep(1); 5 end
  @pure a, b
end
@test fetch(efftask) == (3,5)

using Distributed
efffuture = @spawnat :any @syntax_eff begin
  a = @spawnat :any 3
  b = @spawnat :any begin sleep(1); 5 end
  @pure a, b
end
@test fetch(efffuture) == (3,5)

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
mutable struct Element{dim,N,M} <: JuAFEM.AbstractCell{dim,N,M}
    nodes::Vector{Node}
    boundaries::Vector{Element}

    parent::Union{Nothing,Element}
    children::Vector{Element}

    boundaryof::Union{Nothing,Element}
    attachedto::Vector{Element}
    ismaster::Bool
    function Element(n, p)
        sideorders = [ [] ,[[1], [2]],  [[1,2],[2,3],[3,4],[4,1]],
            [[1,4,3,2], [1,2,6,5],[2,3,7,6], [3,4,8,7], [1,5,8,4], [5,6,7,8]]]
        N = length(n)
        dim = trunc(Int, log2(N))
        M = 2 * dim
        bds = (x -> Element(n[x], nothing)).(sideorders[dim + 1])
        v = Vector{Element}(bds)
        self = new{dim,N,M}(n, v, p, Vector{Element}(), nothing, Vector{Element}(), true)
        for x in self.boundaries
            x.boundaryof = self
        end
        return  self
    end
end

function ultimate_ancestor(e::Element)
    if (e.parent === nothing) return e end
    return ultimate_ancestor(e.parent)
end

function collect_children(e::Element)
    if length(e.children) == 0 return [e] end
    return collect(Iterators.flatten(collect_children.(e.children)))
end

function collect_children(v::Union{AbstractSet{C},AbstractVector{C}}) where C <: Element
    return collect(Iterators.flatten(collect_children.(v)))
end

# 0-dim things
function JuAFEM.vertices(e::Element)
    return tuple(e.nodes...)
end

# 1-codim things
function JuAFEM.faces(e::Element)
    return tuple((vertices.(e.boundaries))...)
end

# meant are the 2-codimensional ones, is there an order convention here?
function JuAFEM.edges(e::Element)
    edg = collect(Iterators.flatten((x -> x.boundaries).(e.boundaries)))
    edg = vertices.(edg)
    # all edges exist (in 3d) twice
    return unique!(x -> Set(x), edg) # exclamation mark or better without
end

JuAFEM.default_interpolation(::Type{Element{2,4,4}})  = Lagrange{2,RefCube,1}()
JuAFEM.default_interpolation(::Type{Element{3,8,6}})  = Lagrange{3,RefCube,1}()

function Base.show(io::IO, e::Element)
    print(io, (n -> n.x).(vertices(e)))
end

function eltocell(e::Element{dim,N,M}, n::Vector{Node{dim,T}}) where {dim,N,M,T} # formerly quadtocell
    nodenumbers = (x -> findfirst(y -> y == x, n)).(e.nodes)
    return JuAFEM.Cell{dim,N,M}(tuple(nodenumbers...))
end

# only works correctly on leaves
function get_master(e::Element) # slaves point to master, master points to parent[s!!!] of a list of slaves if several exists,
    if e.ismaster == true return e end
    i = findfirst(x -> x.ismaster == true, e.attachedto)
    return e.attachedto[i]
end

getdepth(n::Nothing) = -1
getdepth(e::Element) =  1 + getdepth(e.parent)

# makes the first one a master, does not consider children etc.
function attach_trivially!(e::Element, f::Element)
    masterface = nothing
    for i in e.boundaries
        for j in f.boundaries
            if issubset(vertices(i), vertices(j)) # better issetequal?
                i.attachedto = [j]
                j.attachedto = [i]
                i.ismaster = true
                j.ismaster = false
                masterface = i
            end
        end
    end
    return masterface
end

function getcontainednodes(f::Element{1,2,2})
    if !f.ismaster error("only applicable to master elements") end
    leaves = collect_children(union(f.attachedto, [f]))
    ns = unique!(collect(Iterators.flatten(vertices.(leaves))))
    return ns
end

function getinteriornodes(f::Element{1,2,2})
    ns = getcontainednodes(f)
    ns = deleteat!(ns, findall(x -> x ∈ vertices(f), ns))
    return ns
end

function getinteriornodeinfo(f::Element{1,2,2})
    ns = getinteriornodes(f)
    sfaces = collect_children(f.attachedto)
    sfaceinfo = (n -> findall(e -> n ∈ vertices(e), sfaces)).(ns)
    return ns, (x -> sfaces[x]).(sfaceinfo)
end

function getneighbourhood(v::Vector{E}) where E <: Element
    output = copy(v)
    for e in v
        for f in e.boundaries
            neigh = (x -> x.boundaryof).(collect_children(union(get_master(f).attachedto, [get_master(f)])))
            output = vcat(output, neigh)
        end
    end
    return convert(Vector{E}, unique!(output))
end

function getnodeneighborhood(e::Element)
    cells = getneighbourhood(getneighbourhood([e])) # sufficient for up to 3D
    return unique!(collect(Iterators.flatten(vertices.(cells))))
end

# only works in dim 3?
function node_in_el(e::Element{1,2,2}, n::Node)
    v1 = e.nodes[1].x - n.x
    v2 = e.nodes[2].x - n.x
    if length(v1) == 3 && !isapprox(norm(cross(v1, v2)), 0) return false end
    s = dot(v1, v2)
    if s > 0 return false end # isapprox version of that?
    return true
end

function getvec(e::Element{1,2,2})
    return e.nodes[2].x - e.nodes[1].x
end
# only works in dim 3?
function node_in_el(e::Element{2,4,4}, n::Node)
    s = getvec.(e.boundaries)
    v = (y -> y.x - n.x).(e.nodes)
    if !isapprox(dot(cross(s[1], s[2]), v[1]), 0) return false end
    rhs = (i -> dot(s[i], v[i]) / dot(s[i], s[i])).(collect(1:4))
    if (maximum(rhs) <= 1) && (minimum(rhs) >= 0) return true end
    if (maximum(rhs) <= 0) && (minimum(rhs) >= -1) return true end
    return false
end

function containmentinfo(e::Element{2,4,4}, n::Node)
    if !node_in_el(e, n) return [] end
    bds = findall(b -> node_in_el(b, n), e.boundaries)
    if length(bds) == 0 return [e] end
    return e.boundaries[bds]
end


text:
struct ByRef{S,T}
   ref::Ref{T}

   function ByRef(::Type{S}, x::Ref{T}) where {S,T}
       return new{S,T}(x)
   end
end

@inline function ByRef(::Type{S}, ::Type{T}) where {S,T}
    return ByRef(S, Ref(zero(T)))
end

@inline function ByRef(::Type{S}, x::T) where {S,T}
    return ByRef(S, Ref(x))
end

for (S,T) in ((:Float64, :UInt64), (:Float32, :UInt32), (:Float16, :UInt16),
              (:Int64, :UInt64), (:Int32, :UInt32), (:Int16, :UInt16), (:Int8, :UInt8),
              (:UInt64, :UInt64), (:UInt32, :UInt32), (:UInt16, :UInt16), (:UInt8, :UInt8))
  @eval begin
    function ByRef(x::$S)
        u = reinterpret($T, x)
        return ByRef($S, Ref(u))
     end
  end
end

@inline sourcetype(x::ByRef{S,T}) where {S,T} = S
@inline carriertype(x::ByRef{S,T}) where {S,T} = T

@inline ref(x::ByRef{S,T}) where {S,T} = x.ref
@inline val(x::ByRef{S,T}) where {S,T} = x.ref[]

@inline refvalue(x::ByRef{S,T}) where {S,T} = reinterpret(sourcetype(x), val(x))

Base.show(io::IO, x::ByRef{S,T}) where {S,T} = show(io, refvalue(x))

function set!(x::ByRef{S,T}, value::T) where {S,T}
    x.ref[] = value
    return x
end

Base.get(x::ByRef{S,T}) where {S,T} = x.ref[]

Base.zero(::Type{ByRef{S,T}}) where {S,T} = ByRef(S,T)

Base.zero(x::ByRef{S,T}) where {S,T} = set!(x, zero(T))


text:
# * Derivatives

"""
    lagrangeder!(xⁱ, m, L′)

Calculate the derivative of the Lagrange interpolating polynomial
Lⁱₘ(x), given the roots `xⁱ`, *at* the roots, and storing the result
in `L′`.

∂ₓ Lⁱₘ(xⁱₘ,) = (xⁱₘ-xⁱₘ,)⁻¹ ∏(k≠m,m′) (xⁱₘ,-xⁱₖ)/(xⁱₘ-xⁱₖ), m≠m′,
                [δ(m,n) - δ(m,1)]/2wⁱₘ,                    m=m′

Eq. (20) Rescigno2000
"""
function lagrangeder!(xⁱ::AbstractVector, wⁱ::AbstractVector,
                      L′::AbstractMatrix)
    δ(a,b) = a == b ? 1 : 0
    n = length(xⁱ)
    for m in 1:n
        L′[m,m] = (δ(m,n)-δ(m,1))/2wⁱ[m]
        for m′ in 1:n
            m′ == m && continue
            f = 1
            for k = 1:n
                k in [m,m′] && continue
                f *= (xⁱ[m′]-xⁱ[k])/(xⁱ[m]-xⁱ[k])
            end
            L′[m,m′] = f/(xⁱ[m]-xⁱ[m′])
        end
    end
end

function diff!(D, L̃, L′, wⁱ, nⁱ)
    n = size(L′,1)
    for m = 1:n
        for m′ = 1:n
            D[m,m′] = nⁱ[m]*nⁱ[m′]*dot(L̃[m,:],wⁱ.*L′[m′,:])
        end
    end
    D
end

function diff(B::FEDVR{T}, n::Integer, i::Integer) where T
    o = B.order[i]

    # L′ contains derivatives of un-normalized basis functions at
    # the quadrature roots.
    L′ = Matrix{T}(undef, o, o)
    lagrangeder!(@elem(B,x,i), B.wⁱ[i],L′)

    # D contains ⟨ξᵢ|χⱼ′⟩ where ξᵢ = χᵢ⁽ⁿ⁻¹⁾
    D = similar(L′)
    L̃ = n == 1 ? Matrix{T}(I, size(L′)...) : -L′ # ∂ᴴ = -∂
    diff!(D, L̃, L′, B.wⁱ[i], @elem(B,n,i))
    i ≥ B.i₀ && (D ./= √(B.eiϕ)) # TODO Check if correct
    D
end

difffun(B::FEDVR, n::Integer) = i -> diff(B,n,i)
difffun(B::RestrictedFEDVR, n::Integer) = i -> diff(parent(B),n,i)

derop!(A, B::FEDVROrRestricted, n::Integer) =
    set_elements!(difffun(B,n), A, B)

ContinuumArrays.MemoryLayout(::Type{<:BasisOrRestricted{<:FEDVR}}) = ContinuumArrays.BasisLayout()
ContinuumArrays.MemoryLayout(::Type{<:AdjointBasisOrRestricted{<:FEDVR}}) = ContinuumArrays.AdjointBasisLayout{ContinuumArrays.BasisLayout}()

@materialize function *(Ac::AdjointBasisOrRestricted{<:FEDVR},
                        D::Derivative,
                        B::BasisOrRestricted{<:FEDVR})
    T -> begin
        Matrix(undef, B, T)
    end
    dest::AbstractMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        derop!(dest, B, 1)
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:FEDVR},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        D::Derivative,
                        B::BasisOrRestricted{<:FEDVR})
    T -> begin
        Matrix(undef, B, T)
    end
    dest::AbstractMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        derop!(dest, B, 2)
    end
end


text:
using PyCall
using StatsBase

mutable struct PIDForest
    forest

	function PIDForest(parameters)
		py"""
		from pidforest.forest import Forest

		def construct_pidforest(kwargs):
			kwargs['sample_axis'] = 1
			kwargs['threshold'] = 0
			return Forest(**kwargs)
		"""
		new(py"construct_pidforest"(parameters))
	end
end

function StatsBase.fit!(model::PIDForest, X::Array{T, 2}) where T<:Real
	# no need for data transposition as the forest is learned on column major
    model.forest.fit(X)
end

function StatsBase.predict(model::PIDForest, x; pct=50) where T<:Real
	# anomaly scores correspond to percentile `pct` over trees
	# the lower the score the more anomalous the sample should be
    -model.forest.predict(x; pct=pct)[end]
end


text:
##------------------------------##
##实现从下位机读取麦克风阵列数据
##1、调用MATLAB函数读取
##2、c++封装的SO读取
##3、自己生成随机数
##------------------------------##
##name = "线路输入 (LIS_MIIS_8X8_192K)"
export  MicrophoneParameter,
        init_matlab,
        getdata

struct MicrophoneParameter{T1<:String, T2<:Int}
    name::T1 #麦克风阵列名称
    fs::T2  #采样频率
    sample::T2  #一帧采样数
    numchannels::T2    #通道数
end

function init_matlab(parameter::MicrophoneParameter)
    s = MATLAB.MSession(0)
    put_variable(s, :name, parameter.name)
    eval_string(s, "H = audioDeviceReader")
    eval_string(s, "H.Device = name ")
    eval_string(s, "H.NumChannels = $(parameter.numchannels)") #通道数要跟实际设备中搜集到的数据一致
    eval_string(s, "H.SampleRate = $(parameter.fs)")
    eval_string(s, "H.SamplesPerFrame = $(parameter.sample)")
    eval_string(s, "H.BitDepth = '16-bit integer' ")
    eval_string(s, "H.OutputDataType = 'double' ")
    return s
end

function (hardware::MicrophoneParameter)(select::Symbol)
    if select == :MATLAB
        return init_matlab(hardware)
    elseif select == :C++
        error("还没实现C++方法")
    elseif select == :RandomData
        error("还没实现随机数生成方法")
    else
        error("这儿没有实现其他方法")
    end
end

function getdata(s::MSession;spin = false, transform_matrix = [2, 1, 3, 4],
    ishilbert = false)
    eval_string(s, "data = H()")
    data = jarray(get_mvariable(s, :data))
    if spin
        data = view(data, :, transform_matrix)
    end
    if ishilbert
        data = hilbert(data)
    end
    return view(data, :, :)
end


text:
using CrossEntropyMethod
using Distributions, Random
using Plots
# The goal of this example is to show how to estimate the optimal importance sampling distribution

# Actual distibution
base_dist = Categorical([0.74899, 0.2, 0.05, 0.001, 0.00001])

# starting sampling distribution
is_dist_0 = Dict{Symbol, Vector{Sampleable}}(:x => [Categorical(5)])

# Sample a random time series from the starting distribution
rand(is_dist_0)

# Define a loss function that identifies the values of 3,4,5 as those of intereste
function l(d, s)
    v = s[:x][1]
    -(v in [3,4,5])
end

# Define the likelihood ratio weighting function
function w(d, s)
    v = s[:x][1]
    exp(logpdf(base_dist, v) - logpdf(d, s))
end

# show what the weighting function does
s = Dict(:x => [5])
w(is_dist_0, s)

# Run the optimization
is_dist_opt = cross_entropy_method(l, is_dist_0; max_iter = 4, N=10000, weight_fn = w, verbose = true)

# See the final distribution.
is_dist_opt[:x][1].p


text:
function show_volume(vdata::Array{Float64, 3}; cmap::Matrix{Int64}=[(i==1)*(j-1) + (i==3)*(256-j) for i in 1:3, j in 1:256], default_cmap="viridis", default_renderstyle="iso", clim1=0, clim2=1, isothreshold=0.15)
    HTML{String}(
"""<!DOCTYPE html>

	<title>three.js webgl - volume rendering example</title>
	<meta charset="utf-8">
	<meta name="viewport" content="width=device-width, user-scalable=no, minimum-scale=1.0, maximum-scale=1.0">
	<link type="text/css" rel="stylesheet" href="main.css">

""")
end


text:
using Images, ImageFeatures, ImageDraw, ImageMorphology

#### Example 1

img = Gray.(load("sample-images/cat-3417184_640.jpg"))
img_f = Float16.(restrict(img))

new_img = Gray.(hcat(
    Float32.(img_f) .* (~fastcorners(img_f, 12, 0.15)),
    Float32.(img_f) .* (~fastcorners(img_f, 12, 0.25)),
    Float32.(img_f) .* (~fastcorners(img_f, 12, 0.35))
))

new_img = Gray.(hcat(
    Float32.(img_f) .* (~fastcorners(img_f, 12, 0.15)),
    Float32.(img_f) .* (~imcorner(img_f, method=harris)),
    Float32.(img_f) .* (~imcorner(img_f, method=shi_tomasi))
))

imshow(new_img)

### Example 2

img = Gray.(load("sample-images/newspaper-37782_640.png"))
img_f = Float16.(img)

new_img = Gray.(hcat(
    img_f,
    Float16.(img_f) .* (~fastcorners(img_f, 12, 0.15))
))

imshow(new_img)

### Example 3
img = Gray.(restrict(load("sample-images/cat-3417184_640.jpg")))

img = restrict(load("sample-images/board-157165_640.png"))
img_f = Float16.(Gray.(img))

img_harris = copy(img)
img_harris[dilate(imcorner(img_f, method=harris)) .> 0.01] = colorant"yellow"

img_shi = copy(img)
img_shi[dilate(imcorner(img_f, method=shi_tomasi)) .> 0.01] = colorant"yellow"

img_rosenfield = copy(img)
img_rosenfield[dilate(imcorner(img_f, method=kitchen_rosenfeld)) .> 0.01] = colorant"yellow"

img_fast = copy(img)
img_fast[dilate(fastcorners(img_f, 12, 0.05)) .> 0.01] = colorant"yellow"

new_img = vcat(
    hcat(
        img_harris,
        img_shi),
    hcat(
        img_rosenfield,
        img_fast)
)

new_img[Int(size(new_img, 1) / 2), :] = colorant"yellow"
new_img[:, Int(size(new_img, 2) / 2)] = colorant"yellow"

imshow(new_img)

### Example 4
img = Gray.(restrict(load("sample-images/cat-3417184_640.jpg")))
img_f = Float16.(Gray.(img))

img_harris = copy(img)
img_harris[dilate(imcorner(img_f, method=harris)) .> 0.01] = colorant"yellow"

img_shi = copy(img)
img_shi[dilate(imcorner(img_f, method=shi_tomasi)) .> 0.01] = colorant"yellow"

img_rosenfield = copy(img)
img_rosenfield[dilate(imcorner(img_f, method=kitchen_rosenfeld)) .> 0.01] = colorant"yellow"

img_fast = copy(img)
img_fast[dilate(fastcorners(img_f, 12, 0.05)) .> 0.01] = colorant"yellow"

new_img = vcat(
    hcat(
        img_harris,
        img_shi),
    hcat(
        img_rosenfield,
        img_fast)
)

imshow(new_img)

### Example 5

img = Gray.(restrict(load("sample-images/cat-3417184_640.jpg")))
img_f = Float16.(Gray.(img))

img_harris = copy(img)
img_harris[dilate(imcorner(img_f, Percentile(95), method=harris)) .> 0.01] = colorant"yellow"

img_shi = copy(img)
img_shi[dilate(imcorner(img_f, Percentile(95), method=shi_tomasi)) .> 0.01] = colorant"yellow"

img_rosenfield = copy(img)
img_rosenfield[dilate(imcorner(img_f, Percentile(95), method=kitchen_rosenfeld)) .> 0.01] = colorant"yellow"

img_fast = copy(img)
img_fast[dilate(fastcorners(img_f, 12, 0.05)) .> 0.01] = colorant"yellow"

new_img = vcat(
    hcat(
        img_harris,
        img_shi),
    hcat(
        img_rosenfield,
        img_fast)
)

imshow(new_img)

### Performance comparison

img = restrict(load("sample-images/cat-3417184_640.jpg"))
img_f = Float16.(Gray.(img))

@btime fastcorners(img_f, 12, 0.15);
@btime fastcorners(img_f, 12, 0.05);
@btime imcorner(img_f, method=harris);
@btime imcorner(img_f, Percentile(95), method=harris);
@btime imcorner(img_f, method=shi_tomasi);
@btime imcorner(img_f, method=kitchen_rosenfeld);


text:
module significance_testing_module

using Distributions
using Statistics
using RCall

"""
# ____________________________
# Effects distribution fitting

`best_fitting_distribution(data::Array{Float64,1})`

Determine the best fitting distribution to model the data and heuristically estimate the p-values.

# Distributions
- [Bernoulli](https://en.wikipedia.org/wiki/Bernoulli_distribution)
- [Beta](https://en.wikipedia.org/wiki/Beta_distribution)
- [Binomial](https://en.wikipedia.org/wiki/Binomial_distribution)
- [Categorical](https://en.wikipedia.org/wiki/Categorical_distribution)
- [Dicscrete Uniform](https://en.wikipedia.org/wiki/Discrete_uniform_distribution)
- [Exponential](https://en.wikipedia.org/wiki/Exponential_distribution)
- [Gaussian](https://en.wikipedia.org/wiki/Normal_distribution)
- [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution)
- [Geometric](https://en.wikipedia.org/wiki/Geometric_distribution)
- [Laplace](https://en.wikipedia.org/wiki/Laplace_distribution)
- [Pareto](https://en.wikipedia.org/wiki/Pareto_distribution)
- [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution)
- [Inverse Gaussian or Wald](https://en.wikipedia.org/wiki/Inverse_Gaussian_distribution)
- [Uniform](https://en.wikipedia.org/wiki/Uniform_distribution_(continuous))

# Examples
```
using Distributions
using UnicodePlots
x1 = Distributions.rand(Distributions.Normal(0, 1), 100)
UnicodePlots.histogram(x1)
significance_testing_module.best_fitting_distribution(x1)
x2 = Distributions.rand(Distributions.Laplace(0, 1), 100)
UnicodePlots.histogram(x2)
significance_testing_module.best_fitting_distribution(x2)
```
"""
function best_fitting_distribution(data::Array{Float64,1})
    DIST_NAMES =   [Distributions.Bernoulli, Distributions.Beta, Distributions.Binomial, Distributions.Categorical,
                    Distributions.DiscreteUniform, Distributions.Exponential, Distributions.Normal, Distributions.Gamma,
                    Distributions.Geometric, Distributions.Laplace, Distributions.Pareto, Distributions.Poisson,
                    Distributions.InverseGaussian, Distributions.Uniform]
    DIST_INSTANCES = [try Distributions.fit_mle(D, data); catch nothing; end for D in DIST_NAMES]
    NEG_LOGLIK = [try -sum(Distributions.logpdf.(D, data)); catch nothing; end for D in DIST_INSTANCES]
    # NEG_LOGLIK = [try sum(Distributions.pdf.(D, data)); catch nothing; end for D in DIST_INSTANCES]
    DISTRIBUTIONS_DF = hcat((DIST_NAMES[NEG_LOGLIK .!= nothing],
                            DIST_INSTANCES[NEG_LOGLIK .!= nothing],
                            NEG_LOGLIK[NEG_LOGLIK .!= nothing])...)
    D = try
        (DISTRIBUTIONS_DF[argmin(DISTRIBUTIONS_DF[:,3]), 2], DISTRIBUTIONS_DF[argmin(DISTRIBUTIONS_DF[:,3]), 1])
    catch
        (nothing, "Failure to fit into any of the distributions tested.")
    end
    return(D)
end

"""
# ________________________________________________________________________________________
# Heuristic estimates of p-values and -log10(p-values) using the best-fitting distribution

`estimate_pval_lod(data::Array{Float64,1})`

# Output
1. p-values
2. -log10(p-values)

# Examples
```
using Distributions
using UnicodePlots
x1 = Distributions.rand(Distributions.Normal(0, 1), 100)
UnicodePlots.histogram(x1)
PVAL, LOD = significance_testing_module.estimate_pval_lod(x1)
UnicodePlots.scatterplot(LOD)
x2 = Distributions.rand(Distributions.Laplace(0, 1), 100)
UnicodePlots.histogram(x2)
PVAL, LOD = significance_testing_module.estimate_pval_lod(x2)
UnicodePlots.scatterplot(LOD)
```
"""
function estimate_pval_lod(data::Array{Float64,1})
    D, D_name = best_fitting_distribution(data)
    println(string("Distribution used: ", D_name))
    if (D == nothing) | (std(data) == 0)
        PVAL = repeat([1.0], inner=length(data))
        LOD = PVAL .- 1.0
    else
        PVAL = [Distributions.cdf(D, x) <= Distributions.ccdf(D, x) ? 2*Distributions.cdf(D, x) : 2*Distributions.ccdf(D, x) for x in data]
        LOD = -log.(10, PVAL)
        # UnicodePlots.scatterplot(LOD)
    end
    return(PVAL, LOD)
end

"""
# _______________________________________________________________________________________________
# Plot -log10(p-values) with a given Bonferroni threshold at user-defined false positive rate (α)

`plot_manhattan(;chrom::Array{Any,1}, pos::Array{Int64,1}, pval::Array{Float64,1}, lod::Array{Float64,1}, fpr::Float64=0.01, png_fname::String=string("Manhattan_plot-", hash(rand()), ".jpeg"))`

# Inputs
1. *chrom* [Array{Any,1}]: vector of chromosome names
2. *pos* [Array{Int64,1}]: vector of allele/SNP positions in base-pairs
3. *pval* [Array{Float64,1}]: vector of p-values
4. *lod* [Array{Float64,1}]: vector of -log10(p-values)
5. *fpr* [Float64=0.01]: false positive rate threshold or the alpha (default=0.01)
6. *png_fname* [String]: output filename of the Manhattan plot in portable network format (png) (default=string("Manhattan_plot-", hash(rand()), ".jpeg"))

# Output
1. Manhattan plot in portable network graphics (png) format

# Examples
```
using Distributions
chrom = repeat([1,2,3,"X", "Z"], inner=20)
pos = repeat(collect(1:20), outer=5)
eff = Distributions.rand(Distributions.Normal(0, 1), 100)
eff[[10,30,60,90]] = [10.0,20.0,30.0,40.0]
pval, lod = significance_testing_module.estimate_pval_lod(eff)
significance_testing_module.plot_manhattan(chrom=chrom, pos=pos, pval=pval, lod=lod)
```
"""
function plot_manhattan(;chrom::Array{Any,1}, pos::Array{Int64,1}, pval::Array{Float64,1}, lod::Array{Float64,1}, fpr::Float64=0.01, png_fname::String=string("Manhattan_plot-", hash(rand()), ".png"))
	### plot in R
    @rput chrom;
    @rput pos;
    @rput pval;
	@rput lod;
	@rput fpr;
	@rput png_fname;
	### remove missing values and convert infinite LOD into the maximum non-infinite LOD
    R"gwas_out = data.frame(CHROM=as.factor(unlist(chrom)), POS=as.numeric(unlist(pos)), PVALUES=as.numeric(unlist(pval)), LOD=as.numeric(unlist(lod)))";
	R"gwas_out = droplevels(gwas_out[gwas_out$CHROM != 'Intercept', ])";
	R"gwas_out$LOD[is.na(gwas_out$LOD)] = 0.0";
	R"gwas_out$LOD[is.infinite(gwas_out$LOD)] = max(gwas_out$LOD[!is.infinite(gwas_out$LOD)])";
    ### arranging chromosomes or scaffold end to end
    R"gwas_out$POS_CONTINUOUS = gwas_out$POS";
    R"contigs = levels(gwas_out$CHROM)";
    R"min_pos = aggregate(POS ~ CHROM, data=gwas_out, FUN=min)"
    R"max_pos = aggregate(POS ~ CHROM, data=gwas_out, FUN=max)"
    R"idx_chrom = max_pos$CHROM[2:nrow(max_pos)]"
    R"add_us = cumsum(max_pos$POS)[1:nrow(max_pos)-1]"
    R"for (i in 1:length(idx_chrom)){
            chr = idx_chrom[i]
            add = add_us[i]
            gwas_out$POS_CONTINUOUS[gwas_out$CHROM == chr] = gwas_out$POS_CONTINUOUS[gwas_out$CHROM == chr] + add
        }"
	### coloring
	R"colours = RColorBrewer::brewer.pal(6, 'Paired')";
	R"colours_points = rep(colours[c(1,3)], times=ceiling(length(contigs)/2))";
	R"colours_labels = rep(colours[c(2,4)], times=ceiling(length(contigs)/2))";
	R"bonferroni_threshold = -log10(fpr/nrow(gwas_out))";
	### plot -log10(pvalue) points
	R"png(png_fname, width=2000, height=800)";
	R"layout(matrix(c(1,1,1,2), nrow=1, byrow=TRUE))";
	R"par(cex=1.5)";
	R"plot(x=gwas_out$POS_CONTINUOUS, y=gwas_out$LOD, col=colours_points[as.numeric(gwas_out$CHROM)], ylim=c(0, max(c(gwas_out$LOD, bonferroni_threshold))), pch=20, main=paste0(' Bonferroni Threshold at ', round(bonferroni_threshold, 2)), xaxt='n', xlab='Chromosome or Scaffold', ylab=expression(log[10]~(pvalue)), las=1)";
	R"xaxis_pos = aggregate(gwas_out$POS_CONTINUOUS ~ gwas_out$CHROM, FUN=median)";
	R"mtext(side=1, at=xaxis_pos[,2], text=xaxis_pos[,1], cex=1.5)";
	R"abline(h=bonferroni_threshold, lty=2, lwd=2, col=colours[5])";
	R"grid()";
	# QQ plot (p-values are assumed to be uniformly distributed)
	# observed_pval = runif(100); observed_pval = observed_pval[order(observed_pval, decreasing=FALSE)] ### TEST
	# sort the observed p-values
	R"observed_pval = gwas_out$PVALUES[order(gwas_out$PVALUES, decreasing=FALSE)]";
	# calculate the pdf of the observed p-values, i.e. the probability density for a p-value interval which corresponds to each observed p-value
	R"observed_pval_density = density(observed_pval, n=length(observed_pval), from=0, to=1) ### calculate the density of the observed p-values between 0 and 1 because we are dealing with p-values which are probabilities!";
	# calculate the cummulative probabilities of the observed p-values based on the pdf: where Prob(pvalue) = Density(pvalue | pvalue_interval) * pvalue_interval
	R"observed_pval_cumprob = cumsum(observed_pval_density$y * (observed_pval_density$x[2]-observed_pval_density$x[1]))";
	# calculate the expected quantiles based on the cummulative probabilities of the observed p-values
	R"expected_pval = qunif(p=observed_pval_cumprob, min=0, max=1) ### calculate the expected quantiles based on the observed cummulative probabilities";
	# transform into -log10 scale
	R"observed_lod = -log10(observed_pval + 1e-200)";
	R"expected_lod = -log10(expected_pval + 1e-200)";
	# plot
	R"plot(x=c(min(observed_lod), max(observed_lod)), y=c(min(observed_lod), max(observed_lod)), type='n', , main='QQ Plot', xlab=expression(Expected~~log[10]~(pvalue)), ylab=expression(Observed~~log[10]~(pvalue)))";
	R"points(x=expected_lod, y=observed_lod, type='p', pch=20, col=colours[5])";
	R"lines(x=c(0, max(observed_lod)), y=c(0, max(observed_lod)), lty=2, lwd=2, col='gray')";
	R"grid()";
	R"dev.off()";
	return(0)
end

end # end of significance_testing_module

max_stars_repo_name: jeffersonfparil/GWAlpha


text:
### Factor Definition
struct ExamplePrior{T} <: IncrementalInference.FunctorSingleton where T <: SamplableBelief
  Z::T
end

### Sampler and Residual Functions
getSample(s::ExamplePrior, N::Int=1) = (reshape(rand(s.Z,N),:,N), )

### Packed Converters
struct PackedExamplePrior <: PackedInferenceType
  Z::String
  PackedExamplePrior() = new()
  PackedExamplePrior(z::AS) where {AS <: AbstractString} = new(z)
end
function convert(::Type{PackedExamplePrior}, d::ExamplePrior)
  return PackedExamplePrior(string(d.Z))
end
function convert(::Type{ExamplePrior}, d::PackedExamplePrior)
  return ExamplePrior(extractdistribution(d.Z))
end

### JSON Converters
"""
Converter: ExamplePrior -> PackedExamplePrior::Dict{String, Any}
"""
function convert(::Type{Dict{String, Any}}, prior::ExamplePrior)
    z = convert(Dict{String, Any}, prior.Z)
    return JSON.parse(JSON.json(Packed_Factor([z], "ExamplePrior")))
end

"""
Converter: PackedPrior::Dict{String, Any} -> ExamplePrior
"""
function convert(::Type{ExamplePrior}, prior::Dict{String, Any})
    # Genericize to any packed type next.
    z = prior["measurement"][1]
    z = convert(_evalType(z["distType"]), z)
    return ExamplePrior(z)
end


text:
import Requires

function __init__()::Nothing
    Requires.@require MLJ="add582a8-e3ab-11e8-2d5e-e98b27df1bc7" begin
        Requires.@require MLJBase="a7f614a8-145f-11e9-1d2a-a57a1082229d" begin
            Requires.@require MLJModels="d491faf4-2d78-11e9-2867-c94bc002c0b7" begin
                import .MLJ
                import .MLJBase
                import .MLJModels
                include("toplevel/conditionally-loaded/PredictMD_MLJ/PredictMD_MLJ.jl")
                import .PredictMD_MLJ
            end
        end
    end

    print_welcome_message()
    return nothing
end


text:
using Graphs, SparseArrays

struct SparseMatrixDiGraph{T,TW} <: AbstractGraph{T}
	A::SparseMatrixCSC{Bool, T}
	X::SparseMatrixCSC{T,T}
	W::TW
end

struct IndexedEdge{T}
	src::T
	dst::T
	idx::T
end

function SparseMatrixDiGraph(A::AbstractMatrix{Bool}, W)
	size(A,1) != size(A,2) && throw(ArgumentError("Matrix should be square"))
	any(A[i,i] for i=1:size(A,1)) && throw(ArgumentError("Self loops are not allowed"))
	A = sparse(A)
	length(W) != nnz(A) && throw(ArgumentError("Number of properties different from number of edges"))
	X = sparse(SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, collect(1:nnz(A)))')
	SparseMatrixDiGraph(A, X, W)
end

function SparseMatrixDiGraph(A::AbstractMatrix)
	SparseMatrixDiGraph(SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, fill(true, length(A.nzval))), A.nzval)
end

Graphs.edges(g::SparseMatrixDiGraph) = (IndexedEdge{Int}(i,g.A.rowval[k],k) for i=1:size(g.A,2) for k=nzrange(g.A,i))

Base.eltype(g::SparseMatrixDiGraph{T}) where T = T

Graphs.edgetype(g::SparseMatrixDiGraph{T}) where T = IndexedEdge{T}

Graphs.has_edge(g::SparseMatrixDiGraph, i, j) = g.A[i,j]

Graphs.has_vertex(g::SparseMatrixDiGraph, i) = i ∈ 1:size(g.A, 2)

Graphs.inneighbors(g::SparseMatrixDiGraph, i) = @view g.A.rowval[nzrange(g.A,i)]

Graphs.outneighbors(g::SparseMatrixDiGraph, i) = @view g.X.rowval[nzrange(g.X,i)]

inedges(g::SparseMatrixDiGraph, i) = (IndexedEdge{Int}(i,g.A.rowval[k],k) for k=nzrange(g.A,i))

outedges(g::SparseMatrixDiGraph, i) = (IndexedEdge{Int}(g.X.rowval[k],i,g.X.nzval[k]) for k=nzrange(g.X,i))

property(g::SparseMatrixDiGraph, e::IndexedEdge) = g.W[e.idx]

Graphs.ne(g::SparseMatrixDiGraph) = nnz(g.A)

Graphs.nv(g::SparseMatrixDiGraph) = size(g.A,2)

Graphs.vertices(g::SparseMatrixDiGraph) = 1:size(g.A,2)

Graphs.is_directed(::Type{SparseMatrixDiGraph{T,TW}}) where {T,TW} = true

Graphs.is_directed(g::SparseMatrixDiGraph) = true

Base.zero(g::SparseMatrixDiGraph) = SparseMatrixDiGraph(zero(g.A), similar(g.W, 0))


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 1f65fc80-a83b-11eb-2583-89bff1ea8372
html"""

Section 3.3

 Why we can't predict the weather

Lecture Video

"""

# ╔═╡ 88b46d2e-220e-11eb-0f7f-b3f523f0214e
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="LaTeXStrings", version="1"),
        Pkg.PackageSpec(name="Roots", version="1"),
		Pkg.PackageSpec(name="DifferentialEquations", version="6"),
    ])

    using Plots, PlutoUI, LaTeXStrings, Roots, DifferentialEquations
	using LinearAlgebra
end

# ╔═╡ f9dbf2e8-574d-4846-af48-f7e5a82a1afc
TableOfContents(aside=true)

# ╔═╡ a2b2eae2-762e-41c3-a546-0caedc79db7d
md"""
# Why we can predict the climate but can't predict the weather
"""

# ╔═╡ 91b34e59-e4d2-45b7-9e2b-1b95748e5d6a
md"""
## Weather vs. climate
"""

# ╔═╡ 22bddfca-fd1c-4dc7-85f3-f864c3b53407
md"""
Will it rain tomorrow? We should be able to predict that using observations about the current state of the atmosphere and mathematical models in the form of partial differential equations that predict how the wind will move the air around, what the temperature will be, and how much water can precipitate out.

What about predicting rain next week? Next year?

We know from experience that we can only predict the weather in the short term, up to a week or so (and still only with a certain probability). Predicting even one month out doesn't seem possible.

On the other hand, we can talk about the *climate* and how it behaves on much longer time scales, e.g. "the average temperature in in Boston in April is 55°F". The short-term weather has a significant random (or apparently random) component, but becomes much more predictable when we look at statistical properties like averages and standard deviations; this is then called **climate**.
"""

# ╔═╡ 6f554dc0-220c-11eb-341f-1b66de841c86
md"""
# Nonlinear dynamics: stability and bifurcations
"""

# ╔═╡ 2c6f4b31-2f08-4e24-8e31-ad80d2700fa8
md"""
In this notebook we will see that the simplest possible model of the weather, a set of 3 coupled, nonlinear ODEs known as the **Lorenz equations**, can have unpredictable behaviour, known as
**(deterministic) chaos**. We will build up to that by looking at two simpler models of dynamics, in 1D and 2D.
"""

# ╔═╡ 822c015a-8bb3-4432-a2d4-2e2ba4f906a6
md"""
## Reminder: Ordinary Differential Equations (ODEs)
"""

# ╔═╡ 816de40c-220c-11eb-2d3a-23e6267cd529
md"""
Recall that we are looking at **differential equations**, with the goal of modelling the evolution in time of the climate.

The simplest such models are **ordinary differential equations** (ODEs), where one or a few continuous variables evolve continuously in time, with a model that specifies their instantaneous rate of change as a function of their current values.

The general form of an ODE (that does not depend explicitly on time, i.e. that is **autonomous**), is

$$\frac{dx(t)}{dt} = f(x(t)),$$

or

$$\dot{x}(t) = f(x(t)),$$

with an initial condition $x(t=0) = x_0$.

Here $\dot{x}(t)$ denotes the derivative of the function $t \mapsto x(t)$ at time $t$.
"""

# ╔═╡ d21f6358-220c-11eb-00b8-03fe0746fdc9
md"""
We have seen that the simplest numerical method to solve such an equation is the **(forward) Euler method**, in which we convert this equation into an explicit time-stepping routine: we take a small time step of length $h$ and approximate the derivative as

$$\frac{dx(t)}{dt} \simeq \frac{x(t + h) - x(t)}{h},$$

giving the numerical method

$$x_{n+1} = x_{n} + h f(x_n),$$

where $x_n$ is the approximation of the true solution $x(t_n)$ at the $n$th time step $t_n$.

"""

# ╔═╡ 735bb47e-220d-11eb-1ba6-8d591a935857
md"""
# 1D: Modelling bacterial growth
"""

# ╔═╡ 79b148b6-220d-11eb-2785-25d05068aeed
md"""
Let's use this to simulate a simple nonlinear ODE describing the dynamics of a population of bacteria. The bacteria reproduce at a rate $\lambda$ *provided* that there is sufficient food, in which case we would have $\dot{x} = \lambda x$. But the need for enough food will put a *limit* on the sustainable population to a value $K$, sometimes called the **carrying capacity**.

The simplest model for the combined effect of growth and saturation is as follows:

$$\dot{x} = \lambda \, x \, (K - x).$$

When $x$ is close to $0$, the growth rate is $\lambda$, but that rate decreases as $x$ increases.

(This is sometimes called the [**logistic**](https://en.wikipedia.org/wiki/Logistic_function#Logistic_differential_equation) differential equation, although that name does not seem particularly helpful.)
"""

# ╔═╡ d43a869e-220d-11eb-009c-a119de57e7da
md"""
Our goal is to use computational thinking, but we will actually not be interested so much in the exact dynamics in time, but rather in the **qualitative** features of the behaviour of the system. For example, at long times (formally $t \to \infty$) does the population get arbitrarily large? Or does it, for example, oscillate around a particular value? Or does it converge to a particular size? This forms the subject of **nonlinear dynamics** or **dynamical systems** theory.
"""

# ╔═╡ 088eabb2-220e-11eb-1ac0-df419b87f39a
md"""
Let's simulate the system using the Euler method to try to guess the answer to this question. Note that there are many much more sophisticated methods for solving ODEs collected in the DifferentialEquations.jl package. We should never use the Euler method in practice, but should use a tried and tested library instead, and algorithms that provide much better accuracy in the solutions, if we are interested in faithful numerical results.
"""

# ╔═╡ a44cc180-22d4-11eb-2129-211d024c4921
md"""
It's useful to rescale the variables to the simplest form:

$$\dot{x} = x \, (1 - x).$$

[We can do so by defining new adimensional space and time variables $x'$ and $t'$ as follows: $x = K \, x'$, giving $K \dot{x'} = \lambda K \, x' (K - K x')$, and then $t' = \lambda K \, t$, giving $dx' / dt' = x' (1 - x')$.]
"""

# ╔═╡ 57566646-1b4c-4098-abf7-eb242d01ca81
md"""
We define a function representing the right-hand side of the ODE:
"""

# ╔═╡ 6aa1d2f4-220e-11eb-06e0-c346f74aa018
logistic(x) = x * (1 - x)

# ╔═╡ b6d0dd18-22d4-11eb-2083-0b02de030579
md"""
Let's simulate this with the Euler method and plot the trajectory $x(t)$ as a function of time $t$:
"""

# ╔═╡ 6e90ea26-220e-11eb-0c65-bf52b3d2e195
results = euler(logistic, 0.5, 0.01, 20.0)

# ╔═╡ 2d1c6abf-3b15-4638-88c5-89b5d0585c98
md"""
Normally we would *not* choose the Euler method to calculate actual trajectories, but it is enough for our purposes of understanding the *qualitative* behaviour of the equations.
"""

# ╔═╡ f055be76-22d4-11eb-268a-bf70c7d8f1a1
md"""
We see that for this particular initial condition, the solution seems to settle down to a fixed value after some time, and then remains at that value thereafter.

Such a value is called a **fixed point**, a **stationary point**, or a **steady state** of the ODE.
"""

# ╔═╡ 96e22792-220e-11eb-2729-63964507b5f2
begin
	plot(results.ts, results.xs, size=(400, 300), leg=false, xlabel=L"t", ylabel=L"x(t)", lw=3)
	scatter!([(results.ts[1], results.xs[1])])
	ylims!(0.4, 1.1)
end

# ╔═╡ 90ccb392-2216-11eb-1fd8-83b7d7c16b54
md"""
## Qualitative behaviour: Fixed points and their stability
"""

# ╔═╡ 0bf302de-4a0f-4341-b010-d3606e237acc
# gr(fmt=:png, dpi=150, size=(300, 200))
# gr(fmt=:svg)

gr(fmt=:png, dpi=300, size=(400, 300))

# ╔═╡ bef9c2e4-220e-11eb-24d8-bd618d2985ea
md"""
 Let's see what happens for other initial conditions:
"""

# ╔═╡ 17317eb8-22d5-11eb-0f46-2bf1bddd36bb
md"""
x₀ = $(@bind x0 Slider(-0.9:0.001:3.0, default=0.5, show_value=true))
"""

# ╔═╡ 2300b29a-22d5-11eb-3c99-bdec0e5a2685
let
	p = plot(xlabel=L"t", ylabel=L"x(t)",  leg=false, ylim=(-1, 2))

	results = euler(logistic, x0, 0.01, 5.0)

	plot!(results.ts, results.xs, alpha=1, lw=3)
	scatter!([0.0], [x0])

	hline!([0.0], ls=:dash)
	# as_svg(p)
end

# ╔═╡ ae862e10-22d5-11eb-3d75-1d748cf86944
md"""
To get an overview of the behaviour we can draw all the results on a single graph:
"""

# ╔═╡ c5bf15b1-a540-478b-acec-1aa47aad13d1
let
	p = plot(xlabel=L"t", ylabel=L"x(t)", leg=false, ylim=(-1, 2))

	# for x0 in -0.5:0.05:2.0
#for x0 in 0.0:0.1:2.0
	for x0 in -0.5:0.1:2.0
	results = euler(logistic, x0, 0.05, 5.0)

# 	for x0 in -0.5:0.05:2.0
# 		results = euler(logistic, x0, 0.01, 10.0)

		plot!(p, results.ts, results.xs, alpha=0.8, lw=1, arrow=true)
	end

	md"""
	#### Trajectories for $ẋ = x(1-x)$

	$p

	"""
end

# ╔═╡ 446c564e-220f-11eb-191a-6b419e790f3f
md"""
We see that all the curves starting near to $x_0=1.0$ seem to converge to 1 at long times. If the system starts *exactly* at 0 then it stays there forever. However, if it starts close to 0, on either side, then it moves *away* from 0 (on that same side of 0) -- starting from a negative value $x$ becomes ever more negative. (Even though negative populations have no meaning in the original interpretation as the dynamics of a population, we can still ask study the dynamics of the equation with negative initial conditions, since it may model other systems too.)

The special values $x^*_1=1$ and $x^*_2=0$ are called **stationary points** or **fixed points** of the differential equation. If we start at $x^*_i$, then the derivative there is $f'(x^*_i) = 0$, and hence we cannot move away from $x^*_i$! The fixed points can be found as zeros or **roots** of the function $f$, i.e. values $x^*$ such that $f(x^*) = 0$.
"""

# ╔═╡ 38894c32-2210-11eb-26ff-d9796edd7871
md"""
We see, though, that the two types of fixed points are **qualitatively different**: trajectories that start close to $x^*_1 = 1$ move *towards* $x^*_1$, whereas trajectories that start close to $x^*_2 = 0$ move *away* from it. We say that $x^*_1$ is a **stable fixed point** and $x^*_2$ is an **unstable fixed point**.

In general it is not possible to find analytical formulas for the position and stability of fixed points; instead, we can use numerical **root-finding algorithms**, for example the Newton method. There are various Julia packages with implementations of these algorithms; we will use `Roots.jl`, which works well for scalar-valued functions.
"""

# ╔═╡ a196a20a-2216-11eb-11f9-6789bb2a3ece
md"""
## State space: Vector field and phase portrait
"""

# ╔═╡ c15a53e8-2216-11eb-1460-8510aea0f570
md"""
If we want to find the whole trajectory for a given initial condition then we need to solve the equations, either numerically or analytically.

However, we may want less information about the system, for example the **long-time** or **asymptotic** dynamics. It turns out that we can obtain some information about that *without* explicitly solving the ODE! This is the **qualitative approach** to studying nonlinear systems.

Instead of drawing trajectories $x(t)$ as a function of time $t$, as we did above, let's use a different graphical representation, where we draw **state space** or **phase space**: This is the set ("space") of all possible values of the dependent variables ("states"). For the above ODE there is only a single dependent variable, $x$, so the state space is the real line, $\mathbb{R}$.

At each possible value of $x$, the ODE gives us information about the rate of change of $x(t)$ at that point. Let's draw an **arrow** at that point, pointing in the direction that a particle placed at that point would move: to the right if $\dot{x} > 0$ and to the left if $\dot{x} < 0$.
"""

# ╔═╡ 9a833edc-2217-11eb-0701-99862b410bfa
vector_field(logistic)

# ╔═╡ 42c28b16-2218-11eb-304a-e534353fa12b
md"""
This vector field indeed gives us a *qualitative* picture of the dynamics. It does not tell us how fast the dynamics will occur in each region, but it indicates what the *tendency* is. We have coded the fixed points according to their stability; this may be calculated using the derivative evaluated at the fixed point, $f'(x^*)$, since this derivative controls the behaviour of nearby initial conditions $x^* + \delta x$.
The unstable fixed point is shown as a green square, and the stable fixed point as a grey circle; we will use this convention throughout the notebook.
"""

# ╔═╡ e6ed3366-2246-11eb-3a86-f375162c746d
md"""
## Bifurcations
"""

# ╔═╡ ec3c5d9c-2246-11eb-1ec5-e3d9c8a7fa23
md"""
Now suppose that there is a **parameter** $\mu$ in the system that can be varied. For each value of $\mu$ we have a *different* ODE

$$\dot{x} = f_\mu(x).$$

For example,
$$\dot{x} = \mu + x^2.$$

Let's draw the state space for each different value of $\mu$:
"""

# ╔═╡ 3fa21738-2247-11eb-09d5-4dcf0e7377fc
g(μ, x) = μ + x^2

# ╔═╡ 19eb987e-2248-11eb-2982-1930b5840aed
md"""
μ = $(@bind λ Slider(-1.0:0.05:1, show_value=true))
"""

# ╔═╡ 25578fa6-2248-11eb-2838-57fbfb928fc4
begin
	vector_field(x -> g(λ, x))
end

# ╔═╡ 49ec803a-22d7-11eb-3f04-db8f8c119a7d
md"""
Now let's collect all the vector fields into a single plot. The horizontal axis now represents the different possible values of the parameter $\mu$:
"""

# ╔═╡ b9ee4db4-22d7-11eb-38ef-511e30df16cc
bifurcation_diagram(g)

# ╔═╡ 27190ac0-230b-11eb-3367-af2bbbf57e5e
md"""
We see that at the **critical value** $\mu_c = 0$ there is a **qualitative change in behaviour** of the system: for $\mu_c < 0$ there are two fixed points, whereas for $\mu_c > 0$ there are no fixed points at all. Such a qualitative change is called a **bifurcation**. In this particular case the two fixed points collide in a **saddle--node** or **fold** bifurcation.
"""

# ╔═╡ 48ef665e-22d0-11eb-2c72-b722fafd296e
md"""
# 1D: Bistability and hysteresis
"""

# ╔═╡ 51eeb46c-22d0-11eb-1643-c92e725d9da0
md"""
Now let's look at the dynamics of the following system:

$$\dot{x} = \mu + x - x^3.$$
"""

# ╔═╡ 44c1738a-2247-11eb-3232-f30b46fd7ce6
h(μ, x) = μ + x - x^3

# ╔═╡ 7a94642e-39d8-480d-9516-d11e6c1b078b
md"""
Let's plot the bifurcation diagram again:
"""

# ╔═╡ f040d1c4-2247-11eb-0bf0-090e90a86a85
bifurcation_diagram(h)

# ╔═╡ a56ad0bc-22d7-11eb-0418-d36f3ef14109
md"""
We see that there is a range of values of $\mu$ for which there are *three coexisting fixed points*, two stable and one unstable. Since there are two stable fixed points in which the system can remain, we say that the system is **bistable**.
"""

# ╔═╡ 5b9c66ce-224b-11eb-02a4-b9a5a173ee6a
md"""
Now that we understand what the plots mean and the dynamics, let's plot just the fixed points $x^*(\mu)$ as a function of $\mu$. Such a plot is called a **bifurcation diagram**:
"""

# ╔═╡ 0c3c88fe-2249-11eb-03af-bf9c393fadd4
fixed_points(h)

# ╔═╡ aa81663c-230b-11eb-3c0e-39759f36deb6
md"""
The pieces of curve are called **branches**.
"""

# ╔═╡ 55d4d188-2249-11eb-0c47-dfddce9b2027
md"""
## Hysteresis
"""

# ╔═╡ 5ac33518-2249-11eb-2d57-7de185d43e5f
md"""
Suppose we now think about slowly varying the parameter $\mu$. If we change the parameter $\mu$ by a little, the system is no longer at a fixed point, since the position of the fixed point moves when $\mu$ changes. However, the system will then **relax**: it will follow the dynamics at the new value of $\mu$, and will rapidly converge to the new fixed point nearby for that new value of $\mu$.

For example, starting at $\mu=-2$, the system will stay on the lower black (stable) **branch** until $\mu=0.4$ or so. At that point, two fixed points collide and annihilate each other! After that there is no longer a fixed point nearby. However, there is another fixed point much further up that will now attract all trajectories, so the system rapidly transitions to that fixed point.

Now suppose we decrease the parameter again. The system will now track the *upper* branch until $\mu=-0.4$ or so, when again it will jump back down.

For each parameter value $\mu$ in the interval $[-0.4, 0.4]$ there is **bistability**, i.e. **coexistence** of *two* fixed points with the same value of $\mu$ (together with a third, unstable fixed point that is not observable).

The fact that the system tracks different stable branches depending on where we started, i.e. on the history of the dynamics, is known as **hysteresis**.
"""

# ╔═╡ 4c73705e-230c-11eb-3c90-b14536d78808
md"""
Hysteretic behaviour like this is found in many scientific and engineering contexts,  including switches in biology, for example genetic switches, and in the historical dynamics of the earth's climate.
"""

# ╔═╡ 2ced9e26-22d0-11eb-34a5-bf80822f5c43
md"""
## Slow--fast systems
"""

# ╔═╡ 3213d406-22d0-11eb-0e89-71b42ea5cada
md"""
What are we actually doing when we let the parameter $\mu$ vary? Effectively we now have a system of *two* coupled ODEs, for example

$$\dot{x} = \mu + x - x^3;$$
$$\dot{\mu} = \epsilon,$$

where $\mu$ varies at a slow speed $\epsilon$. On a time scale much shorter than $1 / \epsilon$, the dynamics of $x$ "does not know" that $\mu$ is changing, so it will converge to a fixed point $x^*(\mu)$ for the current value of $\mu$. [An associated term is **adiabatic approximation**.] However, in fact $\mu$ does gradually change, so the value of $x$ will effectively "slide along" the curve $x(t) \simeq x^*(\mu(t))$, tracking the curve of fixed points as $\mu$ changes.

Once $\mu$ reaches a critical value $\mu_c$, however, there is no longer a nearby fixed point to move to, and the vector field will move the system to rapidly transition to the far away alternative fixed point.

If we now reverse the dynamics of $\mu$, we slide back along the upper branch.
"""

# ╔═╡ bcce1286-2213-4f25-a947-7a429bd96594
md"""
# 2D: Oscillations in chemical reactions -- the Brusselator model
"""

# ╔═╡ 77ae220e-08bd-4808-b1f5-f80397afe9bc
md"""
Bifurcations are not restricted to 1D systems; indeed, the dynamics can be much richer in higher dimensions.

For example, let's look at the [Brusselator model](https://en.wikipedia.org/wiki/Brusselator). This models a chemical reaction that **oscillates**, known as a [chemical clock](https://en.wikipedia.org/wiki/Chemical_clock):
"""

# ╔═╡ c48b8097-d81e-43a1-a61f-b2d01edbc8d5
md"""
$$\begin{aligned}
\dot{x} &= a + x^2 y - bx - x \\
\dot{y} &= b x - x^2 y
\end{aligned}$$
"""

# ╔═╡ 6b857cb0-182f-4451-9f5a-e170dae185c1
function brusselator(xx, p, t)
	x, y = xx
	a, b = p

	return [a + x^2 * y - b*x - x,
			b * x - x^2 * y]
end

# ╔═╡ e2953d11-6880-45cb-8000-3f6903a85c7d
begin
	u0 = [1, 1]
	tspan = (0.0, 50.0)
	params = [a, b]
end

# ╔═╡ cbb9e351-b496-45cc-8816-519ada1dc5d5
begin
	prob = ODEProblem(brusselator, u0, tspan, params)
	soln = solve(prob)
end;

# ╔═╡ 77e5ff47-1e95-4552-a3d2-2fdf7bbea7cb
md"""
a = $(@bind a Slider(0.0:0.1:5.0, show_value=true, default=1.0))

b = $(@bind b Slider(0.0:0.1:5.0, show_value=true, default=1.5))
"""

# ╔═╡ 92aca04d-376f-4ac7-b4c3-5b0652974157
gr(dpi=300)

# ╔═╡ 34556546-0cab-43e2-bc2c-f7487288392c
let
	tspan = (0.0, 10.0)
	params = [a, b]

	p1 = plot(leg=false, background_color_inside=:black)

	for x in 0:1.0:5
		for y in 0:1.0:5
			u0 = [x, y]

			prob = ODEProblem(brusselator, u0, tspan, params)
			soln = solve(prob)

			plot!(p1, soln, vars=(1, 2), xlims=(0, 5), ylims=(0, 5), ratio=1, lw=1.5)
#	p2 = plot(soln)

#	plot(p1, p2, ylims=(0, 5))
		end
	end

	# plot direction field:
	xs = Float64[]
	ys = Float64[]

	for x in 0:0.1:5
		for y in 0:0.1:5
			v = brusselator([x, y], params, 0)
			v ./= (norm(v) * 30)
			# plot!([x, x + v[1]], [y, y + v[2]], alpha=0.5, c=:gray)

			push!(xs, x - v[1], x + v[1], NaN)
			push!(ys, y - v[2], y + v[2], NaN)

		end
	end

	plot!(xs, ys, alpha=0.7, c=:gray)

	# as_svg(p1)

	md"""
	#### Trajectories in state space for the Brusselator model

	$(p1)

	"""

end

# ╔═╡ c90e3885-33bd-40d2-8a64-946fd7676656
md"""
We see that there is a critical value of $b$ above which the fixed point becomes unstable and gives rise to an attracting periodic orbit. This is a [**Hopf bifurcation**](https://en.wikipedia.org/wiki/Hopf_bifurcation).
"""

# ╔═╡ e01fa79b-e92b-4b77-8234-108839ed2598
md"""
# 3D: Chaos in the Lorenz equations
"""

# ╔═╡ 1d91975b-2644-41da-a2a4-667ae6ed5ca1
md"""
The [Lorenz equations](https://en.wikipedia.org/wiki/Lorenz_system) form a (very) simplified model of convection in a layer of fluid representing the atmosphere, first investigated in a famous paper by Edward Lorenz, published in 1963, about work done at MIT.  They represent a pioneering numerical investigation of **chaotic behaviour**.

The Lorenz equations are a set of three coupled ODEs. They have an apparently simple form, with only two nonlinear terms:
"""

# ╔═╡ d99002e1-02ca-45bd-ba26-581f7f1b4fe8
md"""
$$\begin{align}
\dot{x} &= \sigma (y - x) \\[6pt]
\dot{y} &= x (\rho - z) - y \\[6pt]
\dot{z} &= x y - \beta z
\end{align}$$
"""

# ╔═╡ 64642233-3331-40ea-929b-be5daa933393
md"""
Nonetheless, we will see that they can exhibit remarkably complex behaviour. Indeed, most nonlinear ODEs with at least three variables can be expected to exhibit similarly complicated behaviour.
"""

# ╔═╡ 2228c7f6-8692-487e-bf8a-6600e7e93a06
md"""
Let's solve the system using `DifferentialEquations.jl`. We will take the classical parameter values $\sigma = 10$ and $\beta = 8/3$, but we will allow $\rho$ to vary; its classical value is $28$.
"""

# ╔═╡ cfa7b192-f1a0-4cd1-ae34-cec07e7b47f6
function lorenz(u, p, t)
	x, y, z = u
	σ, ρ, β = p

	dx = σ * (y - x)
	dy = x * (ρ - z) - y
	dz = x * y - β * z

	return [dx, dy, dz]
end

# ╔═╡ 2e61c129-7fd8-46b6-8480-40f0244aab47
begin
	lorenz_prob = ODEProblem(lorenz, [0.01, 0.01, 0.01], (0.0, 100.0), lorenz_params)

	lorenz_soln = solve(lorenz_prob)
end;

# ╔═╡ 58dbd153-09e2-4a38-94c1-9f69acf515ae
lorenz_params = (σ = 10.0, ρ = ρ, β = 8/3)

# ╔═╡ e068c809-34c6-4723-8ad5-a55918cfed87
md"""
ρ = $(@bind ρ Slider(0.0:0.1:100.0, show_value=true, default=10.0))
"""

# ╔═╡ 5ffeba9b-8487-4be1-bb64-8d2330824ee5
plot(lorenz_soln, vars=(1, 2, 3), xlabel="x", ylabel="y", zlabel="z", xlims=(-25, 25), ylims=(-25, 25), zlims=(0, 60))

# ╔═╡ a7f067e4-d088-400b-bba4-d58f9cfb87a8
md"""
As $\rho$ increases we see a sequence of bifurcations. Above a critical value, trajectories converge to a fractal **strange attractor**, on which the dynamics is **chaotic**.
"""

# ╔═╡ cc7ee950-6b86-4d43-8524-86254978bd1b
md"""
Deterministic chaos occurs when nearby initial conditions separate exponentially fast in state space. This has been given the name **butterly effect**: the perturbation to the atmosphere's state caused by a butterfly flapping its wings could end up being magnified to modify the direction in which a tornado moves.

We can see this in a simple way by perturbing the initial condition slightly and calculating the distance between the two systems as a function of $t$:
"""

# ╔═╡ aebdd083-4f70-4490-9949-8fe5c1cb2e1e
begin
	ϵ = 1e-10

	lorenz_prob2 = ODEProblem(lorenz, [0.01 + ϵ, 0.01 + ϵ, 0.01 + ϵ], (0.0, 50.0), lorenz_params)

	lorenz_soln2 = solve(lorenz_prob2)
end;

# ╔═╡ f44c9f92-aacd-4556-b494-8ee874387e5d
begin
	ts = 0:0.01:50
	distances = [norm(lorenz_soln2(t) - lorenz_soln(t)) for t in ts]
end

# ╔═╡ c2615e16-93d1-4db0-af18-0c1a758cd58a
plot(ts, distances, yscale=:log10, label="distance", xlabel="t", leg=:topleft)

# ╔═╡ 0b1b4e0e-d146-4c94-bfee-3d0a50323d0b
md"""
We indeed see a significant window over which the distance grows exponentially (straight line on the semi-logarithmic graph) for large enough values of $\rho$, before it saturates at a value of the order of the diameter of the attractor.

The rate of the exponential growth is known as the **Lyapunov exponent**; there are better ways to measure it more accurately.
"""

# ╔═╡ b8423240-1761-490d-a856-f722714cc763
md"""
Since the Lorenz equations model the atmosphere, we could expect that the dynamics of the atmosphere is at least as complicated.
"""

# ╔═╡ 8eb1c83e-18f2-497a-94fd-a4abe11beed3
md"""
We can also plot the $x$ coordinate of each trajectory for comparison:
"""

# ╔═╡ 3188c2ca-ab8f-45a9-bcbc-10159659a9d4
begin
	plot(lorenz_soln, vars=1, label="original", size=(500, 300), leg=:topleft)
	plot!(lorenz_soln2, vars=1, label="perturbed")
	ylabel!("x(t)")
end

# ╔═╡ 75f3d63d-aa5d-46ce-824c-72153623d2c5
md"""
At long enough times the trajectories separate and behave very differently.
"""

# ╔═╡ 019be011-51ed-4997-8e8e-1035c109efe7
md"""
Now let's look at the "climate" in this model, namely the average value of each coordinate over a long time:
"""

# ╔═╡ edd470f6-828e-47f4-b2d1-cb16f6103158
mean(v) = sum(v) / length(v)

# ╔═╡ 5c27ade2-5983-4332-9097-20efb20504bc
begin
	T = 1000.0

	lorenz_prob3 = ODEProblem(lorenz, [1 + ϵ, 1 + ϵ, 1 + ϵ], (0.0, T), lorenz_params)

	lorenz_soln3 = solve(lorenz_prob3)

	lorenz_prob4 = ODEProblem(lorenz, [1, 1, 1], (0.0, T), lorenz_params)

	lorenz_soln4 = solve(lorenz_prob4)
end;

# ╔═╡ 043c35e5-bf4b-445a-b746-9030b96033ca
mean(abs.(lorenz_soln3(t)) for t in T/2:T)

# ╔═╡ d9b34169-85a2-4f97-a7e3-a2ed67ee205a
mean(abs.(lorenz_soln4(t)) for t in T/2:T)

# ╔═╡ 02f295ba-6bb1-41ca-9581-af33766241a4
md"""
We see that the average mean value of each component is approximately the same, even though the individual trajectories are wildly different. This is an example of how statistical properties can be the same, even if individual behaviour is very different, and motivates the idea that climate -- i.e. "average weather" -- can be stable, even if day-to-day variations differ a lot.
"""

# ╔═╡ d0918f0c-22d4-11eb-383a-553e33114c52
md"""
# Function library
"""

# ╔═╡ 319c2ba8-220e-11eb-1ee0-4dd51ff4cd25
euler_step(f, x, h) = x + h * f(x)

# ╔═╡ 42de6dd6-220e-11eb-2a48-a3c69b10a03b
function euler(f, x0, h, t_final)
	ts = [0.0]
	xs = [x0]

	x = x0
	t = 0.0

	while t < t_final
		x = euler_step(f, x, h)
		t += h

		push!(xs, x)
		push!(ts, t)
	end

	return (ts=ts, xs=xs)  # a named tuple
end

# ╔═╡ 546168ce-2218-11eb-198b-1da9f8a9a242
derivative(f, x, h=0.001) = ( f(x + h) - f(x - h) ) / (2h)

# ╔═╡ 46f1e3a2-2217-11eb-184f-f5649b892004
"Draw 1D vector field using centred arrows"
function horiz_vector_field(f)

	xlo, xhi = -2, 2.5
	arrow_size = 0.07
	tol=1e-5  # tolerance to check for fixed point

	p = plot(size=(400, 100), xlim=(xlo, xhi), ylim=(-0.5, 0.5), leg=false, yticks=[])

	for x in xlo:0.2:xhi
		d = f(x)  # derivative

		if d > tol
			plot!([x-arrow_size, x+arrow_size], [0, 0], arrow=true, c=:blue, alpha=0.5, lw=1.5)

		elseif d < -tol
			plot!([x+arrow_size, x-arrow_size], [0, 0], arrow=true, c=:red, alpha=0.5, lw=1.5)

		end

	end

	roots = find_zeros(f, -10, 10)

	for root in roots
		stability = sign(derivative(f, root))

		if stability > 0
			scatter!([root], [0], c=:green, alpha=0.5, m=:square, ms=5)  # unstable
		else
			scatter!([root], [0], c=:black, alpha=0.5, ms=6)  # stable
		end

	end

	p
end

# ╔═╡ 2b1b49b0-2247-11eb-1d25-57f3b4f46d04
"Draw vertical vector field of 1D ODE on plot p with parameter μ"
function vector_field!(p, μ, f)

	xlo, xhi = -2, 2
	arrow_size = 0.07
	tol=1e-5  # tolerance to check for fixed point

	for x in xlo:0.2:xhi
		d = f(x)  # derivative

		if d > tol

			plot!([μ, μ], [x-arrow_size, x+arrow_size], arrow=(3.0, 2.0), arrowstyle=:triangle, c=:blue, alpha=0.4, lw=1.5)

		elseif d < -tol
			plot!([μ, μ], [x+arrow_size, x-arrow_size], arrow=true, c=:red, alpha=0.4, lw=1.5)
		end

	end

	roots = find_zeros(f, -10, 10)

	for root in roots
		stability = sign(derivative(f, root))

		if stability > 0
			scatter!([μ], [root], c=:green, alpha=0.5, m=:square, ms=2)  # unstable
		else
			scatter!([μ], [root], c=:black, alpha=0.5, ms=3)  # stable
		end

	end

	p
end

# ╔═╡ d4fa7bec-2371-11eb-3e9f-69d552546331
function vector_field(f)
	p = plot(leg=false, xticks=[], xrange=(-1, 1), size=(100, 200))
	vector_field!(p, 0, f)
end

# ╔═╡ e55d2780-2247-11eb-01e0-fbdc94bba264
function bifurcation_diagram(h)
	p = plot(leg=false, ratio=1)

	for μ in -2:0.15:2
		vector_field!(p, μ, x -> h(μ, x))
	end

	xlabel!("μ")
	ylabel!("fixed points and dynamics with given μ")

	return p
end

# ╔═╡ e27dace4-2248-11eb-2ae1-953639d8c944
function fixed_points(f)
	p = plot(leg=false, ratio=1)

	for μ in -2:0.01:2

		roots = find_zeros(x -> f(μ, x), -10, 10)

		for root in roots
			stability = sign(derivative(x -> f(μ, x), root))

			if stability > 0
				scatter!([μ], [root], c=:green, alpha=0.5, m=:square)  # unstable
			else
				scatter!([μ], [root], c=:black, alpha=0.5)  # stable
			end
		end

	end

	xlabel!(L"\mu")
	ylabel!(L"x^*(\mu)")

	return p
end

# ╔═╡ Cell order:
# ╟─1f65fc80-a83b-11eb-2583-89bff1ea8372
# ╠═88b46d2e-220e-11eb-0f7f-b3f523f0214e
# ╠═f9dbf2e8-574d-4846-af48-f7e5a82a1afc
# ╟─a2b2eae2-762e-41c3-a546-0caedc79db7d
# ╟─91b34e59-e4d2-45b7-9e2b-1b95748e5d6a
# ╟─22bddfca-fd1c-4dc7-85f3-f864c3b53407
# ╟─6f554dc0-220c-11eb-341f-1b66de841c86
# ╟─2c6f4b31-2f08-4e24-8e31-ad80d2700fa8
# ╟─822c015a-8bb3-4432-a2d4-2e2ba4f906a6
# ╟─816de40c-220c-11eb-2d3a-23e6267cd529
# ╟─d21f6358-220c-11eb-00b8-03fe0746fdc9
# ╟─735bb47e-220d-11eb-1ba6-8d591a935857
# ╟─79b148b6-220d-11eb-2785-25d05068aeed
# ╟─d43a869e-220d-11eb-009c-a119de57e7da
# ╟─088eabb2-220e-11eb-1ac0-df419b87f39a
# ╟─a44cc180-22d4-11eb-2129-211d024c4921
# ╟─57566646-1b4c-4098-abf7-eb242d01ca81
# ╠═6aa1d2f4-220e-11eb-06e0-c346f74aa018
# ╟─b6d0dd18-22d4-11eb-2083-0b02de030579
# ╠═6e90ea26-220e-11eb-0c65-bf52b3d2e195
# ╟─2d1c6abf-3b15-4638-88c5-89b5d0585c98
# ╟─f055be76-22d4-11eb-268a-bf70c7d8f1a1
# ╟─96e22792-220e-11eb-2729-63964507b5f2
# ╟─90ccb392-2216-11eb-1fd8-83b7d7c16b54
# ╠═0bf302de-4a0f-4341-b010-d3606e237acc
# ╟─bef9c2e4-220e-11eb-24d8-bd618d2985ea
# ╟─17317eb8-22d5-11eb-0f46-2bf1bddd36bb
# ╟─2300b29a-22d5-11eb-3c99-bdec0e5a2685
# ╟─ae862e10-22d5-11eb-3d75-1d748cf86944
# ╟─c5bf15b1-a540-478b-acec-1aa47aad13d1
# ╟─446c564e-220f-11eb-191a-6b419e790f3f
# ╟─38894c32-2210-11eb-26ff-d9796edd7871
# ╟─a196a20a-2216-11eb-11f9-6789bb2a3ece
# ╟─c15a53e8-2216-11eb-1460-8510aea0f570
# ╠═9a833edc-2217-11eb-0701-99862b410bfa
# ╟─42c28b16-2218-11eb-304a-e534353fa12b
# ╟─e6ed3366-2246-11eb-3a86-f375162c746d
# ╟─ec3c5d9c-2246-11eb-1ec5-e3d9c8a7fa23
# ╠═3fa21738-2247-11eb-09d5-4dcf0e7377fc
# ╟─19eb987e-2248-11eb-2982-1930b5840aed
# ╠═25578fa6-2248-11eb-2838-57fbfb928fc4
# ╟─49ec803a-22d7-11eb-3f04-db8f8c119a7d
# ╠═b9ee4db4-22d7-11eb-38ef-511e30df16cc
# ╟─27190ac0-230b-11eb-3367-af2bbbf57e5e
# ╟─48ef665e-22d0-11eb-2c72-b722fafd296e
# ╟─51eeb46c-22d0-11eb-1643-c92e725d9da0
# ╠═44c1738a-2247-11eb-3232-f30b46fd7ce6
# ╟─7a94642e-39d8-480d-9516-d11e6c1b078b
# ╠═f040d1c4-2247-11eb-0bf0-090e90a86a85
# ╟─a56ad0bc-22d7-11eb-0418-d36f3ef14109
# ╟─5b9c66ce-224b-11eb-02a4-b9a5a173ee6a
# ╠═0c3c88fe-2249-11eb-03af-bf9c393fadd4
# ╟─aa81663c-230b-11eb-3c0e-39759f36deb6
# ╟─55d4d188-2249-11eb-0c47-dfddce9b2027
# ╟─5ac33518-2249-11eb-2d57-7de185d43e5f
# ╟─4c73705e-230c-11eb-3c90-b14536d78808
# ╟─2ced9e26-22d0-11eb-34a5-bf80822f5c43
# ╟─3213d406-22d0-11eb-0e89-71b42ea5cada
# ╟─bcce1286-2213-4f25-a947-7a429bd96594
# ╟─77ae220e-08bd-4808-b1f5-f80397afe9bc
# ╟─c48b8097-d81e-43a1-a61f-b2d01edbc8d5
# ╠═6b857cb0-182f-4451-9f5a-e170dae185c1
# ╠═e2953d11-6880-45cb-8000-3f6903a85c7d
# ╠═cbb9e351-b496-45cc-8816-519ada1dc5d5
# ╟─77e5ff47-1e95-4552-a3d2-2fdf7bbea7cb
# ╠═92aca04d-376f-4ac7-b4c3-5b0652974157
# ╠═34556546-0cab-43e2-bc2c-f7487288392c
# ╟─c90e3885-33bd-40d2-8a64-946fd7676656
# ╟─e01fa79b-e92b-4b77-8234-108839ed2598
# ╟─1d91975b-2644-41da-a2a4-667ae6ed5ca1
# ╟─d99002e1-02ca-45bd-ba26-581f7f1b4fe8
# ╟─64642233-3331-40ea-929b-be5daa933393
# ╟─2228c7f6-8692-487e-bf8a-6600e7e93a06
# ╠═cfa7b192-f1a0-4cd1-ae34-cec07e7b47f6
# ╠═2e61c129-7fd8-46b6-8480-40f0244aab47
# ╠═58dbd153-09e2-4a38-94c1-9f69acf515ae
# ╟─e068c809-34c6-4723-8ad5-a55918cfed87
# ╠═5ffeba9b-8487-4be1-bb64-8d2330824ee5
# ╟─a7f067e4-d088-400b-bba4-d58f9cfb87a8
# ╟─cc7ee950-6b86-4d43-8524-86254978bd1b
# ╠═aebdd083-4f70-4490-9949-8fe5c1cb2e1e
# ╠═f44c9f92-aacd-4556-b494-8ee874387e5d
# ╠═c2615e16-93d1-4db0-af18-0c1a758cd58a
# ╟─0b1b4e0e-d146-4c94-bfee-3d0a50323d0b
# ╟─b8423240-1761-490d-a856-f722714cc763
# ╟─8eb1c83e-18f2-497a-94fd-a4abe11beed3
# ╟─3188c2ca-ab8f-45a9-bcbc-10159659a9d4
# ╟─75f3d63d-aa5d-46ce-824c-72153623d2c5
# ╟─019be011-51ed-4997-8e8e-1035c109efe7
# ╠═edd470f6-828e-47f4-b2d1-cb16f6103158
# ╠═5c27ade2-5983-4332-9097-20efb20504bc
# ╠═043c35e5-bf4b-445a-b746-9030b96033ca
# ╠═d9b34169-85a2-4f97-a7e3-a2ed67ee205a
# ╟─02f295ba-6bb1-41ca-9581-af33766241a4
# ╟─d0918f0c-22d4-11eb-383a-553e33114c52
# ╠═319c2ba8-220e-11eb-1ee0-4dd51ff4cd25
# ╠═42de6dd6-220e-11eb-2a48-a3c69b10a03b
# ╠═546168ce-2218-11eb-198b-1da9f8a9a242
# ╠═46f1e3a2-2217-11eb-184f-f5649b892004
# ╠═2b1b49b0-2247-11eb-1d25-57f3b4f46d04
# ╠═d4fa7bec-2371-11eb-3e9f-69d552546331
# ╠═e55d2780-2247-11eb-01e0-fbdc94bba264
# ╟─e27dace4-2248-11eb-2ae1-953639d8c944


text:
module nevanlinna_ac
__precompile__()

#setprecision(BigFloat, 128)
#BigFloat defalt precision 256

using LinearAlgebra
using Random; Random.seed!()
using DelimitedFiles

export eye, delta, gaussian, multi_gaussian
export Masubara_freq, Masubara_GF
export h, invh, h1, invh1

export Giωn, MasubaraGF, make_input
export pick_matrix, ispossemidef, isNevanlinnasolvable
export core, evaluation, spectrum_density

#input file format: ωn  Re[G(iωn)]  Im[G(iωn)]
include("utilities.jl")
include("interpolate.jl")

end


text:
"""
    linesfill(xs, ys, lower, upper; kwargs...)

Line plot with a shaded area between `lower` and `upper`. If `lower` and `upper`
are not given, shaded area is between `0` and `ys`.
## Attributes
$(ATTRIBUTES)
"""
@recipe(LinesFill) do scene
    l_theme = default_theme(scene, Lines)
    Attributes(
        color = l_theme.color,
        colormap = l_theme.colormap,
        colorrange = get(l_theme.attributes, :colorrange, automatic),
        linestyle = l_theme.linestyle,
        linewidth = l_theme.linewidth,
        fillalpha = 0.15,
    )
end

function AbstractPlotting.plot!(p::LinesFill)
    lines!(p, p[1:2]...;
        color = p.color,
        linestyle = p.linestyle,
        linewidth = p.linewidth,
        colormap = p.colormap,
        colorrange = p.colorrange,
    )
    if length(p) == 2
        lower, upper = lift(zero, p[2]), p[2]
    else
        lower, upper = p[3], p[4]
    end
    meshcolor = lift(to_color∘tuple, p.color, p.fillalpha)
    band!(p, p[1], lower, upper; color = meshcolor)
end


text:
"""
    AbstractDiffBackend

An abstract type for diff backends. See [`FiniteDifferencesBackend`](@ref) for
an example.
"""
abstract type AbstractDiffBackend end

struct NoneDiffBackend <: AbstractDiffBackend end

"""
    _derivative(f, t[, backend::AbstractDiffBackend])

Compute the derivative of a callable `f` at time `t` computed using the given `backend`,
an object of type [`Manifolds.AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`default_differential_backend`](@ref).

This function calculates plain Euclidean derivatives, for Riemannian differentiation see
for example [`differential`](@ref Manifolds.differential(::AbstractManifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _derivative end

_derivative(f, t) = _derivative(f, t, default_differential_backend())

function _derivative!(f, X, t, backend::AbstractDiffBackend=default_differential_backend())
    return copyto!(X, _derivative(f, t, backend))
end

"""
    _gradient(f, p[, backend::AbstractDiffBackend])

Compute the gradient of a callable `f` at point `p` computed using the given `backend`,
an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`default_differential_backend`](@ref).

This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _gradient end

_gradient(f, p) = _gradient(f, p, default_differential_backend())

function _gradient!(f, X, p, backend::AbstractDiffBackend=default_differential_backend())
    return copyto!(X, _gradient(f, p, backend))
end

"""
    _jacobian(f, p[, backend::AbstractDiffBackend])

Compute the jacobian of a callable `f` at point `p` computed using the given `backend`,
an object of type [`AbstractDiffBackend`](@ref). If the backend is not explicitly
specified, it is obtained using the function [`default_differential_backend`](@ref).

This function calculates plain Euclidean gradients, for Riemannian gradient calculation see
for example [`gradient`](@ref Manifolds.gradient(::AbstractManifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)).

!!! note

    Not specifying the backend explicitly will usually result in a type instability
    and decreased performance.
"""
function _jacobian end

_jacobian(f, p) = _jacobian(f, p, default_differential_backend())

function _jacobian!(f, X, p, backend::AbstractDiffBackend=default_differential_backend())
    return copyto!(X, _jacobian(f, p, backend))
end

"""
    CurrentDiffBackend(backend::AbstractDiffBackend)

A mutable struct for storing the current differentiation backend in a global
constant [`_current_default_differential_backend`](@ref).

# See also

[`AbstractDiffBackend`](@ref), [`default_differential_backend`](@ref), [`set_default_differential_backend!`](@ref)
"""
mutable struct CurrentDiffBackend
    backend::AbstractDiffBackend
end

"""
    _current_default_differential_backend

The instance of [`Manifolds.CurrentDiffBackend`](@ref) that stores the globally default
differentiation backend.
"""
const _current_default_differential_backend = CurrentDiffBackend(NoneDiffBackend())
"""
    default_differential_backend() -> AbstractDiffBackend

Get the default differentiation backend.
"""
default_differential_backend() = _current_default_differential_backend.backend

"""
    set_default_differential_backend!(backend::AbstractDiffBackend)

Set current backend for differentiation to `backend`.
"""
function set_default_differential_backend!(backend::AbstractDiffBackend)
    _current_default_differential_backend.backend = backend
    return backend
end

@doc raw"""
    ODEExponentialRetraction{T<:AbstractRetractionMethod, B<:AbstractBasis} <: AbstractRetractionMethod

Approximate the exponential map on the manifold by evaluating the ODE descripting the geodesic at 1,
assuming the default connection of the given manifold by solving the ordinary differential
equation

```math
\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0,
```

where ``Γ^k_{ij}`` are the Christoffel symbols of the second kind, and
the Einstein summation convention is assumed.

See [`solve_exp_ode`](@ref) for further details.

# Constructor

    ODEExponentialRetraction(
        r::AbstractRetractionMethod,
        b::AbstractBasis=DefaultOrthogonalBasis(),
    )

Generate the retraction with a retraction to use internally (for some approaches)
and a basis for the tangent space(s).
"""
struct ODEExponentialRetraction{T<:AbstractRetractionMethod,B<:AbstractBasis} <:
       AbstractRetractionMethod
    retraction::T
    basis::B
end
function ODEExponentialRetraction(r::T) where {T<:AbstractRetractionMethod}
    return ODEExponentialRetraction(r, DefaultOrthonormalBasis())
end
function ODEExponentialRetraction(::T, b::CachedBasis) where {T<:AbstractRetractionMethod}
    return throw(
        DomainError(
            b,
            "Cached Bases are currently not supported, since the basis has to be implemented in a surrounding of the start point as well.",
        ),
    )
end
function ODEExponentialRetraction(r::ExponentialRetraction, ::AbstractBasis)
    return throw(
        DomainError(
            r,
            "You can not use the exponential map as an inner method to solve the ode for the exponential map.",
        ),
    )
end
function ODEExponentialRetraction(r::ExponentialRetraction, ::CachedBasis)
    return throw(
        DomainError(
            r,
            "Neither the exponential map nor a Cached Basis can be used with this retraction type.",
        ),
    )
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 61530684-3bb5-11ec-01b7-7bbc2d8f97f1
md"""
# Welcome to Genie

## The Highly Productive Web Framework for Julia

Genie is a full stack web framework for the Julia programming language. Genie's goals are: unparalleled developer productivity, excellent run-time performance, and security by default.

The Genie web framework follows in the footsteps of mainstream full stack web frameworks like Ruby on Rails and Django, while staying 100% true to its Julia roots. Genie's architecture and development is driven by features present in other frameworks, but not by their design. Genie takes a no-magic-no-nonsense approach by doing things the Julia way: `Controllers` are plain Julia modules, `Models` leverage types and multiple dispatch, Genie apps are nothing but Julia projects, and versioning and dependency management is provided by Julia's own `Pkg`.

Genie also takes inspiration from Julia's "start simple, grow as needed" philosophy, by allowing developers to bootstrap an app in the REPL or in a Jupyter notebook, or easily create web services and APIs with just a few lines of code.

As the projects grow more complex, Genie allows adding progressively more structure, by exposing a micro-framework which offers features like powerful routing, flexible logging, support for environments, code reload, etc.

If database persistence is needed, support for Genie's ORM, SearchLight, can be added at any time. Finally, full MVC structure can be used in order to drive the development and maintenance of complex, end-to-end, web applications.
"""

# ╔═╡ Cell order:
# ╟─61530684-3bb5-11ec-01b7-7bbc2d8f97f1


text:
using Base.Test
import WaterData

@testset "No config tests necessary" begin
    nothing
end


text:
struct Primals{T<:Real,N,M}
    Z::Vector{T}
    xinds::Vector{SVector{N,Int}}
    uinds::Vector{SVector{M,Int}}
    equal::Bool
end

function Primals(n::Int, m::Int, N::Int, equal=false)
    NN = n*N + m*(N-1) + equal*m
    Z = zeros(NN)
    uN = N-1 + equal

    xinds = [SVector{n}((n+m)*(k-1) .+ (1:n)) for k = 1:N]
    uinds = [SVector{m}(n + (n+m)*(k-1) .+ (1:m)) for k = 1:N]
    Primals(Z,xinds,uinds,equal)
end

function Base.copy(P::Primals)
    Primals(copy(P.Z),P.xinds,P.uinds,P.equal)
end

function Base.copyto!(P::Primals, Z::Traj)
    uN = P.equal ? length(Z) : length(Z)-1
    for k in 1:uN
        inds = [P.xinds[k]; P.uinds[k]]
        P.Z[inds] = Z[k].z
    end
    if !P.equal
        P.Z[P.xinds[end]] = state(Z[end])
    end
    return nothing
end

function Base.copyto!(V::AbstractVector{<:Real}, Z::Traj,
        xinds::Vector{<:AbstractVector}, uinds::Vector{<:AbstractVector})
    n,m,N = traj_size(Z)
    equal = (n+m)*N == length(V)

    uN = equal ? N : N-1
    for k in 1:uN
        inds = [xinds[k]; uinds[k]]
        V[inds] = Z[k].z
    end
    if !equal
        V[xinds[end]] = state(Z[end])
    end
    return nothing
end

function Base.copyto!(Z::Traj, P::Primals)
    uN = P.equal ? length(Z) : length(Z)-1
    for k in 1:uN
        inds = [P.xinds[k]; P.uinds[k]]
        Z[k].z = P.Z[inds]
    end
    if !P.equal
        xN = P.Z[P.xinds[end]]
        Z[end].z = [xN; control(Z[end])]
    end
    return nothing
end

function Base.copyto!(Z::Traj, V::Vector{<:Real},
        xinds::Vector{<:AbstractVector}, uinds::Vector{<:AbstractVector})
    n,m,N = traj_size(Z)
    equal = (n+m)*N == length(V)

    uN = equal ? N : N-1
    for k in 1:uN
        inds = [xinds[k]; uinds[k]]
        Z[k].z = V[inds]
    end
    if !equal
        xN = V[xinds[end]]
        Z[end].z = [xN; control(Z[end])]
    end
    return nothing
end


text:
# SPDX-License-Identifier: X11
# 2020-11-21
# management (300pt)

function main()
  N = parse(Int, readline())
  A = parse.(Int, (split ∘ readline)())
  subs = zeros(Int, N)

  for x ∈ A
    subs[x] += 1
  end

  for x ∈ subs
    println(x)
  end
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module FlexUI

import Blink
import ..JuliaJSBridge
import ..UIVariables
import ..UIControls
using ..Twinkle
import WebIO
import UUIDs

export  App,
        win,
        jjs,
        variables, variables!,
        controls, controls!,
        viewerUrl, viewerUrl!,
        renderFunction, renderFunction!,
        addVariable!,
        run

mutable struct App <: Twinkle.AbstractUIApp
    _win::Union{Nothing, Blink.Window}
    _jjs::Union{Nothing, JuliaJSBridge.JuliaJS}

    _props::Dict{Symbol, Any}
    _variables::Vector{UIVariables.AbstractUIVariable}
    _controls::Vector{UIControls.AbstractUIControl}

    _render_function::Union{Nothing, Function}

    __insideRender::Bool
end

function App(;
    variables::Vector = [],
    controls::Vector = [],
    render_function::Union{Nothing, Function} = nothing
)
    app = App(
        nothing,
        nothing,
        Dict{Symbol, Any}(),
        Vector{UIVariables.AbstractUIVariable}[],
        Vector{UIControls.AbstractUIControl}[],
        nothing,
        false
    )

    variables!(app, variables)
    controls!(app, controls)
    # viewerUrl!(app, viewer_url)
    renderFunction!(app, render_function)

    return app
end

Twinkle.win(app::App) = app._win
function win!(app::App, win::Blink.Window)
    app._win = win
end

jjs(app::App) = app._jjs
function jjs!(app::App, jjs::JuliaJSBridge.JuliaJS)
    app._jjs = jjs
end

props(app::App) = app._props
function Twinkle.prop(app::App, p::Symbol, default_val = nothing)
    return get(props(app), p, default_val)
end
function Twinkle.prop!(app::App, p::Symbol, val::Any)
    props(app)[p] = val
end

variables(app::App) = app._variables
function variables!(app::App, variables::Vector)
    app._variables = convert(Vector{UIVariables.AbstractUIVariable}, variables)
end

controls(app::App) = app._controls
function controls!(app::App, controls::Vector)
    app._controls = convert(Vector{UIControls.AbstractUIControl}, controls)
    fix_controls!(app, app._controls)
end
function controls!(app::App, control::AbstractUIControl)
    app._controls = [control]
    fix_controls!(app, app._controls)
end

# viewerUrl(app::App) = app._viewer_url
# function viewerUrl!(app::App, url::String)
#     app._viewer_url = url
# end

renderFunction(app::App) = app._render_function
function renderFunction!(app::App, render_func::Union{Nothing, Function})
    app._render_function = render_func
end

function findVariable(app::App, var_name::String)
    return findfirst(x->x.name == var_name, variables(app))
end

function var(app::App, var_name::String)
    index = findVariable(app, var_name)
    if (index === nothing)
        error("can't find variable [$var_name]");
    end
    return variables(app)[index]
end

function addVariable!(app::App, var::AbstractUIVariable)
    if (findVariable(app, var.name) !== nothing)
        error("Application already contains a variable [$(var.name)]")
    end
    var._app = app
    push!(variables(app), var)
    return var
end

function prepareApp!(app::App)
    # set the app on all variables and controls
    for var in variables(app)
        var._app = app
    end

    for c in controls(app)
        c._app = app
    end

end

function fix_controls!(app::App, controls::Vector{UIControls.AbstractUIControl})
    for control in controls
        # generate dummy variable for buttons is nececery
        # @info control.type
        if (control.type == "button" && control.variable === nothing)
            unique_var_name = isempty(control._variable_name_to_create) ? "auto_generate_for_button_$(UUIDs.uuid1())" : control._variable_name_to_create
            @info "Adding automatic variable for Button [$unique_var_name]"
            auto_var = UIVariables.Variable(
                name=unique_var_name,
                type="number",
                value=0,
            )
            push!(variables(app), auto_var)
            control.variable = unique_var_name
        end

        # TODO: Should replace these statements with a scan of all fields and run the fix on each which is a UIControls.AbstractUIControl
        if (hasproperty(control, :children))
            fix_controls!(app, control.children)
        end
        if (hasproperty(control, :content))
            if (control.content isa UIControls.AbstractUIControl)
                fix_controls!(app, Vector{UIControls.AbstractUIControl}([control.content]))
            end
        end
        if (hasproperty(control, :tabs))
            fix_controls!(app, Vector{UIControls.AbstractUIControl}(control.tabs))
        end
    end
end

function prepareJSStructure(app::App)
    res = Dict{Symbol, Any}(
        :app => Dict{Symbol, Any}(),
        :variables => Dict{Symbol, Any}(),
        :ui => [],
    )

    res[:app][:title] = prop(app, :title, "?? app title ??")
    @info app._props
    @info prop(app, :title, "?? app title ??")

    # add variables
    for var in variables(app)
        res[:variables][Symbol(var.name)] = UIVariables.typedict(var)
    end
    # @show variables(app)
    # @show res[:variables]

    res[:ui] = [UIControls.typedict(c) for c in controls(app)]

    return res;
end

function getUIFolder()

    res = joinpath(dirname(@__FILE__), "..", "..", "deps", "FlexUI");
    if (isdir(res))
        return res
    end

    res = joinpath(dirname(@__FILE__), "..", "Data", "dist");
    if (isdir(res))
        return res
    end

    res = raw"D:\Projects\Rays\Github\FlexUI\dist\FlexUI"

    return res;
end

function Base.run(app::App)
    prepareApp!(app)

    window_defaults = Blink.@d(
        :title => prop(app, :title, "?? app title ??"),
        :width => prop(app, :winInitWidth, 1600),
        :height => prop(app, :winInitHeight, 1200),
        # this will allow us to load local file which is a security risk
        :webPreferences => Blink.@d(
            :webSecurity => false,
            :experimentalFeatures => true,
        ),
    )
    win =Blink. Window(window_defaults)

    win!(app, win)

    bridge = JuliaJSBridge.JuliaJS(win; update_func=onBlinkUpdate, update_func_tag=app)
    jjs!(app, bridge)

    app_dir = getUIFolder() # raw"D:\Projects\Rays\Tests\UI\FlexUI\dist\FlexUI"
    app_node = JuliaJSBridge.application_node(jjs(app), app_dir)

    ui = WebIO.Node(:dom,
        app_node,
        WebIO.Node(Symbol("app-root"), ""),
    )

    # show the blink window
    Blink.body!(win, ui, async=false)
    # Blink.AtomShell.opentools(win)

    # set the controls through JavaScript and also the viewer url
    # js_controls = [UIControls.typedict(c) for c in controls(app)]
    data = prepareJSStructure(app)
    # viewer_url = viewerUrl(app)
    Blink.@js_ win begin
        # JuliaJS.SetHTMLFromJulia($(guid), $(html))
        # console.log("Controls", $(js_controls))
        window.fireAngularEvent("initializeApplication", [$(data)])
        # window.fireAngularEvent("setViewerUrl", [$(viewer_url)])
    end

    # call the general update function once after initialization
    render!(app)

end

function onBlinkUpdate(args::Dict{Any, Any}, app::App)
    # @info args
    meta_data = args["meta_data"]
    data = args["data"]

    if get(meta_data, "source", "") == "variable-update"
        var_name = get(data, "name", "")
        index = findfirst(x -> (x.name == var_name), variables(app))
        if (index == nothing)
            @error "Failed to locate a variable [$var_name]"
        else
            var = variables(app)[index]
            val = get(data, "value", "?")
            var.value = val;
        end

        # call variable's onChange function
        if (var._on_change_func !== nothing)
            var._on_change_func(val)
        end

        render!(app)

    elseif get(meta_data, "source", "") == "dev-tools-update"
        state = get(data, "state", "")
        if (state)
            Blink.AtomShell.opentools(Twinkle.win(app))
        else
            Blink.AtomShell.closetools(Twinkle.win(app))
        end

    elseif get(meta_data, "source", "") == "set-zoom-level"
        level = get(data, "level", "")
        if (endswith(level, '%'))
            level = level[1:end-1]
        end
        level = parse(Float64, level) / 100.0;

        w = Twinkle.win(app)
        id = w.id
        @info "zoom level", level
        command = """windows["$(id)"].webContents.setZoomFactor($level);"""
        Blink.js(w.shell, Blink.JSString(command), callback=false)
        @info "zoom level done"

    end

    # id = data["id"]
    # index = findfirst(x -> x.id == id, controls(app))
    # c = controls(app)[index]

    # # get the new value
    # val = nothing
    # if (data["type"] == "slider")
    #     val = get(data, "value", 0.0)
    #     c.value = val
    # elseif (data["type"] == "button")
    #     val = get(data, "text", "")
    # elseif (data["type"] == "field")
    #     val = get(data, "value", "")
    #     c.value = string(val)
    # else
    #     @error "unrecognized control type"
    # end

    # if (c._func !== nothing)
    #     c._func(val)
    # end

    # if (renderFunction(app) !== nothing)
    #     renderFunction(app)()
    # end

end

function Twinkle.updateVariable!(app::App, var::UIVariables.AbstractUIVariable)
    if (app !== var._app)
        # need to add variable as new in this app - not supported yet
        @error "unsupported behaviour - please check"
    end

    # send update to julia
    data = UIVariables.typedict(var)
    Blink.@js_ Twinkle.win(app) begin
        window.fireAngularEvent("setVariableFromJulia", [$(data)])
    end
    render!(app)
end

function render!(app::App)
    # guard agaist stack overflow when updating variables inside the render function
    if (app.__insideRender)
        return
    end
    app.__insideRender = true;
    try
    # call the general update function
    if (renderFunction(app) !== nothing)
            renderFunction(app)()
        end
    catch e
        @error "Encountered some errror [$e]" exception=(e, catch_backtrace())
    end
    app.__insideRender = false;
end

end # modeule FlexUI


text:
using Getopt
using Test

@testset "Getopt" begin
	args = ["-xxy1", "-y", "22", "--foo", "arg1", "-", "--wrong", "--bar=2", "--", "-x"]
	count, missing, n_x, l_y = 0, 0, 0, 0
	for (opt, arg) in Getopt.getopt(args, "xy:", ["foo", "bar="])
		if opt != "?" count += 1
		else missing += 1
		end
		if opt == "-x" n_x += 1
		elseif opt == "-y" l_y += length(arg)
		end
	end
	@test count == 6
	@test length(args) == 3
	@test missing == 1
	@test n_x == 2
	@test l_y == 3
end