text:
using DelimitedFiles, PyPlot

res, header = readdlm("res/exper4.csv", ',', Int64, '\n'; header=true)

for m in [25, 50, 100]
    fig = figure(figsize=(3.14961, 3.14961), dpi=1000)

    power = res[findall((res[:, 1] .== m) .& (res[:, 2] .== 1)), 4] ./ 1000

    scatter(0.0:0.05:0.5, power, s=5, marker="o")

    power = res[findall((res[:, 1] .== m) .& (res[:, 2] .== 3)), 4] ./ 1000

    scatter(0.0:0.05:0.5, power, s=5, marker="v")

    power = res[findall((res[:, 1] .== m) .& (res[:, 2] .== 5)), 4] ./ 1000

    scatter(0.0:0.05:0.5, power, s=5, marker="s")

    minorticks_on()
    grid()
    grid(which="minor", ls="dotted", lw=".5")

    title("without Studentization", size="xx-small")
    xlabel("δ", size="xx-small")
    ylabel("proportion of rejections", size="xx-small")

    ylim(0, 1.05)

    tick_params("both", labelsize="xx-small")

    tight_layout()

    savefig("res/exper4_$(m).pdf")

    close(fig)

    fig = figure(figsize=(3.14961, 3.14961), dpi=1000)

    power = res[findall((res[:, 1] .== m) .& (res[:, 2] .== 1)), 5] ./ 1000

    scatter(0.0:0.05:0.5, power, s=5, marker="o")

    power = res[findall((res[:, 1] .== m) .& (res[:, 2] .== 3)), 5] ./ 1000

    scatter(0.0:0.05:0.5, power, s=5, marker="v")

    power = res[findall((res[:, 1] .== m) .& (res[:, 2] .== 5)), 5] ./ 1000

    scatter(0.0:0.05:0.5, power, s=5, marker="s")

    minorticks_on()
    grid()
    grid(which="minor", ls="dotted", lw=".5")

    title("Studentized", size="xx-small")
    xlabel("δ", size="xx-small")
    ylabel("proportion of rejections", size="xx-small")

    ylim(0, 1.05)

    tick_params("both", labelsize="xx-small")

    tight_layout()

    savefig("res/exper4_$(m)_studentized.pdf")

    close(fig)
end


text:
module Rackauckaspackage

# Write your package code here.

end

ext: jl


text:
import Printf:@sprintf
const REGRESS_MARK = ":x:"
const IMPROVE_MARK = ":white_check_mark:"

function create_benchmark(detector)
    SUITE["$detector fit"] = BenchmarkGroup()
    SUITE["$detector score"] = BenchmarkGroup()
    for dimension ∈ (1,10,100)
        for points ∈ (100,1000)
            data = rand(MersenneTwister(1), dimension, points)
            model, scores = fit(detector, data)
            SUITE["$detector fit"]["dim=$dimension, points=$points"] =
                @benchmarkable fit($detector, $data)
            SUITE["$detector score"]["dim=$dimension, points=$points"] =
                @benchmarkable score($detector, $model, $data)
        end
    end
end

function run_benchmarks(name)
    paramspath = joinpath(dirname(@__FILE__), "benchmark.jld")
    if !isfile(paramspath)
        println("Tuning benchmarks...")
        tune!(SUITE)
        JLD.save(paramspath, "SUITE", params(SUITE))
    end
    loadparams!(SUITE, JLD.load(paramspath, "SUITE"), :evals, :samples)
    println("Running benchmarks...")
    results = run(SUITE, verbose = true, seconds = 2)
    JLD.save(joinpath(dirname(@__FILE__), name * ".jld"), "results", results)
end

function generate_report_single(name)
    result = load(joinpath(dirname(@__FILE__), name * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_single.md"), "w") do f
        printreport(f, minimum(result); iscomparisonjob = false)
    end
end

function generate_report_comparison(v1, v2)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    v2_res = load(joinpath(dirname(@__FILE__), v2 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_compare.md"), "w") do f
        printreport(f, judge(minimum(v1_res), minimum(v2_res)); iscomparisonjob = true)
    end
end

function printreport(io::IO, results; iscomparisonjob::Bool=false)

    if iscomparisonjob
        print(io, """
                  A ratio greater than `1.0` denotes a possible regression (marked with $(REGRESS_MARK)), while a ratio less
                  than `1.0` denotes a possible improvement (marked with $(IMPROVE_MARK)). Only significant results - results
                  that indicate possible regressions or improvements - are shown below (thus, an empty table means that all
                  benchmark results remained invariant between builds).
                  | ID | time ratio | memory ratio |
                  |----|------------|--------------|
                  """)
    else
        print(io, """
                  | ID | time | GC time | memory | allocations |
                  |----|------|---------|--------|-------------|
                  """)
    end

    entries = BenchmarkTools.leaves(results)

    try
        entries = entries[sortperm(map(x -> string(first(x)), entries))]
    catch
    end

    for (ids, t) in entries
        if !(iscomparisonjob) || BenchmarkTools.isregression(t) || BenchmarkTools.isimprovement(t)
            println(io, resultrow(ids, t))
        end
    end

    return nothing
end

idrepr(id) = (str = repr(id); str[first(something(findfirst('[', str))):end])
intpercent(p) = string(ceil(Int, p * 100), "%")
resultrow(ids, t::BenchmarkTools.Trial) = resultrow(ids, minimum(t))

function resultrow(ids, t::BenchmarkTools.TrialEstimate)
    t_tol = intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = intpercent(BenchmarkTools.params(t).memory_tolerance)
    timestr = string(BenchmarkTools.prettytime(BenchmarkTools.time(t)), " (", t_tol, ")")
    memstr = string(BenchmarkTools.prettymemory(BenchmarkTools.memory(t)), " (", m_tol, ")")
    gcstr = BenchmarkTools.prettytime(BenchmarkTools.gctime(t))
    allocstr = string(BenchmarkTools.allocs(t))
    return "| `$(idrepr(ids))` | $(timestr) | $(gcstr) | $(memstr) | $(allocstr) |"
end

function resultrow(ids, t::BenchmarkTools.TrialJudgement)
    t_tol = intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = intpercent(BenchmarkTools.params(t).memory_tolerance)
    t_ratio = @sprintf("%.2f", BenchmarkTools.time(BenchmarkTools.ratio(t)))
    m_ratio =  @sprintf("%.2f", BenchmarkTools.memory(BenchmarkTools.ratio(t)))
    t_mark = resultmark(BenchmarkTools.time(t))
    m_mark = resultmark(BenchmarkTools.memory(t))
    timestr = "$(t_ratio) ($(t_tol)) $(t_mark)"
    memstr = "$(m_ratio) ($(m_tol)) $(m_mark)"
    return "| `$(idrepr(ids))` | $(timestr) | $(memstr) |"
end

resultmark(sym::Symbol) = sym == :regression ? REGRESS_MARK : (sym == :improvement ? IMPROVE_MARK : "")


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions related with True Equator Mean Equinox (TEME) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Remarks
#
# As mentioned in [1, p. 233], there is not an official definition of the TEME
# frame. Hence, in this package, it is considered the definition presented in
# [1, p. 233] in which the complete form of the Equation of Equinoxes is used.
# This seems to be the case when comparing the values shown in Table 3-6 [1, p.
# 232].
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export rTEMEtoTOD,  rTODtoTEME
export rTEMEtoMOD,  rMODtoTEME
export rTEMEtoGCRF, rGCRFtoTEME
export rTEMEtoPEF,  rPEFtoTEME

#                                 TEME <=> TOD
# ==============================================================================

"""
    rTEMEtoTOD([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the True of Date (TOD) frame at the Julian Day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p. 233]. Notice
that one can provide corrections for the nutation in obliquity (`δΔϵ_1980`)
\\[rad] and in longitude (`δΔψ_1980`) \\[rad] that are usually obtained from
IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the TOD frame. The rotation
representation is selected by the optional parameter `T`.

"""
rTEMEtoTOD(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rTEMEtoTOD(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

function rTEMEtoTOD(T::Type, JD_TT::Number, δΔϵ_1980::Number = 0,
                    δΔψ_1980::Number = 0)

    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δϵ_1980 += δΔϵ_1980
    Δψ_1980 += δΔψ_1980

    # Compute the obliquity.
    ϵ_1980 = mϵ_1980 + Δϵ_1980

    # Evaluate the Delaunay parameters associated with the Moon in the interval
    # [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    T_TT = (JD_TT - JD_J2000)/36525
    r    = 360
    Ω_m  = 125.04452222 - (5r + 134.1362608)*T_TT +
                          0.0020708*T_TT^2 +
                          2.2e-6*T_TT^3
    Ω_m  = mod(Ω_m, 360)*pi/180

    # Compute the equation of Equinoxes.
    #
    # According to [2], the constant unit before `sin(2Ω_m)` is also in [rad].
    Eq_equinox1982 = Δψ_1980*cos(mϵ_1980) +
        ( 0.002640*sin(1Ω_m) + 0.000063*sin(2Ω_m) )*pi/648000

    # Compute the rotation.
    return angle_to_rot(T, -Eq_equinox1982, 0, 0, :ZYX)
end

"""
    rTODtoTEME([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True of Date (TOD) frame with the True
Equator Mean Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p. 233].
Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) \\[rad] and in longitude (`δΔψ_1980`) \\[rad] that are usually
obtained from IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TOD frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

"""
rTODtoTEME(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rTODtoTEME(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

rTODtoTEME(T::Type, JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    inv_rotation(rTEMEtoTOD(T, JD_TT, δΔϵ_1980, δΔψ_1980))

#                                 TEME <=> MOD
# ==============================================================================

"""
    rTEMEtoMOD([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the Mean of Date (MOD) frame at the Julian Day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p. 233]. Notice
that one can provide corrections for the nutation in obliquity (`δΔϵ_1980`)
\\[rad] and in longitude (`δΔψ_1980`) \\[rad] that are usually obtained from
IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

"""
rTEMEtoMOD(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rTEMEtoMOD(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

function rTEMEtoMOD(T::Type, JD_TT::Number, δΔϵ_1980::Number = 0,
                    δΔψ_1980::Number = 0)

    # Notice that, in this case, we will not use `rTEMEtoTOD`, and `rTODtoMOD`
    # because this would call the function `nutation` twice, leading to a huge
    # performance drop. Hence, the code of those two functions is almost
    # entirely rewritten here.

    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δϵ_1980 += δΔϵ_1980
    Δψ_1980 += δΔψ_1980

    # Compute the obliquity.
    ϵ_1980 = mϵ_1980 + Δϵ_1980

    # Evaluate the Delaunay parameters associated with the Moon in the interval
    # [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    T_TT = (JD_TT - JD_J2000)/36525
    r    = 360
    Ω_m  = 125.04452222 - (5r + 134.1362608)*T_TT +
                          0.0020708*T_TT^2 +
                          2.2e-6*T_TT^3
    Ω_m  = mod(Ω_m, 360)*pi/180

    # Compute the equation of Equinoxes.
    #
    # According to [2], the constant unit before `sin(2Ω_m)` is also in [rad].
    Eq_equinox1982 = Δψ_1980*cos(mϵ_1980) +
        ( 0.002640*sin(1Ω_m) + 0.000063*sin(2Ω_m) )*pi/648000

    # Compute the rotation TEME => TOD.
    r_TOD_TEME = angle_to_rot(T, -Eq_equinox1982, 0, 0, :ZYX)

    # Compute the rotation TOD => MOD.
    r_MOD_TOD = angle_to_rot(T, ϵ_1980, Δψ_1980, -mϵ_1980, :XZX)

    # Return the full rotation.
    return compose_rotation(r_TOD_TEME, r_MOD_TOD)
end

"""
    rMODtoTEME([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the Mean of Date (MOD) frame with the True
Equator Mean Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p. 233].
Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) \\[rad] and in longitude (`δΔψ_1980`) \\[rad] that are usually
obtained from IERS EOP Data (see `get_iers_eop`).  .

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

"""
rMODtoTEME(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rMODtoTEME(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

rMODtoTEME(T::Type, JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    inv_rotation(rTEMEtoMOD(T, JD_TT, δΔϵ_1980, δΔψ_1980))

#                                TEME <=> GCRF
# ==============================================================================

"""
    rTEMEtoGCRF([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the Geocentric Celestial Reference Frame (GCRF) at the Julian Day `JD_TT`
[Terrestrial Time]. This algorithm uses the IAU-76/FK5 theory and TEME
definition in [1, p. 233]. Notice that one can provide corrections for the
nutation in obliquity (`δΔϵ_1980`) \\[rad] and in longitude (`δΔψ_1980`) \\[rad]
that are usually obtained from IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the GCRF frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The EOP data related to the nutation of the obliquity (`δΔϵ_1980`) and the
nutation of the longitude (`δΔψ_1980`) can be omitted. In this case, the GCRF
frame is what is usually called J2000 reference frame.

"""
rTEMEtoGCRF(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rTEMEtoGCRF(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

function rTEMEtoGCRF(T::Type, JD_TT::Number, δΔϵ_1980::Number = 0,
                     δΔψ_1980::Number = 0)

    # Compute the rotation TEME => MOD.
    r_MOD_TEME  = rTEMEtoMOD(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compute the rotation MOD => GCRF.
    r_GCRF_MOD = rMODtoGCRF_fk5(T, JD_TT)

    # Return the full rotation.
    return compose_rotation(r_MOD_TEME, r_GCRF_MOD)
end

"""
    rGCRFtoTEME([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the GCRF frame with the True Equator Mean
Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p. 233]. Notice
that one can provide corrections for the nutation in obliquity (`δΔϵ_1980`)
\\[rad] and in longitude (`δΔψ_1980`) \\[rad] that are usually obtained from
IERS EOP Data (see `get_iers_eop`).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the GCRF frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

# Remarks

The EOP data related to the nutation of the obliquity (`δΔϵ_1980`) and the
nutation of the longitude (`δΔψ_1980`) can be omitted. In this case, the GCRF
frame is what is usually called J2000 reference frame.

"""
rGCRFtoTEME(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
    rGCRFtoTEME(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)

rGCRFtoTEME(T::Type, JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0) =
   inv_rotation(rTEMEtoGCRF(T, JD_TT, δΔϵ_1980, δΔψ_1980))

#                                 TEME <=> PEF
# ==============================================================================

"""
    rTEMEtoPEF([T,] JD_TT::Number)

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the Pseudo-Earth Fixed (PEF) frame at the Julian Day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p. 233].

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the PEF frame. The rotation
representation is selected by the optional parameter `T`.

"""
rTEMEtoPEF(JD_UT1::Number) = rTEMEtoPEF(DCM, JD_UT1)

function rTEMEtoPEF(T::Type, JD_UT1::Number)
    # Compute the Greenwich Mean Sidereal Time.
    θ_gmst = JDtoGMST(JD_UT1)

    # Compute the rotation.
    return angle_to_rot(T, θ_gmst, 0, 0, :ZYX)
end

"""
    rPEFtoTEME([T,] JD_TT::Number)

Compute the rotation that aligns the Pseudo-Earth Fixed (PEF) frame with the
True Equator Mean Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial
Time]. This algorithm uses the IAU-76/FK5 theory and TEME definition in [1, p.
233].

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the PEF frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

"""
rPEFtoTEME(JD_UT1::Number)          = rPEFtoTEME(DCM, JD_UT1)
rPEFtoTEME(T::Type, JD_UT1::Number) = inv_rotation(rTEMEtoPEF(T, JD_UT1))


text:
#
# InterProcessCommunication.jl --
#
# Inter-Process Communication for Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of InterProcessCommunication.jl released under the MIT
# "expat" license.
#
# Copyright (C) 2016-2019, Éric Thiébaut
# (https://github.com/emmt/InterProcessCommunication.jl).
#

module InterProcessCommunication

export
    CLOCK_MONOTONIC,
    CLOCK_REALTIME,
    DynamicMemory,
    FileDescriptor,
    FileStat,
    IPC,
    Semaphore,
    SharedMemory,
    ShmArray,
    ShmId,
    ShmInfo,
    ShmMatrix,
    ShmVector,
    SigAction,
    SigInfo,
    SigSet,
    TimeSpec,
    TimeVal,
    TimeoutError,
    WrappedArray,
    clock_getres,
    clock_gettime,
    clock_settime,
    gettimeofday,
    nanosleep,
    now,
    post,
    shmat,
    shmcfg,
    shmctl,
    shmdt,
    shmget,
    shmid,
    shminfo!,
    shminfo,
    shmrm,
    sigaction!,
    sigaction,
    signal,
    sigpending!,
    sigpending,
    sigprocmask!,
    sigprocmask,
    sigqueue,
    sigsuspend,
    sigwait!,
    sigwait,
    timedlock,
    trywait

using Printf
using Dates
import Dates: now

using Base: elsize, tail, OneTo, throw_boundserror, @propagate_inbounds

import Base: convert, unsafe_convert,
    islocked, lock, unlock, trylock, timedwait, wait, broadcast

# The following is an exported shortcut to the package.
const IPC = InterProcessCommunication
const PARANOID = true

isfile(joinpath(@__DIR__, "..", "deps", "deps.jl")) ||
    error("InterProcessCommunication not properly installed.  Please run Pkg.build(\"InterProcessCommunication\")")
include(joinpath("..", "deps", "deps.jl"))
include("types.jl")
include("wrappedarrays.jl")
include("unix.jl")
include("utils.jl")
include("shm.jl")
include("semaphores.jl")
include("signals.jl")
include("locks.jl")

@deprecate IPC_NEW IPC.PRIVATE

end # module InterProcessCommunication


text:
abstract AbstractModel

### Factory function
model(s::Symbol,args...;keyargs...) = _model(Val{s},args...;keyargs...)

### Initialize is the same for all models
initialize!(i::InitializeFrom,m::AbstractModel,s::AbstractSample) = (initvalues!(i,parameters(m),s.values) ; s)

### For each model, an evaluation function needs to be defined
evaluate!(m::AbstractModel,v::AbstractVector) = throw(MethodError(evaluate!, (m,v)))

### Common utility functions
numparas(m::AbstractModel) = length(m.parameters)
parameters(m::AbstractModel) = m.parameters

### Utility functions that need to be implemented for each model
@inline dataindex(m::AbstractModel) = throw(MethodError(dataindex, (m,)))
@inline measurements(m::AbstractModel) = throw(MethodError(measurements, (m,)))
@inline noisemodel(m::AbstractModel) = throw(MethodError(noisemodel, (m,)))


text:
############# Distribution types vectorizations   ################
# single parameter distributions
import Distributions: Bernoulli, TDist, Exponential, Poisson

for d in [:Bernoulli, :TDist, :Exponential, :Poisson]
	@eval begin
		function ($d)(ps::Array)
			ds = Array($d, size(ps))
			for i in 1:length(ds)
				ds[i] = ($d)(ps[i])
			end
			ds
		end
	end
end

# two parameter distributions
import Distributions: Normal, Uniform, Weibull, Gamma, Cauchy, LogNormal, Binomial, Beta, Laplace

for d in [:Normal, :Uniform, :Weibull, :Gamma, :Cauchy, :LogNormal, :Binomial, :Beta, :Laplace]
	@eval begin
		function ($d)(p1::Array, p2::Array)
			ds = Array($d, size(p1))
			for i in 1:length(ds)
				ds[i] = ($d)(p1[i], p2[i])
			end
			ds
		end

		function ($d)(p1::Array, p2::Real)
			ds = Array($d, size(p1))
			for i in 1:length(ds)
				ds[i] = ($d)(p1[i], p2)
			end
			ds
		end

		function ($d)(p1::Real, p2::Array)
			ds = Array($d, size(p2))
			for i in 1:length(ds)
				ds[i] = ($d)(p1, p2[i])
			end
			ds
		end
	end
end

############# logpdf vectorization on the distribution argument   ################
import Distributions: logpdf, logcdf, logccdf

function logpdf{T<:Distribution}(ds::Array{T}, x::AbstractArray)
	res = Array(Float64, size(ds))
	size(ds) == size(x) || error("x and distributions sizes do not match")
	for i in 1:length(x)
		res[i] = logpdf(ds[i], x[i])
	end
	res
end

function logcdf{T<:Distribution}(ds::Array{T}, x::AbstractArray)
	res = Array(Float64, size(ds))
	size(ds) == size(x) || error("x and distributions sizes do not match")
	for i in 1:length(x)
		res[i] = logcdf(ds[i], x[i])
	end
	res
end

function logccdf{T<:Distribution}(ds::Array{T}, x::AbstractArray)
	res = Array(Float64, size(ds))
	size(ds) == size(x) || error("x and distributions sizes do not match")
	for i in 1:length(x)
		res[i] = logccdf(ds[i], x[i])
	end
	res
end


text:
# Astronomy example from Hairer et al 1992 p.245. As Hairer writes, "a
# celestial mechanics problem( which we call 'the Pleiades'): seven stars in
# the plan with coordinates x_i,y_i and masses m_i = i". Time span from
#t =0 to t= 3. This is a non-stiff problem.
#
# Reference solution and a more in depth description can be found in
# https://archimede.dm.uniba.it/~testset/report/plei.pdf

export plei

# create an instance of TestCaseExplicit by filling in the ...:
plei = let
    tcname = :plei ::Symbol # name should be same as variable name (except for upper/lower case)

    T = Float64 # the datatype used, probably Float64
    Tarr = [Matrix, SparseMatrixCSC][2] # the container datatype used
                                        # for the Jacobian and mass
                                        # matrix, probably Matrix
                                        # (default) or
                                        # SparseMatrixCSC.
    const N = 7
    const dof = 4N::Int # degrees of freedom
    dae = 0 ::Int  # index of DAE, ==0 for ODE
    # stiffness of system, one of the three constants
    stiffness = [nonstiff, mildlystiff, stiff][1]

    ## the problem function
    function fn!(t::T,y::Vector{T},dydt::Vector{T})
        # The ode function dydt = f(t,y) modifying dydt in-place.
        # Note that all indices of dydt need to be written!

        # Problem is of the form z = (x,y) and z'' = f(z), (z(0),z'(0)) = (z_0,z'_0)

        #input (z,w), where w = z'
        #output dydt = (z,w)' = (w, f(z)) , where f(z) = z''

        for i=1:2N
            dydt[i] = y[2N+i]
        end

        dydt[2N+1:end] = 0

        #parameters
        for i = 1:N
          for j =1:N
            if i != j
              r = ((y[i] -y[j])^2 + (y[i+N] -y[j+N])^2)^(3/2)
              #Note, for what follows, dydt[2N+i] = x″[i] and dydt[3N+i] = y″[i]
              dydt[2N+i] += j*(y[j] - y[i])/r  #j factor is mass of particles,
              dydt[3N+i] += j*(y[j+N] - y[i+N])/r
            end
          end
        end

        return nothing
    end
    # initializes storage for y:
    fn!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_)

    # NOTE, using keyword arguments more intuitively like:
    # fn!( ; T::Type=T, dof=dof) = zeros(T,dof)
    # does not work: https://github.com/JuliaLang/julia/issues/9948

    jac! = nothing
    mass! = nothing

    ic =  T[3, 3, -1, -3, 2, -2, 2,  #x(0)
            3, -3, 2, 0, 0, -4, 4,    #y(0)
            0, 0, 0, 0, 0, 1.75, -1.5, #x'(0)
            0, 0, 0, -1.25, 1, 0, 0]  #y'(0)
              # vector of initial conditions
    tspan = T[0, 3] # integration interval

    refsol = T[0.3706139143970502,3.237284092057233,-3.222559032418324,0.6597091455775310, 0.3425581707156584, 1.562172101400631,-0.7003092922212495,
                -3.943437585517392,-3.271380973972550, 5.225081843456543, -2.590612434977470, 1.198213693392275, -0.2429682344935824,1.091449240428980,
                3.417003806314313, 1.354584501625501, -2.590065597810775, 2.025053734714242, -1.155815100160448, -0.8072988170223021, 0.5952396354208710,
                -3.741244961234010, 0.3773459685750630, 0.9386858869551073, 0.3667922227200571, -0.3474046353808490, 2.344915448180937, -1.947020434263292]
               # reference solution at tspan[2]

    refsolinds = trues(dof)   # if refsol does not contain all indices
                              # then specify which indices of the
                              # solution are to be compared to refsol

    scd_absinds = Int[]       # set where refsol is very small or
                              # zero. To avoid it dominating the
                              # relative error

    tc = TestCaseExplicit{tcname, T, Tarr}(
                             stiffness,
                             dae,
                             dof,
                             fn!,
                             jac!,
                             mass!,
                             ic,
                             tspan,
                             refsol,
                             refsolinds,
                             scd_absinds)

    # put into the right buckets
    tc_all[tcname] = tc
    tc_nonstiff[tcname] = tc
    # Note that the return-result of the last statement of the let
    # block needs to be tc!
    tc
end


text:
module Evolutionary

using Random, Base.Threads
using Distributed
using DistributedArrays, DistributedArrays.SPMD
using Printf
using Dates

export
    # Optimization methods
    es, cmaes, ga,
    # Constants
    GAVector, Individual, AbstractGene

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

"""
Abstract Type that represents all types of genes supported.
"""
abstract type AbstractGene end

const Strategy = Dict{Symbol,Any}

const Individual = Vector{AbstractGene}

const GAVector = Union{T, BitVector} where T <: Vector

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

# Wrapping function for strategy
function strategy(; kwargs...)
    result = Dict{Symbol,Any}()
    for (k, v) in kwargs
        result[k] = v
    end
    return result
end

# Inverse function for reversing optimization direction
function inverseFunc(f::Function)
    function fitnessFunc(x::T) where {T <: Vector}
        return 1.0/(f(x)+eps())
    end
    return fitnessFunc
end

# Obtain individual
function getIndividual(init::Individual, N::Int)
    if isa(init, Vector)
        @assert length(init) == N "Dimensionality of initial population must be $(N)"
        individual = init
    elseif isa(init, Matrix)
        @assert size(init, 1) == N "Dimensionality of initial population must be $(N)"
        populationSize = size(init, 2)
        individual = init[:, 1]
    elseif isa(init, Function) # Creation function
        individual = init(N)
    else
        individual = ones(N)
    end
    return individual
end

# Collecting interim values
function keep(interim, v, vv, col)
    if interim
        if !haskey(col, v)
            col[v] = typeof(vv)[]
        end
        push!(col[v], vv)
    end
end

# General Structures
include("structs.jl")

# ES & GA recombination functions
include("recombinations.jl")

# ES & GA mutation functions
include("mutations.jl")

# GA selection functions
include("selections.jl")

# Evolution Strategy
include("es.jl")
include("cmaes.jl")

# Backup functions
include("backup.jl")

# Genetic Algorithms
include("ga.jl")

end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: finspace_data
using AWS.Compat
using AWS.UUIDs

"""
    create_changeset(change_type, dataset_id, source_params, source_type)
    create_changeset(change_type, dataset_id, source_params, source_type, params::Dict{String,<:Any})

Creates a new changeset in a FinSpace dataset.

# Arguments
- `change_type`: Option to indicate how a changeset will be applied to a dataset.
  REPLACE - Changeset will be considered as a replacement to all prior loaded changesets.
  APPEND - Changeset will be considered as an addition to the end of all prior loaded
  changesets.
- `dataset_id`: The unique identifier for the FinSpace dataset in which the changeset will
  be created.
- `source_params`: Source path from which the files to create the changeset will be sourced.
- `source_type`: Type of the data source from which the files to create the changeset will
  be sourced.    S3 - Amazon S3.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"formatParams"`: Options that define the structure of the source file(s).
- `"formatType"`: Format type of the input files being loaded into the changeset.
- `"tags"`: Metadata tags to apply to this changeset.
"""
create_changeset(changeType, datasetId, sourceParams, sourceType; aws_config::AbstractAWSConfig=global_aws_config()) = finspace_data("POST", "/datasets/$(datasetId)/changesets", Dict{String, Any}("changeType"=>changeType, "sourceParams"=>sourceParams, "sourceType"=>sourceType); aws_config=aws_config)
create_changeset(changeType, datasetId, sourceParams, sourceType, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace_data("POST", "/datasets/$(datasetId)/changesets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("changeType"=>changeType, "sourceParams"=>sourceParams, "sourceType"=>sourceType), params)); aws_config=aws_config)

"""
    get_programmatic_access_credentials(environment_id)
    get_programmatic_access_credentials(environment_id, params::Dict{String,<:Any})

Request programmatic credentials to use with Habanero SDK.

# Arguments
- `environment_id`: The habanero environment identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"durationInMinutes"`: The time duration in which the credentials remain valid.
"""
get_programmatic_access_credentials(environmentId; aws_config::AbstractAWSConfig=global_aws_config()) = finspace_data("GET", "/credentials/programmatic", Dict{String, Any}("environmentId"=>environmentId); aws_config=aws_config)
get_programmatic_access_credentials(environmentId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace_data("GET", "/credentials/programmatic", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("environmentId"=>environmentId), params)); aws_config=aws_config)

"""
    get_working_location()
    get_working_location(params::Dict{String,<:Any})

A temporary Amazon S3 location to copy your files from a source location to stage or use as
a scratch space in Habanero notebook.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"locationType"`: Specify the type of the working location.    SAGEMAKER - Use the Amazon
  S3 location as a temporary location to store data content when working with FinSpace
  Notebooks that run on SageMaker studio.    INGESTION - Use the Amazon S3 location as a
  staging location to copy your data content and then use the location with the changeset
  creation operation.
"""
get_working_location(; aws_config::AbstractAWSConfig=global_aws_config()) = finspace_data("POST", "/workingLocationV1"; aws_config=aws_config)
get_working_location(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = finspace_data("POST", "/workingLocationV1", params; aws_config=aws_config)


text:
type XPCallbacks
    start_cdata::Function
    end_cdata::Function
    comment::Function
    character_data::Function
    default::Function
    default_expand::Function
    start_element::Function
    end_element::Function
    start_namespace::Function
    end_namespace::Function
end

type XPStreamHandler{D}
    cbs::XPCallbacks
    parser::XML_Parser
    data::D
end

# Create an XPCallbacks instance filled with nop callbacks
XPCallbacks() = XPCallbacks(
    (handler::XPStreamHandler) -> nothing,
    (handler::XPStreamHandler) -> nothing,
    (handler::XPStreamHandler, txt::AbstractString) -> nothing,
    (handler::XPStreamHandler, txt::AbstractString) -> nothing,
    (handler::XPStreamHandler, txt::AbstractString) -> nothing,
    (handler::XPStreamHandler, txt::AbstractString) -> nothing,
    (handler::XPStreamHandler, name::AbstractString, attrs_in::Dict{AbstractString,AbstractString}) -> nothing,
    (handler::XPStreamHandler, name::AbstractString) -> nothing,
    (handler::XPStreamHandler, prefix::AbstractString, uri::AbstractString) -> nothing,
    (handler::XPStreamHandler, prefix::AbstractString) -> nothing
)

function streaming_start_cdata(p_cbs::Ptr{Void})
    @DBG_PRINT("Found StartCdata")
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler

    h.cbs.start_cdata(h)
    return
end
cb_streaming_start_cdata = cfunction(streaming_start_cdata, Void, (Ptr{Void},))

function streaming_end_cdata(p_cbs::Ptr{Void})
    @DBG_PRINT("Found EndCdata")
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler

    h.cbs.end_cdata(h)
    return;
end
cb_streaming_end_cdata = cfunction(streaming_end_cdata, Void, (Ptr{Void},))

function streaming_cdata(p_cbs::Ptr{Void}, s::Ptr{UInt8}, len::Cint)
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler

    txt = unsafe_string(s, @compat(Int(len)))

    @DBG_PRINT("Found CData : " * txt)
    h.cbs.character_data(h, txt)

    return;
end
cb_streaming_cdata = cfunction(streaming_cdata, Void, (Ptr{Void},Ptr{UInt8}, Cint))

function streaming_start_element(p_cbs::Ptr{Void}, name::Ptr{UInt8}, attrs_in::Ptr{Ptr{UInt8}})
    h = unsafe_pointer_to_objref(p_cbs)::XPStreamHandler
    txt::AbstractString = unsafe_string(name)
    attrs::Dict{AbstractString,AbstractString} = attrs_in_to_dict(attrs_in)

    h.cbs.start_element(h, txt, attrs)

    return
end
cb_streaming_start_element = cfunction(streaming_start_element, Void, (Ptr{Void},Ptr{UInt8}, Ptr{Ptr{UInt8}}))

function streaming_end_element(p_h::Ptr{Void}, name::Ptr{UInt8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt::AbstractString = unsafe_string(name)
    @DBG_PRINT("End element: $txt, current element: $(xph.pdata.name) ")

    h.cbs.end_element(h, txt)

    return
end
cb_streaming_end_element = cfunction(streaming_end_element, Void, (Ptr{Void},Ptr{UInt8}))

function streaming_comment(p_h::Ptr{Void}, data::Ptr{UInt8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt = unsafe_string(data)
    @DBG_PRINT("Found comment : " * txt)

    h.cbs.comment(h, txt)

    return
end
cb_streaming_comment = cfunction(streaming_comment, Void, (Ptr{Void},Ptr{UInt8}))

function streaming_default(p_h::Ptr{Void}, data::Ptr{UInt8}, len::Cint)
    xph = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt = unsafe_string(data)
    @DBG_PRINT("Default : " * txt)

    h.cbs.default(h, txt)

    return;
end
cb_streaming_default = cfunction(streaming_default, Void, (Ptr{Void},Ptr{UInt8}, Cint))

function streaming_default_expand(p_h::Ptr{Void}, data::Ptr{UInt8}, len::Cint)
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    txt = unsafe_string(data)
    @DBG_PRINT("Default Expand : " * txt)

    h.cbs.default_expand(h, txt)

    return;
end
cb_streaming_default_expand = cfunction(streaming_default_expand, Void, (Ptr{Void},Ptr{UInt8}, Cint))

function streaming_start_namespace(p_h::Ptr{Void}, prefix::Ptr{UInt8}, uri::Ptr{UInt8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    prefix = unsafe_string(prefix)
    uri = unsafe_string(uri)
    @DBG_PRINT("start namespace prefix : $prefix, uri: $uri")

    h.cbs.start_namespace(h, prefix, uri)

    return;
end
cb_streaming_start_namespace = cfunction(streaming_start_namespace, Void, (Ptr{Void},Ptr{UInt8}, Ptr{UInt8}))

function streaming_end_namespace(p_h::Ptr{Void}, prefix::Ptr{UInt8})
    h = unsafe_pointer_to_objref(p_h)::XPStreamHandler
    prefix = unsafe_string(prefix)
    @DBG_PRINT("end namespace prefix : $prefix")

    h.cbs.end_namespace(h, prefix)

    return;
end
cb_streaming_end_namespace = cfunction(streaming_end_namespace, Void, (Ptr{Void},Ptr{UInt8}))

# Unsupported callbacks: External Entity, NotationDecl, Not Stand Alone, Processing, UnparsedEntityDecl, StartDocType
# SetBase and GetBase

function make_parser(cbs::XPCallbacks,data=nothing,sep='\0')
    p::XML_Parser = (sep == '\0') ? XML_ParserCreate(C_NULL) : XML_ParserCreateNS(C_NULL, sep);
    if (p == C_NULL) error("XML_ParserCreate failed") end

    h = XPStreamHandler(cbs, p, data)
    p_h = pointer_from_objref(h)
    XML_SetUserData(p, p_h);

    XML_SetCdataSectionHandler(p, cb_streaming_start_cdata, cb_streaming_end_cdata)
    XML_SetCharacterDataHandler(p, cb_streaming_cdata)
    XML_SetCommentHandler(p, cb_streaming_comment)
    XML_SetDefaultHandler(p, cb_streaming_default)
    XML_SetDefaultHandlerExpand(p, cb_streaming_default_expand)
    XML_SetElementHandler(p, cb_streaming_start_element, cb_streaming_end_element)
#    XML_SetExternalEntityRefHandler(p, f_ExternaEntity)
    XML_SetNamespaceDeclHandler(p, cb_streaming_start_namespace, cb_streaming_end_namespace)
#    XML_SetNotationDeclHandler(p, f_NotationDecl)
#    XML_SetNotStandaloneHandler(p, f_NotStandalone)
#    XML_SetProcessingInstructionHandler(p, f_ProcessingInstruction)
#    XML_SetUnparsedEntityDeclHandler(p, f_UnparsedEntityDecl)
#    XML_SetStartDoctypeDeclHandler(p, f_StartDoctypeDecl)

    return h
end

function stop(h::XPStreamHandler)
    XML_StopParser(h.parser, XML_FALSE)
#    XML_ParserFree(h.parser)
end

function pause(h::XPStreamHandler)
    XML_StopParser(h.parser, XML_TRUE)
end

function resume(h::XPStreamHandler)
    XML_ResumeParser(h.parser)
end

function free(h::XPStreamHandler)
    XML_ParserFree(h.parser)
end

function parsefile(filename::AbstractString,callbacks::XPCallbacks; bufferlines=1024, data=nothing)
    h = make_parser(callbacks, data)
    # TODO: Support suspending for files too
    suspended = false
    file = open(filename, "r")
    try
        io = IOBuffer()
        while !eof(file)
            i::Int = 0
            truncate(io, 0)
            while i < bufferlines && !eof(file)
                write(io, readline(file))
                i += 1
            end
            txt = Compat.String(io)
            rc = XML_Parse(h.parser, txt, length(txt.data), 0)
            if (rc != XML_STATUS_OK) && (XML_GetErrorCode(h.parser) != XML_ERROR_ABORTED)
                # Do not fail if the user aborted the parsing
                error("Error parsing document : $rc")
            end
            if XML_GetErrorCode(h.parser) == XML_ERROR_ABORTED
                break
            end
        end
        rc = XML_Parse(h.parser, "", length(""), 1)
        #if (rc == XML_STATUS_SUSPENDED)
        #    suspended = true
        #    return XPStreamHandler(callbacks,  parser)
        #end
        if (rc != XML_STATUS_OK) && (XML_GetErrorCode(h.parser) != XML_ERROR_ABORTED)
            # Do not fail if the user aborted the parsing
            error("Error parsing document : $rc")
        end
    catch e
        stre = string(e)
        (err, line, column, pos) = xp_geterror(h.parser)
        @DBG_PRINT("$e, $err, $line, $column, $pos")
        rethrow("$e, $err, $line, $column, $pos")
    finally
        if !suspended
            XML_ParserFree(h.parser)
        end
        close(file)
    end
end

function parse(txt::AbstractString,callbacks::XPCallbacks; data=nothing)
    h = make_parser(callbacks, data)
    suspended = false

    try
        rc = XML_Parse(h.parser, txt, length(txt.data), 1)
        if (rc == XML_STATUS_SUSPENDED)
            suspended = true
            return h
        end
        if (rc != XML_STATUS_OK) && (XML_GetErrorCode(h.parser) != XML_ERROR_ABORTED)
            # Do not fail if the user aborted the parsing
            error("Error parsing document : $rc")
        end
    catch e
        stre = string(e)
        (err, line, column, pos) = xp_geterror(h.parser)
        @DBG_PRINT("$e, $err, $line, $column, $pos")
        rethrow("$e, $err, $line, $column, $pos")

    finally
        if !suspended
            XML_ParserFree(h.parser)
        end
    end
end


text:
module Blobs

export Blob

using ..Elements
import ..Motions: induce_velocity, mutually_induce_velocity!, self_induce_velocity!, advect

#== Type definition ==#

"""
    Blob <: Elements.Element

An immutable structure representing a regularized point source/vortex

## Fields
- `z`: position
- `S`: strength/circulation
- `δ`: regularization radius

## Constructors

"""
struct Blob{T <: Number} <: Element
    z::ComplexF64
    S::T
    δ::Float64
    Blob{T}(z, s::Real, δ) where T <: Complex = new(z, im*s, δ)
    Blob{T}(z, s, δ)       where T <: Complex = new(z, s, δ)
    Blob{T}(z, s::Real, δ) where T <: Real    = new(z, s, δ)
end
Elements.kind(::Blob) = Singleton
Elements.kind(::Type{Blob{T}}) where T = Singleton

(b::Blob{T})(; z = b.z, S = b.S, δ = b.δ) where T = Blob{T}(z, S, δ)

#== Methods to be extended ==#

Elements.position(b::Blob) = b.z
Elements.streamfunction(z::ComplexF64, b::Blob) = real(-0.5b.S*log(z - b.z)/π)

Elements.complexpotential(z::ComplexF64, b::Blob) = -0.5im*b.S*log(z - b.z)/π

blob_kernel(z, δ) = 0.5im*z/(π*(abs2(z) + δ^2))

function induce_velocity(z::ComplexF64, b::Blob, t)
    b.S*blob_kernel(z - b.z, b.δ)
end

function induce_velocity(target::Blob, source::Blob, t)
    δ = √(0.5(target.δ^2 + source.δ^2))
    source.S'*blob_kernel(target.z - source.z, δ)
end

function mutually_induce_velocity!(ws₁, ws₂,
                                   blobs₁::Vector{Blob{T₁}},
                                   blobs₂::Vector{Blob{T₂}}, t) where {T₁, T₂}
    for (s, source) in enumerate(blobs₁)
        for (t, target) in enumerate(blobs₂)
            δ = √(0.5(target.δ^2 + source.δ^2))
            K = blob_kernel(target.z - source.z, δ)
            ws₂[t] += source.S'*K
            ws₁[s] -= target.S'*K
        end
    end
    nothing
end

function self_induce_velocity!(ws, blobs::Vector{Blob{T}}, t) where T
    N = length(blobs)

    for s in 1:N, t in s+1:N
        δ = √(0.5(blobs[t].δ^2 + blobs[s].δ^2))
        K = blob_kernel(blobs[t].z -blobs[s].z, δ)
        ws[t] += blobs[s].S'*K
        ws[s] -= blobs[t].S'*K
    end
    ws
end

function advect(b::Blob{T}, w::ComplexF64, Δt::Float64) where T
    Blob{T}(b.z + w*Δt, b.S, b.δ)
end
#
#function Base.show(io::IO, b::Blob)
#    print(io, "Vortex Blob: z = $(round(b.z, digits=3)), Γ = $(round(b.Γ, digits=3)), δ = $(round(b.δ, digits=3))")
#end

end


text:
immutable Logistic <: ContinuousUnivariateDistribution
    location::Real
    scale::Real
    function Logistic(l::Real, s::Real)
    	s > zero(s) || error("scale must be positive")
    	new(float64(l), float64(s))
    end
    Logistic(l::Real) = new(float64(l), 1.0)
    Logistic() = new(0.0, 1.0)
end

## Support
@distr_support Logistic -Inf Inf

## Properties
mean(d::Logistic) = d.location
median(d::Logistic) = d.location
mode(d::Logistic) = d.location

std(d::Logistic) = pi * d.scale / sqrt(3.0)
var(d::Logistic) = (pi * d.scale)^2 / 3.0

skewness(d::Logistic) = 0.0
kurtosis(d::Logistic) = 1.2

entropy(d::Logistic) = log(d.scale) + 2.0

## Functions
function pdf(d::Logistic, x::Float64)
    a = exp(-abs((x-d.location)/d.scale))
    a / (d.scale * (1+a)^2)
end
function logpdf(d::Logistic, x::Float64)
    u = -abs((x-d.location)/d.scale)
    u - 2*log1pexp(u) - log(d.scale)
end

cdf(d::Logistic, x::Float64) = logistic((x-d.location)/d.scale)
ccdf(d::Logistic, x::Float64) = logistic((d.location-x)/d.scale)
logcdf(d::Logistic, x::Float64) = -log1pexp((d.location-x)/d.scale)
logccdf(d::Logistic, x::Float64) = -log1pexp((x-d.location)/d.scale)

quantile(d::Logistic, p::Float64) = @checkquantile p d.location + d.scale*logit(p)
cquantile(d::Logistic, p::Float64) = @checkquantile p d.location - d.scale*logit(p)
invlogcdf(d::Logistic, lp::Float64) = d.location - d.scale*logexpm1(-lp)
invlogccdf(d::Logistic, lp::Float64) = d.location + d.scale*logexpm1(-lp)

function gradlogpdf(d::Logistic, x::Float64)
  expterm = exp((d.location - x) / d.scale)
  ((2 * expterm) / (1 + expterm) - 1) / d.scale
end

mgf(d::Logistic, t::Real) = exp(t*d.location)/sinc(d.scale*t)

function cf(d::Logistic, t::Real)
    a = (pi*t)*d.scale
    a == zero(a) ? complex(one(a)) : exp(im*t*d.location) * a / sinh(a)
end


text:
            Angel News

            new

            show

            ask

            submit

            login

{{ user.username }} ({{ user.karma }})

logout


text:
@doc raw"""
Simple perturbation in electric field in (x,y,z) coordinates:
```math
E(x,y,z) = E_0 \, \begin{pmatrix}
0 \\
0 \\
\cos (2\pi z) \\
\end{pmatrix}
```

Parameters: `E₀`
"""
struct EzCosZ{T <: Number} <: CartesianPerturbation
    name::String
    E₀::T
    EzCosZ{T}(E₀::T) where T <: Number = new("EzCosZ", E₀)
end

EzCosZ(E₀::T) where T <: Number = EzCosZ{T}(E₀)

function Base.show(io::IO, equ::EzCosZ)
    print(io, "Simple perturbation in electric field")
end

function φ(x::AbstractArray{T,1}, equ::EzCosZ) where {T <: Number}
   equ.E₀ / (2π) * sin(2π * Z(x,equ))
end

macro ezcosz_perturbation(E₀=1.)
    generate_equilibrium_code(EzCosZ(E₀); output=false)
end


text:
"""
    AbstractBasis

Abstract supertype for sets of basis functions.
"""
abstract type AbstractBasis end

Base.eltype(B::AbstractBasis,x) = eltype(collect(B|x))

"""
    |(B::AbstractBasis, x)

Iterate over the functions in `B` evaluated at point `x`.

# Examples
```
collect( Monomials(4) | 2 ) == [1,2,4,8]
```
"""
Base.:|(B::AbstractBasis,x) = BasisValues(B,wrap(x))

function Base.getindex(B::AbstractBasis, i)
    @assert i in 1:length(B)
    return BasisFunction(B,i)
end

"""
    Matrix(B::AbstractBasis, x) = [ B[j](x[i]) for i = 1:length(x), j = 1:length(B) ]
"""
function Base.Matrix(B::AbstractBasis, x)
    M = Matrix{eltype(B,x[1])}(undef, length(x),length(B))
    for i = 1:length(x)
        copyto!(@view(M[i,:]), B|x[i])
    end
    return M
end

"""
    evaltransform(B,x,c)

Evaluate `Matrix(B,x)*c`.
"""
evaltransform(B::AbstractBasis, x, c) = default_evaltransform(B,x,c)
default_evaltransform(B::AbstractBasis, x, c) = Matrix(B,x)*c
default_evaltransform(B::AbstractBasis, x, c::AbstractVector) = mapreduce(((ci,bix),) -> ci*bix, +, zip(c,B|x))
default_evaltransform(B::AbstractBasis, x::AbstractVector{<:Number}, c::AbstractVector) = evaltransform.(Ref(B),x,Ref(c))

"""
    BasisValues(B,x)

Auxiliary type returned by `B | x`.
"""
struct BasisValues{B<:AbstractBasis,X}
    basis::B
    point::X
end

Base.length(Bx::BasisValues) = length(Bx.basis)
Base.IteratorEltype(::Type{<:BasisValues}) = Base.EltypeUnknown()
Base.iterate(Bx::BasisValues, args...) = iterate_basis(Bx.basis, Bx.point, args...)

"""
    iterate_basis(B::AbstractBasis, x [, state])

Iterate over the functions in `B` evaluated at point `x`.
Analogous to `Base.iterate`.
"""
function iterate_basis(B::AbstractBasis, x, i=1)
    i > length(B) && return nothing
    return B[i](x), i+1
end

"""
    BasisFunction(B,i)

Auxiliary type returned by `B[i]`.
"""
struct BasisFunction{B<:AbstractBasis}
    basis::B
    i::Int
end
(bf::BasisFunction)(x...) = bf(x)
(bf::BasisFunction)(x) = evaluate_basis(bf.basis, bf.i, wrap(x))
evaluate_basis(B::AbstractBasis, i, x) = nth(B|x,i)

struct LinearCombination{N,C<:AbstractArray{<:Number,N},B<:NTuple{N,AbstractBasis}}
    coeffs::C
    basis::B
    function LinearCombination{N,C,B}(coeffs,basis) where {N,C,B}
        @assert size(coeffs) == length.(basis)
        return new{N,C,B}(coeffs,basis)
    end
end

"""
    LinearCombination(c::AbstractVector, B::AbstractBasis) -> p
    LinearCombination(c::AbstractArray{N}, B::NTuple{N,AbstractBasis}) -> p

Linear combination of basis functions.

# Examples
```
julia> x = LinRange(-1,1,11)
       p = LinearCombination([1,0,2], Monomials(3))
       p.(x) ≈ @.( 1 + 2*x^2 )
true
```
"""
function LinearCombination end
LinearCombination(c::AbstractArray{<:Number,N}, B::NTuple{N,AbstractBasis}) where {N} = LinearCombination{N,typeof(c),typeof(B)}(c,B)
LinearCombination(c::AbstractArray{<:Number,N}, B::Vararg{AbstractBasis,N}) where {N} = LinearCombination(c,B)
LinearCombination(c::AbstractArray{<:Number,N}, B::AbstractBasis) where {N} = LinearCombination(c,ntuple(i->B,Val(N)))

coeffs(c::LinearCombination) = c.coeffs
basis(c::LinearCombination) = c.basis
Base.ndims(c::LinearCombination{N}) where {N} = N
Base.ndims(::Type{<:LinearCombination{N}}) where {N} = N

function evaluate_linear_combination(
    c::AbstractArray{<:Any,N},
    B::NTuple{N,AbstractBasis},
    x::NTuple{N,Any}
) where {N}
    @assert size(c) == length.(B)
    tucker(c, map((B,x)->(c->evaltransform(B,wrap(x),c)), B,x))
end

(p::LinearCombination{N})(x::Vararg{Any,N}) where {N} = evaluate_linear_combination(coeffs(p), basis(p), x)
(p::LinearCombination{N})(x::Vararg{Number,N}) where {N} = first(evaluate_linear_combination(coeffs(p), basis(p), x))

grideval(p::LinearCombination{N}, x::NTuple{N,Any}) where {N} = p(x...)

using MacroTools

"""
    @evaluate expr

Currently supported expressions:
 - `@evaluate p(M)*v`: Apply `p(M)` to `v` efficiently.
"""
macro evaluate(expr)
     @capture(expr, p_(M_)*v_) && return :(evaluate_pmtv($(esc(p)),$(esc(M)),$(esc(v))))
end

evaluate_pmtv(p::LinearCombination,M,v) = p(MatrixVectorWrapper(M,v))


text:
using Documenter, ObservationSchemes

makedocs(;
    modules=[ObservationSchemes],
    format=Documenter.HTML(
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                )
            )
        ),
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "Get started" => joinpath("get_started", "overview.md"),
        "User manual" => Any[
            "A single observation" => joinpath("manual", "single_observation.md"),
            "Priors over starting points" => joinpath("manual", "start_pt_prior.md"),
            "Multiple observations" => joinpath("manual", "multiple_observations.md"),
            "Observation schemes" => joinpath("manual", "observation_schemes.md"),
            "Utility fuctions" => joinpath("manual", "utility_functions.md"),
        ],
        "How to..." => Any[
            "... ergonomically deal with simple observation schemes?" => joinpath("how_to_guides", "simple_os.md"),
            "... define custom priors over starting points?" => joinpath("how_to_guides", "custom_priors.md"),
            "... define custom observation schemes?" => joinpath("how_to_guides", "custom_os.md"),
            "... work with dataframes?" => joinpath("how_to_guides", "dataframes.md"),
            "... read data directly from a `.csv` file?" => joinpath("how_to_guides", "csv_files.md"),
            "... save decorated data to a file?" => joinpath("how_to_guides", "save_to_files.md"),
            "... deal with stochastic processes other than diffusions?" => joinpath("how_to_guides", "other_processes.md"),
            "... leverage `AllObservations` to efficiently update parameters?" => joinpath("how_to_guides", "update_params.md"),
        ],
        "Tutorials" => Any[
            "Defining an observation scheme for a single trajectory of a FitzHugh–Nagumo model" => joinpath("tutorials", "single_path.md"),
            #"(TODO) Multiple paths" => joinpath("tutorials", "multi_paths.md"),
            #"(TODO) Multiple paths differently observed" => joinpath("tutorials", "multi_paths_diff_observ.md"),
            #"(TODO) Mixed effect models" => joinpath("tutorials", "mixed_effect_models.md"),
            #"(TODO) Observations of HMM driven by ARCH(1)" => joinpath("tutorials", "hmm.md"),
        ],
        "Index" => "module_index.md",
    ],
    repo="https://github.com/JuliaDiffusionBayes/ObservationSchemes.jl/blob/{commit}{path}#L{line}",
    sitename="ObservationSchemes.jl",
    authors="Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer",
    #assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/ObservationSchemes.jl",
)


text:
using StochasticDelayDiffEq
using Test

#Hayes Equation
begin
    function hayes_modelf(du,u,h,p,t)
        τ,a,b,c,α,β,γ = p
        du .= a.*u .+ b .* h(p,t-τ) .+ c
    end
    function hayes_modelg(du,u,h,p,t)
        τ,a,b,c,α,β,γ = p
        du .= α.*u .+ β.*h(p,t-τ) .+ γ
    end
    h(p,t) = (ones(1) .+ t);
    tspan = (0.,10.)
end

pmul = [1.0,-4.,-2.,10.,-1.3,-1.2, 1.1]
padd = [1.0,-4.,-2.,10.,-0.0,-0.0, 0.1]

prob = @test_nowarn SDDEProblem(hayes_modelf, hayes_modelg, [1.], h, tspan, pmul; constant_lags = (pmul[1],));
@test_nowarn SDDEProblem(hayes_modelf, hayes_modelg, h, tspan, padd; constant_lags = (padd[1],));

sol = @test_nowarn solve(prob,EM(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,LambaEM(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,EulerHeun(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,LambaEulerHeun(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,RKMil(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,RKMil(interpretation=:Stratonovich),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,RKMilCommute(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,RKMilCommute(interpretation=:Stratonovich),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,WangLi3SMil_A(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,WangLi3SMil_B(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,WangLi3SMil_C(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,WangLi3SMil_D(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,WangLi3SMil_E(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,WangLi3SMil_F(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SRI(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SRIW1(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SRIW2(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SOSRI(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SOSRI2(),dt=0.01)
@test sol.u[end] != zeros(1)

# Additive problems
println("additive problems")
prob.p .= padd;
sol = @test_nowarn solve(prob,SRA(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SRA1(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SRA2(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SRA3(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SOSRA(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SOSRA2(),dt=0.01)
@test sol.u[end] != zeros(1)
# @test_nowarn solve(prob,SKenCarp(),dt=0.01) # Not working

# Test SROCK methods
println("SROCK methods")
prob.p .= pmul;
sol = @test_nowarn solve(prob,SROCK1(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SROCK1(interpretation=:Stratonovich),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SROCKEM(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SROCKEM(strong_order_1=false),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SROCK2(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SKSROCK(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SKSROCK(;post_processing=true),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,TangXiaoSROCK2(),dt=0.01)
@test sol.u[end] != zeros(1)
for i in 1:5
    sol = @test_nowarn solve(prob,TangXiaoSROCK2(version_num=i),dt=0.01)
    @test sol.u[end] != zeros(1)
end

# Test Implicit methods
println("implicit methods")
sol = @test_nowarn solve(prob,ImplicitEM(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ImplicitEM(symplectic=true, theta = 1/2),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ImplicitEulerHeun(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ImplicitEulerHeun(symplectic=true, theta = 1/2),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ImplicitRKMil(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ImplicitRKMil(symplectic=true, theta = 1/2),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ImplicitRKMil(interpretation=:Stratonovich, symplectic = true, theta = 1/2),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ISSEM(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ISSEM(symplectic=true, theta = 1/2),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ISSEulerHeun(),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,ISSEulerHeun(symplectic=true, theta = 1/2),dt=0.01)
@test sol.u[end] != zeros(1)
sol = @test_nowarn solve(prob,SKenCarp(),dt=0.01)
@test sol.u[end] != zeros(1)


text:
using UnifiedPseudopotentialFormat.PSlibrary

@testset "Test `potential_exists`" begin
    @test potential_exists("S.rel-pz-n-rrkjus_psl.0.1.UPF")
    @test potential_exists("S.pbe-nl-rrkjus_psl.1.0.0.UPF")
    @test potential_exists("S.rel-pbesol-n-kjpaw_psl.0.1.UPF")
end

@testset "Test `download_potential`" begin
    download_potential("S.rel-pz-n-rrkjus_psl.0.1.UPF", "./S.rel-pz-n-rrkjus_psl.0.1.UPF")
    @test isfile("./S.rel-pz-n-rrkjus_psl.0.1.UPF")
end


text:
"""
    redirect_output(fun)

Create an `IOBuffer`, pass it to `fun`, and return the captured output
as a string.

# Examples

```jldoctest
julia> redirect_output() do io
           println(io, "Hello")
           println(io, "World")
       end
"Hello\\nWorld\\n"
```
"""
function redirect_output(fun::Function)
    buf = IOBuffer()
    fun(IOContext(buf, :color => true))
    String(take!(buf))
end

"""
    horizontal_line([io=stdout; char="━", color=:light_black])

Draw a `color`ed horizontal line of `char`s across the whole screen.

# Example

```jldoctest
julia> horizontal_line()
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

julia> horizontal_line(char="─")
────────────────────────────────────────────────────────────────────────────────

```
"""
function horizontal_line(io::IO=stdout; char="━", color=:light_black)
    ncol = displaysize(io)[2]
    printstyled(io, repeat(char, ncol), color=color)
    println(io)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
   PoissonProcess(λ)

A Poisson process with intensity `λ`.
"""
struct PoissonProcess{L<:Union{Real,Function}} <: PointProcess
  λ::L
end

Base.union(p₁::PoissonProcess{<:Real}, p₂::PoissonProcess{<:Real}) =
  PoissonProcess(p₁.λ + p₂.λ)

Base.union(p₁::PoissonProcess{<:Function}, p₂::PoissonProcess{<:Function}) =
  PoissonProcess(x -> p₁.λ(x) + p₂.λ(x))

Base.union(p₁::PoissonProcess{<:Real}, p₂::PoissonProcess{<:Function}) =
  PoissonProcess(x -> p₁.λ + p₂.λ(x))

Base.union(p₁::PoissonProcess{<:Function}, p₂::PoissonProcess{<:Real}) =
  PoissonProcess(x -> p₁.λ(x) + p₂.λ)

ishomogeneous(p::PoissonProcess{<:Real}) = true
ishomogeneous(p::PoissonProcess{<:Function}) = false

default_sampling_algorithm(::PoissonProcess, ::Geometry) = DiscretizedSampling()
default_sampling_algorithm(::PoissonProcess{<:Real}, ::Box) = ProductSampling()

#------------------
# HOMOGENEOUS CASE
#------------------

function rand_single(rng::Random.AbstractRNG,
                     p::PoissonProcess{<:Real}, g::Geometry,
                     ::DiscretizedSampling)
  # simulate number of points
  λ = p.λ; V = measure(g)
  n = rand(rng, Poisson(λ*V))

  # simulate homogeneous process
  pts = sample(g, HomogeneousSampling(n))

  # return point pattern
  PointSet(collect(pts))
end

function rand_single(rng::Random.AbstractRNG,
                     p::PoissonProcess{<:Real}, b::Box,
                     ::ProductSampling)
  # region configuration
  lo, up = coordinates.(extrema(b))

  # simulate number of points
  λ = p.λ; V = measure(b)
  n = rand(rng, Poisson(λ*V))

  # product of uniform distributions
  U = product_distribution([Uniform(lo[i], up[i]) for i in 1:embeddim(b)])

  # return point pattern
  PointSet(rand(rng, U, n))
end

#--------------------
# INHOMOGENEOUS CASE
#--------------------

function rand_single(rng::Random.AbstractRNG,
                     p::PoissonProcess{<:Function}, b::Box,
                     algo::DiscretizedSampling)
  # discretize region
  # TODO

  # discretize retention
  # TODO

  # sample each element
  # TODO
end

function rand_single(rng::Random.AbstractRNG,
                     p::PoissonProcess{<:Function}, b::Box,
                     algo::ThinnedSampling)
  # Lewis-Shedler algorithm
  # TODO
end


text:
using Base.Test
using SparseRF
using ScikitLearnBase

srand(2)
n,m = 10^3, 5 ;
features = rand(n,m);
weights = rand(-1:1,m);
labels = _int(features * weights);

# I wish we could use ScikitLearn.jl's cross-validation, but that'd require
# installing it on Travis
model = fit!(DecisionTreeClassifier(pruning_purity_threshold=0.9), features, labels)
@test mean(predict(model, features) .== labels) > 0.8

model = fit!(RandomForestClassifier(), features, labels)
@test mean(predict(model, features) .== labels) > 0.8

model = fit!(AdaBoostStumpClassifier(), features, labels)
# Adaboost isn't so hot on this task, disabled for now
mean(predict(model, features) .== labels)

srand(2)
N = 3000
X = randn(N, 10)
# TODO: we should probably support fit!(::DecisionTreeClassifier, ::BitArray)
y = convert(Vector{Bool}, randn(N) .< 0)
maxdepth = 5
model = fit!(DecisionTreeClassifier(maxdepth=maxdepth), X, y)
@test depth(model) == maxdepth

## Test that the RNG arguments work as expected
srand(2)
X = randn(100, 10)
y = rand(Bool, 100);

@test predict_proba(fit!(RandomForestClassifier(; rng=10), X, y), X) ==
    predict_proba(fit!(RandomForestClassifier(; rng=10), X, y), X)
@test predict_proba(fit!(RandomForestClassifier(; rng=10), X, y), X) !=
    predict_proba(fit!(RandomForestClassifier(; rng=12), X, y), X)


text:
module GraphingUtils
using PlotlyJS

export graph_all

function graph_all(x, y, z, Fx::Function, Fy::Function, Fz::Function,
                x_min::Number, x_max::Number,
                y_min::Number, y_max::Number,
                z_min::Number, z_max::Number,
                density::Int64, z₀::AbstractArray = [NaN])::Plot

    x_span = range(x_min, stop=x_max, length=density)
    y_span = range(y_min, stop=y_max, length=density)
    z_span = range(z_min, stop=z_max, length=density)
    𝓍 = repeat(x_span, density^3)
    𝓎 = repeat(y_span, inner=density, outer=density^2)
    𝓏 = repeat(z_span, inner=density^3)

    trace₁ = surface(;x=x, y=y, z=z)
    trace₂ = cone(;x=𝓍, y=𝓎, z=𝓏, u=Fx.(𝓍, 𝓎, 𝓏), v=Fy.(𝓍, 𝓎, 𝓏), w=Fz.(𝓍, 𝓎, 𝓏), showscale=false)
    layout = Layout(autosize=false, width=600, height=600)
    if !all(isnan.(z₀))
        trace₀ = surface(;x=x, y=y, z=z₀, showscale = false)
        return Plot([trace₀, trace₁, trace₂], layout)
    else
        return Plot([trace₁, trace₂], layout)
    end
end

# check(u, v, f) = v_min(u) <= v && v_max(u) >= v ? f(u, v) : NaN

end


text:
@testset "HashIndices" begin
    @test HashIndices() isa HashIndices{Any}

    h = HashIndices{Int64}()

    @test isinsertable(h)
    @test length(h) == 0
    @test keys(h) === h
    @test h == h
    @test h == copy(h)
    @test isempty(h)
    @test isequal(copy(h), h)
    @test_throws IndexError h[10]
    @test length(unset!(h, 10)) == 0
    io = IOBuffer(); print(io, h); @test String(take!(io)) == "{}"
    io = IOBuffer(); show(io, MIME"text/plain"(), h); @test String(take!(io)) == "0-element HashIndices{Int64}"
    @test_throws IndexError delete!(h, 10)

    insert!(h, 10)

    @test length(h) == 1
    @test keys(h) === h
    @test unique(h) === h
    @test h == h
    @test h == copy(h)
    @test !isempty(h)
    @test isequal(copy(h), h)
    @test h[10] == 10
    @test_throws IndexError insert!(h, 10)
    @test length(set!(h, 10)) == 1
    @test_throws IndexError insert!(h, 10)
    io = IOBuffer(); print(io, h); @test String(take!(io)) == "{10}"
    io = IOBuffer(); show(io, MIME"text/plain"(), h); @test String(take!(io)) == "1-element HashIndices{Int64}\n 10"
    @test !isequal(h, empty(h))
    @test isequal(h, copy(h))
    @test isempty(empty(h))

    delete!(h, 10)

    @test isequal(h, HashIndices{Int64}())

    for i = 2:2:1000
        insert!(h, i)
    end
    @test issetequal(h, HashIndices(2:2:1000))
    @test all(in(i, h) == iseven(i) for i in 2:1000)
    @test isempty(empty!(h))

    # TODO: token interface
end


text:
### wiener.jl ---  Wiener deconvolution
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: wiener, deconvolution, signal processing
#
# This file is a part of Deconvolution.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file provides the `wiener' function to perform Wiener deconvolution.  It
# has two modes: one for the case of presence of a known blurring, the other one
# without blurring.
#
### Code:

export wiener

### The implementation of the deconvolution algorithm

function _wiener_no_blur(Y::AbstractArray, # Fourier transform of the input
                         S::AbstractArray, # Power spectrum of the signal
                         N::AbstractArray) # Power spectrum of the noise
    # Without blurring, the filter is:
    #   |S|² / (|S|² + |N|²)
    return real(ifft(Y .* S ./ (S .+ N)))
end

function _wiener_with_blur(Y::AbstractArray, # Fourier transform of the input
                           S::AbstractArray, # Power spectrum of the signal
                           N::AbstractArray, # Power spectrum of the noise
                           H::AbstractArray) # Fourier transform of the blurring
    # With blurring, the filter is:
    #   H* / (|H|² + |N/S|²)
    return real(ifft(Y .* conj(H) ./ (abs2.(H) .+ N ./ S)))
end

### User interface

## Without blurring
function wiener(input::AbstractArray, signal::AbstractArray,
                noise::AbstractArray)
    @assert size(input) == size(signal) == size(noise)
    input_ft = fft(input)
    signal_power_spectrum = abs2.(fft(signal))
    noise_power_spectrum = abs2.(fft(noise))
    return _wiener_no_blur(input_ft, signal_power_spectrum,
                           noise_power_spectrum)
end

# Noise as a real (it's converted to an array)
wiener(input::AbstractArray, signal::AbstractArray, noise::Real) =
    wiener(input, signal, fill!(similar(input), one(eltype(input)))*noise)

## With blurring
function wiener(input::AbstractArray, signal::AbstractArray,
                noise::AbstractArray, blurring::AbstractArray)
    @assert size(input) == size(signal) == size(noise) == size(blurring)
    input_ft = fft(input)
    signal_power_spectrum = abs2.(fft(signal))
    noise_power_spectrum = abs2.(fft(noise))
    blurring_ft = fft(blurring)
    return _wiener_with_blur(input_ft, signal_power_spectrum,
                             noise_power_spectrum, blurring_ft)
end

# Noise as a real (it's converted to an array)
wiener(input::AbstractArray, signal::AbstractArray,
       noise::Real, blurring::AbstractArray) =
           wiener(input, signal, fill!(similar(input), one(eltype(input)))*noise, blurring)

"""
    wiener(input, signal, noise[, blurring])

Return the Wiener deconvolution of `input`, using the power spectrum of `signal`
and `noise`.  If the `input` was blurred with a known blurring kernel, pass it
as fourth argument, `blurring`.

All arguments must be arrays in the time/space domain and all of the same size,
they will be converted into the frequency domain internally using the `fft`
function.
"""
wiener


text:
function model_3Cat()

    data = getData(data_path, "roentgen.jld")

    x_data = data["x"]
    y_data = data["y"]

    ratios, levels, returnClassImbalance(y_data)

    println("We have classes ", levels)
    println("We have ratios ", ratios)

    y_data = Flux.onehotbatch(vec(y_data), levels)

    x_data = processX(x_data)

    xs, ys = x_data[:, perm[1:Int(round(size(x_data)[2] * 0.9))]], y_data[:, perm[1:Int(round(size(x_data)[2] * 0.9))]]
    x_test, y_test = x_data[: , perm[Int(round(size(x_data)[2] * 0.9)) + 1:end]], Y_Data[: , perm[Int(round(size(x_data)[2] * 0.9)) + 1:end]]

    train_3cat_NN(xs, ys, 0.001)

    @load "NNRoentgen3Way.bson" model
    @load "NNRoentgen3Way.bson" weights

    Flux.loadparams!(model, weights)

    #gt = Int64.(Y_test)
    gt = classify(y_test)
    pred = classify(model(x_test))
    #pred = vec(Int64.(model(X_test)))

    confusmat(3, gt , pred) # not sure about this, we have two classes, but doesnt like it

end

# function getData()

#     fileLocation = "data/"
#     return load(fileLocation * "roentgen.jld")

# end
# #sort out class imbalance

# data = getData()

# X_Data = data["x"]
# Y_Data = data["y"]

# allPositivesInd = findall(x -> x == 1, Y_Data)
# allNegativesInd = findall(x -> x == 0, Y_Data)

# allPositivesX = X_Data[allPositivesInd]
# allNegativesX = X_Data[allNegativesInd]

# allNegativesX_gen = [clamp01nan!(Float64.(imrotate(Gray.(allNegativesX[i] .+ randn()/5), randn()/3, axes(allNegativesX[i])))) .+ rand()/10 for i in 1:length(allNegativesX)];
# allNegativesX = vcat(allNegativesX, allNegativesX_gen)
# furtherNeg = rand(allNegativesX, 824)
# allNegativesX_gen = [clamp01nan!(Float64.(imrotate(Gray.(furtherNeg[i] .+ randn()/5), randn()/3, axes(allNegativesX[i])))) .+ rand()/10 for i in 1:length(furtherNeg)];
# allNegativesX = vcat(allNegativesX, allNegativesX_gen)

# # Put the data toghether
# X_Data = vcat(allNegativesX, allPositivesX)
# Y_Data = Int.(vcat(zeros(length(allNegativesX)), ones(length(allPositivesX))))

# Y_Data = Flux.onehotbatch(vec(Y_Data), [0, 1, 2])

# for i in 1:length(X_Data)

#     X_Data[i] = reshape(X_Data[i], 1, :)

# end

# # create a matrix from an Array of Arrays
# X_Data = vcat(X_Data...)

# # we want one picture to have dimension m features times 1 so m*1 per picture,
# # then the total matrix x should be m*n where m is the length of features and n is the number of examples
# # so also y is a 1 * n matrix

# X_Data = transpose(X_Data)
# #Y_Data = transpose(Y_Data)

# perm = shuffle(1:size(X_Data)[2])

# X_train, Y_train = X_Data[:, perm[1:Int(round(size(X_Data)[2] * 0.9))]], Y_Data[:, perm[1:Int(round(size(X_Data)[2] * 0.9))]]
# X_test, Y_test = X_Data[: , perm[Int(round(size(X_Data)[2] * 0.9)) + 1:end]], Y_Data[: , perm[Int(round(size(X_Data)[2] * 0.9)) + 1:end]]

# # sum(Y_train[1,:]) / size(Y_train)[2]
# # sum(Y_test[1,:]) / size(Y_test)[2]

# # size(X_train)
# # size(X_test)
# # size(Y_train)
# # size(Y_test)

# data_xr = DataLoader(X_train, Y_train, batchsize=1)

# # data_xr.batchsize
# # data_xr.data
# # data_xr.imax
# # data_xr.indices
# # data_xr.nobs
# # data_xr.partial
# # data_xr.shuffle

# #method from baumann, check difference
# #data_xr_batch = (Flux.batch(X_train), Flux.batch(Y_train))

# model = Chain(Dense(size(data_xr.data[1])[1], 32, relu), Dense(32, 3, σ), softmax)

# params(model)

# #check model works:
# model(data_xr.data[1])

# loss(x, y) = Flux.crossentropy(model(x), y) #mse(model(x), y) # loss of the first m=1 examples
# #loss(x, y) = Flux.mse(model(x), y)

# #check loss function works
# loss(data_xr.data[1], data_xr.data[2])

# opt = Descent(0.001)
# @time Flux.train!(loss, params(model), data_xr, opt)

# for _ in 1:100

#     @time Flux.train!(loss, params(model), data_xr, opt)
#     println(loss(data_xr.data[1], data_xr.data[2]))

# end

# @save "NNRoentgen3Way.bson" model

# weights = params(model)

# @save "NNRoentgen3Way.bson" weights

#mean(model() .== Y_test)

# gt = Int64.(Y_test)
# gt = classify(Y_test)
# pred = classify(model(X_test))
# #pred = vec(Int64.(model(X_test)))

# correctrate(gt, pred)
# errorrate(gt, pred)

# confusmat(3, gt , pred) # not sure about this, we have two classes, but doesnt like it


text:
## Mesh class
@fenicsclass Mesh
hmin(mesh::Mesh) = fenicscall(mesh, :hmin)
hmax(mesh::Mesh) = fenicscall(mesh, :hmax)
rmin(mesh::Mesh) = fenicscall(mesh, :rmin)
rmax(mesh::Mesh) = fenicscall(mesh, :rmax)
cells(mesh::Mesh) = fenicscall(mesh, :cells)::Array{UInt32,2}
coordinates(mesh::Mesh) = fenicscall(mesh, :coordinates)::Array{Float64,2}
init(mesh::Mesh) = fenicscall(mesh, :init)
init(mesh::Mesh, dim::Int) = fenicscall(mesh, :init, dim)::Int
size(mesh::Mesh, dim::Int) = fenicscall(mesh, :size, dim)::Int
size(mesh::Mesh) = [size(mesh, dim) for dim = 0:3]
num_cells(mesh::Mesh) = fenicscall(mesh, :num_cells)
num_edges(mesh::Mesh) = fenicscall(mesh, :num_edges)
num_entities(mesh::Mesh, dim::Int) = fenicscall(mesh, :num_entities, dim)
num_faces(mesh::Mesh) = fenicscall(mesh, :num_faces)
num_facets(mesh::Mesh) = fenicscall(mesh, :num_facets)
num_vertices(mesh::Mesh) = fenicscall(mesh, :num_vertices)
# TODO: add more mesh functions
# closest_point(point::Point)::Point, closest_cell(point)::Int

export
  Mesh, hmin, hmax, rmin, rmax, init, coordinates, cells,
  num_cells, num_edges, num_entities, num_faces, num_facets, num_vertices

## Mesh constructors
UnitTriangleMesh() = Mesh(dolfin.UnitTriangleMesh())
UnitSquareMesh(comm::MPI_Comm, nx::Int, ny::Int; diagonal::String="right") = Mesh(dolfin.UnitSquareMesh(comm.pyobject, nx, ny, diagonal))
UnitSquareMesh(nx::Int, ny::Int; diagonal::String="right") = Mesh(dolfin.UnitSquareMesh(nx, ny, diagonal))
# TODO: add more mesh constructors (Box, Rectangle)
# TODO: add mesh loading functions

# TODO: sometime later: add MeshEditor (MeshTopology, MeshGeometry,...), Point
export
  UnitTriangleMesh,
  UnitSquareMesh


text:
struct SplitHLAModelResult
    f::Function
    max_idx::Int

    function SplitHLAModelResult(f, max_idx)
        f(1) isa HLAModelResult || error("f(1) does not return a HLAModelResult.")

        new(f, max_idx)
    end
end

function replacement_summary(
    result::SplitHLAModelResult;
    fdr::Bool = false,fisher_p::Bool = false
)
    df = DataFrame()

    @showprogress for split_idx in 1:result.max_idx
        df = vcat(df, replacement_summary(Base.invokelatest(result.f, split_idx),
            fisher_p = fisher_p))
    end

    sort!(df, :posterior_p, rev = true)

    return df
end

function SplitHLAModelResult(results)
    return SplitHLAModelResult(x -> getindex(results, x), length(results))
end

function Escape.loo(result::SplitHLAModelResult)
    pw = Vector{Loo.PointwiseLoo}()
    draws = 0

    @showprogress for split_idx in 1:result.max_idx
        split_result = Base.invokelatest(result.f, split_idx)
        d = stan_input(split_result)
        split_pw = Vector{Loo.PointwiseLoo}(undef, d["N"])
        sf = stanfit(split_result)
        p = extract(sf)
        draws = sf.iter * sf.chains

        @threads for i in 1:d["N"]
            ll = pointwise_loglikelihoods(sf, p, i)
            split_pw[i] = Loo.pointwise_loo(ll)
        end

        append!(pw, split_pw)
    end

    return Loo.LooResult(pw, (draws, length(pw)))
end

@recipe function f(
    ::Calibration_Plot, result::SplitHLAModelResult
)
    y = Vector{Int}()
    theta_pred = Vector{Float64}()

    for split_idx in 1:result.max_idx
        split_result = Base.invokelatest(result.f, split_idx)
        sf = stanfit(split_result)

        N = sf.data["N"]
        rs = sf.data["rs"]
        idx = sf.data["idx"]
        d = sf.data
        posterior = extract(sf)

        split_y = sf.data["y"]
        split_theta_pred = Vector{Float64}(undef, N)

        @threads for i in 1:N
            theta = mean(theta_i(sf, posterior, i))
            split_theta_pred[i] = theta
        end

        append!(y, split_y)
        append!(theta_pred, split_theta_pred)
    end

    @series begin
        seriestype := :path
        linecolor := "black"
        linestyle := :dash
        [0, 1], [0, 1]
    end

    @series begin
        seriestype := :calibration
        theta_pred, y
    end
end

@recipe function f(
    ::Phylogeny_Calibration, result::SplitHLAModelResult
)

    y = Vector{Int}()
    theta  = Vector{Float64}()

    for split_idx in 1:result.max_idx
        split_result = Base.invokelatest(result.f, split_idx)
        sf = split_result.sf
        N = sf.data["N"]
        rs = sf.data["rs"]
        idx = sf.data["idx"]

        split_theta = [sf.data["phy"][rs[i], idx[i]] for i in 1:N]
        split_y = sf.data["y"]

        append!(y, split_y)
        append!(theta, split_theta)
    end

    @series begin
        seriestype := :calibration
        theta, y
    end

    @series begin
        seriestype := :path
        linecolor := "black"
        linestyle := :dash
        [0, 1], [0, 1]
    end
end