text:
using DifferentialEquations, LinearAlgebra, Printf;

print("FLINT Performance comparison with Julia DiffEq Package\n")
print("\n")

###### @show rtol = 1.0./10.0.^(6:1:13);
@show rtol = 1.0./10.0.^(11);
@show atol = rtol*1e-3;

function cr3bpeom(dX, X, p, t)
    @inbounds begin
        r1 = sqrt((X[1] + mu)^2 + X[2]^2)
        r2 = sqrt((X[1] - 1 + mu)^2 + X[2]^2)

        # @view makes Julia slower!
        #dX[1:3] .= @view X[4:6]
        dX[1] = X[4]
        dX[2] = X[5]
        dX[3] = X[6]
        dX[4] = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX[5] = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3
        dX[6] = 0.0
        #@SVector [X[4],X[5],X[6],dX4,dX5,0.0]
    end
end

function JacobiConstant(odesol)

    X = odesol[1:3,:]
    Xdot = odesol[4:6,:]

    ~,n = size(X)

    # distances
    R1 = copy(X)
    R2 = copy(X)
    R1[1,:] = R1[1,:] .+ mu
    R2[1,:] = R2[1,:] .- 1 .+ mu
    r1 = [norm(R1[:,i]) for i = 1:n]
    r2 = [norm(R2[:,i]) for i = 1:n]

    # radius and velocity nagnitudes
    R =  [sqrt(sum((X.^2)[:,i])) for i = 1:n]
    V = [sqrt(sum((Xdot.^2)[:,i])) for i = 1:n]

    # pseudo-potential
    Omega = 1/2*(R.^2) .+ (1-mu)./r1 .+ mu./r2

    # Jacobi's Constant
    C = (V.^2)/2 .- Omega

    return C
end;

# Event-1: detect a narrow square pulse
# Note that discrete callback wont be able to precisely locate the
# 0-crossings for this scenario
condition1 = function (u, t, integrator)
    if t > 0.5 && t < 0.5005
       1
    else
       -1
    end
end

function affect1(integrator)
    # TBD: how to disable this callback efficiently?
end

function affect1_neg(integrator)
    # TBD: how to disable this callback efficiently?
end

cb1 = ContinuousCallback(condition1, affect1, affect1_neg;
            rootfind=true,save_positions=(true,true),abstol=0.001)

# Event function-2: continuous callback
Mask2 = false
condition2 = function (u, t, integrator)
    if Mask2 == false
        if u[1] < 0
            # TBD: how to safely ignore the event in this region?
            -1
        else
            u[2]
        end
    else
        -1
    end
end

function affect2(integrator)
    # TBD: how to remove the callback?
    Mask2 = true
end

cb2 = ContinuousCallback(condition2, affect2;
            rootfind=true,save_positions=(true,true),abstol=1e-9)

# Event function-3: continuous callback
Mask3 = false
condition3 = function (u, t, integrator)
    if Mask3 == false
        if u[1] > 0
            # TBD: how to safely ignore the event in this region?
            -1
        else
            u[2]
        end
    else
        -1
    end
end

function affect3(integrator)
    # modify solution to set it back to IC
    integrator.u[1:3] = R0
    integrator.u[4:6] = V0
    # TBD: how to remove the callback?
    Mask3 = true
end

cb3 = ContinuousCallback(condition3, affect3;
            rootfind=true,save_positions=(true,true),abstol=1e-9)

# callback set
cb = CallbackSet(cb1,cb2,cb3)

# The 2nd element of each solver sets the dense=true option
# The 3rd element is the name used in results

struct Solver
    solver
    dense::Bool
    name::String
end

ODESolvers = [
                #Solver(dopri5(),false,"DOPRI5"),                              # Hairer's DOPRI5
                Solver(DP5(),true,"DP5"),                                     # Julia's implementation of DOPRI5
                Solver(Tsit5(),true,"TSIT5"),                                 # Tsitouras 5/4 RK
                #Solver(dop853(),false,"DOP853"),                              # Hairer's DOP853
                Solver(DP8(),true,"DP8"),                                     # Julia's implementation of DOP853
                #Solver(ARKODE(Sundials.Explicit(), order = 8),true,"ARKODE8"),# Sundial's explicit RK
                Solver(Vern6(),true,"Vern6"),                                 # Verner's 7/6 RK
                #Solver(Vern7(),true,"Vern7"),                                 # Verner's 7/6 RK
                #Solver(Vern8(),true,"Vern8"),                                 # Verner's 7/6 RK
                Solver(Vern9(),true,"Vern9"),                                 # Verner's 9/8 RK
                #Solver(ddeabm(),false,"DDEABM"),                              # Shampine's DDEABM
                #Solver(VCABM(),true,"JABM"),      # forcing dt_min           # Julia's implementation of Shampine's DDEABM
                #Solver(CVODE_Adams(max_order = 12),true,"SunAM12") # Too much error!   # Sundial's Adams-Moulton solver
            ]

# test result Struct
struct ODETestResult
    FinalStateErr
    IOMErr
    MeanTime
    OdeSol
end

# Integrate ODE
function FireODE(atol, rtol, method, prob, DenseOn)
 if DenseOn == true
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false, dense=false,callback=cb)
    else
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false,callback=cb)
    end
end

# Main ODE Integration Method Test function
function FireODEIntTest(X0, ODE, tspan, OdeMethod,  atol, rtol, FinalState, IOM, IOMFunc::Function, DenseOn)

    n = length(FinalState)

    # run solver bunch of times for benchmarking
    sol = 0
    rtime = 0.0
    X0new = copy(X0)

    rtime = @elapsed for i in 1:Runs
        # randomsize initial conditions
        X0new[1] = X0new[1] + 0.000000000001*X0new[1].*rand()

        # Define ODE problem
        prob = ODEProblem(ODE, X0new, tspan);

        #GC.gc()
        sol = FireODE(atol, rtol, OdeMethod, prob, DenseOn)
    end

    # final state error metric
    serr = norm(sol[1:3,end] - FinalState[1:3])

    # compute IOM variations
    IOMsol = IOMFunc(sol)
    ierr = abs.(IOMsol .- IOM)

    ODETestResult(serr, ierr, rtime, sol)
end

const Runs = 100

# mass-ratio
const mu = 0.012277471::Float64

# Initial conditions
R0 = [0.994::Float64, 0.0, 0.0]
V0 = [0.0, parse(Float64,"-2.00158510637908252240537862224"), 0.0]

# Time period
Period = parse(Float64,"17.0652165601579625588917206249")

# propagate for bunch of orbits
t0 = 0.0::Float64;
tf = 2.5*Period;

tspan = (t0, tf);

X0 = [R0; V0]

op1 = ODEProblem(cr3bpeom,X0,tspan);

nsol = length(ODESolvers)
ntol = length(rtol)

SolResults3B = Matrix(undef, nsol, ntol)

for itr in 1:ntol
    for ctr in 1:nsol
        print(string("Running: ",ODESolvers[ctr].name, "\n"))
        SolResults3B[ctr,itr] = FireODEIntTest(X0, cr3bpeom, tspan, ODESolvers[ctr].solver,  atol[itr], rtol[itr], [R0; V0], 0, JacobiConstant, ODESolvers[ctr].dense);
    end
    print(string("tol: ", rtol[itr],"\n"))
end

# Extract FLINT computation results from results.txt file
cd(@__DIR__)
DOP54_t, DOP853_t,Vern65E_t,Vern98R_t = open("results.txt") do f
    dop54 = 0.0
    dop853 = 0.0
    vern65e = 0.0
    vern98r = 0.0
    for l in eachline(f)
        if occursin("DOP54", l) == true
            strs = split(l,' ',keepempty=false)
            dop54 = parse(Float32,strs[2])
        end
        if occursin("DOP853", l) == true
            strs = split(l,' ',keepempty=false)
            dop853 = parse(Float32,strs[2])
        end
        if occursin("Verner65E", l) == true
            strs = split(l,' ',keepempty=false)
            vern65e = parse(Float32,strs[2])
        end
        if occursin("Verner98R", l) == true
            strs = split(l,' ',keepempty=false)
            vern98r = parse(Float32,strs[2])
        end
        if occursin("Interpolated", l) == true
            break
        end
    end
    (dop54,dop853,vern65e,vern98r)
end

FLINT_res = [DOP54_t,DOP54_t,DOP853_t,Vern65E_t,Vern98R_t]
FLINT_alg = ["DOP54","DOP54","DOP853","Vern65E","Vern98R"]

print("\n")
print("Solver: Julia time, FLINT time, Julia closing errors\n==============================================================\n")
for ctr in 1:nsol
    @Printf.printf("%7s: %.5e, %.5e (FLINT %7s), %.5e\n",
    ODESolvers[ctr].name,SolResults3B[ctr,end].MeanTime,FLINT_res[ctr],FLINT_alg[ctr],SolResults3B[ctr,end].FinalStateErr)
end

max_stars_repo_name: jacobwilliams/FLINT

max_issues_repo_name: jacobwilliams/FLINT

max_forks_repo_name: jacobwilliams/FLINT


text:
struct QCurlGradMonomialBasis{P,V,G,D} <: TensorPolynomialBasis{P,V,G}
  order::Int
  terms::CartesianIndices{D}
  perms::Matrix{Int}
end

function (::Type{QCurlGradMonomialBasis{P,V}})(order::Int) where {P,V}
  D = _length(P)
  @assert _length(V) == D
  G = _gradient_type(V,P)
  n1d = order
  t = fill(n1d,D)
  t[1] = order+1
  terms = CartesianIndices(tuple(t...))
  perms = _prepare_perms(D)
  QCurlGradMonomialBasis{P,V,G,D}(order,terms,perms)
end

# Implementation of the interface

length(b::QCurlGradMonomialBasis{P,V,G,D}) where {P,V,G,D} = D*(b.order+1)*(b.order)^(D-1)

ndims(b::QCurlGradMonomialBasis{P,V,G,D}) where {P,V,G,D} = D

function ScratchData(b::QCurlGradMonomialBasis{P,V}) where {P,V}
  T = eltype(V)
  dim = _length(P)
  n1d = b.order+1
  c = zeros(dim,n1d)
  g = zeros(dim,n1d)
  MonomialBasisCache{T}(c,g)
end

function evaluate!(
  v::AbstractVector{V},
  b::QCurlGradMonomialBasis{P,V,G,D},
  x::P,
  cache::MonomialBasisCache) where {P,V,G,D}
  _evaluate_nd_qgrad!(v,x,b.order,b.terms,b.perms,cache.c)
end

function gradient!(
  v::AbstractVector{G},
  b::QCurlGradMonomialBasis{P,V,G,D},
  x::P,
  cache::MonomialBasisCache) where {P,V,G,D}
  _gradient_nd_qgrad!(v,x,b.order,b.terms,b.perms,cache.c,cache.g,V)
end

ext: jl


text:
"""
    UncertainValueDataset

A dataset of uncertain values.

## Fields

- **`values::AbstractVector{<:AbstractUncertainValue}`**: The uncertain values. Each value is
    represented by an `AbstractUncertainValue`.
"""
struct UncertainValueDataset <: AbstractUncertainValueDataset
    values::AbstractVector{<:AbstractUncertainValue}
end

"""
    ConstrainedUncertainValueDataset

Generic constrained dataset containing uncertain values.

## Fields
- **`values::AbstractVector{<:AbstractUncertainValue}`**: The uncertain values.
"""
struct ConstrainedUncertainValueDataset <: AbstractUncertainValueDataset
    values::AbstractVector{<:AbstractUncertainValue}
end

function UncertainValueDataset(x::AbstractArray{T, 1}) where T
    UncertainValueDataset(CertainValue.(x))
end

export
UncertainValueDataset,
ConstrainedUncertainValueDataset


text:
"""
Implements:
#     r = x-z
#
#     factor = -2*a*exp(-a*dot(r,r))
#
#     return r .* factor
"""
function eval∂SqExp!(   out::Vector{T},
                        x::Vector{T},
                        z::Vector{T},
                        a::T) where T
    #
    @assert length(out) == length(x)
    for d = 1:length(x)
        out[d] = x[d]-z[d]
    end

    factor = -2*a*exp(-a*dot(out,out))

    for d = 1:length(out)
        out[d] = out[d] * factor
    end

    return nothing
end

function eval∂SqExp(x::Vector{T},
                z::Vector{T},
                a::T) where T
    #
    out = Vector{T}(undef, length(x))

    eval∂SqExp!(out, x, z, a)

    return out
end

"""
Implements:
#     r = x-z
#
#     return collect( 2*a*(2*a*r[d]^2-1)*exp(-a*dot(r,r)) for d = 1:length(x) )
# Wolfram Alpha: d^2/dx^2 (exp(-a (x^2 + y^2)))
"""
function eval∂2SqExp!(   out::Vector{T},
                        x::Vector{T},
                        z::Vector{T},
                        a::T)::Nothing where T
    #
    @assert length(out) == length(x)
    for d = 1:length(x)
        out[d] = x[d]-z[d]
    end

    factor = 2*a*exp(-a*dot(out,out))

    for d = 1:length(out)
        out[d] = (2*a*out[d]^2-1)  * factor
    end

    return nothing
end

function eval∂2SqExp(x::Vector{T},
                z::Vector{T},
                a::T) where T
    #
    out = Vector{T}(undef, length(x))

    eval∂2SqExp!(out, x, z, a)

    return out
end

function queryd1(x, c, X, a::T)::Vector{T} where T

    out = zeros(T,length(x))

    queryd1!(out, x, c, X, a)

    return out
end

function queryd1!(out::Vector{T}, x, c, X, a::T)::Nothing where T

    N = length(X)
    @assert length(c) == length(X)

    D = length(x)

    # out = zeros(T, D)
    # for n = 1:N
    #     out += c[n] .* eval∂SqExp(x, X[n], a)
    # end

    tmp = Vector{T}(undef, D)

    resize!(out, D)
    fill!(out, zero(T))
    for n = 1:N
        eval∂SqExp!(tmp, x, X[n], a)

        for d = 1:D
            out[d] += c[n]*tmp[d]
        end
    end

    return nothing
end

function queryd2(x, c, X, a::T)::Vector{T} where T

    out = zeros(T,length(x))

    queryd2!(out, x, c, X, a)

    return out
end

function queryd2!(out::Vector{T}, x, c, X, a::T)::Nothing where T
    N = length(X)
    @assert length(c) == length(X)

    D = length(x)

    # out = zeros(T, D)
    # for n = 1:N
    #     out += c[n] .* eval∂2SqExp(x, X[n], a)
    # end

    tmp = Vector{T}(undef, D)

    #out = zeros(T, D)
    resize!(out, D)
    fill!(out, zero(T))
    for n = 1:N
        eval∂2SqExp!(tmp, x, X[n], a)

        for d = 1:D
            out[d] += c[n]*tmp[d]
        end
    end

    return nothing
end

function querySqExp(x,
                    c::Vector{T},
                    X::Vector{Vector{T}},
                    ϵ_sq::T)::T where T
    #
    D = length(x)
    r = Vector{T}(undef,D)

    out = zero(T)
    for n = 1:length(X)
        for d = 1:D
            r[d] = x[d]-X[n][d]
        end

        out += c[n]*exp(-ϵ_sq*dot(r,r))
    end

    return out
end

##### others.
function sum∂2SqExp(x::Vector{T},
                z::Vector{T},
                a::T)::T where T
    #
    r = x-z

    out = zero(T)
    for d = 1:length(x)
        tmp = 2*a*(2*a*r[d]^2-1)

        out += tmp
    end
    #out = sqrt(out)*exp(-a*dot(r,r))
    out = out*exp(-a*dot(r,r))

    return out
end

function sum∂4SqExp(x::Vector{T},
                z::Vector{T},
                a::T)::T where T
    #
    r = x-z
    out = zero(T)
    for d = 1:length(x)
        tmp = 4*a^2*x[d]^4 - 12*a*x[d]^2 + 3

        out += tmp
    end
    #out = sqrt(out)*4*a^2*exp(-a*dot(r,r))
    out = out*4*a^2*exp(-a*dot(r,r))

    return out
end

# does not reset!
function add∂2diagSqExp!(out::Vector{T},
                w::T,
                x::Vector{T},
                z::Vector{T},
                a::T) where T
    #
    @assert length(out) == length(x) == length(z)
    r = x-z
    factor = w*exp(-a*dot(r,r))

    for d = 1:length(x)
        tmp = 2*a*(2*a*r[d]^2-1)

        out[d] += tmp*factor
    end

    return out
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
# these are all general purpose functions that aren't specifically related to
# radial velocity or GP calculations
using LinearAlgebra
using Random
using Statistics
using Unitful

"""
Nyquist frequency is half of the sampling rate of a discrete signal processing system
(https://en.wikipedia.org/wiki/Nyquist_frequency)
divide by another factor of 4 for uneven spacing
"""
nyquist_frequency(time_span::Union{Real,Quantity}, n_meas::Integer; nyquist_factor::Real=1) = n_meas / time_span / 2 * nyquist_factor
function nyquist_frequency(times::Vector{T}; nyquist_factor::Real=1) where {T<:Union{Real,Quantity}}
    time_span = times[end] - times[1]
    return nyquist_frequency(time_span, length(times), nyquist_factor=nyquist_factor)
end
uneven_nyquist_frequency(times; nyquist_factor=5) = nyquist_frequency(times; nyquist_factor=nyquist_factor)

"""
shamelessly crimped from JuliaAstro.jl
used to calculate range of frequencies to look at in a periodogram
"""
function autofrequency(times::Vector{T} where {T<:Union{Real,Quantity}};
                       samples_per_peak::Integer=5,
                       nyquist_factor::Integer=5,
                       minimum_frequency::Real=NaN,
                       maximum_frequency::Real=NaN)
    time_span = maximum(times) - minimum(times)
    δf = inv(samples_per_peak * time_span)
    f_min = isfinite(minimum_frequency) ? minimum_frequency : (δf / 2)
    if isfinite(maximum_frequency)
        return f_min:δf:maximum_frequency
    else
        return f_min:δf:nyquist_frequency(time_span, length(times); nyquist_factor=nyquist_factor)
    end
end

"assert all passed variables are positive"
function assert_positive(vars...)
    for i in vars
        @assert all(ustrip.(i) .> 0) "passed a negative/0 variable that needs to be positive"
    end
end

import Base.ndims
ndims(x::Cholesky) = 2

"""
Solve a linear system of equations (optionally with variance values at each point or covariance array)
see (https://en.wikipedia.org/wiki/Generalized_least_squares#Method_outline)
"""
function general_lst_sq(
    dm::Matrix{T},
    data::Vector;
    Σ::Union{Cholesky,Diagonal}=Diagonal(ones(length(data))),
    return_ϵ_inv::Bool=false) where {T<:Real}
    @assert ndims(Σ) < 3 "the Σ variable needs to be a 1D or 2D array"

    if return_ϵ_inv
        ϵ_int = Σ \ dm
        ϵ_inv = dm' * ϵ_int
        return ϵ_inv \ (dm' * (Σ \ data)), ϵ_int, ϵ_inv
    else
        return (dm' * (Σ \ dm)) \ (dm' * (Σ \ data))
    end
    # end
end

"Return an amount of indices of local maxima of a data array"
function find_modes(data::Vector{T}; amount::Integer=3) where {T<:Real}

    # creating index list for inds at modes
    mode_inds = [i for i in 2:(length(data)-1) if (data[i]>=data[i-1]) && (data[i]>=data[i+1])]
    if data[1] > data[2]; prepend!(mode_inds, 1) end
    if data[end] > data[end-1]; append!(mode_inds, length(data)) end

    # return highest mode indices
    return mode_inds[partialsortperm(-data[mode_inds], 1:amount)]

end

centered_rand(; rng::AbstractRNG=Random.GLOBAL_RNG, kwargs...) = centered_rand(rng; kwargs...)
centered_rand(rng::AbstractRNG; center::Real=0, scale::Real=1) = (scale * (rand(rng) - 0.5)) + center
centered_rand(d::Integer; rng::AbstractRNG=Random.GLOBAL_RNG, kwargs...) = centered_rand(rng, d; kwargs...)
centered_rand(rng::AbstractRNG, d; center::Real=0, scale::Real=1) = (scale .* (rand(rng, d) .- 0.5)) .+ center

function remove_mean!(vec::AbstractVector)
    vec .-= mean(vec)
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Transformer2W <: ACBranch
        name::String
        available::Bool
        activepower_flow::Float64
        reactivepower_flow::Float64
        arc::Arc
        r::Float64
        x::Float64
        primaryshunt::Float64
        rate::Union{Nothing, Float64}
        _forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

The 2-W transformer model uses an equivalent circuit assuming the impedance is on the High Voltage Side of the transformer. The model allocates the iron losses and magnetizing susceptance to the primary side.

# Arguments
- `name::String`
- `available::Bool`
- `activepower_flow::Float64`
- `reactivepower_flow::Float64`
- `arc::Arc`
- `r::Float64`: System per-unit value
- `x::Float64`: System per-unit value
- `primaryshunt::Float64`: System per-unit value
- `rate::Union{Nothing, Float64}`
- `_forecasts::InfrastructureSystems.Forecasts`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Transformer2W <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    "System per-unit value"
    r::Float64
    "System per-unit value"
    x::Float64
    "System per-unit value"
    primaryshunt::Float64
    rate::Union{Nothing, Float64}
    _forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Transformer2W(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, _forecasts=InfrastructureSystems.Forecasts(), )
    Transformer2W(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, _forecasts, InfrastructureSystemsInternal())
end

function Transformer2W(; name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, _forecasts=InfrastructureSystems.Forecasts(), )
    Transformer2W(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function Transformer2W(::Nothing)
    Transformer2W(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        primaryshunt=0.0,
        rate=nothing,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get Transformer2W name."""
get_name(value::Transformer2W) = value.name
"""Get Transformer2W available."""
get_available(value::Transformer2W) = value.available
"""Get Transformer2W activepower_flow."""
get_activepower_flow(value::Transformer2W) = value.activepower_flow
"""Get Transformer2W reactivepower_flow."""
get_reactivepower_flow(value::Transformer2W) = value.reactivepower_flow
"""Get Transformer2W arc."""
get_arc(value::Transformer2W) = value.arc
"""Get Transformer2W r."""
get_r(value::Transformer2W) = value.r
"""Get Transformer2W x."""
get_x(value::Transformer2W) = value.x
"""Get Transformer2W primaryshunt."""
get_primaryshunt(value::Transformer2W) = value.primaryshunt
"""Get Transformer2W rate."""
get_rate(value::Transformer2W) = value.rate
"""Get Transformer2W _forecasts."""
get__forecasts(value::Transformer2W) = value._forecasts
"""Get Transformer2W internal."""
get_internal(value::Transformer2W) = value.internal


text:
module Pluto
export Notebook, Cell, run

import Pkg

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))
const PLUTO_VERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(PKG_ROOT_DIR, "Project.toml"))["version"])
const PLUTO_VERSION_STR = 'v' * string(PLUTO_VERSION)
const JULIA_VERSION_STR = 'v' * string(VERSION)

@info """\n
    Welcome to Pluto $(PLUTO_VERSION_STR)! ⚡

    Let us know what you think:
    https://github.com/fonsp/Pluto.jl
\n"""

include("./react/ExploreExpression.jl")
using .ExploreExpression

include("./notebookserver/FormatOutput.jl")
include("./notebookserver/Cell.jl")
include("./notebookserver/Notebook.jl")
include("./notebookserver/Client.jl")

include("./react/WorkspaceManager.jl")
include("./react/Errors.jl")
include("./react/React.jl")
include("./react/Run.jl")

include("./webserver/WebServer.jl")
include("./webserver/Static.jl")
include("./webserver/Dynamic.jl")
include("./webserver/Completions.jl")

end


text:
using Plots
using DelimitedFiles
using .GSVD
using LinearAlgebra

function genetest(m, p, n)
    # y = readdlm("m-p-n=$m-$p-$n.txt")
    maxI = 20
    maxJ = 30
    x = zeros(maxI)
    tsvd = zeros(maxI)
    tgsvd = zeros(maxI)

    for i = 1:maxI
        x[i] = m*i*6
        time1 = 0.0
        time2 = 0.0
        for j = 1:maxJ
            A1 = randn(m*i*6, n*i*6)
            B1 = randn(p*i*6, n*i*6)
            A2 = deepcopy(A1)
            B2 = deepcopy(B1)
            ans1 = @timed svd(A1, B1)
            ans2 = @timed gsvd(A2, B2)
            time1 = time1 + ans1[2]
            time2 = time2 + ans2[2]
        end
            tsvd[i] = time1/maxJ
            tgsvd[i] = time2/maxJ

        println("m: ", x[i])
        println("--------------------------------------------------------")
    end

    y = [tsvd tgsvd];

    display(y);

    x = zeros(maxI,1)
    for i=1:maxI
        x[i] = m*6*i
    end
    gr();
    pl = plot(x, y,
    # xaxis=:log,
    yaxis=:log,
    title = "CPU timing of GSVD, m:p:n = $m:$p:$n",
    label = ["svd(A, B)" "gsvd(A, B)"],
    # yformatter = :scientific,
    xlabel = "m: number of rows of A",
    ylabel = "elapsed time in seconds",
    linecolor = [:firebrick2 :steelblue2],
    markershape = [:circle :xcross],
    markercolor = [:firebrick2 :steelblue2],
    markersize = 3.5,
    markeralpha = 0.75,
    markerstrokewidth=0,
    lw = 1.5,
    legend = :topleft)
    # savefig("myplot.png") # Saves the CURRENT_PLOT as a .png
    # savefig(pl, "myplot.pdf") # Saves the plot from p as a .pdf vector graphic
    # return y, pl
    # Check pwd()
    savefig(pl, "m-p-n=$m-$p-$n.pdf")
    savefig("m-p-n=$m-$p-$n.png")

    open("m-p-n=$m-$p-$n.txt", "w") do io
        writedlm(io, y)
    end

    pl = bar(
    string.(6*m:6*m:6*m*maxI),
    y[:,1]./y[:,2],
    xlabel = "m: number of rows of A",
    xticks = :all,
    ylabel = "factors of speedup",
    title = "t(svd(A, B))/t(gsvd(A, B)), m:p:n=$m:$p:$n",
    fillcolor = :royalblue1,
    fillalpha=0.25,
    legend = nothing)
    savefig(pl, "speedup_m-p-n=$m-$p-$n.pdf")
    savefig("speedup_m-p-n=$m-$p-$n.png")
end


text:
using DSGE, Test

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:forecast_horizons, 2)

# Scenario
scen = Scenario(:testscen, "Test Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF")
@testset "Test scenario basic constructor" begin
    @test scen.shock_scaling == 1.0
    @test n_targets(scen) == 2
    @test n_instruments(scen) == 2

    for var in scen.target_names
        scen.targets[!, var] = rand(2)
    end
    @test n_target_horizons(scen) == 2

    global act = targets_to_data(m, scen)
    global expt = copy(scen.targets)
    expt[!, :date] = DSGE.quarter_range(date_forecast_start(m), date_forecast_end(m))
    expt[!, :obs_nominalrate] = fill(missing, 2)
    for col in names(expt)
        global check = false
        if ismissing(act[!, col] == expt[!, col])
            check = (all(ismissing.(act[!, col])) && all(ismissing.(expt[!, col])))
        else
            check = (act[!, col] == expt[!, col])
        end
        @test check
    end
end

# Scenario with shock scaling
@testset "Test scenario constructor with shock scaling" begin
    scale = Scenario(:scaledscen, "Test Shock Scaling Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF",
                     shock_scaling = 2.0)
    @test scale.shock_scaling == 2.0

    # SwitchingScenario
    @test_throws ErrorException SwitchingScenario(:switch, scen, scen, [1.0], [1.0, 1.0])
    @test_throws ErrorException SwitchingScenario(:switch, scen, scen, [-1.0, 1.0], [1.0, 1.0])
    @test_throws ErrorException SwitchingScenario(:switch, scen, scen, [1.0, 1.0], [-1.0, 1.0])

    # ScenarioAggregate
    @test_throws ErrorException ScenarioAggregate(:agg, "Test Scenario Aggregate", AbstractScenario[scen],
                                                  [0.5, 0.5], 10, true, "REF")
    @test_throws ErrorException ScenarioAggregate(:agg, "Test Scenario Aggregate", AbstractScenario[scen],
                                                  [-1.0], 10, true, "REF")
    @test_throws ErrorException ScenarioAggregate(:agg, "Test Scenario Aggregate", AbstractScenario[scen],
                                                  [0.5], 10, true, "REF")
end


text:
# Sparisty-promoting LS-RTM of the 2D Marmousi model with on-the-fly Fourier transforms
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra
using JUDI, SegyIO, HDF5, JOLI, PyPlot, SlimOptim

# Load migration velocity model
if ~isfile("marmousi_model.h5")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_model.h5`)
end
n, d, o, m0 = read(h5open("marmousi_model.h5", "r"), "n", "d", "o", "m0")

# Set up model structure
model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

# Load data
if ~isfile("marmousi_2D.segy")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy`)
end
block = segy_read("marmousi_2D.segy")
d_lin = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key = "source", segy_depth_key = "SourceDepth")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.03)    # 30 Hz wavelet
q = judiVector(src_geometry, wavelet)

# Set up info structure
ntComp = get_computational_nt(q.geometry, d_lin.geometry, model0)  # no. of computational time steps
info = Info(prod(model0.n), d_lin.nsrc, ntComp)

###################################################################################################

# Setup operators
opt = Options(optimal_checkpointing=false, subsampling_factor=3)
M = judiModeling(model0, q.geometry, d_lin.geometry; options=opt)
J = judiJacobian(M, q)

batchsize = 4
# Right-hand preconditioners (model topmute)
Mr = judiTopmute(model0.n, 52, 10)

C = joCurvelet2D(model0.n[1], model0.n[2]; zero_finest = true, DDT = Float32, RDT = Float32)

function breg_obj(x)
    i = randperm(d_lin.nsrc)[1:batchsize]
    Ml = judiMarineTopmute2D(30, d_lin[i].geometry)
    r = Ml*J[i]*Mr*x - Ml*d_lin[i]
    g = adjoint(Mr)*adjoint(J[i])*adjoint(Ml)*r
    return  .5f0*norm(r)^2, g[1:end]
end

opt = bregman_options(maxIter=5, verbose=2, quantile=.9, alpha=1, antichatter=true)#, spg=true)
sol = bregman(breg_obj, 0f0.*vec(m0), opt, C)


text:
function parse_kw(ps::ParseState)
    k = kindof(ps.t)
    if k == Tokens.IF
        return @default ps @closer ps :block parse_if(ps)
    elseif k == Tokens.LET
        return @default ps @closer ps :block parse_let(ps)
    elseif k == Tokens.TRY
        return @default ps @closer ps :block parse_try(ps)
    elseif k == Tokens.FUNCTION
        return @default ps @closer ps :block parse_function(ps)
    elseif k == Tokens.MACRO
        return @default ps @closer ps :block parse_macro(ps)
    elseif k == Tokens.BEGIN
        @static if VERSION < v"1.4"
            return @default ps @closer ps :block parse_begin(ps)
        else
            if ps.closer.inref
                ret = mKEYWORD(ps)
            else
                return @default ps @closer ps :block parse_begin(ps)
            end
        end
    elseif k == Tokens.QUOTE
        return @default ps @closer ps :block parse_quote(ps)
    elseif k == Tokens.FOR
        return @default ps @closer ps :block parse_for(ps)
    elseif k == Tokens.WHILE
        return @default ps @closer ps :block parse_while(ps)
    elseif k == Tokens.BREAK
        return INSTANCE(ps)
    elseif k == Tokens.CONTINUE
        return INSTANCE(ps)
    elseif k == Tokens.IMPORT
        return parse_imports(ps)
    elseif k == Tokens.USING
        return parse_imports(ps)
    elseif k == Tokens.EXPORT
        return parse_export(ps)
    elseif k == Tokens.MODULE || k == Tokens.BAREMODULE
        return @default ps @closer ps :block parse_module(ps)
    elseif k == Tokens.CONST
        return @default ps parse_const(ps)
    elseif k == Tokens.GLOBAL
        return @default ps parse_global(ps)
    elseif k == Tokens.LOCAL
        return @default ps parse_local(ps)
    elseif k == Tokens.RETURN
        return @default ps parse_return(ps)
    elseif k == Tokens.END
        if ps.closer.square
            ret = mKEYWORD(ps)
        else
            ret = mErrorToken(mIDENTIFIER(ps), UnexpectedToken)
            ps.errored = true
        end

        return ret
    elseif k == Tokens.ELSE || k == Tokens.ELSEIF || k == Tokens.CATCH || k == Tokens.FINALLY
        ps.errored = true
        return mErrorToken(mIDENTIFIER(ps), UnexpectedToken)
    elseif k == Tokens.ABSTRACT
        return @default ps parse_abstract(ps)
    elseif k == Tokens.PRIMITIVE
        return @default ps parse_primitive(ps)
    elseif k == Tokens.TYPE
        return mIDENTIFIER(ps)
    elseif k == Tokens.STRUCT
        # return setbinding!(@default ps @closer ps block parse_struct(ps, false))
        return @default ps @closer ps :block parse_struct(ps, false)
    elseif k == Tokens.MUTABLE
        # return setbinding!(@default ps @closer ps block parse_mutable(ps))
        return @default ps @closer ps :block parse_mutable(ps)
    elseif k == Tokens.OUTER
        return mIDENTIFIER(ps)
    else
        ps.errored = true
        return mErrorToken(Unknown)
    end
end
# Prefix

function parse_const(ps::ParseState)
    kw = mKEYWORD(ps)
    arg = parse_expression(ps)
    if !((typof(arg) === BinaryOpCall && kindof(arg.args[2]) === Tokens.EQ) || (typof(arg) === Global && typof(arg.args[2]) === BinaryOpCall && kindof(arg.args[2].args[2]) === Tokens.EQ))
        ps.errored = true
        arg = mErrorToken(arg, ExpectedAssignment)
    end
    ret = EXPR(Const, EXPR[kw, arg])
    return ret
end

function parse_global(ps::ParseState)
    kw = mKEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR(Global, EXPR[kw, arg])
end

function parse_local(ps::ParseState)
    kw = mKEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR(Local, EXPR[kw, arg])
end

function parse_return(ps::ParseState)
    kw = mKEYWORD(ps)
    args = closer(ps) ? NOTHING() : parse_expression(ps)

    return EXPR(Return, EXPR[kw, args])
end

# One line

@addctx :abstract function parse_abstract(ps::ParseState)
    if kindof(ps.nt) == Tokens.TYPE
        kw1 = mKEYWORD(ps)
        kw2 = mKEYWORD(next(ps))
        sig = @closer ps :block parse_expression(ps)
        ret = EXPR(Abstract, EXPR[kw1, kw2, sig, accept_end(ps)])
    else
        ret = mIDENTIFIER(ps)
    end
    return ret
end

@addctx :primitive function parse_primitive(ps::ParseState)
    if kindof(ps.nt) == Tokens.TYPE
        kw1 = mKEYWORD(ps)
        kw2 = mKEYWORD(next(ps))
        sig = @closer ps :ws @closer ps :wsop parse_expression(ps)
        arg = @closer ps :block parse_expression(ps)

        ret = EXPR(Primitive, EXPR[kw1, kw2, sig, arg, accept_end(ps)])
    else
        ret = mIDENTIFIER(ps)
    end
    return ret
end

function parse_imports(ps::ParseState)
    kw = mKEYWORD(ps)
    kwt = is_import(kw) ? Import : Using
    tk = kindof(ps.t)

    arg = parse_dot_mod(ps)

    if kindof(ps.nt) != Tokens.COMMA && kindof(ps.nt) != Tokens.COLON
        ret = EXPR(kwt, vcat(kw, arg))
    elseif kindof(ps.nt) == Tokens.COLON
        ret = EXPR(kwt, vcat(kw, arg))
        push!(ret, mOPERATOR(next(ps)))

        arg = parse_dot_mod(ps, true)
        append!(ret, arg)
        while kindof(ps.nt) == Tokens.COMMA
            accept_comma(ps, ret)
            arg = parse_dot_mod(ps, true)
            append!(ret, arg)
        end
    else
        ret = EXPR(kwt, vcat(kw, arg))
        while kindof(ps.nt) == Tokens.COMMA
            accept_comma(ps, ret)
            arg = parse_dot_mod(ps)
            append!(ret, arg)
        end
    end

    return ret
end

function parse_export(ps::ParseState)
    args = EXPR[mKEYWORD(ps)]
    append!(args, parse_dot_mod(ps))

    while kindof(ps.nt) == Tokens.COMMA
        push!(args, mPUNCTUATION(next(ps)))
        arg = parse_dot_mod(ps)[1]
        push!(args, arg)
    end

    return EXPR(Export, args)
end

# Block

@addctx :begin function parse_begin(ps::ParseState)
    sb = ps.t.startbyte
    kw = mKEYWORD(ps)
    blockargs = parse_block(ps, EXPR[], (Tokens.END,), true)
    if isempty(blockargs)
        block = EXPR(Block, blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb - kw.fullspan
        block = EXPR(Block, blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR(Begin, EXPR[kw, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

@addctx :quote function parse_quote(ps::ParseState)
    kw = mKEYWORD(ps)
    blockargs = parse_block(ps)
    return EXPR(Quote, EXPR[kw, EXPR(Block, blockargs), accept_end(ps)])
end

@addctx :function function parse_function(ps::ParseState)
    kw = mKEYWORD(ps)
    sig = @closer ps :inwhere @closer ps :ws parse_expression(ps)
    if typof(sig) === InvisBrackets && !(typof(sig.args[2]) === TupleH || (typof(sig.args[2]) === Block) || (typof(sig.args[2]) === UnaryOpCall && kindof(sig.args[2].args[2]) === Tokens.DDDOT))
        istuple = true
        sig = EXPR(TupleH, sig.args)
    elseif typof(sig) === TupleH
        istuple = true
    else
        istuple = false
    end

    while kindof(ps.nt) == Tokens.WHERE && kindof(ps.ws) != Tokens.NEWLINE_WS
        sig = @closer ps :inwhere @closer ps :ws parse_operator_where(ps, sig, INSTANCE(next(ps)), false)
    end
    blockargs = parse_block(ps)

    if isempty(blockargs)
        if typof(sig) === Call || typof(sig) === WhereOpCall || (typof(sig) === BinaryOpCall && !is_exor(sig.args[1])) || istuple || (typof(sig) === InvisBrackets && typof(sig.args[2]) === Block)
            args = EXPR[sig, EXPR(Block, blockargs)]
        else
            args = EXPR[sig]
        end
    else
        args = EXPR[sig, EXPR(Block, blockargs)]
    end

    ret = EXPR(FunctionDef, EXPR[kw])
    for a in args
        push!(ret, a)
    end
    accept_end(ps, ret)
    return ret
end

@addctx :macro function parse_macro(ps::ParseState)
    sb  = ps.t.startbyte
    kw = mKEYWORD(ps)
    sig = @closer ps :inwhere @closer ps :ws parse_expression(ps)
    sb1  = ps.nt.startbyte
    blockargs = parse_block(ps)

    if isidentifier(sig)
        ender = accept_end(ps)
        fullspan1 = ps.nt.startbyte - sb
        ret = EXPR(Macro, EXPR[kw, sig, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    elseif isempty(blockargs)
        ender = accept_end(ps)
        fullspan1 = ps.nt.startbyte - sb
        ret = EXPR(Macro, EXPR[kw, sig, EXPR(Block, EXPR[]), ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR(Block, blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
        ender = accept_end(ps)
        fullspan1 = ps.nt.startbyte - sb
        ret = EXPR(Macro, EXPR[kw, sig, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    end
    return ret
end

# loops
@addctx :for function parse_for(ps::ParseState)
    sb  = ps.t.startbyte
    kw = mKEYWORD(ps)
    ranges = parse_iterators(ps)
    sb1  = ps.nt.startbyte
    blockargs = parse_block(ps)

    if isempty(blockargs)
        block = EXPR(Block, blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR(Block, blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR(For, EXPR[kw, ranges, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

@addctx :while function parse_while(ps::ParseState)
    sb = ps.t.startbyte
    kw = mKEYWORD(ps)
    cond = @closer ps :ws parse_expression(ps)
    sb1 = ps.nt.startbyte
    blockargs = parse_block(ps)

    if isempty(blockargs)
        block = EXPR(Block, blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR(Block, blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR(While, EXPR[kw, cond, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

# control flow

"""
    parse_if(ps, ret, nested=false, puncs=[])

Parse an `if` block.
"""
@addctx :if function parse_if(ps::ParseState, nested = false)
    # Parsing
    kw = mKEYWORD(ps)
    if kindof(ps.ws) == NewLineWS || kindof(ps.ws) == SemiColonWS
        ps.errored = true
        cond = mErrorToken(MissingConditional)
    else
        cond = @closer ps :ws parse_expression(ps)
    end
    ifblockargs = parse_block(ps, EXPR[], (Tokens.END, Tokens.ELSE, Tokens.ELSEIF))

    if nested
        ret = EXPR(If, EXPR[cond, EXPR(Block, ifblockargs)])
    else
        ret = EXPR(If, EXPR[kw, cond, EXPR(Block, ifblockargs)])
    end

    elseblockargs = EXPR[]
    if kindof(ps.nt) == Tokens.ELSEIF
        push!(ret, mKEYWORD(next(ps)))
        push!(elseblockargs, parse_if(ps, true))
    end
    elsekw = kindof(ps.nt) == Tokens.ELSE
    if kindof(ps.nt) == Tokens.ELSE
        push!(ret, mKEYWORD(next(ps)))
        parse_block(ps, elseblockargs)
    end

    # Construction
    if !(isempty(elseblockargs) && !elsekw)
        push!(ret, EXPR(Block, elseblockargs))
    end
    !nested && accept_end(ps, ret)

    return ret
end

function is_wrapped_assignment(x::EXPR)
    if is_assignment(x)
        return true
    elseif typof(x) === CSTParser.InvisBrackets && x.args isa Vector{EXPR} && length(x.args) == 3
        return is_wrapped_assignment(x.args[2])
    end
    return false
end

@addctx :let function parse_let(ps::ParseState)
    args = EXPR[mKEYWORD(ps)]
    if !(kindof(ps.ws) == NewLineWS || kindof(ps.ws) == SemiColonWS)
        arg = @closer ps :comma @closer ps :ws  parse_expression(ps)
        if kindof(ps.nt) == Tokens.COMMA || !(is_wrapped_assignment(arg) || isidentifier(arg))
            arg = EXPR(Block, EXPR[arg])
            while kindof(ps.nt) == Tokens.COMMA
                accept_comma(ps, arg)
                startbyte = ps.nt.startbyte
                nextarg = @closer ps :comma @closer ps :ws parse_expression(ps)
                push!(arg, nextarg)
            end
        end
        push!(args, arg)
    end

    blockargs = parse_block(ps)
    push!(args, EXPR(Block, blockargs))
    accept_end(ps, args)

    return EXPR(Let, args)
end

@addctx :try function parse_try(ps::ParseState)
    kw = mKEYWORD(ps)
    ret = EXPR(Try, EXPR[kw])

    tryblockargs = parse_block(ps, EXPR[], (Tokens.END, Tokens.CATCH, Tokens.FINALLY))
    push!(ret, EXPR(Block, tryblockargs))

    #  catch block
    if kindof(ps.nt) == Tokens.CATCH
        next(ps)
        push!(ret, mKEYWORD(ps))
        # catch closing early
        if kindof(ps.nt) == Tokens.FINALLY || kindof(ps.nt) == Tokens.END
            caught = FALSE()
            catchblock = EXPR(Block, EXPR[])
        else
            if kindof(ps.ws) == SemiColonWS || kindof(ps.ws) == NewLineWS
                caught = FALSE()
            else
                caught = @closer ps :ws parse_expression(ps)
            end

            catchblockargs = parse_block(ps, EXPR[], (Tokens.END, Tokens.FINALLY))
            if !(isidentifier(caught) || kindof(caught) == Tokens.FALSE || (typof(caught) === UnaryOpCall && isoperator(caught.args[1]) && kindof(caught.args[1]) == Tokens.EX_OR))
                pushfirst!(catchblockargs, caught)
                caught = FALSE()
            end
            catchblock = EXPR(Block, catchblockargs)
        end
    else
        caught = FALSE()
        catchblock = EXPR(Block, EXPR[])
    end
    push!(ret, caught)
    push!(ret, catchblock)

    # finally block
    if kindof(ps.nt) == Tokens.FINALLY
        if isempty(catchblock.args)
            ret.args[4] = setparent!(FALSE(), ret)
        end
        push!(ret, mKEYWORD(next(ps)))
        finallyblockargs = parse_block(ps)
        push!(ret, EXPR(Block, finallyblockargs))
    end

    push!(ret, accept_end(ps))
    return ret
end

@addctx :do function parse_do(ps::ParseState, ret::EXPR)
    kw = mKEYWORD(next(ps))

    args = EXPR(TupleH, EXPR[])
    @closer ps :comma @closer ps :block while !closer(ps)
        a = parse_expression(ps)
        push!(args, a)
        if kindof(ps.nt) == Tokens.COMMA
            accept_comma(ps, args)
        end
    end

    blockargs = parse_block(ps)

    return EXPR(Do, EXPR[ret, kw, args, EXPR(Block, blockargs), accept_end(ps)])
end

# modules

@addctx :module function parse_module(ps::ParseState)
    sb = ps.t.startbyte
    kw = mKEYWORD(ps)
    @assert kindof(kw) == Tokens.MODULE || kindof(kw) == Tokens.BAREMODULE # work around julia issue #23766
    if isidentifier(ps.nt)
        arg = mIDENTIFIER(next(ps))
    else
        arg = @precedence ps 15 @closer ps :ws parse_expression(ps)
    end
    sb1 = ps.nt.startbyte

    blockargs = parse_block(ps, EXPR[], (Tokens.END,), true)

    if isempty(blockargs)
        block = EXPR(Block, blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR(Block, blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR(is_module(kw) ? ModuleH : BareModule, EXPR[kw, arg, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

function parse_mutable(ps::ParseState)
    if kindof(ps.nt) == Tokens.STRUCT
        kw = mKEYWORD(ps)
        next(ps)
        ret = parse_struct(ps, true)
        pushfirst!(ret, kw)
        update_span!(ret)
    else
        ret = mIDENTIFIER(ps)
    end
    return ret
end

function markparameters!(sig::EXPR)
    signame = rem_where_subtype(sig)
    if typof(signame) === Curly
        for i = 3:length(signame.args) - 1
            if !(typof(signame.args[i]) === PUNCTUATION)
            end
        end
    end
    return sig
end

@addctx :struct function parse_struct(ps::ParseState, mutable::Bool)
    sb = ps.t.startbyte
    kw = mKEYWORD(ps)
    sig = @closer ps :ws parse_expression(ps)

    sb1 = ps.nt.startbyte
    blockargs = parse_block(ps)
    if isempty(blockargs)
        block = EXPR(Block, blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR(Block, blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    if mutable
        ret = EXPR(Mutable, EXPR[kw, sig, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    else
        ret = EXPR(Struct, EXPR[kw, sig, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    end
    return ret
end


text:
"""
    get_morgan_fp(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String

Get Morgan (ECFP like) fingerprints.

```julia
mfp = get_morgan_fp(mol)
```

```julia
fp_details = Dict{String, Any}("nBits" => 512, "radius" => 2)
morgan_fp = get_morgan_fp(mol, fp_details)

```
"""
function get_morgan_fp(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String
    details_json::String = jsonify_details(details)
    val::Cstring = ccall((:get_morgan_fp, librdkitcffi), Cstring, (Cstring, Csize_t, Cstring), mol.mol[], mol.mol_size[], details_json)
    mfp = unsafe_string_and_free(val)
    return mfp
end

"""
    get_morgan_fp_as_bytes(mol::Mol)

Get Morgan (ECFP like) fingerprints as bytes.

```julia
morgan_fp_bytes = get_morgan_fp_as_bytes(mol)
```
"""
function get_morgan_fp_as_bytes(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::Vector{UInt8}
    details_json::String = jsonify_details(details)
    n_bytes = Ref{Csize_t}(0)
    val::Cstring = ccall((:get_morgan_fp_as_bytes, librdkitcffi), Cstring, (Cstring, Csize_t, Ref{Csize_t}, Cstring), mol.mol[], mol.mol_size[], n_bytes, details_json)
    mfp = unsafe_string_and_free(val, n_bytes)
    return Vector{UInt8}(mfp)
end

"""
    get_rdkit_fp(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String

Get RDKit fingerprints.

```julia
rfp = get_rdkit_fp(mol)
```
"""
function get_rdkit_fp(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String
    details_json::String = jsonify_details(details)
    val::Cstring = ccall((:get_rdkit_fp, librdkitcffi), Cstring, (Cstring, Csize_t, Cstring), mol.mol[], mol.mol_size[], details_json)
    rfp = unsafe_string_and_free(val)
    return rfp
end

"""
    get_rdkit_fp_as_bytes(mol::Mol)

Get RDKit fingerprints as bytes.

```julia
rfp_bytes = get_rdkit_fp_as_bytes(mol)
```
"""
function get_rdkit_fp_as_bytes(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::Vector{UInt8}
    details_json::String = jsonify_details(details)
    n_bytes = Ref{Csize_t}(0)
    val::Cstring = ccall((:get_rdkit_fp_as_bytes, librdkitcffi), Cstring, (Cstring, Csize_t, Ref{Csize_t}, Cstring), mol.mol[], mol.mol_size[], n_bytes, details_json)
    rfp = unsafe_string_and_free(val, n_bytes)
    return Vector{UInt8}(rfp)
end

"""
    get_pattern_fp(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String

Get RDKit Pattern fingerprints.

```julia
pfp = get_pattern_fp(mol)
```
"""
function get_pattern_fp(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::String
    details_json::String = jsonify_details(details)
    val::Cstring = ccall((:get_pattern_fp, librdkitcffi), Cstring, (Cstring, Csize_t, Cstring), mol.mol[], mol.mol_size[], details_json)
    pfp = unsafe_string_and_free(val)
    return pfp
end

"""
    get_pattern_fp_as_bytes(mol::Mol)

Get RDKit Pattern fingerprints as bytes.

```julia
pfp_bytes = get_pattern_fp_as_bytes(mol)
```
"""
function get_pattern_fp_as_bytes(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)::Vector{UInt8}
    details_json::String = jsonify_details(details)
    n_bytes = Ref{Csize_t}(0)
    val::Cstring = ccall((:get_pattern_fp_as_bytes, librdkitcffi), Cstring, (Cstring, Csize_t, Ref{Csize_t}, Cstring), mol.mol[], mol.mol_size[], n_bytes, details_json)
    pfp = unsafe_string_and_free(val, n_bytes)
    return Vector{UInt8}(pfp)
end

"""
    get_descriptors(mol::Mol)

Get physico-chemical descriptors.

```julia
descs = get_descriptors(mol)
```
"""
function get_descriptors(mol::Mol)
    val::Cstring = ccall((:get_descriptors, librdkitcffi), Cstring, (Cstring, Csize_t), mol.mol[], mol.mol_size[])
    json = JSON.parse(unsafe_string_and_free(val))
    return json
end


text:
#############################################################################
# contiuous_sets.jl
# Handles generic combinatorial sets and their operations.
#############################################################################

export GenCombiSet, get_elements, cardinality

# Instance of generic Combinatorial Sets
abstract type GenCombiSet <: CombiSet end

function get_elements(x::CombiSet)
  if x.elements == nothing
    error("This set variable has not been assigned elements yet; cannot operate on unassigned set variables.")
  else
    return x.elements
  end
end

function cardinality(x::CombiSet)
  return x.cardinality
end


text:
P = 20000
i = 0.065
# semi-annually means twice a year (every six months)
m = 2
N = 25

F = P * (1 + i / 2)^(N * m)


text:
# [`ComposedProjectiveTransform`](#) implements efficient composition
# of `ProjectiveTransform`s.

"""
    ComposedProjectiveTransform(tfms...)

Wrap multiple projective `tfms` and apply them efficiently.
The projections are fused into a single projection and only
points inside the final crop are evaluated.
"""
struct ComposedProjectiveTransform{T<:Tuple} <: ProjectiveTransform
    tfms::T
end

ComposedProjectiveTransform(tfms...) =
    ComposedProjectiveTransform(tfms)

# Composing `ProjectiveTransform`s should create a `ComposedProjectiveTransform`:

compose(tfm1::ProjectiveTransform, tfm2::ProjectiveTransform) =
    ComposedProjectiveTransform(tfm1, tfm2)

compose(composed::ComposedProjectiveTransform, tfm::ProjectiveTransform) =
    ComposedProjectiveTransform(composed.tfms..., tfm)

compose(tfm::ProjectiveTransform, composed::ComposedProjectiveTransform) =
    ComposedProjectiveTransform(tfm, composed.tfms...)

# The random state is collected from the transformations that make up the
# `ComposedProjectiveTransform`:

getrandstate(tfm::ComposedProjectiveTransform) = getrandstate.(tfm.tfms)

# To combine the concrete projections `P`, we use `ImageTransformations.compose`
# (aliased to `∘`). The bounds are also transformed at each step.

function getprojection(
        composed::ComposedProjectiveTransform,
        bounds;
        randstate = getrandstate(composed))
    P = CoordinateTransformations.IdentityTransformation()
    for (tfm, r) in zip(composed.tfms, randstate)
        P_tfm = getprojection(tfm, bounds; randstate = r)
        bounds = P_tfm.(bounds)
        P = P_tfm ∘ P
    end

    return P
end

# To work nicely with crops, [`cropindices`](#) simply defers to the
# last transform:

function cropindices(composed::ComposedProjectiveTransform, P, bounds; randstate = getrandstate(composed))
    cropindices(composed.tfms[end], P, bounds; randstate = randstate[end])
end


text:
function make_analog_gui(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    push!(b["analog_c_box"],handles.analog.c)
    setproperty!(handles.analog.c,:hexpand,true)
    setproperty!(handles.analog.c,:vexpand,true)
    show(handles.analog.c)

end

function add_analog_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(analog_time_cb,b["analog_x_zoom_minus"],"clicked",Void,(),false,(handles,0))
    signal_connect(analog_time_cb,b["analog_x_zoom_plus"],"clicked",Void,(),false,(handles,1))

    nothing
end

function analog_time_cb(w::Ptr,user_data::Tuple{Tracker_Handles,Int})

    han, time_zoom = user_data

    if time_zoom == 0
        han.analog.t_zoom = han.analog.t_zoom - 1
        if han.analog.t_zoom < 0
            han.analog.t_zoom = 0
        end
    else
        han.analog.t_zoom = han.analog.t_zoom + 1
    end

    update_analog_canvas(han)

    nothing
end

function update_analog_canvas(han::Tracker_Handles)

    frame_range = han.analog.t_zoom
    (lower_id, upper_id) = get_lower_upper_frame(han.displayed_frame,frame_range,han.max_frames)

    ctx=Gtk.getgc(han.analog.c)

    set_source_rgb(ctx,1,1,1)
    paint(ctx)

    w=width(ctx)
    h=height(ctx)

    #center line
    set_source_rgb(ctx,1,0,0)
    cent = (upper_id - han.displayed_frame) / (upper_id - lower_id)
    move_to(ctx,cent * w,0)
    line_to(ctx,cent * w,h)
    stroke(ctx)

    #exclude
    #e_line = make_line(han.man.exclude_block,lower_id,upper_id,w)
    #set_source_rgb(ctx,0,0,0)
    #draw_manual_line(han.analog.c,e_line,2)

    #spikes
    set_source_rgb(ctx,0,0,0)

    #Find points from along signal that are within the frame

    #Get min ind
    a_min=get_min_analog(han,lower_id,1)

    #Get center ind
    get_min_analog(han,han.displayed_frame,1)

    #Get max ind
    a_max=get_max_analog(han,upper_id,1)

    #myline = rand(-100:100,round(Int,w))
    off = h/2

    an_range = range(a_min,a_max,length=round(Int64,w))
    move_to(ctx,0,off+han.analog.var[1][a_min] * han.analog.gains[1])
    aa = a_min + 1
    a_val = 0
    ww = 1
    while aa < a_max

        #=
        if abs(han.analog.var[1][aa]) > a_val
            a_val = han.analog.var[1][aa]
        end

        if aa > an_range[ww]
            line_to(ctx,ww,a_val * han.analog.gains[1] + off)
            a_val = 0
            ww += 1
        end
        =#
        line_to(ctx,ww + aa / an_range[ww],han.analog.var[1][aa] * han.analog.gains[1] + off)

        if aa > an_range[ww]
            ww += 1
        end
        aa += 1
    end
    stroke(ctx)

    reveal(han.analog.c)

    nothing
end

function get_min_analog(han,b1::Int,an::Int)

    findfirst(han.analog.ts[an].>b1)[1]
end

function get_max_analog(han,b1::Int,an::Int)

    findfirst(han.analog.ts[an].>b1)[1]-1
end


text:
include("test_common.jl")
using Dates

@testset "objects" begin
    new_context()

    Object = node"Object"
    node"const a = new Object()"
    a = node"a"
    b = @new Object
    b.a = node"a"
    @test b.a == node"a"
    @test node"(b) => b.a === a"(b)
    b[a] = true
    @test hasproperty(b, :a)
    @test a ∈ b
    @test b[a]
    b[1] = 123
    @test b[1] == 123
    delete!(b, 1)
    delete!(b, "a")
    delete!(b, a)
    @test all((!haskey).([b], ["a", a, 1]))

    node"""class A {
        constructor(a, b) {
            this.c = a * b
        }
    }"""
    ClassA = node"A"
    instance = @new ClassA(4, 5)
    @test instance.c == 20
    @test instanceof(instance, ClassA)
    @test !instanceof(instance, nothing)

    @test node"new Date('2007-01-28')" == DateTime(2007, 1, 28)

    str_dict = Dict("k"=>false, "p"=>123)
    any_dict = Dict(false=>false, 123=>456)
    @test "k" ∈ node_value(str_dict)
    @test Object.keys(str_dict) == ["k", "p"]
    @test Object.keys(any_dict) == ["false", "123"]
    js_map = node"""(() => {
        const m = new Map()
        m.set(1, 2)
        m.set("x", 3)
        return m
    })()"""
    @test value(js_map) == Dict(1=>2, "x"=>3)

    some_set = Set([1, 2])
    mutated_set = node"""(set) => {
        set.add(3)
        return set
    }"""(some_set)
    @test mutated_set == Set([1, 2, 3])
    # `Set` is not mutable.
    @test some_set ≠ mutated_set

    tuple = ("a", "b")
    tuple2 = node"""(tuple) => {
        tuple[0] = "c"
        return tuple
    }"""(tuple)
    @test tuple2 == ("c", "b")
    tuple3 = node"""(tuple) => {
        tuple.push(123)
        return tuple
    }"""(tuple)
    @test tuple3 == ["a", "b", 123]

    named_tuple = (x=1, y=true, z="z")
    @test Object.keys(named_tuple) == ["x", "y", "z"]
    named_tuple2 = node"""(nt) => {
        assert(nt.x === 1 && nt.y && nt.z === 'z')
        nt.x = 5
        return nt
    }"""(named_tuple)
    @test named_tuple2 ≡ (x=5, y=true, z="z")
    named_tuple3 = node"""(nt) => {
        nt.x = 't'
        return nt
    }"""(named_tuple)
    @test named_tuple3 isa JsObject
    @test "x" in named_tuple3
    @test NamedTuple(named_tuple3) ≡ (x="t", y=true, z="z")

    struct Foo
        x::Int
    end
    @test node"(foo) => (foo.x = 3, foo)"(Foo(5)) == Foo(3)
    mutable struct FooM
        x::String
    end
    foo_m = FooM("before")
    node"(foo, s) => (foo.x = s, foo)"(foo_m, "after")
    @test foo_m.x == "after"
    @test node"(foo) => 'x' in foo && foo.x === 'after'"(foo_m)

    @test length(node"{length: 5}") == 5

    delete_context()
end


text:
module TestLinAlg
    using Base.Test
    using DataArrays

    d = @data eye(3, 3)
    d[1, 1] = NA

    svd(d)
end


text:
#=

  Coupon collecter's problem, or card collecter's problem.

  There are N different collecter's cards hidden in a package, but we don't
  know which card there is in the package we buy.
  We want to collect all of them, how many packages must one buy to collect
  all the different cards?

  See https://en.wikipedia.org/wiki/Coupon_collector%27s_problem
  """
  In probability theory, the coupon collector's problem describes 'collect all coupons and win'
  contests. It asks the following question: If each box of a brand of cereals contains a
  coupon, and there are n different types of coupons, what is the probability that more
  than t boxes need to be bought to collect all n coupons?

  An alternative statement is: Given n coupons, how many coupons do you expect you need
  to draw with replacement before having drawn each coupon at least once? The mathematical
  analysis of the problem reveals that the expected number of trials needed grows as
  Θ(n log(n).
  For example, when n = 50 it takes about 225[b] trials on average to collect all 50 coupons.

  ...

  [b]: E(50) = 50(1 + 1/2 + 1/3 + ... + 1/50) = 224.9603, the expected number of trials to
  collect all 50 coupons.

   The approximation n*log(n) + γ*n + 1/2 for this expected number gives in this case
   ≈ 195.6011 + 28.8608 + 0.5 ≈ 224.9619.  [log is the natural logarithm]
  """

  This version draw random integers until all integers 0..N-1 have been selected.
  We only study the length of the generated array.

  NOTE: THIS IS NOT A TURING MODEL.
  Just a plain simulation...

=#

#=
Running n in [10,20,30,50,100,200]
julia> @time include("coupon_collector2.jl")
n:10 theoretical prob:29.289682539682538  s:1000
mean:29.074
Theoretical: 29.289682539682538
  0.029726 seconds (883.17 k allocations: 31.116 MiB)

n:20 theoretical prob:71.95479314287364  s:1000
mean:71.929
Theoretical: 71.95479314287364
  0.241385 seconds (4.69 M allocations: 173.947 MiB, 16.75% gc time)

n:30 theoretical prob:119.84961392761171  s:1000
mean:118.299
Theoretical: 119.84961392761171
  0.499929 seconds (11.57 M allocations: 362.730 MiB, 8.86% gc time)

n:50 theoretical prob:224.96026691647114  s:1000
mean:222.996
Theoretical: 224.96026691647114
  2.083751 seconds (38.59 M allocations: 1.326 GiB, 8.76% gc time)

n:100 theoretical prob:518.737751763962  s:1000
mean:521.378
Theoretical: 518.737751763962
  9.139250 seconds (193.93 M allocations: 5.524 GiB, 8.02% gc time)

n:200 theoretical prob:1175.6061896242888  s:1000
mean:1170.562
Theoretical: 1175.6061896242888
 49.765464 seconds (1.13 G allocations: 32.495 GiB, 7.64% gc time)

 61.848598 seconds (1.38 G allocations: 39.908 GiB, 7.77% gc time)

=#

using Distributions
include("jl_utils.jl")

# Theoretical length of n coupons
# "Exact" probability from https://en.wikipedia.org/wiki/Coupon_collector%27s_problem
# (see footnote [b] above)
function theoretical_prob(n)
    return n*sum([1/i for i in 1:n])
end

function coupon_collect(n,arr)
    if length(arr) > 0 && length(unique(arr))==n
        return arr
    else
        r = rand(1:n) # a little faster
        # return coupon_collect(n,vcat(arr, r))
        return coupon_collect(n,push!(arr, r)) # slightly faster
    end
end

function run_coupon_collect(n=10,s=1000)

    theoretical = theoretical_prob(n)
    println("n:$n theoretical prob:$(theoretical_prob(n))  s:$s")

    lengths = []
    for t in 1:s
        c = coupon_collect(n,[])
        collectLen = length(c)
        # println("c:$c")
        push!(lengths,collectLen)
    end

    println("mean:", mean(lengths))
    println("Theoretical: ", theoretical)
end

n=100
@time run_coupon_collect(n,1000)
#=
for n in [10,20,30,50,100,200]
    @time run_coupon_collect(n,1000)
    println()
end
=#

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@testset verbose = true "HMM" begin

    @testset "Discrete" begin
        tr = DiscreteMarkovChain([0.3, 0.7], [0.9 0.1; 0.2 0.8])
        em1 = Dists.Normal(1, 0.3)
        em2 = Dists.Normal(2, 0.3)
        hmm = HiddenMarkovModel(tr, [em1, em2])

        states, observations = rand(hmm, 1000)

        tr_init = DiscreteMarkovChain(randprobvec(2), randtransmat(2))
        em1_init = Dists.Normal(rand(), 1)
        em2_init = Dists.Normal(rand(), 1)
        hmm_init = HiddenMarkovModel(tr_init, [em1_init, em2_init])

        hmm_est, logL_evolution = baum_welch(hmm_init, observations, iterations = 100)

        error = transition_matrix(hmm_est) - transition_matrix(hmm)

        @test maximum(error) < 0.1
    end

    @testset "Continuous" begin
        tr = ContinuousMarkovChain([0.3, 0.7], [-1.0 1.0; 2.0 -2.0])
        em1 = PoissonProcess(1.0, Dists.Normal(1, 1))
        em2 = PoissonProcess(2.0, Dists.Normal(-1, 1))
        mmpp = MarkovModulatedPoissonProcess(tr, [em1, em2])

        state_history, observations = rand(mmpp, 0.0, 1000.0)
        nb_events(observations)

    end
end


text:
####################################################################################################
# videoio. Media I/O
####################################################################################################

# 1. Reading and Writing Video

# Create the VideoCapture structures
cxx""" cv::VideoCapture VideoCapture(){ cv::VideoCapture capture = cv::VideoCapture(); return(capture); }"""
videoCapture() = @cxx VideoCapture()

cxx""" cv::VideoCapture VideoCapture(const char *filename){ cv::VideoCapture capture = cv::VideoCapture(filename); return(capture); }"""
videoCapture(filename::String) = @cxx VideoCapture(pointer(filename))

cxx""" cv::VideoCapture VideoCapture(int device){ cv::VideoCapture capture = cv::VideoCapture(device); return(capture); }"""
videoCapture(device::Int) = @cxx VideoCapture(device)   # autodetect = 0

# Functions for opening capture and grabbing frames
openVideo(capture,filename::String) = @cxx capture->open(pointer(filename))
openVideo(capture, device::Int) = @cxx capture->open(device)
isOpened(capture) = @cxx capture->isOpened()

# Useful for multi-camera environments
grab(capture) = @cxx capture->grab()
# Decodes and returns the grabbed video frame
retrieve(capture, image, flag=0) = @cxx capture->retrieve(image, flag)
# automatically called by cv::VideoCapture->open
release(capture) = @cxx capture->release()  # automatically called

# Grab, decode and return the next video frame
cxx""" cv::VideoCapture&  videoread(cv::VideoCapture& capture){ cv::Mat frame; capture >> frame; return(capture); } """
videoRead(capture) = @cxx videoread(capture)
videoRead(capture, image) = @cxx capture->read(image)

# Return the specified VideoCapture property
getVideoId(capture, propId::Int) = @cxx capture->get(propId)
# CAP_PROP_POS_MSEC       Current position of the video file (msec or timestamp)
# CAP_PROP_POS_FRAMES     0-based index of the frame to be decoded/captured next
# CAP_PROP_POS_AVI_RATIO  Relative position of the video file: 0 - start of the film, 1 - end of the film
# CAP_PROP_FRAME_WIDTH    Width of the frames in the video stream
# CAP_PROP_FRAME_HEIGHT   Height of the frames in the video stream
# CAP_PROP_FPS            frame rate
# CAP_PROP_FOURCC         4-character code of codec
# CAP_PROP_FRAME_COUNT    Number of frames in the video file
# CAP_PROP_FORMAT         Format of the Mat objects returned by retrieve()
# CAP_PROP_MODE           Backend-specific value indicating the current capture mode
# CAP_PROP_BRIGHTNESS     Brightness of the image (only for cameras)
# CAP_PROP_CONTRAST       Contrast of the image (only for cameras)
# CAP_PROP_SATURATION     Saturation of the image (only for cameras)
# CAP_PROP_HUE            Hue of the image (only for cameras)
# CAP_PROP_GAIN           Gain of the image (only for cameras)
# CAP_PROP_EXPOSURE       Exposure (only for cameras)
# CAP_PROP_CONVERT_RGB    Boolean flags indicating whether images should be converted to RGB
# CAP_PROP_WHITE_BALANCE  Currently not supported
# CAP_PROP_RECTIFICATION  Rectification flag for stereo cameras (note: only supported by DC1394 v 2.x backend currently)

# Sets a property in the VideoCapture
setVideoId(capture, propId::Int, value::Float64) = @cxx capture->set(propId, value)
# same list of prooperties as above

# Create the VideoWriter structures
cxx""" cv::VideoWriter VideoWriter(){ cv::VideoWriter writer = cv::VideoWriter(); return(writer); }"""
videoWriter() = @cxx VideoWriter()

cxx""" cv::VideoWriter VideoWriter(const char *filename, int fourcc, double fps, cv::Size frameSize,
    bool isColor=true){ cv::VideoWriter writer(filename, fourcc, fps, frameSize, isColor); return(writer); }"""
videoWriter(filename::String, fourcc::Int, fps::Float64, frameSize, isColor=true) =
    @cxx VideoWriter(pointer(filename), fourcc, fps, frameSize, isColor)
# Parameters
# filename  – Name of the output video file.
# fourcc    – Fourcc codec, e.g., fourcc('M','J','P','G')
# fps       – Framerate
# frameSize – Size of the video frames
# isColor   –  only supported for Windows

openWriter(writer, filename::String, fourcc::Int, fps::Float64, frameSize, isColor=true) =
    @cxx writer->open(pointer(filename), fourcc, fps, frameSize, isColor)

isOpened(writer) = @cxx writer->isOpened()

# Write the next video frame
writeVideo(writer, image) = @cxx writer->write(image)

# Video write with operator >>
cxx""" cv::VideoWriter&  writeframe(cv::VideoWriter& videowriter){ cv::Mat frame; videowriter << frame; return(videowriter); } """
writeVideo(videowrite) = @cxx writeframe(videowrite)

# Get int for the FourCC codec code
cxx""" int fourcc(const char *cc1, const char *cc2, const char *cc3, const char *cc4) {
           int code = cv::VideoWriter::fourcc (*cc1, *cc2, *cc3, *cc4);
           return(code);
   }
"""
# CV_FOURCC:  const Array(Ptr{UInt8}, 4)
# CV_FOURCC_MPEG
# CV_FOURCC_IYUV  #for yuv420p into an uncompressed AVI
# CV_FOURCC_DIV3  #for DivX MPEG-4 codec
# CV_FOURCC_MP42  #for MPEG-4 codec
# CV_FOURCC_DIVX  #for DivX codec
# CV_FOURCC_PIM1  #for MPEG-1 codec
# CV_FOURCC_I263  #for ITU H.263 codec
# CV_FOURCC_MPEG  #for MPEG-1 codec

fourcc(fcc::Array) = @cxx fourcc(fcc[1],fcc[2], fcc[3], fcc[4])


text:
"""
    Action <: AbstractAction

An Action can be used in the keyword arguments of an [`Object`](@ref) to define small
sub objects on the object function, such as [`appear`](@ref).

An Action should not be created by hand but instead by using one of the constructors.

# Fields
- `frames::Frames`: the frames relative to the parent [`Object`](@ref)
- `anim::Animation`: defines the interpolation function for the transition
- `func::Function`: the function that gets called in each of those frames.
    Takes the following arguments: `video, object, action, rel_frame`
- `transition::Union{Nothing, AbstractTransition}`
- `keep::Bool` defines whether this Action is called even after the last frame it was defined on
- `defs::Dict{Symbol, Any}` any kind of definitions that are relevant for the action.
"""
mutable struct Action <: AbstractAction
    frames::Frames
    anim::Animation
    func::Function
    transition::Union{Nothing,AbstractTransition}
    keep::Bool
    defs::Dict{Symbol,Any}
end

"""
    Action([frames], [Animation], func::Function; keep=true)

An `Action` gives an [`Object`](@ref) the opportunity to move, change color or much more.
It can be defined in many different ways.

# Arguments
- frames can be a `Symbol`, a `UnitRange` or a [`GFrames`](@ref) to define them in a global way.
    - **Default:** If not defined it will be the same as the previous [`Action`](@ref) or
        if it's the first action then it will be applied for the whole length of the object.
    - It defines for which frames the action acts on the object.
    - These are defined in a relative fashion so `1:10` means the first ten frames of the object
        and **not** the first ten frames of the [`Video`](@ref)
- animation can be an easing function or animation which can be defined by Animations.jl
    - **Default:** The default is `linear()`
    - Possible simple easing functions is `sineio()` for more check
        [Animations.jl](https://jkrumbiegel.github.io/Animations.jl/stable/)
- func is the function that describes the actual action
    - It can be a general function which takes in the following four arguments
        - video, object, action, rel_frame
    - If you don't need them you can write `(args...)->your_function(arg1, arg2)`
    - You often don't need an own function and instead can use predefined functions like
        - [`appear`](@ref), [`disappear`](@ref), [`follow_path`](@ref)

# Keywords
- `keep::Bool` defaults to `true` defines whether the [`Action`](@ref) is called
    even for frames after it's last defined.
    In more simple terms: If one has `Action(1:10, anim, translate())`
    It will get translated to the last position on frame `11:END_OF_OBJECT`.
    One can set `; keep = false` to turn off this behavior.

# Example
```julia
function ground(args...)
    background("black")
    sethue("white")
end

video = Video(500, 500)
Background(1:100, ground)
obj = Object((args...)->circle(O, 50, :fill))
act!(obj, Action(1:20, appear(:fade)))
act!(obj, Action(21:50, Translation(50, 50)))
act!(obj, Action(51:80, Translation(-50, -50)))
act!(obj, Action(81:100, disappear(:fade)))
render(video; pathname="test.gif")
```
"""
Action(func::Union{Function,AbstractTransition}; keep = true) =
    Action(:same, func; keep = keep)

Action(frames, easing::Union{ReversedEasing,Easing}, func::Function; keep = true) =
    Action(frames, easing_to_animation(easing), func; keep = keep)

function Action(
    frames,
    easing::Union{ReversedEasing,Easing},
    translation::Translation;
    keep = true,
)
    Action(
        frames,
        easing_to_animation(easing),
        translate();
        transition = translation,
        keep = keep,
    )
end

function Action(
    frames,
    easing::Union{ReversedEasing,Easing},
    rotation::Rotation;
    keep = true,
)
    anim = Animation([0.0, 1.0], [rotation.from, rotation.to], [easing])
    if rotation.center === nothing
        return Action(frames, anim, rotate(), transition = rotation; keep = keep)
    end
    Action(frames, anim, rotate_around(rotation.center); keep = keep)
end

function Action(frames, easing::Union{ReversedEasing,Easing}, scaling::Scaling; keep = true)
    Action(frames, easing_to_animation(easing), scale(); transition = scaling, keep = keep)
end

Action(anim::Animation, func::Function) = Action(:same, anim, func; keep = true)
Action(easing::Union{ReversedEasing,Easing}, func::Function; keep = true) =
    Action(:same, easing_to_animation(easing), func; keep = keep)

Action(frames, func::Union{Function,AbstractTransition}; keep = true) =
    Action(frames, linear(), func; keep = keep)

Action(frames, anim::Animation, func::Function; transition = nothing, keep = true) =
    Action(frames, anim, func, transition, keep, Dict{Symbol,Any}())

Base.copy(a::Action) = Action(copy(a.frames), a.anim, a.func, a.transition, a.keep, a.defs)


text:
using DrWatson
@quickactivate "StatReth"

# %%
using DataFrames
using CSV
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %% 5.45
d = CSV.read(datadir("exp_raw/Howell_1.csv"), DataFrame)

# %% 5.46
μ_female = rand(Normal(178, 20), 10_000)
μ_male = rand(Normal(178, 20), 10_000) .+ rand(Normal(0, 10), 10_000)
DataFrame((; μ_female, μ_male)) |> precis

# %% 5.47
d.sex = ifelse.(d.male .== 1, 2, 1)

# %% 5.48
@model function categ(sex, height)
    σ ~ Uniform(0, 50)
    a ~ filldist(Normal(178, 20), 2)
    μ = a[sex]
    height ~ MvNormal(μ, σ)
end

q5_8 = quap(categ(d.sex, d.height), NelderMead())
# precis

# %% 5.49
post = DataFrame(rand(q5_8.distr, 1000)', q5_8.params)
post.diff_fm = post[:, "a[1]"] .- post[:, "a[2]"]
precis(post)

# %% 5.50
d = CSV.read(datadir("exp_raw/milk.csv"), DataFrame; missingstring = "NA")
d.clade |> unique |> sort

# %% 5.51
# You can't just turn Strings into Integers in Julia but hashing them should give the same
# result
d.clade_id = Int.(indexin(d.clade, unique(d.clade)))

# %% 5.52
d.K = zscore(d.kcal_per_g)

@model function clade(clade_id, K)
    σ ~ Exponential(1)
    α ~ filldist(Normal(0, 0.5), length(unique(clade_id)))
    μ = α[clade_id]
    K ~ MvNormal(μ, σ)
end

q5_9 = quap(clade(d.clade_id, d.K))
post = DataFrame(rand(q5_9.distr, 1000)', q5_9.params)
precis(post)

# TODO: plot

# %% 5.54
d.house = rand(1:4, nrow(d))

@model function clade_house(clade_id, house, K)
    σ ~ Exponential(1)
    a ~ filldist(Normal(0, 0.5), length(unique(clade_id)))
    h ~ filldist(Normal(0, 0.5), length(unique(house)))
    μ = a[clade_id] + h[house]
    K ~ MvNormal(μ, σ)
end

q5_10 = quap(clade_house(d.clade_id, d.house, d.K))
DataFrame(rand(q5_10.distr, 1000)', q5_10.params) |> precis

max_stars_repo_path: scripts/5_3 Categorical variables.jl

max_issues_repo_path: scripts/5_3 Categorical variables.jl

max_forks_repo_path: scripts/5_3 Categorical variables.jl


text:
###################################################################################
# Antony Schutz 2015, ANR - POLCA - 2016
# TOOLS for PAINTER
# Schutz, A., Ferrari, A., Mary, D., Soulez, F., Thiébaut, E., Vannier, M.
# PAINTER: a spatiospectral image reconstruction algorithm for optical interferometry
# J. Opt. Soc. Am. A, Vol. 31, N. 11, pp 2334--2345, 2014.
#
# PAINTER Outpt
#
###################################################################################
# using PyPlot
# verb
function painterplotfct(PDATA::PAINTER_Data,OIDATA::PAINTER_Input)
    x = PDATA.x
    if sum(PDATA.w)==0
      w = ones(PDATA.w)
    else
      w = PDATA.w
    end

    crit1 = PDATA.crit1
    crit2 = PDATA.crit2
    eps1 = OIDATA.eps1
    eps2 = OIDATA.eps2
    wvl = OIDATA.wvl
    nx = OIDATA.nx
    nw = OIDATA.nw
    FOV = OIDATA.FOV
    SubColumn,SubRow = createsubplotindex(nw)

    indpix = linspace(-(FOV / 2), (FOV / 2), nx)
    @printf("stop crit, primal: %e (%e) \t dual %e (%e)\n", crit1[end], eps1, crit2[end], eps2)
###################################################################################
# Reconstructed Cube
###################################################################################
    figure(1)
    clf()
    pos     = round(Int, [1, round(Int, nx / 4), round(Int, nx / 2), round(Int, nx * 3 / 4), nx])
    count_y = 0
    count_x = 0
    for n in 1:nw
        subplot(SubColumn, SubRow, n)
        imshow(x[:, :, n] .* max(0., w[:, :, n]), origin ="lower")
        titlestring = @sprintf("%2.4f µm", wvl[n] * 1e6)
        title(titlestring)
        xticks([])
        yticks([])
        if( n == (nw + 1 - SubRow + count_x) )
            xticks(pos - 1, round(Int, indpix[pos] * 100000) / 100)
            xlabel("FOV (mas)")
            count_x += 1
        end
        if(n == (1 + count_y * SubRow))
            yticks(pos - 1, round(Int, indpix[pos] * 100000) / 100)
            ylabel("FOV (mas)")
            count_y += 1
        end
    end
###################################################################################
# primal and dual residuals plot
###################################################################################
    figure(2)
    clf()
    label1 = "Primal Residual"
    label2 = "Dual Residual"
    label3 = "Primal Threshold"
    label4 = "Dual Threshold"
    figure(2)
    plt_1 = plot(crit1, color="blue", label = label1)
    plt_2 = plot(crit2, color="red" , label = label2)
    plt_3 = plot(eps1 * ones(crit1), linestyle = "--", color = "blue", label = label3)
    plt_4 = plot(eps2 * ones(crit2), linestyle = "--", color = "red" , label = label4)
    ax1 = gca()
    ax1[:set_yscale]("log")
    xlabel("Iteration")
    ylabel("Residual (semilog)")
    title( "Primal and Dual residual")
    legend([plt_1, plt_2, plt_3, plt_4],[label1, label2, label3, label4], loc = "upper right", fancybox = "true")
    grid("on")
    #  show()
end
###################################################################################
# PLOT TOOLS
# ---------------------------------------------------------------------------------
# Automatic Subplot index creation
# (SubColumn,SubRow)= SubIndexPlot(nw::Integer)
function createsubplotindex(nw::Integer)
    if(nw == 1)
        SubColumn = 1
        SubRow = 1
    elseif(nw == 2)
        SubColumn = 1
        SubRow = 2
    elseif(nw == 3)
        SubColumn = 1
        SubRow = 3
    elseif(nw == 4)
        SubColumn = 2
        SubRow = 2
    elseif(nw == 5)
        SubColumn = 2
        SubRow = 3
    elseif(nw == 6)
        SubColumn = 2
        SubRow = 3
    elseif(nw == 7)
        SubColumn = 2
        SubRow = 4
    elseif(nw == 8)
        SubColumn = 2
        SubRow = 4
    elseif(nw == 9)
        SubColumn = 3
        SubRow = 3
    elseif(nw == 10)
        SubColumn = 2
        SubRow = 5
    elseif(nw == 11)
        SubColumn = 3
        SubRow = 4
    elseif(nw == 12)
        SubColumn = 3
        SubRow = 4
    elseif(nw == 13)
        SubColumn = 3
        SubRow = 5
    elseif(nw == 14)
        SubColumn = 3
        SubRow = 5
    elseif(nw == 15)
        SubColumn = 3
        SubRow = 5
    elseif(nw == 16)
        SubColumn = 4
        SubRow = 4
    elseif(nw == 18)
        SubColumn = 3
        SubRow = 6
    else
        println(" Automatic subplot, can be")
        SubColumn = round(Int, sqrt(nw))
        SubRow = round(Int, .5+ sqrt(nw))
    end
    return int(SubColumn), int(SubRow)
end


text:
function klamp(t, a, b)
    if t <= a
        return a
    elseif b <= t
        return b
    else
        da = abs(t - a)
        db = abs(t - b)

        if da < db
            return a
        else
            return b
        end
        # return t
    end
end

function contact_point_box(ps, xc, qc, kx, ky, kz)
    # rotate point into box frame
    ps_box_frame = vector_rotate(ps - xc, inv(qc))

    # find box point in box frame
    px = clamp(ps_box_frame[1], -0.5 * kx, 0.5 * kx)
    py = clamp(ps_box_frame[2], -0.5 * ky, 0.5 * ky)
    pz = clamp(ps_box_frame[3], -0.5 * kz, 0.5 * kz)

    # box point in world frame
    return vector_rotate([px; py; pz], qc) + xc
end

function ∂contact_point_box∂p(ps, xc, qc, kx, ky, kz)

    FiniteDiff.finite_difference_jacobian(p -> contact_point_box(p, xc, qc, kx, ky, kz), ps)
end

function ∂contact_point_box∂x(ps, xc, qc, kx, ky, kz)

    FiniteDiff.finite_difference_jacobian(x -> contact_point_box(ps, x, qc, kx, ky, kz), xc)
end

function ∂contact_point_box∂q(ps, xc, qc, kx, ky, kz)

    FiniteDiff.finite_difference_jacobian(q -> contact_point_box(ps, xc, Quaternion(q..., false), kx, ky, kz), vector(qc))
end

# x_min = -1.0
# x_max = 1.0
# z_min = -1.0
# z_max = 1.0

# p = [-1.0; 1.0] #.- 1.0e-6

# function klamp(t, a, b)
#     if t <= a
#         return a
#     elseif b <= t
#         return b
#     else
#         da = abs(t - a)
#         db = abs(t - b)

#         if da < db
#             return a
#         else
#             return b
#         end
#     end
# end

# function nearest_point(p, x_min, x_max, z_min, z_max)
#     return [klamp(p[1], x_min, x_max), klamp(p[2], z_min, z_max)]
# end

# nearest_point(p, x_min, x_max, z_min, z_max)


text:
@testset "$TEST $G" begin

g10 = StarGraph(10)
c = core_periphery_deg(g10)
@test degree(g10,1) == 9
@test c[1] == 1
for i=2:10
    @test c[i] == 2
end

g10 = StarGraph(10)
g10 = blockdiag(g10,g10)
add_edge!(g10, 1, 11)
c = core_periphery_deg(g10)
@test c[1] == 1
@test c[11] == 1
for i=2:10
    @test c[i] == 2
end
for i=12:20
    @test c[i] == 2
end

end # testset


text:
include("networks/Network-Visualization.jl")
include("ElectionAlgorithm.jl")

function election_visualization!(fgp::FixedGraphPlot{T}, elec_alg::ElectionAlgorithm{T}, process_id_label_map::Dict{T,String}) where {T <: Integer}
    _election_nodecolors!(fgp, elec_alg, process_id_label_map)
    if has_winner(elec_alg)
        _winner_visualization!(fgp, elec_alg.winner.echo_algorithm_instance)
    end
    return fgp
end

function _election_nodecolors!(fgp::FixedGraphPlot{T}, elec_alg::ElectionAlgorithm{T}, process_id_label_map::Dict{T,String}) where {T <: Integer}
    fill!(fgp.nodecolors, colorant"white")
    active_initiator_nodes = Set(echo_alg.initiator for echo_alg ∈ elec_alg.echo_algorithms)
    for initiator_node ∈ values(elec_alg.initiator_id_nodes_map)
        fgp.nodecolors[initiator_node] = initiator_node ∈ active_initiator_nodes ?  colorant"red" : colorant"blue"
    end
    for (i, process_id) ∈ enumerate(elec_alg.process_ids)
        strongest_id = elec_alg.strongest_ids[i]
        strongest_label = strongest_id == 0 ? "" : process_id_label_map[strongest_id]
        fgp.nodelabels[i] = "$(process_id_label_map[process_id]) $strongest_label"
    end
    return fgp
end

function _winner_visualization!(fgp::FixedGraphPlot, echo_alg::EchoAlgorithm)
    fgp.nodecolors[echo_alg.initiator] = colorant"green"
    for (i, j) in enumerate(echo_alg.activation_edges)
        if j != 0
            edge_index = fgp.edgeindices[Edge(i, j)]
            fgp.edgecolors[edge_index] = colorant"green"
        end
    end
    return fgp
end

max_stars_repo_name: tomtuamnuq/dist-algorithms

max_issues_repo_name: tomtuamnuq/dist-algorithms

max_forks_repo_name: tomtuamnuq/dist-algorithms


text:
#mpexperiment.jl
#Algorithm 4.2 of Random Eigenvalues by Alan Edelman

#Experiment:    Marcenko-Pastur for applications
#Plot:          Histogram eigenvalues of the covariance matrix
#Theory:        Marcenko-Pastur as n->Infinity

## Parameters
t = 100         # trials
r = 1           # aspect ratio
n = 100         # matrix column size
dx = 0.1        # binsize

function mp_experiment(n,r,t,dx)
    m = iround(n/r)
    ## Experiment
    v = Float64[]          # eigenvalue samples
    for i = 1:t
        A = randn(m,n) + 4*sqrt(n)*diagm(((1:n).<10))
        A = A + sqrt(n) * diagm((1:n).>(n-1)) * 3 #3+0.1*randn(n,1) #3/sqrt(n)
        append!(v, svd(A)[2]) # eigenvalues
    end
    v = v / sqrt(m) # normalized eigenvalues
    ## Theory
    a = 1-sqrt(r)
    b = 10
    x = a-dx/2:dx:b
    y = real(sqrt((x.^2-a^2).*(2^2-x.^2) + 0im)./(pi*x*r))
    return (hist(v, x), (x, y))
end
((grid, count), (x,y)) = mp_experiment(n,r,t,dx)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count/(t*n*dx), step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(x, y, "linewidth", 2, "color", "blue"))
last = length(count)-2
while count[last] == 0
    last -= 1
end
setattr(p, "xrange", (0, grid[last+2]))
setattr(p, "yrange", (-.1, ceil(max(count/(t*n*dx)))))
if isinteractive()
    Winston.display(p)
else
    file(p, "mpexperiment.png")
end


text:
using egads

width = 5.0
rad   = 1.0

NULL_E = ego(0)

(context, status) = EG_open( )
if (status < 0) error("Can't open, failure code: %i", status) end

# build cylinder 1
datacyl = [0.0, 0.0,    0.0,
           0.0, 0.0,  width,
           rad]
(ebody1, status) = EG_makeSolidBody(context, CYLINDER, datacyl)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# build cylinder 2
datacyl = [0.0, 0.0, 2*width,
           0.0, 0.0, 3*width,
           rad]
(ebody2,status) = EG_makeSolidBody(context, CYLINDER, datacyl)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# Translate cylinder 2
Δx = [0.0,0.0,-width]
matrix = [ 1.0, 0.0, 0.0, Δx[1],
           0.0, 1.0, 0.0, Δx[2],
           0.0, 0.0, 1.0, Δx[3] ]

exform_ptr = Ref{ego}()
status = EG_makeTransform(context, matrix, exform_ptr)
if (status != egads.EGADS_SUCCESS) egads.cleanup(status, context) end
(ebody22,status) = EG_copyObject(ebody2, exform_ptr[])
# join
toler = Cdouble(0.0)
(emodel, status) = EG_join( ebody1, ebody22, toler )
if (status != EGADS_SUCCESS) cleanup(status, context) end

# save model
status = EG_saveModel(emodel, "test.egads")
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_saveModel(emodel, "test.step")
if (status != EGADS_SUCCESS) cleanup(status, context) end

# delete
status = EG_deleteObject(ebody1)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(ebody2)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(ebody22)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(emodel)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# close everything
cleanup( context )


text:
mutable struct GsnTargetLaw{T}
    θ::Vector{Float64}
    P::T

    function GsnTargetLaw(μ, Σ=Matrix{Float64}(I, (length(μ), length(μ))))
        d = length(μ)
        θ = zeros(d*(d+1))
        θ[1:d] .= μ
        θ[d+1:d*d+d] .= vec(Σ)
        P = MvNormal(μ, Σ)
        new{typeof(P)}(θ, P)
    end
end

function set_parameters!(P::GsnTargetLaw, loc2glob_idx, θ)
    P.θ[loc2glob_idx] .= θ
    d = dim(P.P.Σ)
    μ, Σ = P.θ[1:d], reshape(P.θ[(d+1):end], (d,d))
    Σ = Symmetric(triu(Σ))
    P.P = MvNormal(μ, Σ)
end

function Distributions.loglikelihood(P::GsnTargetLaw, observs)
    ll = 0.0
    for obs in observs
        ll += logpdf(P.P, obs)
    end
    ll
end


text:
const utf8_example_text = """
∮ E⋅da = Q,  n → ∞, ∑ f(i) = ∏ g(i), ∀x∈ℝ: ⌈x⌉ = −⌊−x⌋, α ∧ ¬β = ¬(¬α ∨ β),

ℕ ⊆ ℕ₀ ⊂ ℤ ⊂ ℚ ⊂ ℝ ⊂ ℂ, ⊥ < a ≠ b ≡ c ≤ d ≪ ⊤ ⇒ (A ⇔ B),

2H₂ + O₂ ⇌ 2H₂O, R = 4.7 kΩ, ⌀ 200 mm

Linguistics and dictionaries:

  ði ıntəˈnæʃənəl fəˈnɛtık əsoʊsiˈeıʃn
  Y [ˈʏpsilɔn], Yen [jɛn], Yoga [ˈjoːgɑ]

APL:

  ((V⍳V)=⍳⍴V)/V←,V    ⌷←⍳→⍴∆∇⊃‾⍎⍕⌈

Nicer typography in plain text files:

  ╔══════════════════════════════════════════╗
  ║                                          ║
  ║   • ‘single’ and “double” quotes         ║
  ║                                          ║
  ║   • Curly apostrophes: “We’ve been here” ║
  ║                                          ║
  ║   • Latin-1 apostrophe and accents: '´`  ║
  ║                                          ║
  ║   • ‚deutsche‘ „Anführungszeichen“       ║
  ║                                          ║
  ║   • †, ‡, ‰, •, 3–4, —, −5/+5, ™, …      ║
  ║                                          ║
  ║   • ASCII safety test: 1lI|, 0OD, 8B     ║
  ║                      ╭─────────╮         ║
  ║   • the euro symbol: │ 14.95 € │         ║
  ║                      ╰─────────╯         ║
  ╚══════════════════════════════════════════╝

Greek (in Polytonic):

  The Greek anthem:

  Σὲ γνωρίζω ἀπὸ τὴν κόψη
  τοῦ σπαθιοῦ τὴν τρομερή,
  σὲ γνωρίζω ἀπὸ τὴν ὄψη
  ποὺ μὲ βία μετράει τὴ γῆ.

  ᾿Απ᾿ τὰ κόκκαλα βγαλμένη
  τῶν ῾Ελλήνων τὰ ἱερά
  καὶ σὰν πρῶτα ἀνδρειωμένη
  χαῖρε, ὦ χαῖρε, ᾿Ελευθεριά!

  From a speech of Demosthenes in the 4th century BC:

  Οὐχὶ ταὐτὰ παρίσταταί μοι γιγνώσκειν, ὦ ἄνδρες ᾿Αθηναῖοι,
  ὅταν τ᾿ εἰς τὰ πράγματα ἀποβλέψω καὶ ὅταν πρὸς τοὺς
  λόγους οὓς ἀκούω· τοὺς μὲν γὰρ λόγους περὶ τοῦ
  τιμωρήσασθαι Φίλιππον ὁρῶ γιγνομένους, τὰ δὲ πράγματ᾿
  εἰς τοῦτο προήκοντα,  ὥσθ᾿ ὅπως μὴ πεισόμεθ᾿ αὐτοὶ
  πρότερον κακῶς σκέψασθαι δέον. οὐδέν οὖν ἄλλο μοι δοκοῦσιν
  οἱ τὰ τοιαῦτα λέγοντες ἢ τὴν ὑπόθεσιν, περὶ ἧς βουλεύεσθαι,
  οὐχὶ τὴν οὖσαν παριστάντες ὑμῖν ἁμαρτάνειν. ἐγὼ δέ, ὅτι μέν
  ποτ᾿ ἐξῆν τῇ πόλει καὶ τὰ αὑτῆς ἔχειν ἀσφαλῶς καὶ Φίλιππον
  τιμωρήσασθαι, καὶ μάλ᾿ ἀκριβῶς οἶδα· ἐπ᾿ ἐμοῦ γάρ, οὐ πάλαι
  γέγονεν ταῦτ᾿ ἀμφότερα· νῦν μέντοι πέπεισμαι τοῦθ᾿ ἱκανὸν
  προλαβεῖν ἡμῖν εἶναι τὴν πρώτην, ὅπως τοὺς συμμάχους
  σώσομεν. ἐὰν γὰρ τοῦτο βεβαίως ὑπάρξῃ, τότε καὶ περὶ τοῦ
  τίνα τιμωρήσεταί τις καὶ ὃν τρόπον ἐξέσται σκοπεῖν· πρὶν δὲ
  τὴν ἀρχὴν ὀρθῶς ὑποθέσθαι, μάταιον ἡγοῦμαι περὶ τῆς
  τελευτῆς ὁντινοῦν ποιεῖσθαι λόγον.

  Δημοσθένους, Γ´ ᾿Ολυνθιακὸς

Georgian:

  From a Unicode conference invitation:

  გთხოვთ ახლავე გაიაროთ რეგისტრაცია Unicode-ის მეათე საერთაშორისო
  კონფერენციაზე დასასწრებად, რომელიც გაიმართება 10-12 მარტს,
  ქ. მაინცში, გერმანიაში. კონფერენცია შეჰკრებს ერთად მსოფლიოს
  ექსპერტებს ისეთ დარგებში როგორიცაა ინტერნეტი და Unicode-ი,
  ინტერნაციონალიზაცია და ლოკალიზაცია, Unicode-ის გამოყენება
  ოპერაციულ სისტემებსა, და გამოყენებით პროგრამებში, შრიფტებში,
  ტექსტების დამუშავებასა და მრავალენოვან კომპიუტერულ სისტემებში.

Russian:

  From a Unicode conference invitation:

  Зарегистрируйтесь сейчас на Десятую Международную Конференцию по
  Unicode, которая состоится 10-12 марта 1997 года в Майнце в Германии.
  Конференция соберет широкий круг экспертов по  вопросам глобального
  Интернета и Unicode, локализации и интернационализации, воплощению и
  применению Unicode в различных операционных системах и программных
  приложениях, шрифтах, верстке и многоязычных компьютерных системах.

Thai (UCS Level 2):

  Excerpt from a poetry on The Romance of The Three Kingdoms (a Chinese
  classic 'San Gua'):

  [----------------------------|------------------------]
    ๏ แผ่นดินฮั่นเสื่อมโทรมแสนสังเวช  พระปกเกศกองบู๊กู้ขึ้นใหม่
  สิบสองกษัตริย์ก่อนหน้าแลถัดไป       สององค์ไซร้โง่เขลาเบาปัญญา
    ทรงนับถือขันทีเป็นที่พึ่ง           บ้านเมืองจึงวิปริตเป็นนักหนา
  โฮจิ๋นเรียกทัพทั่วหัวเมืองมา         หมายจะฆ่ามดชั่วตัวสำคัญ
    เหมือนขับไสไล่เสือจากเคหา      รับหมาป่าเข้ามาเลยอาสัญ
  ฝ่ายอ้องอุ้นยุแยกให้แตกกัน          ใช้สาวนั้นเป็นชนวนชื่นชวนใจ
    พลันลิฉุยกุยกีกลับก่อเหตุ          ช่างอาเพศจริงหนาฟ้าร้องไห้
  ต้องรบราฆ่าฟันจนบรรลัย           ฤๅหาใครค้ำชูกู้บรรลังก์ ฯ

  (The above is a two-column text. If combining characters are handled
  correctly, the lines of the second column should be aligned with the
  | character above.)

Ethiopian:

  Proverbs in the Amharic language:

  ሰማይ አይታረስ ንጉሥ አይከሰስ።
  ብላ ካለኝ እንደአባቴ በቆመጠኝ።
  ጌጥ ያለቤቱ ቁምጥና ነው።
  ደሀ በሕልሙ ቅቤ ባይጠጣ ንጣት በገደለው።
  የአፍ ወለምታ በቅቤ አይታሽም።
  አይጥ በበላ ዳዋ ተመታ።
  ሲተረጉሙ ይደረግሙ።
  ቀስ በቀስ፥ ዕንቁላል በእግሩ ይሄዳል።
  ድር ቢያብር አንበሳ ያስር።
  ሰው እንደቤቱ እንጅ እንደ ጉረቤቱ አይተዳደርም።
  እግዜር የከፈተውን ጉሮሮ ሳይዘጋው አይድርም።
  የጎረቤት ሌባ፥ ቢያዩት ይስቅ ባያዩት ያጠልቅ።
  ሥራ ከመፍታት ልጄን ላፋታት።
  ዓባይ ማደሪያ የለው፥ ግንድ ይዞ ይዞራል።
  የእስላም አገሩ መካ የአሞራ አገሩ ዋርካ።
  ተንጋሎ ቢተፉ ተመልሶ ባፉ።
  ወዳጅህ ማር ቢሆን ጨርስህ አትላሰው።
  እግርህን በፍራሽህ ልክ ዘርጋ።

Runes:

  ᚻᛖ ᚳᚹᚫᚦ ᚦᚫᛏ ᚻᛖ ᛒᚢᛞᛖ ᚩᚾ ᚦᚫᛗ ᛚᚪᚾᛞᛖ ᚾᚩᚱᚦᚹᛖᚪᚱᛞᚢᛗ ᚹᛁᚦ ᚦᚪ ᚹᛖᛥᚫ

  (Old English, which transcribed into Latin reads 'He cwaeth that he
  bude thaem lande northweardum with tha Westsae.' and means 'He said
  that he lived in the northern land near the Western Sea.')

Braille:

  ⡌⠁⠧⠑ ⠼⠁⠒  ⡍⠜⠇⠑⠹⠰⠎ ⡣⠕⠌

  ⡍⠜⠇⠑⠹ ⠺⠁⠎ ⠙⠑⠁⠙⠒ ⠞⠕ ⠃⠑⠛⠔ ⠺⠊⠹⠲ ⡹⠻⠑ ⠊⠎ ⠝⠕ ⠙⠳⠃⠞
  ⠱⠁⠞⠑⠧⠻ ⠁⠃⠳⠞ ⠹⠁⠞⠲ ⡹⠑ ⠗⠑⠛⠊⠌⠻ ⠕⠋ ⠙⠊⠎ ⠃⠥⠗⠊⠁⠇ ⠺⠁⠎
  ⠎⠊⠛⠝⠫ ⠃⠹ ⠹⠑ ⠊⠇⠻⠛⠹⠍⠁⠝⠂ ⠹⠑ ⠊⠇⠻⠅⠂ ⠹⠑ ⠥⠝⠙⠻⠞⠁⠅⠻⠂
  ⠁⠝⠙ ⠹⠑ ⠡⠊⠑⠋ ⠍⠳⠗⠝⠻⠲ ⡎⠊⠗⠕⠕⠛⠑ ⠎⠊⠛⠝⠫ ⠊⠞⠲ ⡁⠝⠙
  ⡎⠊⠗⠕⠕⠛⠑⠰⠎ ⠝⠁⠍⠑ ⠺⠁⠎ ⠛⠕⠕⠙ ⠥⠏⠕⠝ ⠰⡡⠁⠝⠛⠑⠂ ⠋⠕⠗ ⠁⠝⠹⠹⠔⠛ ⠙⠑
  ⠡⠕⠎⠑ ⠞⠕ ⠏⠥⠞ ⠙⠊⠎ ⠙⠁⠝⠙ ⠞⠕⠲

  ⡕⠇⠙ ⡍⠜⠇⠑⠹ ⠺⠁⠎ ⠁⠎ ⠙⠑⠁⠙ ⠁⠎ ⠁ ⠙⠕⠕⠗⠤⠝⠁⠊⠇⠲

  ⡍⠔⠙⠖ ⡊ ⠙⠕⠝⠰⠞ ⠍⠑⠁⠝ ⠞⠕ ⠎⠁⠹ ⠹⠁⠞ ⡊ ⠅⠝⠪⠂ ⠕⠋ ⠍⠹
  ⠪⠝ ⠅⠝⠪⠇⠫⠛⠑⠂ ⠱⠁⠞ ⠹⠻⠑ ⠊⠎ ⠏⠜⠞⠊⠊⠥⠇⠜⠇⠹ ⠙⠑⠁⠙ ⠁⠃⠳⠞
  ⠁ ⠙⠕⠕⠗⠤⠝⠁⠊⠇⠲ ⡊ ⠍⠊⠣⠞ ⠙⠁⠧⠑ ⠃⠑⠲ ⠔⠊⠇⠔⠫⠂ ⠍⠹⠎⠑⠇⠋⠂ ⠞⠕
  ⠗⠑⠛⠜⠙ ⠁ ⠊⠕⠋⠋⠔⠤⠝⠁⠊⠇ ⠁⠎ ⠹⠑ ⠙⠑⠁⠙⠑⠌ ⠏⠊⠑⠊⠑ ⠕⠋ ⠊⠗⠕⠝⠍⠕⠝⠛⠻⠹
  ⠔ ⠹⠑ ⠞⠗⠁⠙⠑⠲ ⡃⠥⠞ ⠹⠑ ⠺⠊⠎⠙⠕⠍ ⠕⠋ ⠳⠗ ⠁⠝⠊⠑⠌⠕⠗⠎
  ⠊⠎ ⠔ ⠹⠑ ⠎⠊⠍⠊⠇⠑⠆ ⠁⠝⠙ ⠍⠹ ⠥⠝⠙⠁⠇⠇⠪⠫ ⠙⠁⠝⠙⠎
  ⠩⠁⠇⠇ ⠝⠕⠞ ⠙⠊⠌⠥⠗⠃ ⠊⠞⠂ ⠕⠗ ⠹⠑ ⡊⠳⠝⠞⠗⠹⠰⠎ ⠙⠕⠝⠑ ⠋⠕⠗⠲ ⡹⠳
  ⠺⠊⠇⠇ ⠹⠻⠑⠋⠕⠗⠑ ⠏⠻⠍⠊⠞ ⠍⠑ ⠞⠕ ⠗⠑⠏⠑⠁⠞⠂ ⠑⠍⠏⠙⠁⠞⠊⠊⠁⠇⠇⠹⠂ ⠹⠁⠞
  ⡍⠜⠇⠑⠹ ⠺⠁⠎ ⠁⠎ ⠙⠑⠁⠙ ⠁⠎ ⠁ ⠙⠕⠕⠗⠤⠝⠁⠊⠇⠲

  (The first couple of paragraphs of "A Christmas Carol" by Dickens)

Compact font selection example text:

  ABCDEFGHIJKLMNOPQRSTUVWXYZ /0123456789
  abcdefghijklmnopqrstuvwxyz £©µÀÆÖÞßéöÿ
  –—‘“”„†•…‰™œŠŸž€ ΑΒΓΔΩαβγδω АБВГДабвгд
  ∀∂∈ℝ∧∪≡∞ ↑↗↨↻⇣ ┐┼╔╘░►☺♀ ﬁ�⑀₂ἠḂӥẄɐː⍎אԱა

Greetings in various languages:

  Hello world, Καλημέρα κόσμε, コンニチハ

Box drawing alignment tests:                                          █
                                                                      ▉
  ╔══╦══╗  ┌──┬──┐  ╭──┬──╮  ╭──┬──╮  ┏━━┳━━┓  ┎┒┏┑   ╷  ╻ ┏┯┓ ┌┰┐    ▊ ╱╲╱╲╳╳╳
  ║┌─╨─┐║  │╔═╧═╗│  │╒═╪═╕│  │╓─╁─╖│  ┃┌─╂─┐┃  ┗╃╄┙  ╶┼╴╺╋╸┠┼┨ ┝╋┥    ▋ ╲╱╲╱╳╳╳
  ║│╲ ╱│║  │║   ║│  ││ │ ││  │║ ┃ ║│  ┃│ ╿ │┃  ┍╅╆┓   ╵  ╹ ┗┷┛ └┸┘    ▌ ╱╲╱╲╳╳╳
  ╠╡ ╳ ╞╣  ├╢   ╟┤  ├┼─┼─┼┤  ├╫─╂─╫┤  ┣┿╾┼╼┿┫  ┕┛┖┚     ┌┄┄┐ ╎ ┏┅┅┓ ┋ ▍ ╲╱╲╱╳╳╳
  ║│╱ ╲│║  │║   ║│  ││ │ ││  │║ ┃ ║│  ┃│ ╽ │┃  ░░▒▒▓▓██ ┊  ┆ ╎ ╏  ┇ ┋ ▎
  ║└─╥─┘║  │╚═╤═╝│  │╘═╪═╛│  │╙─╀─╜│  ┃└─╂─┘┃  ░░▒▒▓▓██ ┊  ┆ ╎ ╏  ┇ ┋ ▏
  ╚══╩══╝  └──┴──┘  ╰──┴──╯  ╰──┴──╯  ┗━━┻━━┛           └╌╌┘ ╎ ┗╍╍┛ ┋  ▁▂▃▄▅▆▇█
"""


text:
function array end

macro array(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    # array includes iterator by default, as this will almost always be what is expected
    ProxyInterfaces.iterator(p::$ProxyType) = array(p)
    ProxyInterfaces.iterator(p::Type{<:$ProxyType}) = array(p)
    @iterator($ProxyType)

    Base.size(p::$ProxyType) = Base.size(ProxyInterfaces.array(p))
    Base.size(p::$ProxyType, d) = Base.size(ProxyInterfaces.array(p), d)
    Base.getindex(p::$ProxyType, args...)	= Base.getindex(ProxyInterfaces.array(p), args...)
    Base.firstindex(p::$ProxyType) = Base.firstindex(ProxyInterfaces.array(p))
    Base.lastindex(p::$ProxyType)	= Base.lastindex(ProxyInterfaces.array(p))
    # optional
    Base.IndexStyle(T::Type{<:$ProxyType}) = Base.IndexStyle(ProxyInterfaces.array(T))
    Base.length(p::$ProxyType) = Base.length(ProxyInterfaces.array(p))
    Base.similar(p::$ProxyType, args...) = $ProxyType(Base.similar(ProxyInterfaces.array(p), args...))
    # non-traditional
    Base.axes(p::$ProxyType) = Base.axes(ProxyInterfaces.array(p))
    # strides
    Base.strides(p::$ProxyType, args...) = Base.strides(ProxyInterfaces.array(p), args...)
    # extra
    Base.isempty(p::$ProxyType) = Base.isempty(ProxyInterfaces.array(p))
    Base.in(a, p::$ProxyType) = Base.in(a, ProxyInterfaces.array(p))
  end
end

macro array_mutable(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    @array($ProxyType)
    Base.setindex!(p::$ProxyType, args...) = Base.setindex!(ProxyInterfaces.array(p), args...)
    Base.empty!(p::$ProxyType) = Base.empty!(ProxyInterfaces.array(p))
    Base.fill!(p::$ProxyType, x) = $ProxyType(Base.fill!(ProxyInterfaces.array(p), x))
  end
end


text:
function update!(s::Int,smax::Int,data::Array{Float64,1},nodes::Array{Float64,1},w::Float64,t::Float64,r::Float64)
    x = rand(data)
    dist = distance(nodes,x)
    bmu = bestMatchingUnit(dist)
    θ = neighbourhood(bmu,s,smax,nodes,w,t,r)
    for i in 1:length(nodes)
        nodes[i] += θ[i]*dist[i]
    end
end

function SOM(n::Int,smax::Int,data::Array{Float64,1};
            w::Float64=0.1,     # width scale of the Gaussian neighbourhood
            t::Float64=0.001,    # decay rate
            r::Float64=0.001)    # shrink rate

    nodes = maxentropy(n,data)

    w = maxspread(data)/length(nodes)*w

    for s in 1:smax
        update!(s,smax,data,nodes,w,t,r)
    end

    return nodes
end