text:
"""
Analysis data of Monte Carlo simulation
"""
@with_kw mutable struct MCAnalysis
    acc_rate::Float64 = 0.
    prop_local::Int = 0
    acc_local::Int = 0
    acc_rate_global::Float64 = 0.
    prop_global::Int = 0
    acc_global::Int = 0
end

"""
Parameters of Monte Carlo
"""
@with_kw struct MCParameters
    global_moves::Bool = false
    global_rate::Int = 5
    thermalization::Int = 0 # number of thermalization sweeps
    sweeps::Int = 1000 # number of sweeps (after thermalization)
    measure_rate::Int = 1 # measure at every nth sweep (after thermalization)
    print_rate::Int = 1000

    beta::Float64
end

"""
Monte Carlo simulation
"""
mutable struct MC{M<:Model, C} <: MonteCarloFlavor
    model::M
    conf::C

    thermalization_measurements::Dict{Symbol, AbstractMeasurement}
    measurements::Dict{Symbol, AbstractMeasurement}
    p::MCParameters
    a::MCAnalysis

    MC{M,C}() where {M,C} = new()
end

"""
    MC(m::M; kwargs...) where M<:Model

Create a Monte Carlo simulation for model `m` with keyword parameters `kwargs`.
"""
function MC(m::M;
        seed::Int=-1,
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = :default,
        kwargs...
    ) where M<:Model

    conf = rand(MC, m)
    mc = MC{M, typeof(conf)}()
    mc.model = m
    kwdict = Dict(kwargs)
    if :T in keys(kwargs)
        kwdict[:beta] = 1/kwargs[:T]
        delete!(kwdict, :T)
    end
    mc.p = MCParameters(; kwdict...)

    init!(
        mc, seed = seed, conf = conf,
        thermalization_measurements = thermalization_measurements,
        measurements = measurements
    )
    return mc
end

"""
    MC(m::M, params::Dict)
    MC(m::M, params::NamedTuple)

Create a Monte Carlo simulation for model `m` with (keyword) parameters
as specified in the dictionary/named tuple `params`.
"""
MC(m::Model, params::Dict{Symbol, T}) where T = MC(m; params...)
MC(m::Model, params::NamedTuple) = MC(m; params...)

# convenience
@inline beta(mc::MC) = mc.p.beta
@inline model(mc::MC) = mc.model
@inline conf(mc::MC) = mc.conf

# cosmetics
Base.summary(mc::MC) = "MC simulation of $(summary(mc.model))"
function Base.show(io::IO, mc::MC)
    print(io, "Monte Carlo simulation\n")
    print(io, "Model: ", mc.model, "\n")
    print(io, "Beta: ", round(beta(mc), sigdigits=3), " (T ≈ $(round(1/beta(mc), sigdigits=3)))\n")
    N_th_meas = length(mc.thermalization_measurements)
    N_me_meas = length(mc.measurements)
    print(io, "Measurements: ", N_th_meas + N_me_meas, " ($N_th_meas + $N_me_meas)")
end
Base.show(io::IO, m::MIME"text/plain", mc::MC) = print(io, mc)

# implement MonteCarloFlavor interface
"""
    init!(mc::MC[; seed::Real=-1])

Initialize the Monte Carlo simulation `mc`.
If `seed !=- 1` the random generator will be initialized with `Random.seed!(seed)`.
"""
function init!(mc::MC;
        seed::Real=-1,
        conf=rand(MC, mc.model),
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = :default
    )
    seed == -1 || Random.seed!(seed)
    mc.conf = conf
    mc.a = MCAnalysis()

    mc.thermalization_measurements = thermalization_measurements
    if measurements isa Dict{Symbol, AbstractMeasurement}
        mc.measurements = measurements
    elseif measurements == :default
        mc.measurements = default_measurements(mc, mc.model)
    else
        @warn(
            "`measurements` should be of type Dict{Symbol, AbstractMeasurement}, but are " *
            "$(typeof(measurements)). No measurements have been set."
        )
        mc.measurements = Dict{Symbol, AbstractMeasurement}()
    end

    nothing
end

"""
    run!(mc::MC[; verbose::Bool=true, sweeps::Int, thermalization::Int])

Runs the given Monte Carlo simulation `mc`.
Progress will be printed to `stdout` if `verbose=true` (default).
"""
function run!(mc::MC; verbose::Bool=true, sweeps::Int=mc.p.sweeps,
        thermalization=mc.p.thermalization)

    do_th_measurements = !isempty(mc.thermalization_measurements)
    do_me_measurements = !isempty(mc.measurements)
    !do_me_measurements && @warn(
        "There are no measurements set up for this simulation!"
    )
    total_sweeps = sweeps + thermalization

    start_time = now()
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    _time = time()
    do_th_measurements && prepare!(mc.thermalization_measurements, mc, mc.model)
    for i in 1:total_sweeps
        sweep(mc)

        if mc.p.global_moves && mod(i, mc.p.global_rate) == 0
            mc.a.prop_global += 1
            mc.a.acc_global += global_move(mc, mc.model, conf(mc))
        end

        # For optimal performance whatever is most likely to fail should be
        # checked first.
        if i <= thermalization && iszero(mod(i, mc.p.measure_rate)) && do_th_measurements
            measure!(mc.thermalization_measurements, mc, mc.model, i)
        end
        if (i == thermalization+1)
            do_th_measurements && finish!(mc.thermalization_measurements, mc, mc.model)
            do_me_measurements && prepare!(mc.measurements, mc, mc.model)
        end
        if i > thermalization && iszero(mod(i, mc.p.measure_rate)) && do_me_measurements
            measure!(mc.measurements, mc, mc.model, i)
        end

        print_rate = mc.p.print_rate
        if print_rate != 0 && iszero(mod(i, print_rate))
            mc.a.acc_rate /= print_rate
            mc.a.acc_rate_global /= print_rate / mc.p.global_rate
            sweep_dur = (time() - _time)/print_rate
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
                @printf("\t\tacc rate (local) : %.1f%%\n", mc.a.acc_rate*100)
                if mc.p.global_moves
                    @printf("\t\tacc rate (global): %.1f%%\n", mc.a.acc_rate_global*100)
                    @printf("\t\tacc rate (global, overall): %.1f%%\n", mc.a.acc_global/mc.a.prop_global*100)
                end
            end

            mc.a.acc_rate = 0.0
            mc.a.acc_rate_global = 0.0
            flush(stdout)
            _time = time()
        end
    end
    do_me_measurements && finish!(mc.measurements, mc, mc.model)

    mc.a.acc_rate = mc.a.acc_local / mc.a.prop_local
    mc.a.acc_rate_global = mc.a.acc_global / mc.a.prop_global

    end_time = now()
    verbose && println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
    verbose && @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)

    nothing
end

"""
    sweep(mc::MC)

Performs a sweep of local moves.
"""
@inline function sweep(mc::MC)
    c = conf(mc)
    m = model(mc)
    β = beta(mc)
    @inbounds for i in eachindex(c)
        ΔE, Δsite = propose_local(mc, m, i, c)
        mc.a.prop_local += 1
        # Metropolis
        if ΔE <= 0 || rand() < exp(- β*ΔE)
            accept_local!(mc, m, i, c, Δsite, ΔE)
            mc.a.acc_rate += 1/nsites(m)
            mc.a.acc_local += 1
        end
    end

    nothing
end

include("MC_mandatory.jl")
include("MC_optional.jl")


text:
# This file was generated, do not modify it. # hide
train, test = partition(eachindex(yc), 0.8, shuffle=true, rng=111)
Xtrain = selectrows(Xc, train)
Xtest = selectrows(Xc, test)
ytrain = selectrows(yc, train)
ytest = selectrows(yc, test);

ext: jl


text:
using Dierckx
using Plots
using RollingFunctions
using LinearAlgebra
using Dates

include("constants.jl")
include("grid.jl")
include("parameters.jl")
include("time.jl")

include("delta.jl")
include("heaviside.jl")
include("angle_of_incidence.jl")
include("bulk.jl")
include("holtslag.jl")
include("large.jl")
include("orad.jl")
include("initialize.jl")
include("ForwardStep.jl")

mtime=stime # model clock
Q=zeros(n)
sw=zeros(n)
lw=zeros(n)
lh=zeros(n)
sh=zeros(n)
sst=zeros(n)
for i=2:n+1

    global mtime
    println(i,"   ",mtime)
    ForwardStep(mtime)

    if rem(i,360)==0;#3600*6/dt
      title = plot(title = mtime, grid = false, axis = false,
       showaxis = false, bottom_margin = -50Plots.px,yaxis=nothing)
      f1=plot(-ZOC,TO,ylabel=("TO"));#,ylims=(19,21)
      f2=plot(-ZOC,UO,ylabel=("UO"));
      f3=plot(-ZOC,KOm,ylabel=("KOm"));
      f4=plot(ZAC,UA,ylabel=("UA"));
      f5=plot(ZAC,Θv.-273.15,ylabel=("Θv")); #ylim([19 21]);...
      f6=plot(ZAC,qA,ylabel=("q"));
      f7=plot(ZAC,KAm,ylabel=("KAm"));
      plt=plot(title,f1,f2,f3,f4,f5,f6,f7,layout=(8,1),legend=false,titlefontsize=6,
        size = (500, 1000))
      display(plt)
    end
    mtime=mtime+Second(ΔT)
    Q[i-1]=Qnet
    sw[i-1]=sum(SW)
    lw[i-1]=LW
    lh[i-1]=LH
    sh[i-1]=SH
    sst[i-1]=SST
end


text:
"""
Implementations of pure, extended and unscented Kalman filters.

Algorithms based on Ch.5 of Bayesian Filtering & Smoothing (Särkka, 2014)

Wouter Kouw
04-06-2020
"""

using Distributions
using Random

include("util.jl")

function kalman_filter(observations,
                       transition_matrix,
                       emission_matrix,
                       process_noise,
                       measurement_noise,
                       state0)
    """
    Kalman filter (Th. 4.2)

    This filter is built for a linear Gaussian dynamical system with known
    transition coefficients, process and measurement noise.
    """

    # Dimensionality
    Dx = size(process_noise,1)
    Dy = size(measurement_noise,1)

    # Recast process noise to matrix
    if Dx == 1
        if typeof(process_noise) != Array{Float64,2}
            process_noise = reshape([process_noise], 1, 1)
        end
        if typeof(measurement_noise) != Array{Float64,2}
            measurement_noise = reshape([measurement_noise], 1, 1)
        end
    end

    # Time horizon
    time_horizon = length(observations)

    # Initialize estimate arrays
    mx = zeros(Dx, time_horizon)
    Px = zeros(Dx, Dx, time_horizon)

    # Initial state prior
    m_0, P_0 = state0

    # Start previous state variable
    m_tmin = m_0
    P_tmin = P_0

    for t = 1:time_horizon

        # Prediction step
        m_t_pred = transition_matrix*m_tmin
        P_t_pred = transition_matrix*P_tmin*transition_matrix' .+ process_noise

        # Update step
        v_t = observations[:,t] .- emission_matrix*m_t_pred
        S_t = emission_matrix*P_t_pred*emission_matrix' .+ measurement_noise
        K_t = P_t_pred*emission_matrix'*inv(S_t)
        m_t = m_t_pred .+ K_t*v_t
        P_t = P_t_pred .- K_t*S_t*K_t'

        # Store estimates
        mx[:,t] = m_t
        Px[:,:,t] = P_t

        # Update previous state variable
        m_tmin = m_t
        P_tmin = P_t
    end
    return mx, Px
end

function extended_kalman_filter(observations::Array{Float64,2},
                                transition_function::Function,
                                emission_function::Function,
                                process_noise::Array{Float64,2},
                                measurement_noise::Array{Float64,2},
                                state0::Tuple{Array{Float64,1}, Array{Float64,2}};
                                second_order=false)
    """
    Extended Kalman filter for additive noise case
        First-order Taylor approximation = Alg 5.4 (Särkkä - BF&S)
        Second-order Taylor approxumation = Alg 5.6 (Särkkä - BF&S)

    This filter is built for a linear Gaussian dynamical system with known
    transition coefficients, process and measurement noise.
    """

    # Dimensionality
    Dx = size(process_noise,1)
    Dy = size(measurement_noise,1)

    # Recast process noise to matrix
    if Dx == 1
        if typeof(process_noise) != Array{Float64,2}
            process_noise = reshape([process_noise], 1, 1)
        end
        if typeof(measurement_noise) != Array{Float64,2}
            measurement_noise = reshape([measurement_noise], 1, 1)
        end
    end

    # Time horizon
    time_horizon = length(observations)

    # Initialize estimate arrays
    mx = zeros(Dx, time_horizon)
    Px = zeros(Dx, Dx, time_horizon)

    # Initial state prior
    m_0, P_0 = state0

    # Start previous state variable
    m_tmin = m_0
    P_tmin = P_0

    for t = 1:time_horizon

        "Prediction step"

        # Compute Jacobian evaluated at previous mean
        Fx = Jacobian(transition_function, m_tmin)

        if second_order
            # Compute per-element Hessian evaluated at previous mean
            Fxx = Hessian(transition_function, m_tmin)
        end

        # Compute predicted mean
        m_t_pred = transition_function(m_tmin)

        # Add second-order term from Taylor approx
        if second_order
            nn = length(m_t_pred)
            e = zeros(nn,)
            for i = 1:nn
                e[i] = tr(Fxx[i,:,:]*P_tmin)
            end
            m_t_pred += 1/2. *e
        end

        # Compute predicted covariance
        P_t_pred = Fx*P_tmin*Fx' + process_noise

        # Add second-order term from Taylor approx
        if second_order
            nn = length(m_t_pred)
            E = zeros(nn,nn)
            for i = 1:nn
                for j = 1:nn
                    E[i,j] = tr(Fxx[i,:,:]*P_tmin*Fxx[j,:,:]*P_tmin)
                end
            end
            P_t_pred += 1/2. *E
        end

        "Update step"

        # Compute Jacobian evaluated at predicted mean
        Hx = Jacobian(emission_function, m_t_pred)

        if second_order
            # Compute per-element Hessian evaluated at predicted mean
            Hxx = Hessian(emission_function, m_t_pred)
        end

        # Update equations
        v_t = observations[:,t] .- emission_function(m_t_pred)
        S_t = Hx*P_t_pred*Hx' .+ measurement_noise

        if second_order
            nn = length(m_t_pred)
            e = zeros(nn,)
            for i = 1:nn
                e[i] = tr(Hxx[i,:,:]*P_t_pred)
            end
            v_t += -1/2. *e

            E = zeros(nn,nn)
            for i = 1:nn
                for j = 1:nn
                    E[i,j] = tr(Hxx[i,:,:]*P_t_pred*Hxx[j,:,:]*P_t_pred)
                end
            end
            S_t += -1/2. *E
        end

        K_t = P_t_pred*Hx'*inv(S_t)
        m_t = m_t_pred .+ K_t*v_t
        P_t = P_t_pred .- K_t*S_t*K_t'

        # Store estimates
        mx[:,t] = m_t
        Px[:,:,t] = P_t

        # Update previous state variable
        m_tmin = m_t
        P_tmin = P_t
    end
    return mx, Px
end

function unscented_kalman_filter(observations::Array{Float64,2},
                                 transition_function::Function,
                                 emission_function::Function,
                                 process_noise::Array{Float64,2},
                                 measurement_noise::Array{Float64,2},
                                 state0::Tuple{Array{Float64,1}, Array{Float64,2}};
                                 α=1.1, β=2., κ=.1)
    """
    Unscented Kalman filter with additive noise (Alg. 5.14)

    This filter is built for a linear Gaussian dynamical system with unknown
    transition function, but with known process and measurement noise.
    """

    # Dimensionality
    Dx = size(process_noise,1)
    Dy = size(measurement_noise,1)

    # Recast process noise to matrix
    if Dx == 1
        if typeof(process_noise) != Array{Float64,2}
            process_noise = reshape([process_noise], 1, 1)
        end
        if typeof(measurement_noise) != Array{Float64,2}
            measurement_noise = reshape([measurement_noise], 1, 1)
        end
    end

    # Time horizon
    time_horizon = length(observations)

    # Set algorithm parameters
    λ = α^2*(Dx + κ) - Dx

    if λ < 0
        error("Unsuitable algorithm parameters chosen; λ is negative.")
    end

    # Compute constant weights
    weights_m = [[λ/(Dx + λ)]; repeat( [1/(2*(Dx + λ))], 2*Dx)]
    weights_c = [[λ/(Dx + λ) + (1-α^2+β)]; repeat( [1/(2*(Dx + λ))], 2*Dx)]

    if sum(weights_m) != 1.0
        error("Unsuitable algorithm parameters chosen; first-order weights do not sum to 1.")
    end

    # Initialize estimate arrays
    mx = zeros(Dx, time_horizon)
    Px = zeros(Dx, Dx, time_horizon)

    # Initial state prior
    m_0, P_0 = state0

    # Start previous state variable
    m_tmin = m_0
    P_tmin = P_0

    for t = 1:time_horizon

        "Prediction step"

        # Square root of previous covariance matrix
        sP_tmin = sqrt(P_tmin)

        # Initialize sigma point array
        ξ = zeros(Dx, 2*Dx+1)

        # Center point
        ξ[:, 1] = m_tmin

        # Positive and negative spread points
        for i = 1:Dx
            ξ[:, 1+i] = m_tmin + sqrt(Dx + λ)*sP_tmin[:,i]
            ξ[:, 1+Dx+i] = m_tmin - sqrt(Dx + λ)*sP_tmin[:,i]
        end

        # Propagate sigma points through transition function
        ζ = transition_function.(ξ)

        # Compute predicted mean
        m_t_pred = zeros(Dx,)
        for i=1:2Dx+1
            m_t_pred += weights_m[i] .* ζ[:,i]
        end

        # Compute predicted covariance
        P_t_pred = zeros(Dx,Dx)
        for i = 1:2*Dx+1
            P_t_pred += weights_c[i]*(ζ[:,i] .- m_t_pred)*(ζ[:,i] .- m_t_pred)' .+ process_noise
        end

        "Update step"

        # Square root of previous covariance matrix
        sP_t_pred = sqrt(P_t_pred)

        # Initialize sigma point array
        ξ_ = zeros(Dx, 2*Dx+1)

        # Center point
        ξ_[:, 1] = m_t_pred

        # Positive and negative spread points
        for i = 1:Dx
            ξ_[:, 1+i] = m_t_pred + sqrt(Dx + λ)*sP_t_pred[:,i]
            ξ_[:, 1+Dx+i] = m_t_pred - sqrt(Dx + λ)*sP_t_pred[:,i]
        end

        # Propagate sigma points through emission function
        ζ_ = emission_function.(ξ_)

        # Compute mean of inverted likelihood
        μ_t = zeros(Dx,)
        for i = 1:2*Dx+1
            μ_t += weights_m[i] * ζ_[:,i]
        end

        # Compute covariance and cross-covariance of inverted likelihood
        S_t = zeros(Dx,Dx)
        C_t = zeros(Dx,Dx)
        for i = 1:2*Dx+1
            S_t += weights_c[i]*(ζ_[:,i] - μ_t)*(ζ_[:,i] - μ_t)' .+ measurement_noise
            C_t += weights_c[i]*(ξ_[:,i] - m_t_pred)*(ζ_[:,i] - μ_t)'
        end

        # Compute Kalman gain
        K_t = C_t*inv(S_t)

        # Compute filtered state parameters
        m_t = m_t_pred .+ K_t*(observations[:,t] - μ_t)
        P_t = P_t_pred .- K_t*S_t*K_t'

        # Store estimates
        mx[:,t] = m_t
        Px[:,:,t] = P_t

        # Update previous state variable
        m_tmin = m_t
        P_tmin = P_t
    end
    return mx, Px
end

max_stars_repo_name: wmkouw/kalman-filters

max_issues_repo_name: wmkouw/kalman-filters

max_forks_repo_name: wmkouw/kalman-filters


text:
module TestProperties

using Atomix: AtomicRefArray
using Test

function test()
    A = AtomicRefArray(ones(Int, 3))
    @atomic A[1]._ = 123
    @test (@atomic A[1]._) === 123
    @test (@atomic A[1]._ += 1) === 124
    @test (@atomic A[1]._) === 124
    @test (@atomicswap A[1]._ = 567) === 124
    @test (@atomic A[1]._) === 567
    @test (@atomicreplace A[1]._ 567 => 123) === (old = 567, success = true)
    @test (@atomicreplace A[1]._ 567 => 123) === (old = 123, success = false)
end

end  # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "spglib"
version = v"1.16.3"

# Collection of sources required to build spglib
sources = [
    ArchiveSource("https://github.com/atztogo/spglib/archive/v$(version).tar.gz",
                  "1dfe313b460f71de90ee8a01d9f2cd250cd59e16836e1bf64924500dd2aa7dc6"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/spglib-*/
if [[ ${target} == *-mingw32 ]]; then
    sed -i -e 's/LIBRARY/RUNTIME/' CMakeLists.txt
fi
mkdir _build
cd _build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      ..
make -j${nproc}
make install VERBOSE=1
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libsymspg", :libsymspg)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Statistics, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.7"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	dist1 = TriangularDist(0,1,1)
	dist2 = TriangularDist(0,1,0)
	N=10^6

	data1, data2 = rand(dist1,N), rand(dist2,N)
	dataSum = data1 + data2
end

# ╔═╡ 8f60b2f6-20f9-11eb-0a66-771e655035a3
mgf(s) = 4(1+(s-1)*MathConstants.e^s)*(MathConstants.e^s-1-s)/s^4

# ╔═╡ 8f61044a-20f9-11eb-1485-d145145c4c57
mgfPointEst(s) = mean([MathConstants.e^(s*z) for z in
					rand(dist1,20) + rand(dist2,20)])

# ╔═╡ 8f6181a4-20f9-11eb-3eae-6dc331b9a88a
begin
	p1 = histogram(dataSum, bins=80, normed=:true,
		ylims=(0,1.4), xlabel="z", ylabel="PDF")

	sGrid = -1:0.01:1
	p2 = plot(sGrid, mgfPointEst.(sGrid), c=:blue, ylims=(0,3.5))
	p2 = plot!(sGrid, mgf.(sGrid), c=:red)
	p2 = plot!( [minimum(sGrid),maximum(sGrid)],
		[minimum(sGrid),maximum(sGrid)].+1,
		c=:black, xlabel="s", ylabel="MGF")

	plot(p1, p2, legend=:none, size=(800, 400))
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.7"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═8f60b2f6-20f9-11eb-0a66-771e655035a3
# ╠═8f61044a-20f9-11eb-1485-d145145c4c57
# ╠═8f6181a4-20f9-11eb-3eae-6dc331b9a88a
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
n = 2
func = sphere
mu = false
lambda = false
max_fevals = 2000
covarMatrixSampler = SparseCholeskyCovarSampler
utilitiesFunc = linear_utilities
trace = true


text:
"""
    add_nums
Add numbers together.
"""
# Write your package code he
function add_nums end

"""
    add_nums(a, b)
Add two numbers.
"""
function add_nums(a, b)
    a+b
end

"""
    add_nums(a, b ,c)
Add three numbers.
"""
function add_nums(a, b, c)
   a+b+c
end


text:
#=
We shall say that an n-digit number is pandigital if it makes use of all the digits 1 to n exactly once. For example, 2143 is a 4-digit pandigital and is also prime.

What is the largest n-digit pandigital prime that exists?
=#
include("projecteulerutils.jl")
using Combinatorics
using Primes

function calc()
  largestpandigitalprime = 1
  for n in 2:9
    digit_arr = [i for i in 1:n]
    for permutation in permutations(digit_arr)
      x = nd(permutation)
      if isprime(x) && x > largestpandigitalprime
        largestpandigitalprime = x
      end
    end
  end
  largestpandigitalprime
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using DelimitedFiles, DataFrames, CSV
using Random, Printf
using ThreadTools

include("src/parzenwindow.jl")
include("src/gmm.jl")
include("src/eval_metrics.jl")

@doc """
    Main function to load dataset, train Gaussian mixture model and Parzen window estimation, set hyperparametrs on valid data and eval models
    on test data.

    # Examples
    ```jldoctest
    g::Float64, p::Float64 = main(true, "yeast")
    ```
    where the 'g' is GMM roc-AUC on test data and 'p' is PWE roc-AUC on test data.
    The main(show_report::Bool, dataset_name::String) has got two arguments. If 'show_report' is set to true, the whole evaluation report is
    printed, 'dataset_name' is a string name of dataset to be loaded.
    """ ->
function main(show_report::Bool, dataset_name::String)

    # Load all data.
    data_anomal_raw = readdlm(string("data_anomalyproject/", dataset_name, "/anomalous.txt"))'
    data_normal_raw = readdlm(string("data_anomalyproject/", dataset_name, "/normal.txt"))'
    data_anomal::Matrix{Float64} = data_anomal_raw[:, shuffle(1:end)]
    data_normal::Matrix{Float64} = data_normal_raw[:, shuffle(1:end)]

    # Prepare all data.
    N_normal::Int64 = size(data_normal)[2]
    trn_data::Matrix{Float64} = data_normal[:, begin:Int(round(N_normal/2))]
    valid_data::Matrix{Float64} = data_normal[:, Int(round(N_normal/2)):Int(round(3*N_normal/4))]
    test_data_n::Matrix{Float64} = data_normal[:, Int(round(3*N_normal/4)):end]
    N::Int64 = size(data_normal)[1]

    # Define function to get likelihood proxy param.
    lt(x) = x < -10e5 ? -10e5 : x
    likelihood(model, params, data) = mean(tmap(x->lt(log(model(params, Vector{Float64}(x)))), eachcol(data)))
    likelihood(model, data) = mean(tmap(x->lt(log(model(Vector{Float64}(x)))), eachcol(data)))

    @doc """
    Function to choose optimal count of components (hyperparametr) of Gaussian mixture model. The choose is
    performed based on maximization of likelihood (proxy parametr) on given (validation) data.
    Minimum of components is 2, maximum is 10.

    # Examples
    ```jldoctest
    model, Θ = choose_gmm_model(data)
    ```
    where the 'Θ' is dictionary of params of the GMM,
    'data' is set of (validation) data (size(dim, N) ~ (dimension of data, data count),
    'model' is a Gaussian mixture model.
    """ ->
    function choose_gmm_model(data::Matrix{Float64})
        models_dict = Dict{Float64, Vector{Any}}()  # likelihood => [model, params]
        lh::Float64 = 0
        for K::Int64=2:10
            ps::Dict{Symbol, Vector}, gmm_model, gm_model = create_gmm(K, N); # prepare model
            EM!(ps, trn_data, K, gmm_model, gm_model, 60); # learn model params
            lh = likelihood(gmm_model, ps, data)
            models_dict[lh] = [gmm_model, ps]
        end
        return models_dict[maximum(keys(models_dict))]
    end

    @doc """
    Function to choose optimal window-size (hyperparametr) of Parzen window estimation. The choose is
    performed based on maximization of likelihood (proxy parametr) on given (validation) data.
    Minimum of window-size is 0.01, maximum is 10, step is 0.01. The kernel function is Gaussian kernel.

    # Examples
    ```jldoctest
    h = choose_parzenwindow_model(data)
    ```
    where the 'h' is window-size and
    'data' is set of (validation) data (size(dim, N) ~ (dimension of data, data count).
    """ ->
    function choose_parzenwindow_model(data::Matrix{Float64})
        models_dict = Dict{Float64, Float64}()  # likelihood => window-size
        lh::Float64 = 0
        kernel(x) = k(x)
        for step::Float64=0.01:0.01:10
            model(x) = create_parzen_window(step, trn_data, kernel, x) # prepare model
            lh = likelihood(model, data)
            models_dict[lh] = step
        end
        return models_dict[maximum(keys(models_dict))]
    end

    # Learn GMM on train data and choose the best count of components by validation data.
    gmm_model, params::Dict{Symbol, Vector} = choose_gmm_model(valid_data)
    @printf("Best K: %d\n", size(params[:μ])[1])

    # Learn Parzen window estimation on train data and choose the best window-size by validation data.
    h::Float64 = choose_parzenwindow_model(valid_data)
    kernel(x) = k(x)
    parzenwindow(x) = create_parzen_window(h, trn_data, kernel, x)
    @printf("Best window-size: %.3f\n", h)

    # Prepare testing data.
    testing_data::Matrix{Float64} = hcat(test_data_n, data_anomal);
    testing_labels::Vector{Bool} = Vector{Bool}(vcat(ones(size(test_data_n)[2], 1),
                                    zeros(size(data_anomal)[2], 1))[:, 1]);

    # Print evaluation report.
    gmm_auc::Float64 = eval_report(gmm_model, params, testing_data, testing_labels, show_report);
    parzen_auc::Float64 = eval_report(parzenwindow, testing_data, testing_labels, show_report);
    return gmm_auc, parzen_auc
end

@doc """
    Function to learn models on dataset (on each dataset 10-times.) and store roc-AUC into CSV.

    # Examples
    ```jldoctest
    compare_models()
    ```
    """ ->
function compare_models()
    datasets_names::Vector{String} = filter(x->x[1]!='.', readdir("data_anomalyproject/"))
    gmm_auc_arr = Vector{Float64}([])
    parzen_auc_arr = Vector{Float64}([])
    for name in datasets_names
        for iter=1:10
            println(name, " / ", iter)
            g::Float64, p::Float64 = main(false, name)
            @printf("gmm auc: %.3f, parzen auc: %.3f\n", g, p)
            println()
            append!(gmm_auc_arr, g)
            append!(parzen_auc_arr, p)
        end
    end
    df = DataFrame(gmm = gmm_auc_arr, parzen = parzen_auc_arr)
    CSV.write("auc_stat.csv", df)
end

max_stars_repo_name: kozvojtex/AnomalyDetection

max_issues_repo_name: kozvojtex/AnomalyDetection

max_forks_repo_name: kozvojtex/AnomalyDetection


text:
module Eyeball

Base.Experimental.@compiler_options compile=min optimize=1 # infer=false

using Base.Docs

using InteractiveUtils
using REPL: REPL, AbstractTerminal

using REPL.TerminalMenus
import REPL.TerminalMenus: request
using AbstractTrees

import TerminalPager
using Statistics

# TODO: de-vendor once changes are upstreamed.
# include("vendor/FoldingTrees/src/FoldingTrees.jl")
using FoldingTrees

# include("ui.jl")

export eye

default_terminal() = REPL.LineEdit.terminal(Base.active_repl)
"""
Explore objects and types.
```
eye(x = Main, depth = 10; interactive = true)
```
`depth` controls the depth of folding.
`interactive` means browse the object `x` interactively.
With `interactive` set to `false`, `eye` returns the tree as a `FoldingTrees.Node`.

During interactive browsing of the object tree, the following keys are available:

* `↑` `↓` `←` `→` -- Up and down moves through the tree. Left collapses a tree. Right expands a folded tree. Vim movement keys (`h` `j` `k` `l`) are also supported.
* `d` -- Docs. Show documentation on the object.
* `e` -- Expand. Show more subobjects.
* `f` -- Toggle fields. By default, parameters are shown for most objects.
  `f` toggles between the normal view and a view showing the fields of an object.
* `m` -- Methodswith. Show methods available for objects of this type. `M` specifies `supertypes = true`.
* `o` -- Open. Open the object in a new tree view. `O` opens all (mainly useful for modules).
* `r` -- Return tree. Return the tree (a `FoldingTrees.Node`).
* `s` -- Show object.
* `t` -- Typeof. Show the type of the object in a new tree view.
* `z` -- Summarize. Toggle a summary of the object and child objects.
  For arrays, this shows the mean and 0, 25, 50, 75, and 100% quantiles (skipping missings).
* `0`-`9` -- Fold to depth. Also toggles expansion of items normally left folded.
* `enter` -- Return the object.
* `q` -- Quit.
"""
function eye(x = Main, depth = 10; interactive = true, all = false)
    if all
        x = All(x)
    end
    cursor = Ref(1)
    returnfun = x -> x.data.value
    redo = false
    foldctx = (cursoridx = 1, depth = depth, expandall = true)
    function resetterm()
        println(term.out_stream)
        REPL.Terminals.raw!(term, true)
        print(term.out_stream, "\x1b[?25l")  # hide cursor
    end
    function keypress(menu::TreeMenu, i::UInt32)
        if i == Int('j')
            cursor[] = TerminalMenus.move_down!(menu, cursor[])
        elseif i == Int('k')
            cursor[] = TerminalMenus.move_up!(menu, cursor[])
        elseif i == Int('l') || i == Int(TerminalMenus.ARROW_RIGHT)
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            node.foldchildren = false
            menu.pagesize = min(menu.maxsize, count_open_leaves(menu.root))
        elseif i == Int('h') || i == Int(TerminalMenus.ARROW_LEFT)
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            node.foldchildren = has_children(node)
            menu.pagesize = min(menu.maxsize, count_open_leaves(menu.root))
        elseif i == Int('f')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            o = node.data.value
            if isstructtype(typeof(o))
                node.children = Node[]
                node.data.showfields[] = !node.data.showfields[]
                treelist(o, 8, 16, node)
                node.foldchildren = false
                menu.pagesize = min(menu.maxsize, count_open_leaves(menu.root))
            end
        elseif i == Int('d')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            _pager(getdoc(node.data.value))
            resetterm()
        elseif i == Int('e')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            o = node.data.value
            n = length(node.children)
            if length(node.children) < 1 || !(node.children[1].data.value isa ExpandPlaceholder)
                return false
            end
            placeholder = popfirst!(node.children)
            ctx = placeholder.data.value
            shown = length(node.children)  # double what's shown
            _iterate(ctx.itr, ctx.depth, shown, node, ctx.history)
            menu.pagesize = min(menu.maxsize, count_open_leaves(menu.root))
        elseif i == Int('m')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            println(term.out_stream, "\n\nOpening methodswith(`$(node.data.key)`) ...\n")
            o = node.data.value
            newchoice = eye(methodswith(o isa DataType ? o : typeof(o)))
            resetterm()
            if newchoice !== nothing
                returnfun = x -> newchoice
                menu.chosen = true
                node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
                return true
            end
        elseif i == Int('M')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            println(term.out_stream, "\n\nOpening methodswith(`$(node.data.key)`, supertypes = true) ...\n")
            o = node.data.value
            newchoice = eye(methodswith(o isa DataType ? o : typeof(o), supertypes = true))
            resetterm()
            if newchoice !== nothing
                returnfun = x -> newchoice
                menu.chosen = true
                node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
                return true
            end
        elseif i == Int('o')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            println(term.out_stream, "\n\nOpening `$(node.data.key)` ...\n")
            newchoice = eye(node.data.value)
            resetterm()
            if newchoice !== nothing
                returnfun = x -> newchoice
                menu.chosen = true
                node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
                return true
            end
        elseif i == Int('O')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            println(term.out_stream, "\n\nOpening `$(node.data.key)` ...\n")
            newchoice = eye(All(node.data.value))
            resetterm()
            if newchoice !== nothing
                returnfun = x -> newchoice
                menu.chosen = true
                node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
                return true
            end
        elseif i == Int('r')
            returnfun = identity
            menu.chosen = true
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            return true
        elseif i == Int('s')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            io = IOContext(IOBuffer(), :displaysize => displaysize(term), :limit => true, :color => true)
            show(io, MIME"text/plain"(), node.data.value)
            sobj = String(take!(io.io))
            _pager(node.data.value)
            resetterm()
        elseif i == Int('S')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            hasmethod(sort, (typeof(node.data.value),)) || return false
            println(term.out_stream, "\n\nOpening sort(`$(node.data.key)`) ...\n")
            choice = eye(sort(node.data.value))
            resetterm()
        elseif i == Int('t')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            println(term.out_stream, "\n\nOpening typeof(`$(node.data.key)`) ...\n")
            choice = eye(typeof(node.data.value))
            resetterm()
        elseif i == Int('z')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            d = node.data
            d.summarize[] = !d.summarize[]
            if d.summarize[]
                d.string[] = tostring(d.key, Summarize(d.value))
            else
                d.string[] = tostring(d.key, d.value)
            end
            for n in node.children
                dc = n.data
                dc.summarize[] = d.summarize[]
                if dc.summarize[]
                    dc.string[] = tostring(dc.key, Summarize(dc.value))
                else
                    dc.string[] = tostring(dc.key, dc.value)
                end
            end
        elseif i in Int.('0':'9')
            node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
            expandall = foldctx.cursoridx == menu.cursoridx && foldctx.depth == i && foldctx.expandall
            fold!(node, i - 48, expandall)
            foldctx = (cursoridx = menu.cursoridx, depth = i, expandall = !expandall)
            menu.pagesize = min(menu.maxsize, count_open_leaves(menu.root))
        end
        return false
    end
    shown = 50
    root = treelist(x, depth - 1, shown)
    if interactive
        term = default_terminal()
        while true
            menu = TreeMenu(root, pagesize = REPL.displaysize(term)[1] - 2, dynamic = true, keypress = keypress)
            choice = TerminalMenus.request(term, "[f] fields [d] docs [e] expand [m/M] methodswith [o] open [r] tree [s] show [S] Sort [t] typeof [z] summarize [q] quit", menu; cursor=cursor)
            choice !== nothing && return returnfun(choice)
            if redo
                redo = false
                root = treelist(x, depth - 1, shown)
            else
                return
            end
        end
    else
        menu = TreeMenu(root, pagesize = 20, dynamic = true, maxsize = 30, keypress = keypress)
        return root
    end
end

struct ObjectWrapper{K,V}
    key::K
    value::V
    string::Base.RefValue{Union{Nothing,String}}
    showfields::Base.RefValue{Bool}
    summarize::Base.RefValue{Bool}
    ObjectWrapper{K,V}(key, value, string = Ref{Union{Nothing,String}}(nothing), showfields = Ref(false), summarize = Ref(false)) where {K,V} =
        new(key, value, string, showfields, summarize
)
end

Base.show(io::IO, x::ObjectWrapper) = print(io, tostring(x.key, x.value))

struct UNDEFPlaceholder end
const UNDEF = UNDEFPlaceholder()

struct ExpandPlaceholder{T}
    itr::T
    depth::Int
    shown::Int
    history
end

struct All{T}
    value
end
All(x) = All{typeof(x)}(x)

struct Summarize{T}
    value
end
Summarize(x) = Summarize{typeof(x)}(x)
# Don't wrap some types
Summarize(x::ExpandPlaceholder) = x
Summarize(x::UNDEFPlaceholder) = x

function Base.show(io::IO, x::Summarize{<:AbstractArray{T}}) where T<:Union{Real,Missing}
    if eltype(x) >: Missing
        #nm = nmissing(x)
        v = skipmissing(x.value)
    else
        v = x
    end
    print(io, "x̄=", mean(v), ", q=", quantile(v, [0, .25, 0.5, 0.75, 1]))
    # nm > 0 && print(io, ", nm=", nm)
    nothing
end
Base.show(io::IO, x::Summarize) = show(io, x.value)

getdoc(x) = doc(x)
getdoc(x::Method) = doc(x.module.eval(x.name))

has_children(x) = length(x.children) > 0

function fold!(node, depth, expandall)
    if depth == 0
        node.foldchildren = has_children(node)
    else
        if foldobject(node.data.value) && !expandall
            node.foldchildren = has_children(node)
        else
            node.foldchildren = false
            for child in AbstractTrees.children(node)
                fold!(child, depth - 1, expandall)
            end
        end
    end
end

# from https://github.com/MichaelHatherly/InteractiveErrors.jl/blob/5e2e90f9636d748aa3aae0887e18df388829b8e7/src/InteractiveErrors.jl#L52-L61
function style(str; kws...)
    sprint(; context = :color => true) do io
        printstyled(
            io, str;
            bold = get(kws, :bold, false),
            color = get(kws, :color, :normal),
        )
    end
end

struct UndefIterator{T}
    x::T
end
Base.iterate(a::UndefIterator, state = 1) = (isassigned(a.x, state) ? a.x[state] : UNDEF, state + 1)
Base.length(a::UndefIterator) = length(a.x)
undefiterator(x) = x
undefiterator(x::AbstractArray) = UndefIterator(x)

function treelist(x, depth = 0, shown = 10, parent = Node{ObjectWrapper}(ObjectWrapper{String,typeof(x)}("", x)), history = Base.IdSet{Any}((x,)))
    usefields = parent.data.showfields[] && isstructtype(typeof(x)) && !(x isa DataType)
    itr = Iterators.Stateful(usefields ? getfields(x) : getoptions(x))
    _iterate(itr, depth, shown, parent, history)
    # parent.itr[] = itr
    return parent
end

function _iterate(itr, depth, shown, parent, history)
    for (k,v) in Iterators.take(itr, shown)
        node = Node{ObjectWrapper}(ObjectWrapper{typeof(k),typeof(v)}(k, v),
                    parent,
                    foldobject(v) || (depth < 1 && shouldrecurse(v)))
        if depth > -20 && v ∉ history && shouldrecurse(v)
            treelist(v, depth - 1, 10, node, push!(copy(history), v))
        end
        if !has_children(node)
            node.foldchildren = false
        end
    end
    if 0 < length(itr) < 4
        _iterate(itr, depth, shown, parent, history)
    end
    if length(itr) > 0
        marker = Node{ObjectWrapper}(ObjectWrapper{Symbol,ExpandPlaceholder}(:!, ExpandPlaceholder(itr, depth, shown, history)))
        marker.parent = parent
        pushfirst!(parent.children, marker)
    end
end

function FoldingTrees.writeoption(buf::IO, obj::ObjectWrapper, charsused::Int; width::Int=(displaysize(stdout)::Tuple{Int,Int})[2])
    if obj.string[] == nothing
        obj.string[] = tostring(obj.key, obj.value)
    end
    FoldingTrees.writeoption(buf, obj.string[], charsused; width=width)
end

function _pager(object)
     buffer = IOBuffer()
     ctx = IOContext(buffer, :color => true, :compact => true, :limit => true, :displaysize => (3000,500))
     show(ctx, "text/plain", object)
     TerminalPager.pager(String(take!(buffer)))
end

############################
#  API functions
############################

"""
```
tostring(pn, obj)
```
Returns a string with the text representation of `obj` with key `pn`.
"""
function tostring(key, value)
    io = IOContext(IOBuffer(), :compact => true, :limit => true, :color => true, :displaysize => (3,200))
    show(io, value)
    svalue = String(take!(io.io))
    string(style(string(key), color = :cyan), ": ",
           style(string(typeof(value)), color = :green), " ",
           extras(value), " ",
           svalue)
end

function tostring(key, value::UNDEFPlaceholder)
    string(style(string(key), color = :cyan), ": #undef")
end

function tostring(key, value::ExpandPlaceholder)
    N = length(value.itr.itr)
    n = length(value.itr)
    string(style(string(key), color = :red), "   Showing items 1-", N-n, " of ", N, ". Items remaining=", n, ". Hit [e] to expand.")
end

function tostring(key, x::Summarize)
    io = IOContext(IOBuffer(), :compact => true, :limit => true, :color => true, :displaysize => (3,200))
    show(io, x)
    svalue = String(take!(io.io))
    string(style(string(key), color = :cyan), ": ",
           style(string(typeof(x.value)), color = :green), " ",
           extras(x.value), " ",
           svalue)
end

"""
```
extras(x) = ""
```
Returns a string with any extra information about `x`.
For AbstractArrays, this returns size information.
"""
extras(x) = ""
extras(x::AbstractArray) = style(string(size(x)), color = :magenta) * " " * style(string(sizeof(x)), color = :yellow)

"""
```
getfields(x)
```
Return a tuple of two arrays describing the child objects to be shown for `x`.
The first array has the field names of `x`, and the second array has the fields.
Normally, this should not have a custom definition for a type.
Use `getoptions` for that.
"""
function getfields(x::T) where T
    !isstructtype(T) && return nothing
    keys = fieldnames(typeof(x))
    values = (isdefined(x, i) ? getfield(x, i) : UNDEF for i in 1:length(keys))
    return zip(keys, values)
end

"""
```
getoptions(x)
```
Return a tuple of two arrays describing the child objects to be shown for `x`.
The first array has the keys or indexes of the child objects, and the second array is the child objects.
"""
function getoptions(x)
    keys = propertynames(x)
    if keys === fieldnames(typeof(x))
        values = (isdefined(x, pn) ? getproperty(x, pn) : UNDEF for pn in keys)
    else
        values = (getproperty(x, pn) for pn in keys)
    end
    return zip(keys, values)
end
function getoptions(x::All)
    return getoptions(x.value)
end
function getoptions(x::All{Module})
    y = x.value
    keys = names(y, all = true)
    values = (isdefined(y, pn) ? getproperty(y, pn) : UNDEF for pn in keys)
    return zip(keys, values)
end
function getoptions(x::DataType)
    if isabstracttype(x)
        st = filter(t -> t !== Any, subtypes(x))
        return zip(Iterators.repeated(Symbol("")), st)
    end
    x <: Tuple && return nothing
    if x.name === Base.NamedTuple_typename && !(x.parameters[1] isa Tuple)
        # named tuple type with unknown field names
        return nothing
    end
    fields = fieldnames(x)
    fieldtypes = Base.datatype_fieldtypes(x)
    return zip(fields, fieldtypes)
end
function getoptions(x::AbstractArray{T}) where T
    return zip(eachindex(x), UndefIterator(x))
end
function getoptions(x::AbstractDict{<:S, T}) where {S<:Union{AbstractString, Symbol, Number},T}
    return x
end
function getoptions(x::AbstractDict)
    return zip(Iterators.repeated(Symbol("")), (k => v for (k,v) in x))
end
function getoptions(x::AbstractSet{T}) where T
    return zip(Iterators.repeated(Symbol("")), x)
end

iscorejunk(x) = parentmodule(parentmodule(parentmodule(x))) === Core && !isabstracttype(x) && isstructtype(x)

"""
```
shouldrecurse(x)
```
A boolean that controls whether eye recurses into object `x`.
`len` is the length of the object.
This defaults to true.
"""
shouldrecurse(x) = true
shouldrecurse(::Module) = false
shouldrecurse(::Method) = false
shouldrecurse(::Core.MethodInstance) = false
shouldrecurse(::TypeVar) = false
shouldrecurse(x::DataType) = x !== Any && x !== Function && !iscorejunk(x)
shouldrecurse(::Union) = false
shouldrecurse(::UNDEFPlaceholder) = false

"""
```
foldobject(x)
```
A boolean that controls whether `eye` automatically folds `x`.
This is useful for types where the components are usually not needed.
This defaults to false.
"""
foldobject(x) = false
foldobject(x::AbstractArray) = true
foldobject(x::AbstractVector{Any}) = length(x) > 10
foldobject(::UnitRange) = true
foldobject(x::Number) = true
foldobject(x::DataType) = !isabstracttype(x) && isstructtype(x)
foldobject(x::UnionAll) = true
foldobject(x::Pair) = true

end


text:
# Defines a set of standard callback functions

"""
Default callback function that stores x & W of the simulations.
"""
function default_callback(out)
    return function(t,x,W)
    push!(out[1], copy(x))
    push!(out[2], copy(W))
    end
end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
export interactive_billiard, interactive_billiard_bmap,
       billiard_video, billiard_video_timeseries

"""
    interactive_billiard(bd::Billiard [, x, y, φ] [, ω=nothing]; kwargs...)
    interactive_billiard(bd::Billiard, ps::Vector{<:AbstractParticle}; kwargs...)

Launch an interactive application that evolves particles in a dynamical billiard `bd`, using
[DynamicalBilliards.jl](https://juliadynamics.github.io/DynamicalBilliards.jl/dev/).
Requires `DynamicalBilliards`.

You can either specify exactly the particles that will be used `ps` or provide
some initial conditions `x,y,φ,ω`, which by default are random in the billiard.
The particles are evolved in real time instead of being pre-calculated,
so the application can be left to run for infinite time.

See also [`interactive_billiard_bmap`](@ref) and [`billiard_video`](@ref).

## Interaction
Push "play" to start evolving particles in the billiard, and "reset" to restore them
to their (latest) initial condition. The "particles" hides or shows the particles.
The "speed" slider controls the animation speed (in fact, it controls how often
are the plots updated).

Clicking and dragging inside the billiard plot shows a line. When the line is selected,
new particles are created that have the direction of this line, as well as its starting
position, using the function `particlebeam` from `DynamicalBilliards`.

## Further keywords
* `N = 100` : if exact particles are not given, `N` are created. Otherwise it is `length(ps)`.
* `dx = 0.01` : width of the particle beam.
* `dt = 0.001` : time resolution of the animation.
* `tail = 1000` : length of the tail of the particles (multiplies `dt`).
* `colors = JULIADYNAMICS_COLORS` : If a symbol (colormap name) each particle gets
  a color from the map. If Vector of length `N`, each particle gets a color form the vector.
  If Vector with length < `N`, linear interpolation across contained colors is done.
* `tailwidth = 1` : Linewidth of the particle tail.
* `fade = true` : Whether to add fadeout to the particle tail.
* `sleept = nothing` : If the slowest speed of the animation is already too fast,
  give a small number to `sleept`.
* `plot_particles = true` : If false, the particles are not plotted (as balls and arrows).
  This makes the application faster (you cannot show them again with the button).
* `particle_size = 1.0` A multiplier for the size of the particles.
"""
interactive_billiard(bd::Billiard, ω::Union{Nothing, Real} = nothing; kwargs...) =
interactive_billiard(bd::Billiard, DynamicalBilliards.randominside_xyφ(bd)..., ω; kwargs...)

function interactive_billiard(bd::Billiard, x::Real, y::Real, φ::Real, ω = nothing;
    kwargs...)
    N = get(kwargs, :N, 100)
    dx = get(kwargs, :dx, 0.01)
    ps = DynamicalBilliards.particlebeam(x, y, φ, N, dx, ω, eltype(bd))
    interactive_billiard(bd::Billiard, ps; kwargs...)
end

function interactive_billiard(bd::Billiard, ps::Vector{<:AbstractParticle};
        dt = 0.001, tail = 1000, dx = 0.01, colors = JULIADYNAMICS_COLORS,
        plot_particles = true, α = 1.0, N = 100, res = (800, 800),
        intervals = nothing, sleept = nothing, fade = true,
        backgroundcolor = DEFAULT_BG,
        vr = _estimate_vr(bd),
        add_controls = true,
        displayfigure = true,
        tailwidth = 1,
        particle_size = 1.0,
    )

    N = length(ps)
    p0s = deepcopy(ps) # deep is necessary because vector of mutables
    ω0 = DynamicalBilliards.ismagnetic(ps[1]) ? ps[1].ω : nothing

    # Initialized inside process
    cs = if !(colors isa Vector) || length(colors) ≠ N
        colors_from_map(colors, α, N)
    else
        to_color.(colors)
    end
    figure = Figure(resolution = res, backgroundcolor = backgroundcolor)
    ax = figure[1, 1] = Axis(figure, backgroundcolor = backgroundcolor)
    tight_ticklabel_spacing!(ax)
    ax.autolimitaspect = 1
    allparobs = [ParObs(p, bd, tail) for p in ps]

    # Plot tails
    for (i, p) in enumerate(allparobs)
        x = to_color(cs[i])
        if fade
            x = [RGBAf0(x.r, x.g, x.b, i/tail) for i in 1:tail]
        end
        lines!(ax, p.tail; color = x, linewidth = tailwidth)
    end

    if plot_particles # plot ball and arrow as a particle
        balls = Observable([Point2f0(p.p.pos) for p in allparobs])
        vels = Observable([particle_size * vr * Point2f0(p.p.vel) for p in allparobs])
        particle_plots = (
            scatter!(
                ax, balls; color = darken_color.(cs),
                marker = MARKER, markersize = 8*particle_size*AbstractPlotting.px
            ),
            arrows!(
                ax, balls, vels; arrowcolor = darken_color.(cs),
                linecolor = darken_color.(cs),
                normalize = false, arrowsize = particle_size*vr/3,
                linewidth  = particle_size*4,
            )
        )
    end

    # Plot billiard (after particles, so that obstacles hide the collision points)
    bdplot!(ax, bd)

    # Controls
    if add_controls
        resetbutton = Button(figure;
            label = "reset", buttoncolor = RGBf0(0.8, 0.8, 0.8),
            height = 40, width = 80
        )
        runbutton = Button(figure; label = Observable("run"),
            buttoncolor = Observable(RGBf0(0.8, 0.8, 0.8)),
            buttoncolor_hover = Observable(RGBf0(0.7, 0.7, 0.9)),
            buttoncolor_active = Observable(RGBf0(0.6, 0.6, 1.0)),
            labelcolor = Observable((RGBf0(0,0,0))),
            labelcolor_active = Observable((RGBf0(1,1,1))),
            height = 40, width = 70,
        )
        nslider = Slider(figure, range = 0:50, startvalue=0)
        controls = [resetbutton, runbutton, Label(figure, "speed:"), nslider]
        if plot_particles
            particlebutton = Button(figure, label = "particles",
                buttoncolor = RGBf0(0.8, 0.8, 0.8),
                height = 40, width = 100
            )
            pushfirst!(controls, particlebutton)
        end

        # Functionality that CREATES the play/stop button
        # TODO: will be deleted once MakieLayout has proper togglable button
        on(runbutton.clicks) do n
            t = runbutton.label[] == "run" ? "stop" : "run"
            runbutton.label[] = t
            for (s1, s2) in ((:buttoncolor, :buttoncolor_active), (:labelcolor, :labelcolor_active))
                getproperty(runbutton, s1)[], getproperty(runbutton, s2)[] =
                    getproperty(runbutton, s2)[], getproperty(runbutton, s1)[]
            end
            runbutton.labelcolor_hover[] = runbutton.labelcolor[]
        end

        # Create the "control panel"
        figure[2, 1] = grid!(hcat(controls...,), tellwidth = false, tellheight = true)

        # Play/stop
        on(runbutton.clicks) do nclicks
            @async while runbutton.label[] == "stop"
                for i in 1:N
                    p = ps[i]
                    parobs = allparobs[i]
                    animstep!(parobs, bd, dt, true, intervals)
                    if plot_particles
                        balls[][i] = parobs.p.pos
                        vels[][i] = vr * parobs.p.vel
                    end
                end
                if plot_particles
                    balls[] = balls[]
                    vels[] = vels[]
                end
                for _ in 1:nslider.value[]
                    for i in 1:N
                        p = ps[i]
                        parobs = allparobs[i]
                        animstep!(parobs, bd, dt, false, intervals)
                    end
                end
                if sleept == nothing
                    yield()
                else
                    sleep(sleept)
                end
                isopen(figure.scene) || break
            end
        end

        # Resetting
        on(resetbutton.clicks) do nclicks
            # update parobs
            for i in 1:N
                parobs = allparobs[i]
                rebind_partobs!(parobs, p0s[i], bd)
                if plot_particles
                    balls[][i] = parobs.p.pos
                    vels[][i] = vr * parobs.p.vel
                end
            end
            if plot_particles
                balls[] = balls[]
                vels[] = vels[]
            end
            yield()
        end

        # Show/hide particles
        if plot_particles
            on(particlebutton.clicks) do nclicks
                particle_plots[1].visible[] = !particle_plots[1].visible[]
                for i in 1:2
                particle_plots[2].plots[i].visible[] = !particle_plots[2].plots[i].visible[]
                end
            end
        end

        # Selecting new particles
        MakieLayout.deactivate_interaction!(ax, :rectanglezoom)
        sline = select_line(ax.scene)
        on(sline) do val
            pos = val[1]
            dir = val[2] - val[1]
            φ = atan(dir[2], dir[1])
            p0s .= DynamicalBilliards.particlebeam(pos..., φ, N, dx, ω0, eltype(bd))
            resetbutton.clicks[] += 1
        end
    end # adding controls blocks

    displayfigure && display(figure)
    if add_controls
        return figure, allparobs, resetbutton, p0s, sline
    elseif plot_particles
        return figure, allparobs, balls, vels, vr
    else
        return figure, allparobs, nothing, nothing, nothing
    end
end

"""
    billiard_video(file, bd::Billiard [, x, y, φ] [, ω=nothing]; kwargs...)
    billiard_video(file, bd::Billiard, ps::Vector{<:AbstractParticle}; kwargs...)

Perform the same animation like in [`interactive_billiard`](@ref), but there is no
interaction; the result is saved directly as a video in `file` (no buttons are shown).

## Keywords
* `N, dt, tail, dx, colors, plot_particles, fade, tailwidth, backgroundcolor`:
  same as in `interactive_billiard`, but `plot_particles` is `false` by default here.
* `speed = 5`: Animation "speed" (how many `dt` steps are taken before a frame is recorded)
* `frames = 1000`: amount of frames to record.
* `framerate = 60`: exported framerate.
* `res = nothing`: resolution of the frames. If nothing, a resolution matching the
  the billiard aspect ratio is estimated. Otherwise pass a 2-tuple.

Notice that the animation performs an extra step for every `speed` steps and the
first frame saved is always at time 0. Therefore the following holds:
```julia
total_time = (frames-1)*(speed+1)*dt
time_covered_per_frame = (speed+1)*dt
```
"""
billiard_video(file::String, bd::Billiard, ω::Union{Nothing, Real} = nothing; kwargs...) =
billiard_video(file::String, bd::Billiard, DynamicalBilliards.randominside_xyφ(bd)..., ω; kwargs...)

function billiard_video(file::String, bd::Billiard, x::Real, y::Real, φ::Real, ω = nothing;
    kwargs...)
    N = get(kwargs, :N, 500)
    dx = get(kwargs, :dx, 0.01)
    ps = DynamicalBilliards.particlebeam(x, y, φ, N, dx, ω, eltype(bd))
    billiard_video(file, bd, ps; kwargs...)
end

function billiard_video(file::String, bd::Billiard, ps::Vector{<:AbstractParticle};
        plot_particles = false, res = nothing, dt = 0.001,
        speed = 5, frames = 1000, framerate = 60, kwargs...
    )

    dt = eltype(bd)(dt)
    if res == nothing
        xmin, ymin, xmax, ymax = DynamicalBilliards.cellsize(bd)
        aspect = (xmax - xmin)/(ymax-ymin)
        res = (round(Int, aspect*1000), 1000)
    end

    figure, allparobs, balls, vels, vr = interactive_billiard(bd, ps;
        res = res, plot_particles=plot_particles, kwargs..., add_controls = false,
        displayfigure = false
    )
    N = length(ps)

    AbstractPlotting.inline!(true) # to not display figure while recording
    record(figure, file, 1:frames; framerate = framerate) do j
        for i in 1:N
            p = ps[i]
            parobs = allparobs[i]
            animstep!(parobs, bd, dt, true)
            if plot_particles
                balls[][i] = parobs.p.pos
                vels[][i] = vr * parobs.p.vel
            end
        end
        if plot_particles
            balls[] = balls[]
            vels[] = vels[]
        end
        for _ in 1:speed
            for i in 1:N
                p = ps[i]
                parobs = allparobs[i]
                animstep!(parobs, bd, dt, false)
            end
        end
    end
    AbstractPlotting.inline!(false)
    return
end

"""
    billiard_video_timeseries(file, bd::Billiard, ps, f; kwargs...)

Perform the same animation like in [`billiard_video`](@ref), but in addition show the
timeseries of a chosen observable above the billiard. The observable is given using the
**function** `f`, which takes as an input a particle and outputs the observable.
E.g. `f(p) = p.pos[2]` or `f(p) = atan(p.vel[2], p.vel[1])`.
The video is saved directly into `file`.

## Keywords
* `N, dt, tail, dx, colors, plot_particles, fade, tailwidth, backgroundcolor`:
  same as in `interactive_billiard`.
* `speed, frames, framerate, res`: As in `billiard_video`.
* `total_span = 10.0`: Total span of the x-axis of the timeseries plot in real time units.
* `ylim = (0, 1)`: Limits of the y-axis of the timeseries plot.
* `ylabel = "f"`: Label of the y-axis of the timeseries plot.
"""
function billiard_video_timeseries(file::AbstractString, bd::Billiard, ps::Vector{<:AbstractParticle}, f;
        plot_particles = true, dt = 0.001,
        speed = 5, frames = 1000, framerate = 60,
        total_span = 10.0, colors = JULIADYNAMICS_COLORS,
        res = (800, 800), displayfigure = false, ylim = (0, 1),
        ylabel = "f",
        kwargs...
    )

    N = length(ps)
    tdt = total_span/20
    cs = if !(colors isa Vector) || length(colors) ≠ N
        colors_from_map(colors, α, N)
    else
        to_color.(colors)
    end
    dt = eltype(bd)(dt)

    figure, allparobs, balls, vels, vr = interactive_billiard(bd, ps;
        res, plot_particles, kwargs..., add_controls = false,
        displayfigure, colors = cs,
    )

    # Add the axis on the top
    tsax = figure[0, :] = Axis(figure; height = Relative(8/9))
    # Make the axis occupy 1/3 instead of 1/2:
    rowsize!(figure.layout, 1, Relative(1/3))
    tsax.xlabel = "time"
    tsax.ylabel = ylabel

    all_ts = [Observable([Point2f0(0, f(p))]) for p in ps]
    all_balls = [Observable(Point2f0(0, f(p))) for p in ps]

    for i in 1:length(ps)
        lines!(tsax, all_ts[i];
            color = cs[i],
            linewidth = 4,
        )
        scatter!(
            tsax, all_balls[i];
            markersize = 20*AbstractPlotting.px,
            color = InteractiveDynamics.to_alpha(cs[i], 0.75),
        )
    end
    ylims!(tsax, ylim)
    t_current = 0.0
    xlims!(tsax, -tdt, total_span+tdt)

    !displayfigure && AbstractPlotting.inline!(true)
    record(figure, file, 1:frames; framerate = framerate) do j
        # This loop propagates the particles for `speed` steps but doesn't update the plot
        t_current += dt*speed
        for _ in 1:speed
            for i in 1:N
                p = ps[i]
                parobs = allparobs[i]
                InteractiveDynamics.animstep!(parobs, bd, dt, false)
            end
        end

        # This loop propagates the particles for 1 step & updates the plots
        t_current += dt
        for i in 1:N
            p = ps[i]
            parobs = allparobs[i]
            InteractiveDynamics.animstep!(parobs, bd, dt, true)
            # Update timeseries
            current_point = Point2f0(t_current, f(p))
            all_balls[i][] = current_point
            InteractiveDynamics.pushupdate!(all_ts[i], current_point)
            if plot_particles
                balls[][i] = parobs.p.pos
                vels[][i] = vr * parobs.p.vel
            end
        end
        if plot_particles
            balls[] = balls[]
            vels[] = vels[]
        end

        t_prev = max(0, t_current - total_span)
        xlims!(tsax, t_prev-tdt, max(t_current, total_span)+tdt)
    end
    !displayfigure && AbstractPlotting.inline!(false)
    return figure
end

"""
    interactive_billiard_bmap(bd::Billiard, ω=nothing; kwargs...)
Launch an interactive application whose left part is [`interactive_billiard`](@ref)
and whose write part is an interactive boundary map of the billiard (see "Phase spaces"
in DynamicalBilliards.jl).

A particle evolved in the real billiard is also shown on the boundary map.
All interaction of the billiard works as before, but there is also interaction in
the boundary map: clicking on it will generate a particle whose boundary map is
the clicked point.

The mean collision time "m.c.t." of the particle is shown as well.

## Keywords
* `newcolor = randomcolor` A function which takes as input `(pos, vel, ξ, sφ)` and
  outputs a color (for the scatter points in the boundary map).
* `ms = 12` markersize (in pixels).
* `dt, tail, sleept, fade` : propagated to `interactive_billiard`.
"""
function interactive_billiard_bmap(bd::Billiard, ω=nothing;
        newcolor = randomcolor, ms = 10, lock = true,
        kwargs...
    )

    intervals = DynamicalBilliards.arcintervals(bd)
    figure, allparobs, resetbutton, p0s, sline = interactive_billiard(
        bd, ω; kwargs..., N = 1, intervals = intervals, res = (1600, 800)
    )
    parobs = allparobs[1] # only one exists.

    sublayout = GridLayout()
    cleanbutton = Button(figure, label = "clean", tellwidth = false)
    sublayout[2, 1] = cleanbutton
    mct = Observable("m.c.t. = 0.0")
    mcttext = Label(figure, mct, haligh = :left, tellwidth = false)
    sublayout[2, 2] = mcttext
    bmapax = sublayout[1,:] = Axis(figure)
    bmapax.xlabel = "arclength, ξ"
    bmapax.ylabel = "sine of normal angle, sin(θ)"
    bmapax.targetlimits[] = BBox(intervals[1], intervals[end], -1, 1)

    current_color = Observable(newcolor(parobs.p.pos, parobs.p.vel, parobs.ξsin[]...))
    scatter_points = Observable(Point2f0[])
    scatter_colors = Observable(RGBAf0[])

    scatter!(bmapax.scene, scatter_points; color = scatter_colors,
        marker = MARKER, markersize = ms*AbstractPlotting.px
    )

    # Make obstacle axis, add info about where each obstacle terminates
    if length(intervals) > 2 # at least 2 obstacles
        add_obstacle_axis!(figure, sublayout, intervals, bmapax, lock)
    end

    # Obtain new color when selecting line in main plot
    on(sline) do val
        current_color[] = newcolor(parobs.p.pos, parobs.p.vel, parobs.ξsin[]...)
    end
    # Whenever boundary map is updated, plot the update
    on(parobs.ξsin) do v
        pushupdate!(scatter_points, v)
        pushupdate!(scatter_colors, current_color[])
        τ = parobs.T/parobs.n
        mct[] = "m.c.t. = $(rpad(string(τ), 25))"
    end

    # Selecting a point on the boundary map:
    MakieLayout.deactivate_interaction!(bmapax, :rectanglezoom)
    spoint = select_point(bmapax.scene)
    on(spoint) do ξsin
        ξ, sφ = ξsin
        sφ = clamp(sφ, -1, 1)
        ξsin = Point2f0(ξ, sφ)
        if !(intervals[1] ≤ ξ ≤ intervals[end])
            error("selected point is outside the boundary map")
        end
        pos, vel = DynamicalBilliards.from_bcoords(ξsin..., bd, intervals)
        current_color[] = newcolor(pos, vel, ξsin...)
        p0s[1] = ω ≠ nothing ? DynamicalBilliards.MagneticParticle(pos, vel, ω) : DynamicalBilliards.Particle(pos, vel)
        DynamicalBilliards.propagate!(p0s[1], 10eps(eltype(bd))) # ensure no bad behavior with boundary
        rebind_partobs!(parobs, p0s[1], bd, ξsin)
        resetbutton.clicks[] += 1 # reset main billiard
    end

    # Clean button functionality
    on(cleanbutton.clicks) do nclicks
        scatter_points[] = Point2f0[]
        scatter_colors[] = RGBAf0[]
    end

    # Lock zooming
    if lock
        bmapax.xpanlock = true
        bmapax.ypanlock = true
        bmapax.xzoomlock = true
        bmapax.yzoomlock = true
    end

    figure[:, 2] = sublayout
    display(figure)
    return figure, parobs
end

function add_obstacle_axis!(figure, sublayout, intervals, bmapax, lock)
    ticklabels = ["$(round(ξ, sigdigits=4))" for ξ in intervals[2:end-1]]
    bmapax.xticks = (Float32[intervals[2:end-1]...], ticklabels)
    for (i, ξ) in enumerate(intervals[2:end-1])
        lines!(bmapax.scene, [Point2f0(ξ, -1), Point2f0(ξ, 1)], linestyle = :dash, color = :black)
    end
    obstacle_ticklabels = String[]
    obstacle_ticks = Float32[]
    for (i, ξ) in enumerate(intervals[1:end-1])
        push!(obstacle_ticks, ξ + (intervals[i+1] - ξ)/2)
        push!(obstacle_ticklabels, string(i))
    end

    obax = sublayout[1,:] = Axis(figure)
    obax.xticks = (obstacle_ticks, obstacle_ticklabels)
    obax.xaxisposition = :top
    obax.xticklabelalign = (:center, :bottom)
    obax.xlabel = "obstacle index"
    obax.xgridvisible = false
    hideydecorations!(obax)
    hidespines!(obax)
    xlims!(obax, 0, intervals[end])
    if lock
        obax.xpanlock = true
        obax.ypanlock = true
        obax.xzoomlock = true
        obax.yzoomlock = true
    end
    return obax
end

export billiard_bmap_plot

"""
    billiard_bmap_plot(bd::Billiard, ps::Vector{<:AbstractParticle}; kwargs...)
Return a static figure which has particles plotted on both the real billiard as well
the boundary map, each with its own color (keyword `colors`), and the same color is used
for the corresponding scatter points in the boundary map.

All keyword arguments are the same as [`interactive_billiard_bmap`](@ref), besides
the interaction part of course. The additional keyword `steps` counts how many
times to progress the particles (in multiples of `dt`).
"""
function billiard_bmap_plot(bd::Billiard, ps::Vector{<:AbstractParticle};
        ms = 8, plot_particles=true, colors = JULIADYNAMICS_COLORS,
        dt = 0.001, steps = round(Int, 10/dt), kwargs...
    )

    N = length(ps)
    intervals = DynamicalBilliards.arcintervals(bd)
    figure, allparobs, balls, vels, vr = interactive_billiard(bd, ps;
        kwargs..., dt = dt, add_controls =false, plot_particles=plot_particles,
        intervals = intervals, res = (1600, 800), colors = colors
    )
    cs = (!(colors isa Vector) || length(colors) ≠ N) ? InteractiveDynamics.colors_from_map(colors, 1.0, N) : colors

    sublayout = GridLayout()
    bmapax = sublayout[1,1] = Axis(figure)
    bmapax.xlabel = "arclength, ξ"
    bmapax.ylabel = "sine of normal angle, sin(θ)"
    ylims!(bmapax, -1, 1)
    xlims!(bmapax, intervals[1], intervals[end])
    bmapax.xticklabelsize = 28
    bmapax.yticklabelsize = 28
    bmapax.xlabelsize = 36
    bmapax.ylabelsize = 36
    bmapax.ylabelpadding = 20
    if length(intervals) > 2 # at least 2 obstacles
        obax = add_obstacle_axis!(figure, sublayout, intervals, bmapax, false)
        obax.xticklabelsize = 28
        obax.yticklabelsize = 28
        obax.xlabelsize = 36
        obax.ylabelsize = 36
        obax.ylabelpadding = 20
    end

    # create listeners that update boundary map points
    all_bmap_scatters = [Point2f0[] for i in 1:N]
    for i in 1:N
        parobs = allparobs[i]
        vector = all_bmap_scatters[i]
        on(parobs.ξsin) do val
            push!(vector, val)
        end
    end

    # evolve particles until necessary time
    for j in 1:steps-1
        for i in 1:N
            p = ps[i]
            parobs = allparobs[i]
            animstep!(parobs, bd, dt, false, intervals)
            if plot_particles
                balls[][i] = parobs.p.pos
                vels[][i] = vr * parobs.p.vel
            end
        end
    end
    # last step (that actually updates plot)
    for i in 1:N
        p = ps[i]
        parobs = allparobs[i]
        animstep!(parobs, bd, dt, true)
    end
    if plot_particles
        balls[] = balls[]
        vels[] = vels[]
    end

    # Scatterplots:
    for i in 1:N
        c = cs[i]
        vector = all_bmap_scatters[i]
        scatter!(bmapax, vector; color = cs[i],
            marker = MARKER, markersize = ms*AbstractPlotting.px
        )
    end
    ylims!(bmapax, -1, 1)
    xlims!(bmapax, intervals[1], intervals[end])
    figure[:, 2] = sublayout
    return figure, bmapax
end


text:
@safetestset "HTML attributes rendering" begin
  @safetestset "No attributes" begin
    using Genie.Renderer.Html

    r = html("

")
    @test String(r.body) == "

"
  end;

  @safetestset "No attributes force parse" begin
    using Genie.Renderer.Html

    r = html("

", forceparse = true)
    @test String(r.body) == "<!DOCTYPE html>

"
  end;

  @safetestset "Regular attribute" begin
    using Genie.Renderer.Html

    r = html("""

""")
    @test String(r.body) == """

"""
  end;

  @safetestset "Regular attribute force parse" begin
    using Genie.Renderer.Html

    r = html("""

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

"""
  end;

  @safetestset "Dashed attributes" begin
    using Genie.Renderer.Html

    r = html("""

""")
    @test String(r.body) == """

"""
  end;

  @safetestset "Dashed attributes force parse" begin
    using Genie.Renderer.Html

    r = html("""

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

"""
  end;

  @safetestset "Multiple dashed attributes" begin
    using Genie.Renderer.Html

    r = html("""

""")
    @test String(r.body) == """

"""
  end;

  @safetestset "Multiple dashed attributes force parse" begin
    using Genie.Renderer.Html

    r = html("""

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

"""
  end;

  @safetestset "Single quotes" begin
    using Genie.Renderer.Html

    r = html("

")
    @test String(r.body) == """

"""
  end;

  @safetestset "Single quotes force parse" begin
    using Genie.Renderer.Html

    r = html("

", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

"""
  end;

  @safetestset "Vue args force parse" begin
    using Genie
    using Genie.Renderer.Html

    r = html("""
    Hover your mouse over me for a few seconds
    to see my dynamically bound title!
  """, forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>  Hover your mouse over me for a few seconds
    to see my dynamically bound title!
  """

    r = html("""

  Now you see me

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>
Now you see me
"""

    r = html("""

        {{ todo.text }}

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

  {{ todo.text }}

"""

    r = html("""

{{ message }}

    <button v-on:click="reverseMessage">Reverse Message</button>

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

{{ message }}
<button v-on:click="reverseMessage">Reverse Message</button>
"""

    r = html("""

{{ message }}

    <input v-model="message">

""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>

{{ message }}
<input v-model="message"$(Genie.config.html_parser_close_tag)>
"""

    Genie.Renderer.Html.register_element("todo-item")

    r = html("""

    <todo-item></todo-item>

""", forceparse = true)
    @test_broken String(r.body) == """<!DOCTYPE html>
<todo-item></todo-item>
"""

    r = html("""

      <todo-item
        v-for="item in groceryList"
        v-bind:todo="item"
        v-bind:key="item.id"
      ></todo-item>

""", forceparse = true)

    @test_broken String(r.body) == """<!DOCTYPE html>

<todo-item v-for="item in groceryList" v-bind:todo="item" v-bind:key="item.id"></todo-item>

"""

    r = html("""""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>"""

    r = html("""""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>"""

    r = html("""""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>"""

    r = html("""""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>"""
  end;

  @safetestset "Embedded Julia" begin
    using Genie
    using Genie.Renderer.Html

    id = 10
    r = html(raw"""""", id = 10)
    @test String(r.body) == """<!DOCTYPE html>"""

    r = html(raw"""""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>"""

    r = html("""""", forceparse = true)
    @test String(r.body) == """<!DOCTYPE html>"""

    r = html("""""", forceparse = true)
    @test_skip String(r.body) == """<!DOCTYPE html>"""
  end;
end;


text:
using TimeZones: Transition, TZFILE_MAX

abbrs = b"LMT\0WSST\0SDT\0WSDT\0"  # Pacific/Apia
@test TimeZones.abbreviation(abbrs, 5) == "WSST"
@test TimeZones.abbreviation(abbrs, 6) == "SST"
@test TimeZones.abbreviation(abbrs, 10) == "SDT"
@test TimeZones.abbreviation(abbrs, 14) == "WSDT"

# Extracts Transitions such that the two arrays start and stop at the
# same DateTimes.
function overlap(a::Array{Transition}, b::Array{Transition})
    start_dt = max(first(a).utc_datetime, first(b).utc_datetime)
    end_dt = min(last(a).utc_datetime, last(b).utc_datetime)

    within = t -> start_dt <= t.utc_datetime <= end_dt
    return a[findall(within, a)], b[findall(within, b)]
end

function issimilar(x::Transition, y::Transition)
    x == y || x.utc_datetime == y.utc_datetime && x.zone.name == y.zone.name && isequal(x.zone.offset, y.zone.offset)
end

@test_throws AssertionError TimeZones.read_tzfile(IOBuffer(), "")

# Compare tzfile transitions with those we resolved directly from the Olson zones/rules

# Ensure that read_tzfile returns a FixedTimeZone with the right data
utc = FixedTimeZone("UTC", 0)
open(joinpath(TZFILE_DIR, "Etc", "UTC")) do f
    tz = TimeZones.read_tzfile(f, "UTC")
    @test tz == utc
end

# Fixed time zone using version 2 data.
utc_plus_6 = FixedTimeZone("UTC+6", 6 * 3600)
open(joinpath(TZFILE_DIR, "Etc", "GMT-6")) do f
    tz = TimeZones.read_tzfile(f, "UTC+6")
    @test tz == utc_plus_6
end

warsaw = first(compile("Europe/Warsaw", tzdata["europe"]))
open(joinpath(TZFILE_DIR, "Europe", "Warsaw")) do f
    tz = TimeZones.read_tzfile(f, "Europe/Warsaw")
    @test string(tz) == "Europe/Warsaw"
    @test first(tz.transitions).utc_datetime == DateTime(1915,8,4,22,36)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX
    @test ==(overlap(tz.transitions, warsaw.transitions)...)
end

# Read version 1 compatible data
open(joinpath(TZFILE_DIR, "Europe", "Warsaw (Version 2)")) do f
    version, tz = TimeZones._read_tzfile(f, "Europe/Warsaw")
    @test version == '2'
    @test string(tz) == "Europe/Warsaw"
    @test first(tz.transitions).utc_datetime == typemin(DateTime)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX

    # File skips 1879-12-31T22:36:00
    @test ==(overlap(tz.transitions, warsaw.transitions[3:end])...)
end

# Read version 2 data
open(joinpath(TZFILE_DIR, "Europe", "Warsaw (Version 2)")) do f
    tz = TimeZones.read_tzfile(f, "Europe/Warsaw")
    @test string(tz) == "Europe/Warsaw"
    @test first(tz.transitions).utc_datetime == typemin(DateTime)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX
    @test ==(overlap(tz.transitions, warsaw.transitions)...)
end

godthab = first(compile("America/Godthab", tzdata["europe"]))
open(joinpath(TZFILE_DIR, "America", "Godthab")) do f
    tz = TimeZones.read_tzfile(f, "America/Godthab")
    @test string(tz) == "America/Godthab"
    @test first(tz.transitions).utc_datetime == DateTime(1916,7,28,3,26,56)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX
    @test ==(overlap(tz.transitions, godthab.transitions)...)
end

# Read version 1 compatible data
open(joinpath(TZFILE_DIR, "America", "Godthab (Version 3)")) do f
    version, tz = TimeZones._read_tzfile(f, "America/Godthab")
    @test version == '3'
    @test string(tz) == "America/Godthab"
    @test first(tz.transitions).utc_datetime == typemin(DateTime)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX
    @test ==(overlap(tz.transitions, godthab.transitions)...)
end

# Read version 3 data
open(joinpath(TZFILE_DIR, "America", "Godthab (Version 3)")) do f
    tz = TimeZones.read_tzfile(f, "America/Godthab")
    @test string(tz) == "America/Godthab"
    @test first(tz.transitions).utc_datetime == typemin(DateTime)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX
    @test ==(overlap(tz.transitions, godthab.transitions)...)
end

# "Pacific/Apia" was the time zone I was thinking could be an issue for the
# DST calculation. The entire day of 2011/12/30 was skipped when they changed from a
# -11:00 GMT offset to 13:00 GMT offset
apia = first(compile("Pacific/Apia", tzdata["australasia"]))
open(joinpath(TZFILE_DIR, "Pacific", "Apia")) do f
    tz = TimeZones.read_tzfile(f, "Pacific/Apia")
    @test string(tz) == "Pacific/Apia"
    @test first(tz.transitions).utc_datetime == DateTime(1911,1,1,11,26,56)
    @test last(tz.transitions).utc_datetime == DateTime(2037,9,26,14)
    @test tz.cutoff == TZFILE_MAX
    @test ==(overlap(tz.transitions, apia.transitions)...)
end

# Because read_tzfile files only store a single offset if both utc and dst change at the same
# time then the resulting utc and dst might not be quite right. Most notably during
# midsomer back in 1940's there were 2 different dst one after another, we get a
# different utc and dst than Olson.
paris = first(compile("Europe/Paris", tzdata["europe"]))
open(joinpath(TZFILE_DIR, "Europe", "Paris")) do f
    tz = TimeZones.read_tzfile(f, "Europe/Paris")
    @test string(tz) == "Europe/Paris"
    @test first(tz.transitions).utc_datetime == DateTime(1911,3,10,23,51,39)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX

    tz_transitions, paris_transitions = overlap(tz.transitions, paris.transitions)

    # Indices 56:2:58 don't match due to issues with Midsummer time.
    mask = tz_transitions .== paris_transitions
    neg_mask = map(!, mask)
    @test sum(neg_mask) == 2
    @test tz_transitions[mask] == paris_transitions[mask]
    @test all(map(issimilar, tz_transitions[neg_mask], paris_transitions[neg_mask]))
end

madrid = first(compile("Europe/Madrid", tzdata["europe"]))
open(joinpath(TZFILE_DIR, "Europe", "Madrid")) do f
    tz = TimeZones.read_tzfile(f, "Europe/Madrid")
    @test string(tz) == "Europe/Madrid"
    @test first(tz.transitions).utc_datetime == DateTime(1917,5,5,23)
    @test last(tz.transitions).utc_datetime == DateTime(2037,10,25,1)
    @test tz.cutoff == TZFILE_MAX

    tz_transitions, madrid_transitions = overlap(tz.transitions, madrid.transitions)

    # Indices 24:2:32 don't match due to issues with Midsummer time.
    mask = tz_transitions .== madrid_transitions
    neg_mask = map(!, mask)
    @test sum(neg_mask) == 5
    @test tz_transitions[mask] == madrid_transitions[mask]
    @test all(map(issimilar, tz_transitions[neg_mask], madrid_transitions[neg_mask]))
end

# "Australia/Perth" test processing a tzfile that should not contain a cutoff
perth = first(compile("Australia/Perth", tzdata["australasia"]))
open(joinpath(TZFILE_DIR, "Australia", "Perth")) do f
    tz = TimeZones.read_tzfile(f, "Australia/Perth")
    @test string(tz) == "Australia/Perth"
    @test first(tz.transitions).utc_datetime == DateTime(1916,12,31,16,1)
    @test last(tz.transitions).utc_datetime == DateTime(2009,3,28,18)
    @test tz.cutoff === nothing

    tz_transitions, perth_transitions = overlap(tz.transitions, perth.transitions)

    # Index 1 doesn't match up
    mask = tz_transitions .== perth_transitions
    neg_mask = map(!, mask)
    @test sum(neg_mask) == 1
    @test tz_transitions[mask] == perth_transitions[mask]
    @test all(map(issimilar, tz_transitions[neg_mask], perth_transitions[neg_mask]))
end


text:
"""
An `efficientIV` model consists of outcomes, `y` ∈ (0,1),  regressors
`x` and instruments `z`.  The moment condition is
``E[ g(z)(y-x β) ] = 0``

The dimensions of `x`, `y`, and `z` must be such that
`length(y) == size(x,1) == size(z,1)`.
In this version, `size(x,2) = size(z,2)=1`.
"""
mutable struct efficientIV <: GMMModel
  x::Vector{Float64}
  y::Vector{Float64}
  z::Vector{Float64}
end

"""
efficientIV(n::Integer, β::Real, ρ::Real, π::Real, a::AbstractVector)
Simulate an efficientIV model.

# Arguments
- `n` number of observations
- `β` coefficients on `x`
- `π` first stage coefficients `x = z*π + v`
- `ρ` correlation betwwen `v` and `ϵ`
- `a` values that z take
Returns an efficientIV model.
"""

function efficientIV(n::Integer, β::Real, ρ::Real, π::Real, a::AbstractVector)
    # (a) generate z
    w=rand(n)
    w1=(w.<0.2)*1.0
    w2=(w.>=0.2).*(w.<0.4)*1.0
    w3=(w.>=0.4).*(w.<0.6)*1.0
    w4=(w.>=0.6)*1.0
    z=a[1].*w1+a[2].*w2+a[3].*w3+a[4].*w4

    # (b) Matrix of Errors
    A=[1 ρ; ρ 1]
    ev=randn(n,2)*sqrt(A)
    one=ones(n)
    u=(one+z).*ev[:,1]

    # (c) Creating the X, Y vectors
    x = π.*z+ ev[:, 2]
    y = x.*β+u

    #Return the simulated data
    return (efficientIV(x,y,z))
end

number_parameters(model::efficientIV) = size(model.x,2)
number_observations(model::efficientIV) = length(model.y)
number_moments(model::efficientIV) = size(model.x,2)

function gz_hat(model::efficientIV)
    n=number_observations(model)
    β_2sls = (model.z'*model.x)^(-1)*model.z'*model.y
    uhat_2sls = model.y .- model.x .* β_2sls

    zunique=unique(model.z)
    xz4=ones((size(zunique)))
    u2z4=ones((size(zunique)))
    for i=1:size(zunique,1)
        xz4[i]=(sum(model.x.*(model.z.==zunique[i]),dims=1)/sum((model.z.==zunique[i]),dims=1))[1]
    end
    for i=1:size(zunique,1)
        u2z4[i]=(sum(uhat_2sls.*uhat_2sls.*(model.z.==zunique[i]),dims=1)/sum((model.z.==zunique[i]),dims=1))[1]
    end
    ratio4=xz4./u2z4
    gz_hat=ones(n)
    for i=1:n
            for j=1:size(zunique,1)
                if model.z[i]==zunique[j]
                    gz_hat[i]=ratio4[j]
                else
                    gz_hat[i]=gz_hat[i]
                end
            end
        end

    return gz_hat
end

function get_gi(model::efficientIV)
   function(β)
        gz_hat=gz_hat(model::efficientIV)
        ξ=y-model.x.*β
        m=ξ.*gz_hat
   end
end

function estimation(model::efficientIV)
    gzhat=gz_hat(model::efficientIV)
    β_gz_hat=(gzhat'*model.x)^(-1)*gzhat'*model.y
    return (β_gz_hat)
end


text:
include("centrality.jl")

using Statistics

function main()
    m = 3
    dim = 7
    T = SymTensor3([1, 1, 3, 5],
                   [2, 2, 5, 6],
                   [3, 4, 6, 7],
                   ones(Float64, 4),
                   dim)

    z, conv = Z_evec_dynsys(T)
    println("dynamical systems algorithm converged: $conv")
    z2 = z / norm(z, 2)
    λ1 = mean(apply(T, z2) ./ z2)
    println("tensor z-eigenvalue: $(λ1)")
    println("tensor z-eigenvector (2-norm normalized): $(z2)")

    # Check for stability
    Ustar = Array(qr((I - z2 * z2') * randn(dim, dim - 1)).Q)
    maxval = maximum(z2' * Ustar )
    println("orthogonality test (should be close to 0): $(maxval)")
    diff = Ustar' * Ustar - I
    diffnorm = [norm(vec(diff[:, ind]), 2) for ind in 1:size(diff)[2]]
    println("orthonormality test (should be close to 0): $(diffnorm)")

    C = Ustar' * ((m - 1) * reduce(T, z2) - λ1 * I) * Ustar
    U = triu(C, 1)  # for symmetrization
    C_sym = Matrix(U + U' + Diagonal(C))
    spectrum = eigen(C_sym).values
    println("Spectrum of projected Hessian of the Lagrangian: $spectrum")
end
;

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
#=
Description
    Module for UAV type objects.

    state: [ x,   xdot,       // x-direction
             y,   ydot,       // y-direction
             psi, psidot,     // heading
             phi, phidot  ]   // bank angle
    input:   phiddot          // bank angle 'acceleration'
=#

module UAVs

using ODE

export UAV, generate_trajectory, set_state!

const SPEED_DEF = 10.0
const MU = 0.0
const SIGMA = 0.05
const G = 9.8
const NSTATES = 8
const NINPUTS = 1

type UAV
    state       :: Vector{Float64}
    nStates     :: Int64
    nInputs     :: Int64
    speed       :: Float64

    function UAV(speed::Float64 = SPEED_DEF,
                 nStates::Int64 = NSTATES,
                 nInputs::Int64 = NINPUTS)
        return new(zeros(nStates), nStates, nInputs, speed)
    end # function UAV
end # type UAV

function dynamics(uav::UAV, policy::Function, nT::Int64, dT::Float64)
    V = uav.speed
    tStart = 0
    tEnd = tStart + dT * nT
    tVec = tStart:dT:tEnd

    function dxdt(t::Float64, state::Vector{Float64})
        dx = zeros(size(state))
        dx[1] = V * cos(state[5])                   # xdot
        dx[3] = V * sin(state[5])                   # ydot
        dx[5] = G * tan(state[7]) / V               # psidot
        dx[7] = state[8]                            # phidot
        dx[8] = policy(t, state)                    # phiddot
        dx[2] = -V * sin(state[5]) * dx[7]          # xddot
        dx[4] = V * cos(state[5]) * dx[7]           # yddot
        dx[6] = G / V * sec(state[8])^2 * dx[8]     # psiddot
        return dx
    end # function dxdt

    X = zeros(uav.nStates, nT + 1)
    state = uav.state
    X[:, 1] = state
    for ti = 1:length(tVec) - 1
        tSpan = [tVec[ti], tVec[ti + 1]]
        t, xRaw = ode45(dxdt, state, tSpan) # non-stiff
        # t, xRaw = ode23s(dxdt, state, tSpan) # stiff
        state = xRaw[end]
        X[:, ti + 1] = state
    end # for ti
    return tVec, X
end # function discrete_dynamics

#=
Description
    Returns the time-state trajectory pair for the given input
    policy |policy|, number of time steps |nT| and discrete time
    period |dT|.
=#
function generate_trajectory(uav::UAV, policy::Function, nT::Int64, dT::Float64)
    if uav.nInputs != size(policy(0.0, uav.state), 1)
        error("Dimension mismatch between uav.nInputs and policy output")
    end # if
    T, X = dynamics(uav, policy, nT, dT)
    return T, X
end # function generate_trajectory

function set_state!(uav::UAV, newState::Vector{Float64})
    if uav.nStates == size(newState, 1)
        uav.state = newState
    else
        error("Dimension mismatch between uav.nStates and newState")
    end # if
end # function set_state!

end # module UAVs


text:
import Base.rand
import Distributions, StatsBase
import IntervalArithmetic: interval
import Statistics

abstract type TheoreticalDistributionScalarValue <: AbstractUncertainValue end

Base.rand(uv::TheoreticalDistributionScalarValue) = rand(uv.distribution)
Base.rand(uv::TheoreticalDistributionScalarValue, n::Int) = rand(uv.distribution, n)

function Distributions.support(fd::TheoreticalDistributionScalarValue)
    s = support(fd.distribution)
    interval(s.lb, s.ub)
end

Distributions.pdf(fd::TheoreticalDistributionScalarValue, x) = pdf(fd.distribution, x)
Statistics.mean(fd::TheoreticalDistributionScalarValue) = mean(fd.distribution)
Statistics.median(fd::TheoreticalDistributionScalarValue) = median(fd.distribution)
Statistics.middle(fd::TheoreticalDistributionScalarValue) = middle(fd.distribution)
Statistics.quantile(fd::TheoreticalDistributionScalarValue, q) = quantile(fd.distribution, q)
Statistics.std(fd::TheoreticalDistributionScalarValue) = std(fd.distribution)
Statistics.var(fd::TheoreticalDistributionScalarValue) = var(fd.distribution)
StatsBase.mode(fd::TheoreticalDistributionScalarValue) = mode(fd.distribution)

abstract type AbstractUncertainOneParameterScalarValue{S <: ValueSupport, T1} <: TheoreticalDistributionScalarValue end
abstract type AbstractUncertainTwoParameterScalarValue{S <: ValueSupport, T1, T2} <: TheoreticalDistributionScalarValue end
abstract type AbstractUncertainThreeParameterScalarValue{S <: ValueSupport, T1, T2, T3} <: TheoreticalDistributionScalarValue end


text:
using Test
using JVODE
using Sundials

prob = ODEProblem((du, u, p, t)->begin
                      du[1] = u[2]
                      du[2] = (1-u[1]^2)*3 * u[2] - u[1]
                  end, [2, 0.], (0.0, 1.39283880203), 0.9)

integ = @inferred init(prob, Adams(), reltol=0.0, abstol=1e-6)
step!(integ)
integs = init(prob, CVODE_Adams())
step!(integs)
@test integ.dt == integs.dt


text:
using Random

Random.rand!(A::oneWrappedArray) = Random.rand!(GPUArrays.default_rng(oneArray), A)
Random.randn!(A::oneWrappedArray) = Random.randn!(GPUArrays.default_rng(oneArray), A)


text:
using BrainFlow
using BrainFlowViz

function get_some_board_data(board_shim, nsamples, nchannels=4)
    data = BrainFlow.get_current_board_data(nsamples, board_shim)
    data = transpose(data)
    eeg_chans = BrainFlow.get_eeg_channels(board_shim.board_id)
    emg_data = view(data, :, eeg_chans)
    for chan in 1:nchannels
        emg_channel_data = view(emg_data, :, chan)
        BrainFlow.detrend(emg_channel_data, BrainFlow.CONSTANT)
    end
    return emg_data
end

### Start streaming
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)
params = BrainFlowInputParams(serial_port="COM5")
board_shim = BrainFlow.BoardShim(BrainFlow.MUSE_2_BLED_BOARD, params)
BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)

nchannels = length(BrainFlow.get_eeg_channels(board_shim.board_id))
nsamples = 512
data_func = ()->get_some_board_data(board_shim, nsamples, nchannels)

# start a live plotting task
BrainFlowViz.plot_data(
    data_func,
    nsamples,
    nchannels;
    #y_lim = [-200 200],
    theme = :dark,
    color = :lime,
    )

# this should go into a 'finally' of try/catch/finally
"""
BrainFlow.stop_stream(board_shim)
BrainFlow.release_session(board_shim)
"""


text:
using ConvexTests, ProxSDP

@info "Starting ProxSDP tests"
do_tests("ProxSDP", () -> ProxSDP.Optimizer(
    tol_gap = 5e-8, tol_feasibility = 1e-7, max_iter = 100_000_000, time_limit = 100.0
); exclude = [r"mip", r"exp"], description =
"Tests run with `tol_gap = 5e-8`, `tol_feasibility = 1e-7`, `max_iter = 100_000_000` and `time_limit = 100.0`.")

max_stars_repo_name: ericphanson/ConvexTests

max_issues_repo_name: ericphanson/ConvexTests

max_forks_repo_name: ericphanson/ConvexTests


text:
module Born

import GeoPhyInv: AGeom, SrcWav
import GeoPhyInv: Medium
import GeoPhyInv.Data
using SpecialFunctions
using FFTW
using DSP

mutable struct Param
end

"""
k = wavenumber = 2pif/v0
"""
function G0_homo_acou(x::T, z::T, k::T, rho0::T) where {T<:AbstractFloat}
	sqrt(x*x + z*z) == 0.0 ? error("distance cannot be zero") : nothing
	k == 0.0 ? error("wavenumber cannot be zero") : nothing
	G0 = -0.25 * rho0 * im * complex.(besselj0(k*sqrt(x*x + z*z)), -bessely0(k*sqrt(x*x + z*z)))
	return G0
end

function mod(;
	     vp0::Float64=2000.0,
	     rho0::Float64=2000.0,
	     medium_pert::Medium=nothing,
             born_flag::Bool=false,
	     tgridmod::StepRangeLen=nothing,
	     tgrid::StepRangeLen = tgridmod,
	     ageom::AGeom=nothing,
	     srcwav::SrcWav=nothing,
	     src_flag::Int64=2,
		  )
	(length(ageom) != length(srcwav))  ? error("different supersources") : nothing
	@info "fix this"
#	(length(ageom) != getfield(srcwav,:ns))  ? error("different sources") : nothing

	nt = (length(tgridmod) == length(srcwav[1].grid)) ? length(tgrid) : error("srcwav tgrid")
	np2 = nextpow(2, 2*nt);

	if(born_flag)
		mesh_x=medium_pert.mgrid[2]
		mesh_z=medium_pert.mgrid[1]
		nz=length(mesh_z)
		nx=length(mesh_x)
		δx = step(mesh_x)
		δz = step(mesh_z)

		δmodtt = medium_pert[:KI] - (vp0 * vp0 * rho0)^(-1)
		δmodrr = medium_pert[:rhoI] - (rho0)^(-1)
	end

	fnpow2grid = FFTW.fftfreq(np2,inv(step(tgridmod)));

	nss = length(ageom)
	data=Data(tgridmod, ageom, [:P])

	for iss in 1:nss, ifield in 1:length(data[iss].d)
		sx = ageom[iss].sx
		rx = ageom[iss].rx
		sz = ageom[iss].sz
		rz = ageom[iss].rz

		for ir = 1:ageom[iss].nr

		dtemp=zeros(nt)
		dpow2all=complex.(zeros(np2), zeros(np2));
		wpow2=complex.(zeros(np2), zeros(np2));

		for is=1:ageom[iss].ns
			# zero pad wavelet
			for it in 1:nt
				wpow2[it]=complex(srcwav[iss].d[ifield][it,is])
			end
			FFTW.fft!(wpow2) # source wavelet in the frequency domain

			x = sx[is] - rx[ir]
			z = sz[is] - rz[ir]
			# analytical expression for every frequency, except zero
			dpow2 = complex.(zeros(np2), zeros(np2));
			dpow2[1] = complex.(0.0, 0.0)
			for iω in 2:np2
				ω = 2. * pi * abs(fnpow2grid[iω])
				k = ω / vp0

				if(born_flag)
					term = complex(0., 0.)
					for ix=1:nx
						@simd for iz=1:nz
							if(δmodtt[iz,ix] ≠ 0.0)
								term += (G0_homo_acou(sx[is]-mesh_x[ix], sz[is]-mesh_z[iz], k, rho0)[1]
									 * G0_homo_acou(rx[ir]-mesh_x[ix], rz[ir]-mesh_z[iz], k, rho0)[1] .* ω .* ω .* δmodtt[iz,ix]) * δx * δz
								# factor due to integration
							end
						end
					end

				else
					if(src_flag == 2)
						term = G0_homo_acou(x, z, k, rho0);
					elseif(src_flag == 1)
						dpow2[iω] = G0_homo_acou(x, z, k, rho0) * im * abs(fnpow2grid[iω])
					else
						error("invalid src_flag")
					end
				end
				if(fnpow2grid[iω] > 0)
					dpow2[iω] = term;
				elseif(fnpow2grid[iω] < 0)
					dpow2[iω] = conj(term);
				end

			end
			# convolution and stack over simultaneous sources
			dpow2all += dpow2 .* wpow2;
		end

		# back to time domain
		FFTW.ifft!(dpow2all)

		# truncate
		for it in 1:nt
			data[iss].d[ifield][it,ir] = real(dpow2all[it])
		end
		end
	end
	return data
end

end # module