text:
using Test
using Stripeline
using Dates
using Healpix
const Sl = Stripeline

eps = 3e-4 # Corresponds to pointing precision of 1.08 arcseconds

t_start = DateTime(2018, 01, 01, 0, 0, 0)

crab_az_stellarium_rad = 3.168070267969909
crab_alt_stellarium_rad = 1.4612458881566615

dirs = [π / 2 - crab_alt_stellarium_rad 2π - crab_az_stellarium_rad]

crab_ra_astropy_rad = 1.4596726619436968
crab_dec_astropy_rad = 0.3842255081802917
crab_position = sqrt(crab_ra_astropy_rad^2 + crab_dec_astropy_rad^2)

# Invert Crab coordinates into telescope pointing directions
groundq = telescopetoground(_->(0, deg2rad(20), 0), 0)
rotmatr = rotationmatrix_normalized(groundq)
vector = Healpix.ang2vec(dirs...)
dir = inv(rotmatr) * vector

# Compute skydirs
(skydirs, skyψ) = genpointings(_->(0, deg2rad(20), 0),
                               dir,
                               [0],
                               t_start,
                               latitude_deg = TENERIFE_LATITUDE_DEG,
                               longitude_deg = TENERIFE_LONGITUDE_DEG,
                               height_m = TENERIFE_HEIGHT_M)
crab_position_skydirs = sqrt(skydirs[1]^2 + skydirs[2]^2)
@test skydirs[1] ≈ crab_dec_astropy_rad atol = eps
@test skydirs[2] ≈ crab_ra_astropy_rad atol = eps
@test crab_position_skydirs ≈ crab_position atol = eps

t_1 = DateTime(2019, 02, 01, 2, 0, 0)
t_2 = DateTime(2019, 02, 02, 2, 0, 0)
t_3 = DateTime(2019, 02, 03, 2, 0, 0)

days = [t_1, t_2, t_3]

dirs = [π / 2 - 0.6174703397894518 2π - 4.852881197778698 ;
        π / 2 - 0.6024799007695448 2π - 4.859519751514131 ;
        π / 2 - 0.5875049757874335 2π - 4.866139397516001]

crab_ra_astropy_rad = 1.4596726619436968
crab_dec_astropy_rad = 0.3842255081802917
crab_position = sqrt(crab_ra_astropy_rad^2 + crab_dec_astropy_rad^2)

# Invert crab coordinates into telescope pointing directions
skydirs = Array{Float64}(undef, 3, 2)
for (idx, day) in enumerate(days)
    vector = Healpix.ang2vec(dirs[idx, 1], dirs[idx, 2])
    dir = inv(rotmatr) * vector

    (skydirections, skyψ) = genpointings(_->(0, deg2rad(20), 0),
                                         dir,
                                         [0],
                                         day,
                                         latitude_deg = TENERIFE_LATITUDE_DEG,
                                         longitude_deg = TENERIFE_LONGITUDE_DEG,
                                         height_m = TENERIFE_HEIGHT_M)

    skydirs[idx, 1] = skydirections[1]
    skydirs[idx, 2] = skydirections[2]
end

crab_position_skydirs = sqrt.(skydirs[:, 1].^2 + skydirs[:, 2].^2)

@test skydirs[1, 1] ≈ crab_dec_astropy_rad atol = eps
@test skydirs[1, 2] ≈ crab_ra_astropy_rad atol = eps
@test skydirs[2, 1] ≈ crab_dec_astropy_rad atol = eps
@test skydirs[2, 2] ≈ crab_ra_astropy_rad atol = eps
@test skydirs[3, 1] ≈ crab_dec_astropy_rad atol = eps
@test skydirs[3, 2] ≈ crab_ra_astropy_rad atol = eps
@test crab_position_skydirs[1] ≈ crab_position atol = eps
@test crab_position_skydirs[2] ≈ crab_position atol = eps
@test crab_position_skydirs[3] ≈ crab_position atol = eps

# Test flag `ground`
db = Sl.InstrumentDB()

time_duration = 86400
sampling_rate = 50.0
spin_velocity = 1

τ_s = 1 / sampling_rate
times = 0:τ_s:time_duration

(dirs, ψ) = genpointings(db.focalplane["I0"].orientation,
                         times;
                         latitude_deg = TENERIFE_LATITUDE_DEG) do time_s
    (0,
                              deg2rad(20.0),
                              Sl.timetorotang(time_s, spin_velocity))
end

expected_nsamples = convert(Int, time_duration * sampling_rate + 1)
@test size(dirs) == (expected_nsamples, 2)
@test size(ψ) == (expected_nsamples,)

(dirs, ψ) = genpointings(db.focalplane["I0"].orientation,
                         times;
                         ground = true,
                         latitude_deg = TENERIFE_LATITUDE_DEG) do time_s
    (0,
                              deg2rad(20.0),
                              Sl.timetorotang(time_s, spin_velocity))
end

expected_nsamples = convert(Int, time_duration * sampling_rate + 1)
@test size(dirs) == (expected_nsamples, 4)
@test size(ψ) == (expected_nsamples, 2)

# Check that all the values in the 3rd column are the same
@test dirs[1:end - 1, 3] ≈ dirs[2:end, 3]

# We strieve for an accuracy of a few arcseconds
@test rad2deg(dirs[1, 3]) ≈ 20.0 atol = 1e-3

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

# Check that the sign of the wheel angles is correct

let defaultdb = InstrumentDB()
    wheelfn(time_s) = (0, deg2rad(20), Stripeline.timetorotang(time_s, 1))

    # We do not start from t = 0s because we want the longitude to be positive
    timerange = 7200:1:(7200 + 60)
    G0_vec = defaultdb.focalplane["G0"].orientation
    V0_vec = defaultdb.focalplane["V0"].orientation

    dirG0, psiG0 = Stripeline.genpointings(wheelfn, G0_vec, timerange)
    dirV0, psiV0 = Stripeline.genpointings(wheelfn, V0_vec, timerange)

    # We expect G0 to draw larger circles in the sky
    @test minimum(dirG0[:, 1]) < minimum(dirV0[:, 1])
    @test minimum(dirG0[:, 1]) < minimum(dirV0[:, 1])
    @test maximum(dirG0[:, 1]) > maximum(dirV0[:, 1])
    @test maximum(dirG0[:, 2]) > maximum(dirV0[:, 2])
end


text:
using CairoMakie
using ElectronDisplay
using FFTW
using LinearAlgebra
using Scapin.Elasticity
using Scapin.Bri17

const T = Float64
const d = 2

C = Hooke{d,T}(1.0, 0.3)
α = (0.25, 0.25) # Fraction of the domained that is polarized
N_coarse = (4, 4)
r_max = 9
N_fine = (2^r_max) .* N_coarse

results = Dict()

for r ∈ 0:r_max
    N = (2^r) .* N_coarse
    𝒩 = CartesianIndices(N)
    h = 1.0 ./ N
    Γ̂ = DiscreteGreenOperatorBri17{d,T}(C, N, h)

    τ = zeros(T, 3, N...)
    τ[3, fill(1:2^r, d)...] .= 1
    τ̂ = fft(τ, 2:(d+1))
    ε̂ = Array{eltype(τ̂)}(undef, size(τ̂)...)

    for n ∈ 𝒩
        apply!(view(ε̂, :, n), Γ̂, n, τ̂[:, n])
    end

    ε = real.(ifft(ε̂, 2:(d+1)))

    ε_fine = zeros(T, 3, N_fine...)
    s = 2^(r_max - r)
    for n ∈ map(Tuple, 𝒩)
        n₁ = CartesianIndex(s .* (n .- 1) .+ 1)
        n₂ = CartesianIndex(s .* n)
        ε_fine[:, n₁:n₂] .= ε[:, n...]
    end

    # fig, ax, hm = heatmap(ε_fine[3, :, :])
    # save("eps_xy-$(N[1])x$(N[2]).png", fig)
    results[r] = ε_fine
end

x = Array{T}(undef, r_max)
y = Array{T}(undef, r_max)
for r ∈ 0:(r_max-1)
    ε₁ = results[r]
    ε₂ = results[r_max]
    x[r+1] = 2^r * N_coarse[1]
    y[r+1] = norm(ε₂-ε₁)
end

C = x[end] * y[end]
fig = scatter(x, y,  axis = (xscale=log10, yscale = log10))
lines!(x, C ./ x)
#save("convergence.png", fig)
electrondisplay(fig)

ext: jl


text:
# Example parser for a very basic Julia-like language of expressions, calls and
# function definitions.

using JuliaSyntax: @K_str, is_literal, is_keyword, is_operator

function parse_toplevel(st)
    mark = position(st)
    while true
        bump_trivia(st, skip_newlines=true)
        if peek(st) == K"EndMarker"
            break
        end
        parse_statement(st)
    end
    emit(st, mark, K"toplevel")
end

function parse_statement(st)
    mark = position(st)
    if peek(st) == K"function"
        parse_function_def(st)
    else
        parse_assignment(st)
    end
end

function parse_function_def(st)
    mark = position(st)
    @assert peek(st) == K"function"
    bump(st, TRIVIA_FLAG)
    parse_call(st)
    parse_block(st, K"end")
    emit(st, mark, K"function")
end

function parse_block(st, closing_kind, mark=position(st))
    while true
        bump_trivia(st, skip_newlines=true)
        if peek(st) == closing_kind
            bump(st, TRIVIA_FLAG)
            break
        elseif peek(st) == K"EndMarker"
            emit_diagnostic(st, error="Unexpected end of input")
            break
        end
        parse_assignment(st)
    end
    emit(st, mark, K"block")
end

function parse_assignment(st)
    mark = position(st)
    parse_expression(st)
    if peek(st) == K"="
        bump(st, TRIVIA_FLAG)
        parse_expression(st)
        emit(st, mark, K"=")
    end
end

function parse_expression(st)
    mark = position(st)
    parse_term(st)
    while peek(st) in (K"+", K"-")
        bump(st)
        parse_term(st)
        emit(st, mark, K"call", INFIX_FLAG)
    end
end

function parse_term(st)
    mark = position(st)
    parse_call(st)
    while peek(st) in (K"*", K"/")
        bump(st)
        parse_call(st)
        emit(st, mark, K"call", INFIX_FLAG)
    end
end

function parse_call(st)
    mark = position(st)
    parse_atom(st)
    if peek(st) == K"("
        bump(st, TRIVIA_FLAG)
        need_comma = false
        while true
            k = peek(st)
            if need_comma && k == K","
                bump(st, TRIVIA_FLAG)
                k = peek(st)
                need_comma = false
            end
            if k == K")"
                bump(st, TRIVIA_FLAG)
                break
            elseif k == K"EndMarker"
                emit_diagnostic(st, error="Unexpected end of input")
                break
            elseif need_comma
                bump_invisible(st, K"error", TRIVIA_TOKEN, error="Expected a `,`")
            end
            parse_expression(st)
            need_comma = true
        end
        emit(st, mark, K"call")
    end
end

function parse_atom(st)
    bump_trivia(st, skip_newlines=true)
    mark = position(st)
    k = peek(st)
    if k == K"Identifier" || is_literal(k)
        bump(st)
    elseif k in (K"-", K"+")
        bump(st)
        parse_atom(st)
        emit(st, mark, K"call")
    elseif k == K"("
        bump(st, TRIVIA_FLAG)
        parse_expression(st)
        if peek(st) == K")"
            bump(st, TRIVIA_FLAG)
            # emit(st, mark, K"(")
        else
            bump_invisible(st, K"error", TRIVIA_FLAG,
                           error="Expected `)` following expression")
        end
    elseif k == K"begin"
        bump(st, TRIVIA_FLAG)
        parse_block(st, K"end", mark)
    else
        bump(st)
        emit(st, mark, K"error",
             error="Expected literal, identifier or opening parenthesis")
    end
end

function parse_and_show(production::Function, code)
    st = ParseStream(code)
    production(st)
    t = JuliaSyntax.build_tree(GreenNode, st, wrap_toplevel_as_kind=K"error")
    show(stdout, MIME"text/plain"(), t, code, show_trivia=true)
    if !isempty(st.diagnostics)
        println()
        for d in st.diagnostics
            JuliaSyntax.show_diagnostic(stdout, d, code)
        end
    end
    t
end

println()
println("Example good parse:")
parse_and_show(parse_toplevel,
               """
               function f(x, y)
                   z = x - y
                   begin
                       a
                       b
                   end
                   z * z
               end

               f(1,2)
               """)

println()
println("Example diagnostics:")
parse_and_show(parse_expression, "(x + a*y) * (b")

nothing


text:
using Distributions
using Statistics

@enum ScenarioLoopPlacement OUTER INNER

"""
    RandomVariable{T}

A RandomVariable can be "assigned" to different model parameters. Its
value can be used directly or applied by adding to, or multiplying by,
reference values for the indicated parameter. Note that the distribution
must be instantiated, e.g., call as RandomVariable(:foo, Normal()), not
RandomVariable(:foo, Normal).
"""
struct RandomVariable{T}
    name::Symbol
    dist::T

    function RandomVariable(name::Symbol, dist::T) where T
        self = new{T}(name, dist)
        return self
    end
end

Base.eltype(rv::RandomVariable) = eltype(rv.dist)

distribution(rv::RandomVariable) = rv.dist

abstract type PseudoDistribution end

"""     ReshapedDistribution
A pseudo-distribution that returns a reshaped array of values from the
stored distribution and dimensions.

Example:
    rd = ReshapedDistribution([5, 5], Dirichlet(25,1))
"""
struct ReshapedDistribution <: PseudoDistribution
    dims::Vector{Int}
    dist::Distribution
end

function Base.rand(rd::ReshapedDistribution, draws::Int=1)
    return [reshape(rand(rd.dist), rd.dims...) for i in 1:draws]
end

# SampleStore is a faux Distribution that implements base.rand()
# to yield stored values.
mutable struct SampleStore{T} <: PseudoDistribution
    values::Vector{T}       # generally Int or Float64
    idx::Int                # index of next value to return
    dist::Union{Nothing, Distribution, PseudoDistribution}  # original distribution, if any

    function SampleStore(values::Vector{T}, dist::Union{Nothing, Distribution, PseudoDistribution}=nothing) where T
        return new{T}(values, 1, dist)
    end
end

Base.eltype(ss::SampleStore) = eltype(ss.values)

#
# TBD: maybe have a different SampleStore subtype for values drawn from a dist
# versus those loaded from a file, which would be treated as immutable?
#

function Statistics.quantile(ss::SampleStore{T}, q::Float64) where T
    return quantile(sort(ss.values), q)
end

# TBD: This interpolates values between those in the vector. Is this reasonable?
# Probably shouldn't use correlation on values loaded from a file rather than
# from a proper distribution.
function Statistics.quantile(ss::SampleStore{T}, probs::AbstractArray) where T
    return quantile(sort(ss.values), probs)
end

Base.length(ss::SampleStore{T}) where T = length(ss.values)

Base.iterate(ss::SampleStore{T}) where T = iterate(ss.values)
Base.iterate(ss::SampleStore{T}, idx) where T = iterate(ss.values, idx)

struct TransformSpec
    paramname::Symbol
    op::Symbol
    rvname::Symbol
    dims::Vector{Any}

    function TransformSpec(paramname::Symbol, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T
        if ! (op in (:(=), :(+=), :(*=)))
            error("Valid operators are =, +=, and *= (got $op)")
        end

        return new(paramname, op, rvname, dims)
    end
end

"""
    Base.rand(s::SampleStore{T}, n::Int=1) where T

Pseudo `rand`` function that just returns the next `n` values from the vector of
stored values. Currently does not support cycling through results multiple times,
but you can `reset` the `SampleStore` to reuse it.
"""
function Base.rand(s::SampleStore{T}, n::Int=1) where T
    idx = s.idx
    s.idx = next_idx = idx + n
    return n > 1 ? s.values[idx:(next_idx - 1)] : s.values[idx]
end

function Base.reset(s::SampleStore{T}) where T
    s.idx = 1
    return nothing
end

abstract type AbstractSimulationData end

"""
    SimulationDef

Holds all the data that defines a simulation.
"""
mutable struct SimulationDef{T}
    rvdict::OrderedDict{Symbol, RandomVariable}
    translist::Vector{TransformSpec}
    savelist::Vector{Tuple{Symbol, Symbol}}
    nt_type::Any                    # a generated NamedTuple type to hold data for a single trial
    data::T                         # data specific to a given sensitivity analysis method
    payload::Any                    # opaque (to Mimi) data the user wants access to in callbacks

    function SimulationDef{T}(rvlist::Vector,
                           translist::Vector{TransformSpec},
                           savelist::Vector{Tuple{Symbol, Symbol}},
                           data::T) where T <: AbstractSimulationData
        self = new()
        self.rvdict = OrderedDict([rv.name => rv for rv in rvlist])
        self.translist = translist
        self.savelist = savelist

        names = (keys(self.rvdict)...,)
        types = [eltype(fld) for fld in values(self.rvdict)]
        self.nt_type = NamedTuple{names, Tuple{types...}}

        self.data = data
        self.payload = nothing

        return self
    end
end

"""
    SimulationInstance{T}

Holds all the data that defines simulation results.
"""
mutable struct SimulationInstance{T}
    trials::Int
    current_trial::Int
    current_data::Any               # holds data for current_trial when current_trial > 0
    sim_def::SimulationDef{T} where T <: AbstractSimulationData
    models::Vector{M} where M <: AbstractModel
    results::Vector{Dict{Tuple, DataFrame}}
    payload::Any

    function SimulationInstance{T}(sim_def::SimulationDef{T}) where T <: AbstractSimulationData
        self = new()
        self.trials = 0
        self.current_trial = 0
        self.current_data = nothing
        self.sim_def = deepcopy(sim_def)
        self.payload = deepcopy(self.sim_def.payload)

        # These are parallel arrays; each model has a corresponding results dict
        self.models = Vector{AbstractModel}(undef, 0)
        self.results = [Dict{Tuple, DataFrame}()]

        return self
    end
end

"""
    SimulationDef{T}() where T <: AbstractSimulationData

Allow creation of an "empty" SimulationDef instance.
"""
function SimulationDef{T}() where T <: AbstractSimulationData
    SimulationDef{T}([], TransformSpec[], Tuple{Symbol, Symbol}[], T())
end

"""
    set_payload!(sim_def::SimulationDef, payload)

Attach a user's `payload` to the `SimulationDef`. A copy of the payload object
will be stored in the `SimulationInstance` at run time so it can be
accessed in scenario and pre-/post-trial callback functions. The value
is not used by Mimi in any way; it can be anything useful to the user.
"""
set_payload!(sim_def::SimulationDef, payload) = (sim_def.payload = payload)

"""
    payload(sim_def::SimulationDef)

Return the `payload` value set by the user via `set_payload!()`.
"""
payload(sim_def::SimulationDef) = sim_def.payload

"""
    payload(sim_inst::SimulationInstance)

Return the copy of the `payload` value stored in the `SimulationInstance` set by the user via `set_payload!()`.
"""
payload(sim_inst::SimulationInstance) = sim_inst.payload

struct MCSData <: AbstractSimulationData end

const MonteCarloSimulationDef = SimulationDef{MCSData}
const MonteCarloSimulationInstance = SimulationInstance{MCSData}

struct SimIterator{NT, T}
    sim_inst::SimulationInstance{T}

    function SimIterator{NT, T}(sim_inst::SimulationInstance{T}) where {NT <: NamedTuple, T <: AbstractSimulationData}
        return new{NT, T}(sim_inst)
    end
end

function getdataframe(sim_inst::SimulationInstance, comp_name::Symbol, datum_name::Symbol; model::Int = 1)
    return sim_inst.results[model][(comp_name, datum_name)]
end

function Base.getindex(sim_inst::SimulationInstance, comp_name::Symbol, datum_name::Symbol; model::Int = 1)
    msg = "getindex method for `SimulationInstance` has been replaced with getdataframe function for consistency of return type, please use getdataframe instead"
    Base.depwarn("$msg, $(stacktrace())", :getindex)
    return sim_inst.results[model][(comp_name, datum_name)]
end


text:
# progress bar

type TwProgressData
    uiChannel::RemoteRef
    statusChannel::RemoteRef
    progress::Float64
    showProgress::Bool
    cursorPos::Int # where is the progress bar cursor
    text::UTF8String
    redrawTime::Float64
    statusTime::Float64
    startTime::Float64
    TwProgressData() = new( RemoteRef(), RemoteRef(), 0.0, true, 1, utf8(""), time(), time(), time() )
end

twGlobProgressData = TwProgressData()

function updateProgressChannel( status::Symbol, v::Any )
    global twGlobProgressData
    if isready( twGlobProgressData.statusChannel )
        take!( twGlobProgressData.statusChannel )
    end
    put!( twGlobProgressData.statusChannel, ( status, v ) )
end

function progressMessage( s::UTF8String )
    global twGlobProgressData
    st = :normal
    val = nothing
    if isready( twGlobProgressData.statusChannel )
        (st, val) = take!( twGlobProgressData.statusChannel )
        if st == :init
            st = :normal
        end
    end
    if typeof( val ) <: Dict && eltype( val ) <: (Symbol, Any)
        val[ :message ] = utf8(s)
        put!( twGlobProgressData.statusChannel, ( st, val ) )
    else
        if st != :error && st != :done
            put!( twGlobProgressData.statusChannel, (st, @compat Dict{Symbol,Any}( :message => utf8(s) ) ) )
        end
    end
end

function progressUpdate( n::Float64 )
    global twGlobProgressData
    st = :normal
    val = nothing
    if isready( twGlobProgressData.statusChannel )
        (st, val) = take!( twGlobProgressData.statusChannel )
        if st == :init
            st = :normal
        end
    end
    if typeof( val ) <: Dict && eltype( val ) <: (Symbol, Any)
        val[ :progress ] = n
        put!( twGlobProgressData.statusChannel, ( st, val ) )
    else
        if st != :error && st != :done
            put!( twGlobProgressData.statusChannel, (st, @compat Dict{Symbol,Any}( :progress => n ) ) )
        end
    end
end

# standalone panel
# as a subwin as part of another widget (see next function)
# w include title width, if it's shown on the left
# the function f takes no argument. It's started right-away
# the function can call
# * TermWin.progressMessage( s::UTF8String ) # make sure height can accommodate the content
# * TermWin.progressUpdate( n::Float64 ) # 0.0 <= n <= 1.0
function newTwProgress( scr::TwObj; height::Real=5, width::Real=40, posy::Any=:center,posx::Any=:center, box=true, title = utf8("") )
    global twGlobProgressData
    obj = TwObj( TwProgressData(), Val{ :Progress } )
    obj.data = twGlobProgressData
    obj.box = box
    obj.title = title
    obj.borderSizeV= box ? 1 : 0
    obj.borderSizeH= box ? 1 : 0
    link_parent_child( scr, obj, height,width,posy,posx )
    obj
end

function draw( o::TwObj{TwProgressData} )
    werase( o.window )
    if o.box
        box( o.window, 0,0 )
    end
    if !isempty( o.title ) && o.box
        mvwprintw( o.window, 0, (@compat round(Int, ( o.width - length(o.title) )/2 )), "%s", o.title )
    end
    starty = o.borderSizeV
    startx = o.borderSizeH
    viewContentHeight = o.height - o.borderSizeV * 2
    viewContentWidth  = o.width  - o.borderSizeH * 2

    wattron( o.window, COLOR_PAIR(15)) #white on blue for progress bar

    p = max( 0.0, min( 1.0, twGlobProgressData.progress ) )
    left = max( 0, (@compat round(Int, viewContentWidth * p )) - 1 )
    bar = repeat( string( '\U2592' ), left+1 ) * repeat( " ", viewContentWidth - left - 1 )
    mvwprintw( o.window, starty, startx, "%s", bar )
    wattroff( o.window, COLOR_PAIR(15))
    if twGlobProgressData.text != ""
        mvwprintw( o.window, starty+1,startx, "%s", twGlobProgressData.text )
    end
end

function inject( o::TwObj{TwProgressData}, token::Any )
    global twGlobProgressData
    t = time()
    dorefresh = (t - twGlobProgressData.redrawTime) > 1.0
    retcode = :got_it # default behavior is that we know what to do with it

    if token == :progressupdate
        (st, val) = take!( twGlobProgressData.statusChannel )
        twGlobProgressData.statusTime = t
        if st == :init
            twGlobProgressData.startTime = t
            twGlobProgressData.progress = 0.0
            twGlobProgressData.text = utf8( "init" * strftime( " %H:%M:%S", time() )  )
        elseif st == :error
            o.value = val
            retcode = :exit_nothing
        elseif st == :done
            o.value = val
            retcode = :exit_ok
        else
            if typeof( val ) <: Dict && eltype( val ) <: (Symbol,Any)
                if haskey( val, :message ) && typeof( val[:message] ) <: AbstractString
                    twGlobProgressData.text = utf8( val[ :message ] )
                end
                if haskey( val, :progress ) && typeof( val[ :progress ] ) == Float64
                    twGlobProgressData.progress = val[ :progress ]
                    #twGlobProgressData.text = utf8( string( st ) * strftime( " %H:%M:%S", time() ) )
                end
            end
            dorefresh = true
        end
    elseif token == :ctrl_k # kill, preemptively
    elseif token == :ctrl_p # pause, cooperatively
    else
        retcode = :pass # I don't know what to do with it
    end

    if dorefresh
        twGlobProgressData.redrawTime = t
        refresh(o)
    end

    return retcode
end


text:
import DataFrames
using Gtk.ShortNames
using Gtk
using Plots, LaTeXStrings
using CSV
pyplot(dpi=80)

# Create window
win = Window()
# Position of win [3 == center]
set_gtk_property!(win, :title,"PSD Analizer 1.0")
set_gtk_property!(win, :window_position,3)
#set_gtk_property!(win, :height_request,300)
set_gtk_property!(win, :accept_focus,true)

################################################################################
# Menu settings for file
menu1 = MenuItem("_File")
filemenu1 = Menu(menu1)
new = MenuItem("New")
push!(filemenu1, new)
open = MenuItem("Open")
push!(filemenu1, open)
save = MenuItem("Save")
push!(filemenu1, save)
push!(filemenu1, SeparatorMenuItem())
exit = MenuItem("Exit")
idnew = signal_connect(exit, :activate) do widget
  destroy(win)
end
push!(filemenu1, exit)
mb = MenuBar()
push!(mb, menu1)

# Menu setting for edit
menu2 = MenuItem("_Edit")
filemenu2 = Menu(menu2)
clearmenu = MenuItem("Clear")
idnew = signal_connect(clearmenu, :activate) do widget
  set_gtk_property!(setFactor, :label, "Set")
  set_gtk_property!(factor, :sensitive, true)
  set_gtk_property!(factor, :text, "conversion factor")

  set_gtk_property!(setNumFiles, :label, "Set")
  set_gtk_property!(numFiles, :sensitive, true)
  set_gtk_property!(numFiles, :text,"# of databases")
end
push!(filemenu2, clearmenu)
push!(mb, menu2)

# Menu setting for help
menu3 = MenuItem("_Help")
filemenu3 = Menu(menu3)
help = MenuItem("Help...")
push!(filemenu3, help)

about = MenuItem("About...")
idnew = signal_connect(about, :activate) do widget
  win3 = Window()
  gabout = Grid()
  push!(win3, gabout)
  set_gtk_property!(gabout, :column_spacing, 10)
  set_gtk_property!(gabout, :row_spacing, 10)
  set_gtk_property!(gabout, :margin_bottom, 10)
  set_gtk_property!(gabout, :margin_top, 30)
  set_gtk_property!(gabout, :column_homogeneous,true)
  # Position of win [3 == center]
  set_gtk_property!(win3, :title,"About")
  set_gtk_property!(win3, :window_position,3)
  set_gtk_property!(win3, :width_request,390)

  okabout = Button("Ok")
  signal_connect(okabout, :clicked) do widget
    destroy(win3)
  end
  label = Label("Hola")
  GAccessor.markup(label,"""PhD Kelvyn B. Sánchez Sánchez\nInstituto Tecnológico de Celaya\nkelvyn.baruc@gmail.com\n
  PhD Eusebio Bolanos Reynoso\nInstituto Tecnológico de Orizaba\neusebio.itorizaba@gmail.com""")
  GAccessor.justify(label,Gtk.GConstants.GtkJustification.CENTER)
  gabout[1:3,1] = label
  gabout[2,2] = okabout

  showall(win3)
end
push!(filemenu3, about)
push!(mb, menu3)

################################################################################
# General grid
g = Grid()
set_gtk_property!(g, :column_spacing, 10)
set_gtk_property!(g, :row_spacing, 10)
set_gtk_property!(g, :margin_bottom, 10)
set_gtk_property!(g, :margin_left, 10)
set_gtk_property!(g, :margin_right, 10)
#set_gtk_property!(g, :column_homogeneous, true)

g[1:3,1] = mb

################################################################################
# Frame for app
load = Frame("Load Data")
set_gtk_property!(load, :width_request,300)
set_gtk_property!(load, :height_request,400)
set_gtk_property!(load, :label_xalign, 0.50)

gload = Grid()
gscroll = Grid()
scrollwin = ScrolledWindow(gscroll)
set_gtk_property!(gload, :margin_top,10)
set_gtk_property!(gload, :margin_bottom,10)
set_gtk_property!(gload, :margin_left,10)
set_gtk_property!(gload, :margin_right,10)
set_gtk_property!(gload, :column_spacing, 10)
set_gtk_property!(gload, :row_spacing, 10)

# Table for data
datalist = ListStore(String, Int, Bool)

# Visual table
dataview = TreeView(TreeModel(datalist))
set_gtk_property!(dataview, :width_request, 280)
set_gtk_property!(dataview, :reorderable, true)

# Set selectable
selmodel = G_.selection(dataview)
set_gtk_property!(dataview, :width_request,280)
set_gtk_property!(dataview, :height_request,380)

# Define type of cell
rTxt = CellRendererText()
rTog = CellRendererToggle()

# Define the source of data
c1 = TreeViewColumn("Label", rTxt, Dict([("text",0)]))
c2 = TreeViewColumn("Data", rTxt, Dict([("text",1)]))
c3 = TreeViewColumn("Status", rTog, Dict([("active",2)]))

# Allows to select rows
for c in [c1, c2, c3]
  GAccessor.resizable(c, true)
end
set_gtk_property!(dataview, :enable_grid_lines, 3)
set_gtk_property!(dataview, :expand, true)

push!(dataview, c1, c2, c3)

gscroll[1,1] = dataview
gload[1:3,1] = scrollwin

# Defining bottons to manipulate the table
del = Button("Delete")

signal_connect(del, "clicked") do widget
  if hasselection(selmodel)
    global currentIt
    currentIt = selected(selmodel)
    index = Gtk.index_from_iter(datalist, currentIt)
    deleteat!(datalist, currentIt)
    DataFrames.deletecols!(datadiam,index)
  end
end

cleardata = Button("Clear")

signal_connect(cleardata, "clicked") do widget
  while length(datalist) > 0
    deleteat!(datalist,1)
    DataFrames.deletecols!(datadiam,1)
  end
end

add = Button("Add")
datadiam = DataFrames.DataFrame()
signal_connect(add, "clicked") do widget
  # Create window
  win2 = Window()
  # Position of win [3 == center]
  set_gtk_property!(win2, :title,"Add data")
  set_gtk_property!(win2, :window_position,3)
  set_gtk_property!(win2, :height_request,190)
  #set_gtk_property!(win2, :height_request,200)
  #set_gtk_property!(win2, :accept_focus,true)
  loadframe = Frame()

  gadd = Grid()
  set_gtk_property!(gadd, :margin_top,50)
  set_gtk_property!(gadd, :margin_left,10)
  set_gtk_property!(gadd, :margin_right,10)
  #set_gtk_property!(gadd, :column_homogeneous, true)
  set_gtk_property!(gadd, :column_spacing, 10)
  set_gtk_property!(gadd, :row_spacing, 10)

  labeldata = Entry()
  set_gtk_property!(labeldata, :tooltip_markup, "Enter the label")
  set_gtk_property!(labeldata, :width_request, 200)
  set_gtk_property!(labeldata, :text,"label")

  canceldata = Button("Cancel")
  signal_connect(canceldata, :clicked) do widget
    destroy(win2)
  end

  browsedata = Button("Browse")
  signal_connect(browsedata, :clicked) do widget
    global dlg
    dlg = open_dialog("Choose file...", win2, ("*.txt, *.csv",), select_multiple=false)
  end

  addtodata = Button("Add")
  signal_connect(addtodata, :clicked) do widget
    global dlg
    global datadiam
    label1 = get_gtk_property(labeldata, :text, String)
    data = CSV.read(dlg, datarow=1)
    datadiam[Symbol(label1)] = data[1]
    # save to global data
    ll = size(datadiam[Symbol(label1)],1)
    status = datadiam[Symbol(label1)] != String
    push!(datalist,(label1,ll,status))
    destroy(win2)
  end

  gadd2 = Grid()
  set_gtk_property!(gadd2, :column_spacing, 10)

  gadd[1:2,1] = gadd2
  gadd2[1,1] = labeldata
  gadd2[2,1] = browsedata
  gadd[1,2] = canceldata
  gadd[2,2] = addtodata

  push!(loadframe,gadd)
  push!(win2,loadframe)
  showall(win2)
end

gload[1,3] = del
gload[2,3] = cleardata
gload[3,3] = add
push!(load,gload)

################################################################################
setting = Frame("Setting Analysis")
set_gtk_property!(setting, :width_request,300)
set_gtk_property!(setting, :height_request,400)
set_gtk_property!(setting, :label_xalign, 0.50)

gsetting = Grid()
set_gtk_property!(gsetting, :margin_top,10)
set_gtk_property!(gsetting, :margin_bottom,10)
set_gtk_property!(gsetting, :margin_left,10)
set_gtk_property!(gsetting, :margin_right,10)
set_gtk_property!(gsetting, :column_spacing, 10)
set_gtk_property!(gsetting, :column_homogeneous, true)
set_gtk_property!(gsetting, :row_spacing, 30)

factor = Entry()
set_gtk_property!(factor, :tooltip_markup, "Enter the conversion factor")
set_gtk_property!(factor, :width_request, 100)

setFactor = Button("Set")
set_gtk_property!(setFactor, :width_request, 30)
set_gtk_property!(setFactor, :xalign, 0.50)
idnew = signal_connect(setFactor, :clicked) do widget
  msgFactor = get_gtk_property(factor, :text, String)

  # Check for non a number
  try
    global numFactor
    numFactor = parse(Float64,msgFactor)

    # Check for number > 0
    try
      if numFactor ≤ 0
        lolol
      end
      set_gtk_property!(setFactor, :label, "Set ✔")
      set_gtk_property!(factor, :sensitive, false)
    catch
      warn_dialog("Oops!... Factor must be higher than 0", win)
    end
  catch
    warn_dialog("Oops!... Please write a number", win)
  end
end
signal_connect(win, "key-press-event") do widget, event
  if event.keyval == 65293
    msgFactor = get_gtk_property(factor, :text, String)
    # Check for non a number
    try
      global numFactor
      numFactor = parse(Float64,msgFactor)

      # Check for number > 0
      try
        if numFactor ≤ 0
          lolol
        end
        set_gtk_property!(setFactor, :label, "Set ✔")
        set_gtk_property!(factor, :sensitive, false)
      catch
        warn_dialog("Oops!... Factor must be higher than 0", win)
      end
    catch
      warn_dialog("Oops!... Please write a number", win)
    end
  end
end

distr = Frame()
choices = ["% Number", "% Length", "% Surface", "% Volume"]
f = Gtk.GtkBox(:v)
r = Vector{RadioButton}(undef, 4)
r[1] = RadioButton(choices[1])
push!(f,r[1])
r[2] = RadioButton(r[1],choices[2])
push!(f,r[2])
r[3] = RadioButton(r[2],choices[3])
push!(f,r[3])
r[4] = RadioButton(r[3],choices[4])
push!(f,r[4])
push!(distr,f)

labeldistr = Label("Select distribution:")
labelfactor = Label("Enter the conversion factor:")

gfactor = Grid()
set_gtk_property!(gfactor, :column_spacing, 10)
set_gtk_property!(gfactor, :row_spacing, 10)

over = Frame()
gover = Grid()
set_gtk_property!(gover, :margin_top,6)
set_gtk_property!(gover, :margin_bottom,6)
set_gtk_property!(gover, :margin_left,6)
set_gtk_property!(gover, :margin_right,6)
set_gtk_property!(gover, :column_homogeneous, true)
oversize = CheckButton("Oversize")
undersize = CheckButton("Undersize")
set_gtk_property!(oversize,:active,false)
set_gtk_property!(undersize,:active,false)
gover[1,1] = oversize
gover[2,1] = undersize
push!(over, gover)

gfactor[1:2,1] = labelfactor
gfactor[1,2] = factor
gfactor[2,2] = setFactor
gsetting[1:2,1] = gfactor
gsetting[1,3] = labeldistr
gsetting[2,3] = distr
gsetting[1:2,4] = over
push!(setting,gsetting)

###############################################################################
# Results Section
plots = Frame("Results & Plots")
set_gtk_property!(plots, :width_request,300)
set_gtk_property!(plots, :height_request,400)
set_gtk_property!(plots, :label_xalign, 0.50)

gplot = Grid()
set_gtk_property!(gplot, :column_spacing, 10)  # introduce a 15-pixel gap between columns
set_gtk_property!(gplot, :row_spacing, 10)  # introduce a 15-pixel gap between columns
set_gtk_property!(gplot, :column_homogeneous, true)

set_gtk_property!(gplot, :margin_top,10)
set_gtk_property!(gplot, :margin_bottom,10)
set_gtk_property!(gplot, :margin_left,10)
set_gtk_property!(gplot, :margin_right,10)

gtkimg = Image()
set_gtk_property!(gtkimg, :width_request,250)
set_gtk_property!(gtkimg, :height_request,200)
set_gtk_property!(gtkimg, :xpad,10)
set_gtk_property!(gtkimg, :ypad,10)
#push!(frameimg, gtkimg)
gplot[1:3,2] = gtkimg

exit = Button("Clear")
clear = Button("Drawn")
run = Button("Export")
gplot[1,3] = exit
gplot[2,3] = clear
gplot[3,3] = run

push!(plots, gplot)

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

g[1,2] = load
g[2,2] = setting
g[3,2] = plots

################################################################################
# Main bottons
exit = Button("Exit")
#set_gtk_property!(exit, :width_request, 300)
idnew = signal_connect(exit, :clicked) do widget
  destroy(win)
end

clear = Button("Clear")
#set_gtk_property!(clear, :width_request, 300)
idnew = signal_connect(clear, :clicked) do widget
  set_gtk_property!(setFactor, :label, "Set")
  set_gtk_property!(factor, :sensitive, true)
  set_gtk_property!(factor, :text, "")
  while length(datalist) > 0
    deleteat!(datalist,1)
    DataFrames.deletecols!(datadiam,1)
  end
  DataFrames.deletecols!(datadiam, :)
  empty!(gtkimg)
end

run = Button("Run")
idnew = signal_connect(run, :clicked) do widget
  empty!(gtkimg)
  function px2dtc(diam; factor = 0.4355)
    diameter = sqrt.(diam) / factor
    num = zeros(length(diameter),1)
    den = zeros(length(diameter),1)

    sum_num_D = 0
    sum_den_D = 0
    sum_num_S = 0
    sum_den_S = 0

    DTC = zeros(4,2)
    for m = 1:4
      # Cálculo de D[m,m-1]
      sum_num = 0
      sum_den = 0

      for i = 1:length(diameter)
        sum_num = diameter[i]^m + sum_num
      end

      for i = 1:length(diameter)
        sum_den = diameter[i]^(m-1) + sum_den
      end

      # Se guarda el valor del diametro
      DTC[m,1] = sum_num / sum_den

      # Cálculo de S[m,m-1]
      sum_num = 0
      sum_den = 0

      for i = 1:length(diameter)
        sum_num = (diameter[i]^(m-1) * (diameter[i] - DTC[m,1])^2) + sum_num
      end

      for i = 1:length(diameter)
        sum_den = diameter[i]^(m-1) + sum_den
      end

      # Se guarda el valor de la desviación
      DTC[m,2] = sqrt(sum_num / sum_den)

    end
    return DTC
  end

  global numFactor
  global datadiam
  global PSD
  PSD = []
  for i=1:size(datadiam,2)
    psd = px2dtc(datadiam[i], factor = numFactor)
    push!(PSD,psd)
  end

  global PSDTotal

  PSDTotal = []
  for j=1:size(datadiam,2)
    PSDInd = zeros(4,2)
    for i=1:4
      # Solve for lognorm parameters
      x0 = 6 # Initial value for D43
      y0 = 0.6 # Initial value for S43
      e = 1 # Error
      while e > 0.0001
        df1x = exp(x0+0.5*y0^2)
        df1y = y0*exp(x0+0.5*y0^2)
        df2x = 2.0*exp(2*x0+y0^2)*(exp(y0^2)-1.0)
        df2y = 2*y0*exp(2*x0+y0^2)*(exp(y0^2)-1)+2*exp(2.0*x0+y0^2)*y0*exp(y0^2)
        f1 = exp(x0+0.5*y0^2) - PSD[j][i,1]
        f2 = (exp(2*x0+y0^2))*(exp(y0^2)-1.0)- PSD[j][i,2]
        A = [df1x df1y; df2x df2y]
        b = [-f1; -f2]
        hj = A\b
        x1 = x0 + hj[1]
        y1 = y0 + hj[2]
        e = ((x1-x0)^2+(y1-y0)^2)^0.5
        x0 = x1
        y0 = y1

      end
      PSDInd[i,1] = x0
      PSDInd[i,2] = y0
    end
    push!(PSDTotal, PSDInd)
  end

  function lognorm(mean, sd; x = 0:1.04:1000)
    fd = zeros(length(x), 1)

    for i=1:length(x)
      fd[i, 1] = (1/(sd*sqrt(2*pi)))*exp(-((log(x[i]))-mean)^2/(2*sd^2))
    end

    density = (fd.*100)./sum(fd);
    return x, density
  end

  plt = []
  set_gtk_property!(gtkimg, :file, "")

  global choicedistr
  choisedistr = []
  for i=1:4
    if get_gtk_property(r[i], :active, Bool) == true
      global choicedistr
      choicedistr = i
    end
  end

  global y, x
  x = []
  y = []
  for i=1:length(datalist)
    global y, x, yy
    x, yy = lognorm(PSDTotal[i][choicedistr,1], PSDTotal[i][choicedistr,2], x = 0:10:1000)
    push!(y,yy)

    plt = plot!(x,y[i], label = datalist[i,1], lw=1.5,
    xlabel = "Size (μm)", ylabel=choices[choicedistr],
    xtickfont=font(12),ytickfont=font(12),guidefont=font(12), legendfontsize=10,
    framestyle=:box, xlim=(0,maximum(x)),ylim=(0,100))

    if length(datalist) == 1
      #################################################################
      # Initialization for vectors
      global y_under, y_over, yy
      y_under = zeros(length(yy),1)
      y_over = zeros(length(yy),1)
      y_under[1,1] = yy[1,1]
      y_over[1,1] = 100 - y_under[1,1]

      #################################################################
      # Oversize & Undersize calculation
      for l = 2:length(yy)
        # Undersize
        y_under[l,1] = y_under[l-1,1] + yy[l]
        # Oversize
        y_over[l,1] = 100 - y_under[l,1]
      end
      overs = get_gtk_property(oversize,:active,Bool)
      unders = get_gtk_property(undersize,:active,Bool)

      if overs == true
        plot!(x,y_over,ls=:dash, lw=1.5, label = "Oversize")
      end

      if unders == true
        plot!(x,y_under,ls=:dot, lw=1.5, label = "Undersize")
      end
    end
  end

  savefig(plt,"PSD.png")
  set_gtk_property!(gtkimg, :file, "PSD.png")
  rm("PSD.png")

end
#set_gtk_property!(run, :width_request, 300)
g[1,3] = exit
g[2,3] = clear
g[3,3] = run

#g[1:3,3] = g2
push!(win,g)
showall(win)


text:
using PlotThemes
using Colors
import PlotThemes._themes
# See:
#    https://medialab.github.io/iwanthue/
# and:
#    https://chir.ag/projects/name-that-color/
const cornflowerblue = colorant"cornflowerblue"; export cornflowerblue
const crimson = colorant"crimson"; export crimson
const cucumber = colorant"#77ab58"; export cucumber
const copper = colorant"#c37940"; export copper
const fuchsia = colorant"#be5ea7"; export fuschia
const dark = colorant"#282C34"; export dark

dark_palette = [
    cornflowerblue,
    crimson,
    cucumber,
    copper
]

# Not actually transparent, for some reason that default doesn't propagate
_themes[:transparentdark] = PlotThemes.PlotTheme(
    bg = colorant"#282C34",
    bginside = colorant"#282C34",
    fg = colorant"#FFFFFF",
    fg_axis = colorant"#FFFFFF",
    foreground_color_border = colorant"#FFFFFF",
    foreground_color_text = colorant"#FFFFFF",
    fgtext = colorant"#FFFFFF",
    fgguide = colorant"#FFFFFF",
    fglegend = colorant"#FFFFFF",
    legendfontcolor = colorant"#FFFFFF",
    legendtitlefontcolor = colorant"#FFFFFF",
    titlefontcolor = colorant"#FFFFFF",
    framestyle = :box,
    palette = PlotThemes.expand_palette(RGB(0, 0, 0), dark_palette; lchoices = [57], cchoices = [100]),
    colorgradient = :viridis
)


text:
using Test
using Liblet

# Operators

test_grammar_eq() = Grammar("S -> A B | B\nA -> a\nB -> b") == Grammar("S -> B | A B\nB -> b\nA -> a")

test_grammar_hash() = begin
    S = Dict(
        Grammar("S -> A B | B\nA -> a\nB -> b") => 1,
        Grammar("S -> B | A B\nB -> b\nA -> a") => 2
    )
    length(S) == 1
end

test_grammar_nondisjoint() = Grammar(Set(["S", "A"]),Set(["A", "a"]), [], "S")

test_grammar_wrongstart() = Grammar(Set(["T", "A"]),Set(["a"]), [], "S")

test_grammar_cf() = Grammar("S -> T U\nT -> t\nU -> u").iscontextfree

test_grammar_not_cf() = Grammar("S -> T U\nT x -> t\nT U -> u", false).iscontextfree

test_grammar_wrong_cf() = Grammar(Set(["S"]), Set(["s"]), [Production("S", ["s"]), Production("T", ["s"])], "S")

test_grammar_extrasymbol() = Grammar(Set(["S"]), Set(["s"]), [Production("S", ["s"]), Production(["S"], ["s", "t"])], "S")

test_grammar_from_to_string() = begin
    G = Grammar("""
        Z -> E
        E -> T | E + T
        T -> i | ( E )
    """)
    s = "" # to finish
    s == string(G)
end

# restrict

test_grammar_restrict_to() = begin
    G = Grammar("""
        Z -> E
        E -> T | E + T | X
        X -> x
        T -> i | ( E )
    """)
    Gr = Grammar("""
        Z -> E
        E -> T | E + T
        T -> i | ( E )
    """)
    symbols = G.T ∪ G.N - Set(["X", "x"])
    Gt = restrict(G, symbols)
    Gt == Gr
end

test_grammar_restrict_to_no_start() = begin
    G = Grammar("""
    Z -> E
    E -> T | E + T | X
    X -> x
    T -> i | ( E )
    """)
    restrict(G, (G.T ∪ G.N) - Set(["Z"]))
end

# alternative

test_alternatives() = begin
    G = Grammar("""
    Z -> E k
    E -> T | E + T
    T -> i | ( E )
    """)
    expected = Set([["T"], ["E", "+", "T"]])
    Set(alternatives(G, "E")) == expected
end

# Tests

function rungrammartests()
    @testset "grammar_tests" begin
        @test test_grammar_eq()
        #@test test_grammar_hash()
        @test_throws ArgumentError test_grammar_nondisjoint()
        @test_throws ArgumentError test_grammar_wrongstart()
        @test test_grammar_cf()
        @test !test_grammar_not_cf()
        #@test test_grammar_from_to_string()
        @test_throws ArgumentError test_grammar_wrong_cf()
        @test_throws ArgumentError test_grammar_extrasymbol()

        # restrict
        @test test_grammar_restrict_to()
        @test_throws ArgumentError test_grammar_restrict_to_no_start()

        # alternative
        @test test_alternatives()
    end
end

rungrammartests()


text:
struct Senseval3{S}
    filepaths::Vector{S}
end

function init_datadeps(::Type{Senseval3})
    register(DataDep("Senseval 3",
        """
        Website: http://web.eecs.umich.edu/~mihalcea/senseval/senseval3/index.html
        Data Website: http://web.eecs.umich.edu/~mihalcea/downloads.html#sensevalsemcor
        Prepared by: Rada Mihalcea, Timothy Chklovski, Adam Kilgarriff
        Current Maintainer: Rada Mihalcea

        Senseval-3 took place in March-April 2004, followed by a workshop held in July 2004 in Barcelona, in conjunction with ACL 2004.
        This is corpus was updated by Rada Mihalcea with the word senses according to WordNet 1.7.1, 2.0, 2.1, 3.0.

        The details of the task are in the following paper, please cite the following publication if you are using the corpora.

        Mihalcea, Rada, Timothy Chklovski, and Adam Kilgarriff. "The Senseval-3 English lexical sample task." In Proceedings of SENSEVAL-3, the third international workshop on the evaluation of systems for the semantic analysis of text. 2004.
        http://www.aclweb.org/anthology/W04-0807
        """,
        "http://web.eecs.umich.edu/~mihalcea/downloads/senseval.semcor/senseval3.semcor.tar.gz",
        "ab5a1230d1109ac7847b86babf442b2006276d73e6ee47d67f7f7c372c50d9eb";
        post_fetch_method = fn -> begin
            unpack(fn)
            innerdir = "senseval3.semcor"
            innerfiles = readdir(innerdir)
            # Move everything to current directory, under same name
            mv.(joinpath.(innerdir, innerfiles), innerfiles)
            rm(innerdir)
        end
    ))
end

function Senseval3(dirpath, ver=3.0)
    @assert(isdir(dirpath), dirpath)

    if "wordnet$ver" ∈ readdir(dirpath)
        paths = joinpath.(dirpath,"wordnet$ver", readdir(joinpath(dirpath, "wordnet$ver")))
    else
        paths = joinpath.(dirpath, readdir(dirpath))
    end
    Senseval3(paths)
end

Senseval3() = Senseval3(datadep"Senseval 3")

MultiResolutionIterators.levelname_map(::Type{Senseval3}) = [
    :doc=>1, :contextfile=>1, :context=>1, :document=>1,
    :sent=>2, :sentence=>2,
    :word=>3, :token=>3,
    :char=>4, :character=>4
    ]

function parse_sense_annotated_word_senseval3(line::AbstractString)
    captures = match(r"<wf cmd=done.* pos=(.*) lemma=(.*) wnsn=(.*) lexsn=(\d.*:\d*|[U]|NOANSWER|NOWN).*>(.*).*</wf>", line).captures
    pos, lemma, wnsn, lexsn, word = captures
    if ';' in wnsn
        # Discard Extra Senses
        wnsn = split(wnsn, ';') |> first
        lexsn = split(lexsn, ';') |> first
    end

    SenseAnnotatedWord(pos, lemma, parse(Int, wnsn), lexsn, word)
end

function parse_senseval3file(filename)
    local sent
    lines = @NestedVector(TaggedWord,2)()
    context = Document(intern(basename(filename)), lines)

    # structure

    function new_sentence(line)
        sent = @NestedVector(TaggedWord,1)()
        push!(lines, sent)
    end

    # words
    get_tagged(line) = push!(sent, parse_tagged_word(line))
    get_sense_annotated(line) = push!(sent, parse_sense_annotated_word_senseval3(line))
    get_punc(line) = push!(sent, parse_punc(line))

    # Parse
    subparsers = [
        "<wf cmd=tag"=>          get_tagged,
        "<wf cmd=ignore"=>       get_tagged,
        "<wf cmd=done"=>         line -> occursin("lemma=", line) ? get_sense_annotated(line) : get_tagged(line),
        "<punc>"=>               get_punc,
        "<context"=>             ignore,
        "</context" =>           ignore,
        "<s" =>                  new_sentence,
        "</s" =>                 ignore
    ]

    apply_subparsers(filename,subparsers)

    return context
end

function load(corpus::Senseval3, doc_buffersize=16)
    Channel(;ctype=Document{@NestedVector(TaggedWord, 2), String}, csize=doc_buffersize) do ch
        for fn in corpus.filepaths
            doc = parse_senseval3file(fn)
            put!(ch, doc)
        end
    end
end


text:
using KernelAbstractions, Test

function asynccopy_testsuite(backend, ArrayT)
    M = 1024

    A = ArrayT(rand(Float64, M))
    B = ArrayT(rand(Float64, M))

    a = Array{Float64}(undef, M)
    event = async_copy!(backend(), a, B, dependencies=Event(CPU()))
    event = async_copy!(backend(), A, a, dependencies=event)
    wait(event)

    @test isapprox(a, Array(A))
    @test isapprox(a, Array(B))
end


text:
import GreyDecision.Utility
import GreyDecision.GreyNumbers: GreyNumber

@testset "Utility" begin

    @testset "zeros" begin
        @testset "zero" begin
            @test zero(GreyNumber{Int64}) == GreyNumber(0, 0)
            @test zero(GreyNumber(0, 0)) == GreyNumber(0, 0)

            @test zero(GreyNumber{Float64}) == GreyNumber(0.0, 0.0)
            @test zero(GreyNumber(0.0, 0.0)) == GreyNumber(0.0, 0.0)
        end

        @testset "zeros vector" begin
            z = zeros(GreyNumber, 10)
            g = z[1]
            @test length(z) == 10
            @test g isa GreyNumber
            @test g.a isa Float64
            @test g.b isa Float64
        end

        @testset "zeros matrix" begin
            z = zeros(GreyNumber, 10, 3)
            @test size(z) == (10, 3)
            @test z isa Matrix
            @test z[1, 1] == GreyNumber(0.0, 0.0)
        end

        @testset "zeros matrix with tuple" begin
            z = zeros(GreyNumber, (20, 2))
            @test size(z) == (20, 2)
            @test z isa Matrix
            @test z[1, 1] == GreyNumber(0.0, 0.0)
        end

    end

    @testset "ones" begin
        @testset "ones" begin
            @test one(GreyNumber{Int64}) == GreyNumber(1, 1)
            @test one(GreyNumber(1, 1)) == GreyNumber(1, 1)

            @test one(GreyNumber{Float64}) == GreyNumber(1.0, 1.0)
            @test one(GreyNumber(1.0, 1.0)) == GreyNumber(1.0, 1.0)
        end

        @testset "ones vector" begin
            z = ones(GreyNumber, 10)
            g = z[1]
            @test length(z) == 10
            @test g isa GreyNumber
            @test g.a isa Float64
            @test g.b isa Float64
        end

        @testset "ones matrix" begin
            z = ones(GreyNumber, 10, 3)
            @test size(z) == (10, 3)
            @test z isa Matrix
            @test z[1, 1] == GreyNumber(1.0, 1.0)
        end

        @testset "ones matrix with tuple" begin
            z = ones(GreyNumber, (20, 2))
            @test size(z) == (20, 2)
            @test z isa Matrix
            @test z[1, 1] == GreyNumber(1.0, 1.0)
        end

    end

    @testset "Euclidean distance" begin
        @testset "Distance to origin" begin
            @test Utility.euclidean([GreyNumber(), GreyNumber()]) == GreyNumber(0.0, 0.0)
        end
        @testset "Distance to two grey numbers" begin
            @test Utility.euclidean(
                [GreyNumber(0.0, 0.0),GreyNumber(0.0, 1.0)],
                [GreyNumber(0.0, 0.0),GreyNumber(0.0, 2.0)]
                ) == GreyNumber(1.0, 2.0)
        end
    end

    @testset "Apply columns" begin
        @testset "Apply sum to columns" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_columns(sum, greys)
            @test length(result) == 2
            @test result[1] == GreyNumber(8, 11)
            @test result[2] == GreyNumber(80, 110)
            @test result == Utility.colsums(greys)
        end

        @testset "Apply minimum to columns" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_columns(minimum, greys)
            @test length(result) == 2
            @test result[1] == GreyNumber(1, 2)
            @test result[2] == GreyNumber(10, 20)
            @test result == Utility.colmins(greys)
        end

        @testset "Apply maximum to columns" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_columns(maximum, greys)
            @test length(result) == 2
            @test result[1] == GreyNumber(4, 5)
            @test result[2] == GreyNumber(40, 50)
            @test result == Utility.colmaxs(greys)
        end

        @testset "Apply mean to columns" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_columns(Utility.mean, greys)
            @test length(result) == 2
            @test result[1] == GreyNumber(8/3, 11/3)
            @test result[2] == GreyNumber(80/3, 110/3)
            @test result == Utility.colmeans(greys)
        end

    end

    @testset "Apply rows" begin
        @testset "Apply sum to rows" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_rows(sum, greys)
            @test length(result) == 3
            @test result[1] == GreyNumber(11, 22)
            @test result[2] == GreyNumber(33, 44)
            @test result[3] == GreyNumber(44, 55)
            @test result == Utility.rowsums(greys)
        end

        @testset "Apply minimum to rows" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_rows(minimum, greys)
            @test length(result) == 3
            @test result[1] == GreyNumber(1, 2)
            @test result[2] == GreyNumber(3, 4)
            @test result[3] == GreyNumber(4, 5)
            @test result == Utility.rowmins(greys)
        end

        @testset "Apply maximum to rows" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_rows(maximum, greys)
            @test length(result) == 3
            @test result[1] == GreyNumber(10, 20)
            @test result[2] == GreyNumber(30, 40)
            @test result[3] == GreyNumber(40, 50)
            @test result == Utility.rowmaxs(greys)
        end

        @testset "Apply mean to rows" begin
            greys = [
                GreyNumber(1, 2) GreyNumber(10, 20);
                GreyNumber(3, 4) GreyNumber(30, 40);
                GreyNumber(4, 5) GreyNumber(40, 50);
            ]
            result = Utility.apply_rows(Utility.mean, greys)
            @test length(result) == 3
            @test result[1] == GreyNumber(11/2, 22/2)
            @test result[2] == GreyNumber(33/2, 44/2)
            @test result[3] == GreyNumber(44/2, 55/2)
            @test result == Utility.rowmeans(greys)
        end
    end

    @testset "isvalid" begin
        @test isvalid(GreyNumber())
        @test isvalid(GreyNumber(7, 8))
        # if lower bound is not less or equal
        # than the upper bound, the type
        # constructor automatically replace
        # them to make the number valid.
        @test isvalid(GreyNumber(8, 7))
    end

    @testset "convert" begin
        g = GreyNumber(2, 5)
        arr = convert(Array{Float64, 1}, g)
        @test arr isa Array{Float64, 1}
        @test arr[1] == 2.0
        @test arr[2] == 5.0
    end

    @testset "Reverse minimum & maximum array" begin
        fns = [minimum, maximum, maximum, minimum, maximum]
        revfns = [maximum, minimum, minimum, maximum, minimum]
        @test Utility.reverseminmax(fns) == revfns
        @test Utility.reverseminmax(revfns) == fns
    end

    @testset "Make Array of minimum and maximum" begin
        result1 = Utility.makeminmax([maximum, maximum, maximum, maximum])
        @test typeof(result1) == Array{Function,1}
        @test typeof(result1[1]([1.0, 2.0, 3.0])) == Float64
        @test result1[1]([1.0, 2.0, 3.0]) == 3.0
    end

    @testset "Product weight vector with columns of matrix" begin
        data = [
            GreyNumber(1.0, 2.0) GreyNumber(2.0, 3.0) GreyNumber(3.0, 4.0);
            GreyNumber(1.0, 2.0) GreyNumber(2.0, 3.0) GreyNumber(3.0, 4.0);
            GreyNumber(1.0, 2.0) GreyNumber(2.0, 3.0) GreyNumber(3.0, 4.0)
        ]
        w = [0.5, 10, 100]
        result = Utility.weightize(data,  w)

        @test size(result) == (3, 3)
        @test typeof(result) == Array{GreyNumber{Float64}, 2}

        @test result[:,1] == [
            GreyNumber(0.5, 1.0),
            GreyNumber(0.5, 1.0),
            GreyNumber(0.5, 1.0)
        ]

        @test result[:,2] == [
            GreyNumber(2.0 * 10, 3.0 * 10),
            GreyNumber(2.0 * 10, 3.0 * 10),
            GreyNumber(2.0 * 10, 3.0 * 10)
        ]

        @test result[:,3] == [
            GreyNumber(3.0 * 100, 4.0 * 100),
            GreyNumber(3.0 * 100, 4.0 * 100),
            GreyNumber(3.0 * 100, 4.0 * 100)
        ]

    end

end # End of test Utility


text:
module Tilde

import Base.rand
using Random
using Reexport: @reexport

@reexport using StatsFuns
@reexport using MeasureTheory
using MeasureBase: productmeasure, Returns

import DensityInterface: logdensityof
import DensityInterface: densityof
import DensityInterface: DensityKind
using DensityInterface

using NamedTupleTools
using SampleChains
# using SymbolicCodegen

# using SymbolicUtils: Symbolic
# const MaybeSym{T} = Union{T, Symbolic{T}}

# MeasureTheory.For(f, dims::MaybeSym{<: Integer}...) = ProductMeasure(mappedarray(i -> f(Tuple(i)...), CartesianIndices(dims)))
# MeasureTheory.For(f, dims::MaybeSym{<: AbstractUnitRange}...) = ProductMeasure(mappedarray(i -> f(Tuple(i)...), CartesianIndices(dims)))

import MacroTools: prewalk, postwalk, @q, striplines, replace, @capture
import MacroTools
import MLStyle
# import MonteCarloMeasurements
# using MonteCarloMeasurements: Particles, StaticParticles, AbstractParticles

using Requires
using ArrayInterface: StaticInt
using Static

using IfElse: ifelse
using TransformVariables: as, asℝ, as𝕀, asℝ₊
import TransformVariables as TV

# using SimplePosets: SimplePoset
# import SimplePosets

using GeneralizedGenerated
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)
using MeasureBase: AbstractKleisli

using NestedTuples: TypelevelExpr

using MeasureTheory: ∞
import MeasureTheory: xform

"""
we use this to avoid introduce static type parameters
for generated functions
"""
_unwrap_type(a::Type{<:Type}) = a.parameters[1]

export model, Model, tilde, @model

using MLStyle
include("callify.jl")
import GeneralizedGenerated as GG

@generated function MeasureTheory.For(f::GG.Closure{F,Free}, inds::I) where {F,Free,I<:Tuple}
    freetypes = Free.types
    eltypes = eltype.(I.types)
    T = Core.Compiler.return_type(F, Tuple{freetypes..., eltypes...})
    quote
        $(Expr(:meta, :inline))
        For{$T,GG.Closure{F,Free},I}(f, inds)
    end
end

include("optics.jl")
include("maybe.jl")
include("core/models/abstractmodel.jl")
include("core/models/astmodel/astmodel.jl")
include("core/models/model.jl")
include("core/dependencies.jl")
include("core/utils.jl")
include("core/models/closure.jl")
include("core/models/posterior.jl")
include("primitives/interpret.jl")
include("distributions/iid.jl")

include("primitives/rand.jl")
include("primitives/logdensity.jl")
include("primitives/insupport.jl")

# include("primitives/basemeasure.jl")
include("primitives/testvalue.jl")
include("primitives/testparams.jl")

include("transforms/utils.jl")

function __init__()
    @require SampleChainsDynamicHMC = "6d9fd711-e8b2-4778-9c70-c1dfb499d4c4" begin
        include("samplechains/dynamichmc.jl")
    end
end

end # module


text:
"""
`concat(list1,list2)` concatenates two lists (Julia type `Vector`)
returning a new list. Designed for use with `@distributed`.

## Example
```julia
julia> using Distributed

julia> addprocs(5)

julia> @everywhere include("concat.jl")

julia> combo = @distributed (concat) for j=1:100
       rand(100)
       end;

julia> length(combo)
10000
```
"""
function concat(x::Vector{T}, y::Vector{T}) where T
    nx = length(x)
    ny = length(y)
    result = Vector{T}(undef,nx + ny)
    for i = 1:nx
        @inbounds result[i] = x[i]
    end
    for i = 1:ny
        @inbounds result[i+nx] = y[i]
    end
    return result
end


text:
include("src/unified_real_ctrl_fvecs_bettis_PIs_new.jl") ## load the data analysis utils

## Usage demonstration
data_path = "your_data_path.jld" ## data source path
d_up = 5 ## dimension upper limit for persistent homology computation
cutoff_size = 1 ## row removal cutoff threshold
cutoff_ind_ratio = 0.50 ## filtration index cutoff

real_data_fvec_bettis_PIs(data_path, d_up, cutoff_size, cutoff_ind_ratio)


text:
##
using CSV, DataFrames, DataFramesMeta, StatsPlots, Dates, DelimitedFiles
cd(@__DIR__)
# download and unzip the data file, from
# https://realized.oxford-man.ox.ac.uk/images/oxfordmanrealizedvolatilityindices.zip

## read the data into a DataFrame
data = CSV.read("oxfordmanrealizedvolatilityindices.csv",DataFrame)

## select the desired columns
data = data[:,[:Column1, :Symbol, :close_price, :rv10, :bv]]
# set desired names
rename!(data, :Column1 => :Date)
rename!(data, :rv10 => :rv)
# select desired index
sp500 = @subset(data, :Symbol .== ".SPX")

## drop times from date by taking the leading 10 characters (yyyy-mm-dd)
sp500.Date = [sp500.Date[i][1:10] for i = 1:size(sp500,1)]

## compute percentage returns
sp500.rets = vcat(0.0,[100.0 *(log(sp500.close_price[i]) - log(sp500.close_price[i-1])) for i=2:size(sp500,1)])
# scale volatility measures to be compatible with percentage returns (log price mult. by 100)
sp500.rv = 10000.0 .* sp500.rv
sp500.bv = 10000.0 .* sp500.bv

## select the date range
sp500 = @subset(sp500, :Date .>= "2013-12-17") # these dates give 1000 obs.
sp500 = @subset(sp500, :Date .<= "2017-12-05")

## write out variables to plain text file
data = sp500[:,[:Date, :rets, :rv, :bv]]
#CSV.write("sp500.csv", data) # write as CSV

##
# some plots
plot(sp500.Date, sp500.rets, tickfontsize=5, legend=false)

##
plot(sp500.Date, [sp500.rv, sp500.bv], tickfontsize=5, label=["rv" "bv"])

##
# the next plot shows that volatility is higher when
# returns are in the tails of distribution, and that
# the overall correlation is negative
marginalkde(sp500.rets, sp500.bv)

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
convert(::Type{Sparams}, raw::TouchstoneData) = _Raw_to_S(raw)

"""
Conversion from `TouchstoneData` to `NetworkData{Sparams}`
"""
function _Raw_to_S(touchstone::TouchstoneData)
    nPort = touchstone.nPort
    nPoint = touchstone.nPoint
    Z0 = touchstone.impedance
    freq_unit = touchstone.freq_unit
    data_type = touchstone.data_type
    format_type = touchstone.format_type
    data = touchstone.data

    ports = [Port(Z0) for n in 1:nPort]
    if data_type != "S"
        error("Data type not a scattering parameter")
    else
        # extracting frequency
        if freq_unit == "HZ"
            freq_multiplier = 1.0
        elseif freq_unit == "KHZ"
            freq_multiplier = 1.0e3
        elseif freq_unit == "MHZ"
            freq_multiplier = 1.0e6
        elseif freq_unit == "GHZ"
            freq_multiplier = 1.0e9
        elseif freq_unit == "THZ"
            freq_multiplier = 1.0e12
        end
        freq = data[1, :] * freq_multiplier

        # extracting scattering matrix

        if format_type == "RI"
            params = __touchstone_sparams_ri(nPort, nPoint, data)
        elseif format_type == "MA"
            params = __touchstone_sparams_ma(nPort, nPoint, data)
        elseif format_type == "DB"
            params = __touchstone_sparams_db(nPort, nPoint, data)
        end
    end
    return NetworkData(ports, freq, params)
end

function __touchstone_sparams_ri(nPort, nPoint, data)
    params = [[begin
            idx1 = 1 + 2 * (nPort * (j - 1) + i - 1) + 1
            idx2 = 1 + 2 * (nPort * (j - 1) + i - 1) + 2
            data[idx1, n] + im * data[idx2, n]
        end for i in 1:nPort, j in 1:nPort] for n in 1:nPoint]
    return Sparams(params)
end

function __touchstone_sparams_ma(nPort, nPoint, data)
    params = [[begin
            idx1 = 1 + 2 * (nPort * (j - 1) + i - 1) + 1
            idx2 = 1 + 2 * (nPort * (j - 1) + i - 1) + 2
            data[idx1, n] * exp(im * π / 180 * data[idx2, n])
        end for i in 1:nPort, j in 1:nPort] for n in 1:nPoint]
    return Sparams(params)
end

function __touchstone_sparams_db(nPort, nPoint, data)
    params = [[begin
            idx1 = 1 + 2 * (nPort * (j - 1) + i - 1) + 1
            idx2 = 1 + 2 * (nPort * (j - 1) + i - 1) + 2
            10^(data[idx1, n]/20) * exp(im * π/180 * data[idx2, n])
        end for i in 1:nPort, j in 1:nPort] for n in 1:nPoint]
    return Sparams(params)
end


text:
using TestImages, ImageCore, ImageShow
using Test, Suppressor
using FileIO
using ImageDistances

@testset "GIF" begin
    toucan = testimage("toucan") # 150×162 RGBA image
    moon = testimage("moon") # 256×256 Gray image

    # This might be a non-error in future
    @test_throws ArgumentError ImageShow.gif(toucan, moon)
    warning_msg = Suppressor.@capture_err ImageShow.gif([toucan, moon]; fps=0.5)
    @test occursin("FPS should be larger than 1", warning_msg)

    @testset "show" begin
        # TODO: the `fps` keyword is unable to test

        img = RGB.(testimage("mri-stack")) # ImageMagick saves GIF as RGB/ARGB array
        ref_gif = ImageShow.gif(img)

        @test !showable(MIME("image/gif"), img)
        @test showable(MIME("image/gif"), ref_gif)

        workdir = "tmp"
        mktempdir() do workdir
            fn = joinpath(workdir, "writemime.gif")
            open(fn, "w") do file
                show(file, MIME("image/gif"), ImageShow.gif(img))
            end
            @test load(fn) == ref_gif

            framestack = [img[:, :, i] for i in axes(img, 3)]
            fn = joinpath(workdir, "writemime.gif")
            open(fn, "w") do file
                show(file, MIME("image/gif"), ImageShow.gif(framestack))
            end
            @test load(fn) == ref_gif

            framestack = [toucan, moon]
            fn = joinpath(workdir, "writemime.gif")
            open(fn, "w") do file
                show(file, MIME("image/gif"), ImageShow.gif(framestack))
            end
            @test euclidean(RGB.(load(fn)), RGB.(ImageShow.gif(framestack))) <= 5.5

            sizes = 16:4:64
            values = range(0, stop=1, length=length(sizes))
            gif = ImageShow.gif(values, sizes) do v, x
                fill(v, ntuple(_->x, 2)...)
            end
            fn = joinpath(workdir, "writemime.gif")
            open(fn, "w") do file
                show(file, MIME("image/gif"), gif)
            end
            @test euclidean(RGB.(load(fn)), RGB.(gif)) <= 0.5
        end
    end

    @testset "gif as a readonly array" begin
        A = fill(1.0, 4, 8)
        B = fill(2.0, 2, 2)
        gif = ImageShow.gif([A, B])
        @test axes(gif) == (1:4, 1:8, 1:2)
        @test size(gif) == (4, 8, 2)
        cv = channelview(gif)
        @test all(cv[:, :, 1] .== 1.0) && all(cv[2:3, 4:5, 2] .== 2.0)
    end
end


text:
using PyPlot, PyCall, IJulia

"""
    t = delaunay(p)

Delaunay triangulation `t` of N x 2 node array `p`.
"""
function delaunay(p)
    tri = pyimport("matplotlib.tri")
    t = tri[:Triangulation](p[:,1], p[:,2])
    return Int64.(t[:triangles] .+ 1)
end

"""
    edges, boundary_indices = all_edges(t)

Find all unique edges in the triangulation `t` (ne x 2 array)
Second output is indices to the boundary edges.
"""
function all_edges(t)
    etag = vcat(t[:,[1,2]], t[:,[2,3]], t[:,[3,1]])
    etag = hcat(sort(etag, dims=2), 1:3*size(t,1))
    etag = sortslices(etag, dims=1)
    dup = all(etag[2:end,1:2] - etag[1:end-1,1:2] .== 0, dims=2)[:]
    keep = .![false;dup]
    edges = etag[keep,1:2]
    emap = cumsum(keep)
    invpermute!(emap, etag[:,3])
    emap = reshape(emap,:,3)
    dup = [dup;false]
    dup = dup[keep]
    bndix = findall(.!dup)
    return edges, bndix, emap
end

"""
    e = boundary_nodes(t)

Find all boundary nodes in the triangulation `t`.
"""
function boundary_nodes(t)
    edges, boundary_indices = all_edges(t)
    return unique(edges[boundary_indices,:][:])
end

"""
    tplot(p, t, u=nothing)

If `u` == nothing: Plot triangular mesh with nodes `p` and triangles `t`.
If `u` == solution vector: Plot filled contour color plot of solution `u`.
"""
function tplot(p, t, u=nothing)
    clf()
    axis("equal")
    if u == nothing
        tripcolor(p[:,1], p[:,2], t .- 1, 0*t[:,1],
                  cmap="Set3", edgecolors="k", linewidth=1)
    else
        tricontourf(p[:,1], p[:,2], t .- 1, u, 20)
    end
    draw()
end

"""
    inside = inpolygon(p, pv)

Determine if each point in the N x 2 node array `p` is inside the polygon
described by the NE x 2 node array `pv`.
"""
function inpolygon(p::Array{Float64,2}, pv::Array{Float64,2})
    path = pyimport("matplotlib.path")
    poly = path[:Path](pv)
    inside = [poly[:contains_point](p[ip,:]) for ip = 1:size(p,1)]
end

function remove_outside_tris(p, t, pv)
    pmid = dropdims(sum(p[t,:], dims=2), dims=2) / 3
    is_inside = inpolygon(pmid, pv)
    t = t[is_inside,:]
end

function triarea(p, t)
    d12 = @. p[t[:,2],:] - p[t[:,1],:]
    d13 = @. p[t[:,3],:] - p[t[:,1],:]
    @. abs(d12[:,1] * d13[:,2] - d12[:,2] * d13[:,1]) / 2
end

function remove_tiny_tris(p, t)
    t = t[triarea(p,t) .> 1e-14,:]
end

function circumcenter(p)
    dp1 = @. p[2,:] - p[1,:]
    dp2 = @. p[3,:] - p[1,:]

    mid1 = @. (p[1,:] + p[2,:]) / 2
    mid2 = @. (p[1,:] + p[3,:]) / 2

    s = [ -dp1[2] dp2[2]; dp1[1] -dp2[1]] \ (-mid1 .+ mid2)
    pc = @. mid1' + s[1] * [-dp1[2] dp1[1]]
end

function edge_midpoints(p, t)
    pmid = reshape(p[t,:] + p[t[:,[2,3,1]],:], :, 2) / 2
    pmid = unique(pmid, dims=1)
end

function pmesh(pv, hmax, nref)
    p = zeros(Float64, 0, 2)
    for i = 1:size(pv,1) - 1
        pp = pv[i:i+1,:]
        L = sqrt(sum(diff(pp, dims=1).^2, dims=2))[1]
        if L > hmax
            n = ceil(Int, L / hmax)
            ss = (0:n) / n
            pp = [1 .- ss ss] * pp
        end
        p = [p; pp[1:end-1,:]]
    end

    t = zeros(Int64, 0, 3)
    while true
        t = delaunay(p)
        t = remove_tiny_tris(p, t)
        t = remove_outside_tris(p, t, pv)
        # tplot(p,t), pause(1e-3)

        area = triarea(p, t)
        maxarea, ix = findmax(area)
        if maxarea < hmax^2 / 2
            break
        end
        pc = circumcenter(p[t[ix,:],:])
        p = [p; pc]
    end

    for iref = 1:nref
        p = [p; edge_midpoints(p, t)]
        t = delaunay(p)
        t = remove_tiny_tris(p, t)
        t = remove_outside_tris(p, t, pv)
        # tplot(p, t), pause(1e-3)
    end

    e = boundary_nodes(t)
    p, t, e
end

using SparseArrays
using LinearAlgebra

function waveguide_edges(p,t)
    elist, bedgeindex = all_edges(t)
    bnode = boundary_nodes(t)
    ein = []
    eout = []
    ewall = []
    for i = 1:size(bedgeindex)[1]
        if p[elist[bedgeindex[i],1],1] == 0 && p[elist[bedgeindex[i],2],1] == 0
            ein = [ein; bedgeindex[i]]
        elseif p[elist[bedgeindex[i],1],1] == 5 && p[elist[bedgeindex[i],2],1] == 5
            eout = [eout; bedgeindex[i]]
        else
            ewall = [ewall; bedgeindex[i]]
        end

    end
    ein = elist[ein,:]
    eout = elist[eout,:]
    ewall = elist[ewall,:]
    return ein, eout, ewall

end

using SymPy

function coeffmx(RHS, p, t, k)
    nm = p[t[k,1:3],:]
    LHS = zeros(3,3)
    for n = 1:3
        LHS[n,:] = [1 nm[n,1] nm[n,2]]
    end
    return (LHS\RHS)
end

function kentry(cmx, i, j, ak)
    return ak*(cmx[2,i]*cmx[2,j] + cmx[3,i]*cmx[3,j])
end

function mentry(cmx, i, j, ak, p, t, k)
    f(x,y) = (cmx[1,i] + cmx[2,i]*x + cmx[3,i]*y)*(cmx[1,j] + cmx[2,j]*x + cmx[3,j]*y)
    result = 0.0
    for n = 1:3
        base = (p[t[k,1],:] + p[t[k,2],:] + p[t[k,3],:])/6 + p[t[k,n],:]/2
        result = result + f(base[1],base[2])
    end
    result = (result*ak)/3
    return result
end

function find_bnd(e,t)
  ind=intersect(hcat(getindex.(findall(x -> x==e[1], t),1)),hcat(getindex.(findall(x -> x==e[2],t),1)))
  i=0
  j=0
  for k=1:3
    if t[ind[1],k]==e[1]
      i=k
    elseif t[ind[1],k]==e[2]
      j=k
    end
  end
  return ind[1],i,j
end

function boundsint(p, t, e)
    k, i, j = find_bnd(e, t)
    cmx = coeffmx(Matrix{Float64}(I, 3, 3), p, t, k)
    cross = 0
    f(x,y) = cmx[1,i]*cmx[1,j]*y + x*y*(cmx[1,j]*cmx[2,i] + cmx[1,i]*cmx[2,j]) + (y^2*(cmx[3,i]*cmx[1,j] + cmx[3,j]*cmx[1,i]))/2 + (y^2*x*(cmx[2,i]*cmx[3,j] + cmx[2,j]*cmx[3,i]))/2 + cmx[2,i]*cmx[2,j]*x^2*y + (cmx[3,i]*cmx[3,j]*y^3)/3
    h(x,y) = cmx[1,i]*cmx[1,i]*y + x*y*(cmx[1,i]*cmx[2,i] + cmx[1,i]*cmx[2,i]) + (y^2*(cmx[3,i]*cmx[1,i] + cmx[3,i]*cmx[1,i]))/2 + (y^2*x*(cmx[2,i]*cmx[3,i] + cmx[2,i]*cmx[3,i]))/2 + cmx[2,i]*cmx[2,i]*x^2*y + (cmx[3,i]*cmx[3,i]*y^3)/3
    g(x,y) = cmx[1,j]*cmx[1,j]*y + x*y*(cmx[1,j]*cmx[2,j] + cmx[1,j]*cmx[2,j]) + (y^2*(cmx[3,j]*cmx[1,j] + cmx[3,j]*cmx[1,j]))/2 + (y^2*x*(cmx[2,j]*cmx[3,j] + cmx[2,j]*cmx[3,j]))/2 + cmx[2,j]*cmx[2,j]*x^2*y + (cmx[3,j]*cmx[3,j]*y^3)/3
    if p[t[k,i],2] > p[t[k,j],2]
        cross = f(p[t[k,i],1],p[t[k,i],2]) - f(p[t[k,i],1],p[t[k,j],2])
        id = h(p[t[k,i],1],p[t[k,i],2]) - h(p[t[k,i],1],p[t[k,j],2])
        jd = g(p[t[k,i],1],p[t[k,i],2]) - g(p[t[k,i],1],p[t[k,j],2])
    elseif p[t[k,i],2] < p[t[k,j],2]
        cross = f(p[t[k,i],1],p[t[k,j],2]) - f(p[t[k,i],1],p[t[k,i],2])
        id = h(p[t[k,i],1],p[t[k,j],2]) - h(p[t[k,i],1],p[t[k,i],2])
        jd = g(p[t[k,i],1],p[t[k,j],2]) - g(p[t[k,i],1],p[t[k,i],2])
    end
    return [id cross; cross jd]
end

function loadvec(p, t, e)
    k, i, j = find_bnd(e, t)
    cmx = coeffmx(Matrix{Float64}(I, 3, 3), p, t, k)

    f(x,y) = (cmx[1,i] + cmx[2,i]*x)*y + (cmx[3,i]*y^2)/2
    g(x,y) = (cmx[1,j] + cmx[2,j]*x)*y + (cmx[3,j]*y^2)/2

    if p[t[k,i],2] > p[t[k,j],2]
        id = f(p[t[k,i],1],p[t[k,i],2]) - f(p[t[k,i],1],p[t[k,j],2])
        jd = g(p[t[k,i],1],p[t[k,i],2]) - g(p[t[k,i],1],p[t[k,j],2])
    else
        id = f(p[t[k,i],1],p[t[k,j],2]) - f(p[t[k,i],1],p[t[k,i],2])
        jd = g(p[t[k,i],1],p[t[k,j],2]) - g(p[t[k,i],1],p[t[k,i],2])
    end

    return [id jd]
end

function femhelmholtz(p, t, ein, eout)
    areas = triarea(p,t)
    n = size(p)[1]
    K = spzeros(n,n); M = spzeros(n,n); Bin = spzeros(n,n); Bout = spzeros(n,n); bin = zeros(n);

    for kn = 1:size(t)[1]
        kinsert = zeros(3,3)
        minsert = zeros(3,3)
        cmx = coeffmx(Matrix{Float64}(I, 3, 3), p, t, kn)
        for i = 1:3, j = 1:3
            kinsert[i,j] = kentry(cmx,i,j,areas[kn])
            minsert[i,j] = mentry(cmx, i, j, areas[kn], p, t, kn)
        end
        K[t[kn,:],t[kn,:]] = K[t[kn,:],t[kn,:]]+kinsert
        M[t[kn,:],t[kn,:]] = M[t[kn,:],t[kn,:]]+minsert
    end

    for m = 1:size(ein)[1]
        einsert = boundsint(p,t,ein[m,:])
        Bin[ein[m,:],ein[m,:]] = Bin[ein[m,:],ein[m,:]] + einsert
    end

    for m = 1:size(eout)[1]
        einsert = boundsint(p,t,eout[m,:])
        Bout[eout[m,:],eout[m,:]] = Bout[eout[m,:],eout[m,:]] + einsert
    end

    for m = 1:size(ein)[1]
        binsert = loadvec(p,t,ein[m,:])
        bin[ein[m,:]] = bin[ein[m,:]] + binsert'
    end

    return K , M , Bin , Bout , bin
end

function helmholtz(pv, hmax, nref, k)
    p, t, e = pmesh(pv, hmax, nref)
    ein, eout, ewall = waveguide_edges(p,t)
    K, M, Bin, Bout, bin = femhelmholtz(p, t, ein, eout)
    A = K-k^2*M+k*im*(Bout + Bin)
    B = bin*2*k*im
    u = A\B
    ur = real.(u)
    "tplot(p,t,ur)"
    return u, ur, Bout
end

function helmerrors(pv, hmax, nrefmax, k)
   pmax, tmax, emax = pmesh(pv, hmax, nrefmax)
    e=Base.MathConstants.e
    f(x,y) = e.^(-k*x*im)
   umax = f(pmax[:,1],1)
   errors = zeros(nrefmax)
   for ncur = 1:(nrefmax)
       "pprox, tprox, eprox = pmesh(pv, hmax, ncur)"
       uprox = helmholtz(pv, hmax, ncur, k)
       errorarray = maximum(abs.(umax[1:size(uprox)[1]] - uprox))
       #errors = [errors ; errorarray]
       errors[ncur]=errorarray
   end
   return errors
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{LegendreSphereCoeff{Float32},Int64,Int64,Float32})
    Base.precompile(Tuple{LegendreSphereCoeff{Float32},UnitRange{Int64},Int64,Float32})
    Base.precompile(Tuple{LegendreSphereCoeff{Float32},UnitRange{Int64},UnitRange{Int64},Float32})
    Base.precompile(Tuple{LegendreSphereCoeff{Float64},Int64,Int64,Float64})
    Base.precompile(Tuple{LegendreSphereCoeff{Float64},UnitRange{Int64},Int64,Float64})
    Base.precompile(Tuple{LegendreSphereCoeff{Float64},UnitRange{Int64},UnitRange{Int64},Float64})
    Base.precompile(Tuple{LegendreSphereNorm,Int64,Int64,Float32})
    Base.precompile(Tuple{LegendreSphereNorm,Int64,Int64,Float64})
    Base.precompile(Tuple{LegendreSphereNorm,UnitRange{Int64},Int64,Float32})
    Base.precompile(Tuple{LegendreSphereNorm,UnitRange{Int64},Int64,Float64})
    Base.precompile(Tuple{LegendreSphereNorm,UnitRange{Int64},UnitRange{Int64},Float32})
    Base.precompile(Tuple{LegendreSphereNorm,UnitRange{Int64},UnitRange{Int64},Float64})
    Base.precompile(Tuple{LegendreUnitCoeff{Float32},Int64,Int64,Float32})
    Base.precompile(Tuple{LegendreUnitCoeff{Float32},UnitRange{Int64},Int64,Float32})
    Base.precompile(Tuple{LegendreUnitCoeff{Float32},UnitRange{Int64},UnitRange{Int64},Float32})
    Base.precompile(Tuple{LegendreUnitCoeff{Float64},Int64,Int64,Float64})
    Base.precompile(Tuple{LegendreUnitCoeff{Float64},UnitRange{Int64},Int64,Float64})
    Base.precompile(Tuple{LegendreUnitCoeff{Float64},UnitRange{Int64},UnitRange{Int64},Float64})
    Base.precompile(Tuple{LegendreUnitNorm,Int64,Int64,Float32})
    Base.precompile(Tuple{LegendreUnitNorm,Int64,Int64,Float64})
    Base.precompile(Tuple{LegendreUnitNorm,UnitRange{Int64},Int64,Float32})
    Base.precompile(Tuple{LegendreUnitNorm,UnitRange{Int64},Int64,Float64})
    Base.precompile(Tuple{LegendreUnitNorm,UnitRange{Int64},UnitRange{Int64},Float32})
    Base.precompile(Tuple{LegendreUnitNorm,UnitRange{Int64},UnitRange{Int64},Float64})
    Base.precompile(Tuple{Type{LegendreSphereCoeff{Float32}},Int64})
    Base.precompile(Tuple{Type{LegendreSphereCoeff{Float64}},Int64})
    Base.precompile(Tuple{Type{LegendreUnitCoeff{Float32}},Int64})
    Base.precompile(Tuple{Type{LegendreUnitCoeff{Float64}},Int64})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreSphereCoeff{Float32},UnitRange{Int64},Int64,Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreSphereCoeff{Float64},UnitRange{Int64},Int64,Vector{Float64}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreSphereNorm,UnitRange{Int64},Int64,Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreSphereNorm,UnitRange{Int64},Int64,Vector{Float64}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreSphereNorm,UnitRange{Int64},UnitRange{Int64},Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreSphereNorm,UnitRange{Int64},UnitRange{Int64},Vector{Float64}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitCoeff{Float32},UnitRange{Int64},Int64,Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitCoeff{Float32},UnitRange{Int64},UnitRange{Int64},Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitCoeff{Float64},UnitRange{Int64},Int64,Vector{Float64}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitCoeff{Float64},UnitRange{Int64},UnitRange{Int64},Vector{Float64}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitNorm,UnitRange{Int64},Int64,Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitNorm,UnitRange{Int64},Int64,Vector{Float64}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitNorm,UnitRange{Int64},UnitRange{Int64},Float32})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitNorm,UnitRange{Int64},UnitRange{Int64},Vector{Float32}})
    Base.precompile(Tuple{typeof(Base.Broadcast.broadcasted),LegendreUnitNorm,UnitRange{Int64},UnitRange{Int64},Vector{Float64}})
    Base.precompile(Tuple{typeof(show),Base.TTY,MIME{Symbol("text/plain")},LegendreSphereCoeff{Float32}})
    Base.precompile(Tuple{typeof(show),Base.TTY,MIME{Symbol("text/plain")},LegendreSphereCoeff{Float64}})
    Base.precompile(Tuple{typeof(show),Base.TTY,MIME{Symbol("text/plain")},LegendreUnitCoeff{Float32}})
    Base.precompile(Tuple{typeof(show),Base.TTY,MIME{Symbol("text/plain")},LegendreUnitCoeff{Float64}})
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module BlockMatchingSOMs

using Distances

include("initializer.jl")
include("bmsom.jl")

end


text:
# Default RNG type for when nothing is specified
const _BASE_RNG = Random123.Philox2x

"""
    TracedRNG{R,N,T}

Wrapped random number generator from Random123 to keep track of random streams during model evaluation
"""
mutable struct TracedRNG{R,N,T<:Random123.AbstractR123{R}} <: Random.AbstractRNG
    "Model step counter"
    count::Int
    "Inner RNG"
    rng::T
    "Array of keys"
    keys::Array{R,N}
end

"""
    TracedRNG(r::Random123.AbstractR123=AdvancedPS._BASE_RNG())
Create a `TracedRNG` with `r` as the inner RNG.
"""
function TracedRNG(r::Random123.AbstractR123=_BASE_RNG())
    Random123.set_counter!(r, 0)
    return TracedRNG(1, r, typeof(r.key)[])
end

# Connect to the Random API
Random.rng_native_52(rng::TracedRNG) = Random.rng_native_52(rng.rng)
Base.rand(rng::TracedRNG, ::Type{T}) where {T} = Base.rand(rng.rng, T)

"""
    split(key::Integer, n::Integer=1)

Split `key` into `n` new keys
"""
function split(key::Integer, n::Integer=1)
    T = typeof(key) # Make sure the type of `key` is consistent on W32 and W64 systems.
    return T[hash(key, i) for i in UInt(1):UInt(n)]
end

"""
    load_state!(r::TracedRNG)

Load state from current model iteration. Random streams are now replayed
"""
function load_state!(rng::TracedRNG)
    key = rng.keys[rng.count]
    Random.seed!(rng.rng, key)
    return Random123.set_counter!(rng.rng, 0)
end

"""
    update_rng!(rng::TracedRNG)

Set key and counter of inner rng in `rng` to `key` and the running model step to 0
"""
function Random.seed!(rng::TracedRNG, key)
    Random.seed!(rng.rng, key)
    return Random123.set_counter!(rng.rng, 0)
end

"""
    save_state!(r::TracedRNG)

Add current key of the inner rng in `r` to `keys`.
"""
function save_state!(r::TracedRNG)
    return push!(r.keys, r.rng.key)
end

Base.copy(r::TracedRNG) = TracedRNG(r.count, copy(r.rng), deepcopy(r.keys))

"""
    set_counter!(r::TracedRNG, n::Integer)

Set the counter of the inner rng in `r`, used to keep track of the current model step
"""
Random123.set_counter!(r::TracedRNG, n::Integer) = r.count = n

"""
    inc_counter!(r::TracedRNG, n::Integer=1)

Increase the model step counter by `n`
"""
inc_counter!(r::TracedRNG, n::Integer=1) = r.count += n


text:
using Plots

data = readdlm(IOBuffer("""
    1 |             2000 |       16.234234 |     0.328893463
    2 |             2000 |       17.744895 |     0.300894045
    4 |             2000 |        7.110266 |     0.750933013
    8 |             2000 |        4.138830 |     1.290058600
   16 |             2000 |        3.829481 |     1.394270721
   32 |             2000 |        3.991391 |     1.337712400
   64 |             2000 |        4.183819 |     1.276186478
"""), '|')

plot!([1,64], [data[1,4], data[1,4]*64], c=:grey, label="ideal Julia")
plot!(data[:,1], data[:,4], xlabel="# cores", ylabel="Gflops", c=:blue, label="Julia", ylims=(0, 1.5), xlims=(0, 80))
scatter!(data[:,1], data[:,4], c=:blue, label="")
png("../plots/hpl")


text:
using LinearAlgebra

import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

import ClimaCore.Geometry: ⊗

using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

const FT = Float64

# https://github.com/CliMA/CLIMAParameters.jl/blob/master/src/Planet/planet_parameters.jl#L5
const MSLP = 1e5 # mean sea level pressure
const grav = 9.8 # gravitational constant
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const γ = 1.4 # heat capacity ratio
const C_p = R_d * γ / (γ - 1) # heat capacity at constant pressure
const C_v = R_d / (γ - 1) # heat capacit at constant volume
const R_m = R_d # moist R, assumed to be dry

const f = 5e-5
const ν = 0.01
const L = 2e2
const nelems = 30
const Cd = ν / (L / nelems)
const ug = 1.0
const vg = 0.0
const uvg = Geometry.UVVector(ug, vg)
const d = sqrt(2 * ν / f)
domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(0.0),
    Geometry.ZPoint{FT}(L);
    boundary_names = (:bottom, :top),
)
#mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
mesh = Meshes.IntervalMesh(domain; nelems = nelems)

cspace = Spaces.CenterFiniteDifferenceSpace(mesh)
fspace = Spaces.FaceFiniteDifferenceSpace(cspace)

function ρ_from_pθ(p, θ_liq)
    T = θ_liq * (p / MSLP)^(R_d / C_p)
    return p / (R_d * T)
end

# https://github.com/CliMA/Thermodynamics.jl/blob/main/src/TemperatureProfiles.jl#L115-L155
# https://clima.github.io/Thermodynamics.jl/dev/TemperatureProfiles/#DecayingTemperatureProfile
function adiabatic_temperature_profile(z; T_surf = 300.0, T_min_ref = 230.0)
    ## Initial conditions
    # Scale height for surface temperature

    # Temperature
    Γ = grav / C_p
    T = max(T_surf - Γ * z, T_min_ref)

    # Pressure
    p = MSLP * (T / T_surf)^(grav / (R_d * Γ))
    if T == T_min_ref
        z_top = (T_surf - T_min_ref) / Γ
        H_min = R_d * T_min_ref / grav
        p *= exp(-(z - z_top) / H_min)
    end

    θ = FT(T_surf)

    ρ = ρ_from_pθ(p, θ)
    ρθ = ρ * θ

    uv = Geometry.UVVector(FT(1.0), FT(0.0))
    return (ρ = ρ, uv = uv, ρθ = ρθ)
end

Π(ρθ) = C_p * (R_d * ρθ / MSLP)^(R_m / C_v)
Φ(z) = grav * z

zc = Fields.coordinate_field(cspace)
Yc = adiabatic_temperature_profile.(zc.z)
w = Geometry.WVector.(zeros(FT, fspace))

Y_init = copy(Yc)
w_init = copy(w)

# Y = (Yc, w)

function tendency!(dY, Y, _, t)
    Yc = Y.Yc
    w = Y.w

    dYc = dY.Yc
    dw = dY.w

    ρ = Yc.ρ
    uv = Yc.uv
    ρθ = Yc.ρθ

    dρ = dYc.ρ
    duv = dYc.uv
    dρθ = dYc.ρθ

    # density
    If = Operators.InterpolateC2F()
    ∂f = Operators.GradientC2F()
    ∂c = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.WVector(zero(FT))),
        top = Operators.SetValue(Geometry.WVector(zero(FT))),
    )
    @. dρ = -∂c(w * If(ρ))

    # potential temperature
    If = Operators.InterpolateC2F()
    ∂f = Operators.GradientC2F()
    ∂c = Operators.DivergenceF2C(
        bottom = Operators.SetValue(Geometry.WVector(zero(FT))),
        top = Operators.SetValue(Geometry.WVector(zero(FT))),
    )
    # TODO!: Undesirable casting to vector required
    @. dρθ = -∂c(w * If(ρθ)) + ρ * ∂c(ν * ∂f(ρθ / ρ))

    uv_1 = Operators.getidx(uv, Operators.Interior(), 1)
    u_wind = LinearAlgebra.norm(uv_1)

    A = Operators.AdvectionC2C(
        bottom = Operators.SetValue(Geometry.UVVector(0.0, 0.0)),
        top = Operators.SetValue(Geometry.UVVector(0.0, 0.0)),
    )

    # uv
    bcs_bottom = Operators.SetValue(Geometry.WVector(Cd * u_wind) ⊗ uv_1)
    bcs_top = Operators.SetValue(uvg)
    ∂c = Operators.DivergenceF2C(bottom = bcs_bottom)
    ∂f = Operators.GradientC2F(top = bcs_top)
    duv .= (uv .- Ref(uvg)) .× Ref(Geometry.WVector(f))
    @. duv += ∂c(ν * ∂f(uv)) - A(w, uv)

    # w
    If = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    ∂f = Operators.GradientC2F()
    ∂c = Operators.GradientF2C()
    Af = Operators.AdvectionF2F()
    divf = Operators.DivergenceC2F()
    B = Operators.SetBoundaryOperator(
        bottom = Operators.SetValue(Geometry.WVector(zero(FT))),
        top = Operators.SetValue(Geometry.WVector(zero(FT))),
    )
    @. dw = B(
        Geometry.WVector(-(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z))) +
        divf(ν * ∂c(w)) - Af(w, w),
    )

    return dY
end

Y = Fields.FieldVector(Yc = Yc, w = w)
dY = tendency!(similar(Y), Y, nothing, 0.0)

Δt = 1.0 / 100.0
# Solve the ODE operator
prob = ODEProblem(tendency!, Y, (0.0, 60 * 60))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 600, # save every hour
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "hydrostatic_ekman"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

z_centers = parent(Fields.coordinate_field(cspace))
z_faces = parent(Fields.coordinate_field(fspace))

function ekman_plot(u; title = "", size = (1024, 600))
    u_ref =
        ug .-
        exp.(-z_centers / d) .*
        (ug * cos.(z_centers / d) + vg * sin.(z_centers / d))
    sub_plt1 = Plots.plot(
        u_ref,
        z_centers,
        marker = :circle,
        xlabel = "u",
        label = "Ref",
    )
    # get u component of uv vector
    sub_plt1 = Plots.plot!(
        sub_plt1,
        parent(u.Yc.uv.components.data.:1),
        z_centers,
        label = "Comp",
    )

    v_ref =
        vg .+
        exp.(-z_centers / d) .*
        (ug * sin.(z_centers / d) - vg * cos.(z_centers / d))
    sub_plt2 = Plots.plot(
        v_ref,
        z_centers,
        marker = :circle,
        xlabel = "v",
        label = "Ref",
    )
    # get v component of uv vector
    sub_plt2 = Plots.plot!(
        sub_plt2,
        parent(u.Yc.uv.components.data.:2),
        z_centers,
        label = "Comp",
    )
    return Plots.plot(
        sub_plt1,
        sub_plt2,
        title = title,
        layout = (1, 2),
        size = size,
    )
end

anim = Plots.@animate for (i, u) in enumerate(sol.u)
    ekman_plot(u, title = "Hour $(i)")
end
Plots.mp4(anim, joinpath(path, "hydrostatic_ekman.mp4"), fps = 10)

Plots.png(ekman_plot(sol[end]), joinpath(path, "hydrostatic_ekman_end.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    relpath(
        joinpath(path, "hydrostatic_ekman_end.png"),
        joinpath(@__DIR__, "../.."),
    ),
    "ekman end",
)


text:
using DiffEqNoiseProcess, Test, Random
using StochasticDiffEq, LinearAlgebra
@testset "NoiseWrapper" begin
  _W = WienerProcess(0.0,0.0,0.0)

  dt = 0.1
  calculate_step!(_W,dt,nothing,nothing)

  for i in 1:10
    accept_step!(_W,dt,nothing,nothing)
  end

  W2 = NoiseWrapper(_W)

  dt = 0.1
  calculate_step!(W2,dt,nothing,nothing)

  for i in 1:10
    accept_step!(W2,dt,nothing,nothing)
  end

  _W = WienerProcess(0.0,0.0,0.0)

  dt = 0.1
  calculate_step!(_W,dt,nothing,nothing)

  for i in 1:10
    accept_step!(_W,dt,nothing,nothing)
  end

  old_W = copy(_W.W)

  W2 = NoiseWrapper(_W)
  lspace =range(_W.t[1], stop=_W.t[end], length=1000)
  dt = lspace[2]-lspace[1]
  calculate_step!(W2,dt,nothing,nothing)
  for t in lspace
    accept_step!(W2,dt,nothing,nothing)
  end

  @test W2.W[end] ≈ _W(W2.t[end])[1]
  @test W2.Z[end] ≈ _W(W2.t[end])[2]

  # Inplace

  _W = WienerProcess!(0.0,zeros(4),zeros(4))

  dt = 0.1
  calculate_step!(_W,dt,nothing,nothing)

  for i in 1:10
    accept_step!(_W,dt,nothing,nothing)
  end

  W2 = NoiseWrapper(_W)

  dt = 0.1
  calculate_step!(W2,dt,nothing,nothing)

  for i in 1:10
    accept_step!(W2,dt,nothing,nothing)
  end

  _W = WienerProcess!(0.0,zeros(4),zeros(4))

  dt = 0.1
  calculate_step!(_W,dt,nothing,nothing)

  for i in 1:10
    accept_step!(_W,dt,nothing,nothing)
  end

  W2 = NoiseWrapper(_W)
  lspace =range(_W.t[1], stop=_W.t[end], length=20)
  dt = lspace[2]-lspace[1]
  calculate_step!(W2,dt,nothing,nothing)
  for t in lspace
    accept_step!(W2,dt,nothing,nothing)
  end

  @test W2.W[end]≈ _W(W2.t[end])[1]
  @test W2.Z[end]≈ _W(W2.t[end])[2]

  @test W2.W[end] ≈ _W.W[end-1]
end

@testset "NoiseWrapper restart tests to interpolate" begin
  seed = 100
  u₀ = [0.75,0.5]
  p = [-1.5,0.05,0.2, 0.01]
  trange = (0.0,0.1)
  dtmix = trange[2]/1e3

  function f_mixing!(du,u,p,t)
    du[1] = p[1]*u[1] + p[2]*u[2]
    du[2] = p[2]*u[1] + p[1]*u[2]
    nothing
  end

  function g_mixing!(du,u,p,t)
    du[1] = p[3]*u[1] + p[4]*u[2]
    du[2] = p[3]*u[1] + p[4]*u[2]
    nothing
  end

  function f_mixing(u,p,t)
    dx = p[1]*u[1] + p[2]*u[2]
    dy = p[2]*u[1] + p[1]*u[2]
    [dx,dy]
  end

  function g_mixing(u,p,t)
    dx = p[3]*u[1] + p[4]*u[2]
    dy = p[3]*u[1] + p[4]*u[2]
    [dx,dy]
  end

  Random.seed!(seed)
  prob = SDEProblem(f_mixing!,g_mixing!,u₀,trange,p)

  soltsave = collect(trange[1]:dtmix:trange[2])
  sol = solve(prob, EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave )

  Random.seed!(seed)
  proboop = SDEProblem(f_mixing,g_mixing,u₀,trange,p)
  soloop = solve(proboop,EulerHeun(), dt=dtmix, save_noise=true, saveat=soltsave)

  @test soloop.u ≈ sol.u

  interval = (sol.t[end-1], sol.t[end])

  # oop

  _sol = deepcopy(soloop)
  soloop.W.save_everystep = false
  _sol.W.save_everystep = false

  forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=1000)
  checkWrapper = solve(remake(_sol.prob, tspan=interval, u0=_sol(interval[1]), noise=forwardnoise), _sol.alg, save_noise=false; dt=abs(_sol.W.dt))

  @test checkWrapper.u[end-1] ≈ soloop.u[end] rtol=1e-10
  @test checkWrapper.W.W[end] ≈ soloop.W.W[end] rtol=1e-16

  @show checkWrapper.u[end] - soloop.u[end]

  forwardnoise = DiffEqNoiseProcess.NoiseGrid(_sol.W.t[1000:1001], _sol.W.W[1000:1001])
  checkGrid = solve(remake(_sol.prob, tspan=interval, u0=_sol(interval[1]), noise=forwardnoise), _sol.alg, save_noise=false; dt=abs(_sol.W.dt))

  @test checkGrid.u[end-1] ≈ soloop.u[end] rtol=1e-10
  @test checkGrid.W.W[end] ≈ soloop.W.W[end] rtol=1e-16

  @show checkGrid.u[end] - soloop.u[end]

  # inplace

  _sol = deepcopy(sol)
  sol.W.save_everystep = false
  _sol.W.save_everystep = false

  forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=1000)
  checkWrapper = solve(remake(_sol.prob, tspan=interval, u0=_sol(interval[1]), noise=forwardnoise), _sol.alg, save_noise=false; dt=abs(_sol.W.dt))

  @test checkWrapper.u[end-1] ≈ sol.u[end] rtol=1e-10
  @test checkWrapper.W.W[end] ≈ sol.W.W[end] rtol=1e-16

  @show checkWrapper.u[end] - sol.u[end]

  forwardnoise = DiffEqNoiseProcess.NoiseGrid(_sol.W.t[1000:1001], _sol.W.W[1000:1001])
  checkGrid = solve(remake(_sol.prob, tspan=interval, u0=_sol(interval[1]), noise=forwardnoise), _sol.alg, save_noise=false; dt=abs(_sol.W.dt))

  @test checkGrid.u[end-1] ≈ sol.u[end] rtol=1e-10
  @test checkGrid.W.W[end] ≈ sol.W.W[end] rtol=1e-16

  @show checkGrid.u[end] - sol.u[end]
end

@testset "Interpolation with small eps jumps" begin
  seed = 100
  u₀ = [1.0]
  p = [1.1,0.87]
  trange = (0.0,1.0)
  dtmix = trange[2]/1e4

  function f!(du,u,p,t)
    du[1] = p[1]*u[1]
    nothing
  end

  function g!(du,u,p,t)
    du[1] = p[2]*u[1]
    nothing
  end

  function f(u,p,t)
    dx = p[1]*u[1]
    [dx]
  end

  function g(u,p,t)
    dx = p[2]*u[1]
    [dx]
  end

  Random.seed!(seed)
  prob = SDEProblem(f!,g!,u₀,trange,p)
  sol = solve(prob, EM(), dt=dtmix, save_noise=true)

  Random.seed!(seed)
  proboop = SDEProblem(f,g,u₀,trange,p)
  soloop = solve(proboop,EM(), dt=dtmix, save_noise=true,save_end=false)

  @test soloop.u ≈ sol.u

  @test length(sol.u) != 1e4+1

  indx1 = Int(length(sol.t) - 100)
  interval = (sol.t[indx1], sol.t[end])

  #inplace

  _sol = deepcopy(sol)
  sol.W.save_everystep = false
  _sol.W.save_everystep = false

  forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=indx1)

  checkWrapper = solve(remake(_sol.prob, tspan=interval, u0=_sol(interval[1]), noise=forwardnoise), _sol.alg, save_noise=false; dt=abs(_sol.W.t[end-1]-_sol.W.t[end-2]),tstops = sol.t[indx1:end])

  @test length(checkWrapper.u) != length(sol.u)

  @test checkWrapper.u[1] ≈ sol.u[indx1] rtol=1e-10
  @test checkWrapper.u[end-1] ≈ sol.u[end-1] rtol=1e-10
  @test checkWrapper.u[end] ≈ sol.u[end] rtol=1e-10
  @test checkWrapper.W.W[end] ≈ sol.W.W[end] rtol=1e-10
  @test checkWrapper(sol.t[indx1:end]).u ≈ sol.u[indx1:end]  rtol=1e-10

  _sol = deepcopy(soloop)
  sol.W.save_everystep = false
  _sol.W.save_everystep = false

  forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=indx1)

  checkWrapper = solve(remake(_sol.prob, tspan=interval, u0=_sol(interval[1]), noise=forwardnoise), _sol.alg, save_noise=false; dt=abs(_sol.W.t[end-1]-_sol.W.t[end-2]),tstops = sol.t[indx1:end])

  @test length(checkWrapper.u) != length(soloop.u)

  @test checkWrapper.u[1] ≈ sol.u[indx1] rtol=1e-10
  @test checkWrapper.u[end-1] ≈ soloop.u[end-1] rtol=1e-10
  @test checkWrapper.u[end] ≈ soloop.u[end] rtol=1e-10
  @test checkWrapper.W.W[end] ≈ soloop.W.W[end] rtol=1e-10
  @test checkWrapper(soloop.t[indx1:end]).u ≈ soloop.u[indx1:end]  rtol=1e-10
end


text:
module MultiFEFunctionsTests

using Test

using Gridap

order = 1
diritag = "boundary"
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,3))
fespace = H1ConformingFESpace(Float64,model,order,diritag)

ufun1(x) = x[1] + x[2]
U1 = TrialFESpace(fespace,ufun1)

ufun2(x) = x[1] + x[2]
U2 = TrialFESpace(fespace,ufun2)

U = MultiFESpace([U1,U2])

n = num_free_dofs(U)

x = rand(n)

uh = MultiFEFunction(x,U)

@test string(uh) == "MultiFEFunction object with 2 fields"

@test length(uh) == 2
@test free_dofs(uh[1]) == x[1:2]
@test free_dofs(uh[2]) == x[3:4]

uh1, state = iterate(uh)
@test free_dofs(uh1) == x[1:2]

uh2, state = iterate(uh,state)
@test free_dofs(uh2) == x[3:4]

zh = zero(U)
@test isa(zh,MultiFEFunction)
@test free_dofs(zh) == zeros(num_free_dofs(U))

zh = zero([U1,U2])
@test isa(zh,MultiFEFunction)
@test free_dofs(zh) == zeros(num_free_dofs(U))

zh = FEFunction(U,free_dofs(zh))
@test isa(zh,MultiFEFunction)
@test free_dofs(zh) == zeros(num_free_dofs(U))

end # module MultiFEFunctionsTests