text:
# https://github.com/JuliaLang/julia/pull/35973
if v"1.3.0" <= VERSION < v"1.6.0-DEV.84"
    # Declare a new `==` function to avoid type piracy and triggering invalidations
    function == end
    @inline ==(a, b) = Base.:(==)(a, b)

    function ==(a::Union{String, SubString{String}}, b::Union{String, SubString{String}})
        s = sizeof(a)
        s == sizeof(b) && 0 == Base._memcmp(a, b, s)
    end
end

if v"1.3.1-pre.18" <= VERSION < v"1.3.2"
    using Pkg.Artifacts: do_artifact_str, find_artifacts_toml, load_artifacts_toml

    # A copy of `Pkg.Artifacts.@artifact_str` where `name` is properly escaped on
    # Julia 1.3.1
    # https://github.com/JuliaLang/Pkg.jl/issues/1912
    # https://github.com/JuliaLang/Pkg.jl/pull/1580
    macro artifact_str(name)
        # Load Artifacts.toml at compile time, so that we don't have to use `__source__.file`
        # at runtime, which gets stale if the `.ji` file is relocated.
        local artifacts_toml = find_artifacts_toml(string(__source__.file))
        if artifacts_toml === nothing
            error(string(
                "Cannot locate '(Julia)Artifacts.toml' file when attempting to use artifact '",
                name,
                "' in '",
                __module__,
                "'",
            ))
        end

        local artifact_dict = load_artifacts_toml(artifacts_toml)
        return quote
            # Invalidate .ji file if Artifacts.toml file changes
            Base.include_dependency($(artifacts_toml))

            # Use invokelatest() to introduce a compiler barrier, preventing many backedges from being added
            # and slowing down not only compile time, but also `.ji` load time.  This is critical here, as
            # artifact"" is used in other modules, so we don't want to be spreading backedges around everywhere.
            Base.invokelatest(do_artifact_str, $(esc(name)), $(artifact_dict), $(artifacts_toml), $__module__)
        end
    end
elseif VERSION >= v"1.3.0"
    # Issue does not exist in 1.3.0, and >= 1.4.0. Also, assume that the issue would also be
    # fixed on >= 1.3.2
    using Pkg.Artifacts: @artifact_str
end


text:
struct Exponential <: Distribution{Float64} end

"""
    exponential(rate::Real)

Sample a `Float64` from the exponential distribution with rate parameter `rate`.
"""
const exponential = Exponential()

function Gen.logpdf(::Exponential, x::Real, rate::Real)
    scale = 1/rate
    Distributions.logpdf(Distributions.Exponential(scale), x)
end

function Gen.logpdf_grad(::Exponential, x::Real, rate::Real)
    scale = 1/rate
    x_grad = Distributions.gradlogpdf(Distributions.Exponential(scale), x)
    rate_grad = 1/rate - x
    (x_grad, rate_grad)
end

function Gen.random(::Exponential, rate::Real)
    scale = 1/rate
    rand(Distributions.Exponential(scale))
end

is_discrete(::Exponential) = false

(::Exponential)(rate) = random(Exponential(), rate)

Gen.has_output_grad(::Exponential) = true
Gen.has_argument_grads(::Exponential) = (true,)

export exponential

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## thread/task locking abstraction

"""
    AbstractLock

Abstract supertype describing types that
implement the synchronization primitives:
[`lock`](@ref), [`trylock`](@ref), [`unlock`](@ref), and [`islocked`](@ref).
"""
abstract type AbstractLock end
function lock end
function unlock end
function trylock end
function islocked end
unlockall(l::AbstractLock) = unlock(l) # internal function for implementing `wait`
relockall(l::AbstractLock, token::Nothing) = lock(l) # internal function for implementing `wait`
assert_havelock(l::AbstractLock) = assert_havelock(l, Threads.threadid())
assert_havelock(l::AbstractLock, tid::Integer) =
    (islocked(l) && tid == Threads.threadid()) ? nothing : error("concurrency violation detected")
assert_havelock(l::AbstractLock, tid::Task) =
    (islocked(l) && tid === current_task()) ? nothing : error("concurrency violation detected")
assert_havelock(l::AbstractLock, tid::Nothing) = error("concurrency violation detected")

"""
    AlwaysLockedST

This struct does not implement a real lock, but instead
pretends to be always locked on the original thread it was allocated on,
and simply ignores all other interactions.
It also does not synchronize tasks; for that use a real lock such as [`RecursiveLock`](@ref).
This can be used in the place of a real lock to, instead, simply and cheaply assert
that the operation is only occurring on a single cooperatively-scheduled thread.
It is thus functionally equivalent to allocating a real, recursive, task-unaware lock
immediately calling `lock` on it, and then never calling a matching `unlock`,
except that calling `lock` from another thread will throw a concurrency violation exception.
"""
struct AlwaysLockedST <: AbstractLock
    ownertid::Int16
    AlwaysLockedST() = new(Threads.threadid())
end
assert_havelock(l::AlwaysLockedST) = assert_havelock(l, l.ownertid)
lock(l::AlwaysLockedST) = assert_havelock(l)
unlock(l::AlwaysLockedST) = assert_havelock(l)
trylock(l::AlwaysLockedST) = l.ownertid == Threads.threadid()
islocked(::AlwaysLockedST) = true

## condition variables

"""
    GenericCondition

Abstract implementation of a condition object
for synchonizing tasks objects with a given lock.
"""
struct GenericCondition{L<:AbstractLock}
    waitq::InvasiveLinkedList{Task}
    lock::L

    GenericCondition{L}() where {L<:AbstractLock} = new{L}(InvasiveLinkedList{Task}(), L())
    GenericCondition{L}(l::L) where {L<:AbstractLock} = new{L}(InvasiveLinkedList{Task}(), l)
    GenericCondition(l::AbstractLock) = new{typeof(l)}(InvasiveLinkedList{Task}(), l)
end

assert_havelock(c::GenericCondition) = assert_havelock(c.lock)
lock(c::GenericCondition) = lock(c.lock)
unlock(c::GenericCondition) = unlock(c.lock)
trylock(c::GenericCondition) = trylock(c.lock)
islocked(c::GenericCondition) = islocked(c.lock)

"""
    wait([x])

Block the current task until some event occurs, depending on the type of the argument:

* [`Channel`](@ref): Wait for a value to be appended to the channel.
* [`Condition`](@ref): Wait for [`notify`](@ref) on a condition.
* `Process`: Wait for a process or process chain to exit. The `exitcode` field of a process
  can be used to determine success or failure.
* [`Task`](@ref): Wait for a `Task` to finish. If the task fails with an exception, the
  exception is propagated (re-thrown in the task that called `wait`).
* [`RawFD`](@ref): Wait for changes on a file descriptor (see the `FileWatching` package).

If no argument is passed, the task blocks for an undefined period. A task can only be
restarted by an explicit call to [`schedule`](@ref) or [`yieldto`](@ref).

Often `wait` is called within a `while` loop to ensure a waited-for condition is met before
proceeding.
"""
function wait(c::GenericCondition)
    ct = current_task()
    assert_havelock(c)
    push!(c.waitq, ct)
    token = unlockall(c.lock)
    try
        return wait()
    catch
        list_deletefirst!(c.waitq, ct)
        rethrow()
    finally
        relockall(c.lock, token)
    end
end

"""
    notify(condition, val=nothing; all=true, error=false)

Wake up tasks waiting for a condition, passing them `val`. If `all` is `true` (the default),
all waiting tasks are woken, otherwise only one is. If `error` is `true`, the passed value
is raised as an exception in the woken tasks.

Return the count of tasks woken up. Return 0 if no tasks are waiting on `condition`.
"""
notify(c::GenericCondition, @nospecialize(arg = nothing); all=true, error=false) = notify(c, arg, all, error)
function notify(c::GenericCondition, @nospecialize(arg), all, error)
    assert_havelock(c)
    cnt = 0
    while !isempty(c.waitq)
        t = popfirst!(c.waitq)
        schedule(t, arg, error=error)
        cnt += 1
        all || break
    end
    return cnt
end

notify_error(c::GenericCondition, err) = notify(c, err, true, true)

n_waiters(c::GenericCondition) = length(c.waitq)

"""
    isempty(condition)

Return `true` if no tasks are waiting on the condition, `false` otherwise.
"""
isempty(c::GenericCondition) = isempty(c.waitq)

# default (Julia v1.0) is currently single-threaded
# (although it uses MT-safe versions, when possible)
"""
    Condition()

Create an edge-triggered event source that tasks can wait for. Tasks that call [`wait`](@ref) on a
`Condition` are suspended and queued. Tasks are woken up when [`notify`](@ref) is later called on
the `Condition`. Edge triggering means that only tasks waiting at the time [`notify`](@ref) is
called can be woken up. For level-triggered notifications, you must keep extra state to keep
track of whether a notification has happened. The [`Channel`](@ref) and [`Event`](@ref) types do
this, and can be used for level-triggered events.

This object is NOT thread-safe. See [`Threads.Condition`](@ref) for a thread-safe version.
"""
const Condition = GenericCondition{AlwaysLockedST}

## async event notifications

"""
    AsyncCondition()

Create a async condition that wakes up tasks waiting for it
(by calling [`wait`](@ref) on the object)
when notified from C by a call to `uv_async_send`.
Waiting tasks are woken with an error when the object is closed (by [`close`](@ref).
Use [`isopen`](@ref) to check whether it is still active.
"""
mutable struct AsyncCondition
    handle::Ptr{Cvoid}
    cond::Condition
    isopen::Bool

    function AsyncCondition()
        this = new(Libc.malloc(_sizeof_uv_async), Condition(), true)
        associate_julia_struct(this.handle, this)
        finalizer(uvfinalize, this)
        err = ccall(:uv_async_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
            eventloop(), this, uv_jl_asynccb::Ptr{Cvoid})
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_async_init", err))
        end
        return this
    end
end

"""
    AsyncCondition(callback::Function)

Create a async condition that calls the given `callback` function. The `callback` is passed one argument,
the async condition object itself.
"""
function AsyncCondition(cb::Function)
    async = AsyncCondition()
    waiter = Task(function()
        while isopen(async)
            success = try
                wait(async)
                true
            catch exc # ignore possible exception on close()
                isa(exc, EOFError) || rethrow()
            end
            success && cb(async)
        end
    end)
    # must start the task right away so that it can wait for the AsyncCondition before
    # we re-enter the event loop. this avoids a race condition. see issue #12719
    yield(waiter)
    return async
end

## timer-based notifications

"""
    Timer(delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object).

Waiting tasks are woken after an initial delay of `delay` seconds, and then repeating with the given
`interval` in seconds. If `interval` is equal to `0`, the timer is only triggered once. When
the timer is closed (by [`close`](@ref) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.
"""
mutable struct Timer
    handle::Ptr{Cvoid}
    cond::Condition
    isopen::Bool

    function Timer(timeout::Real; interval::Real = 0.0)
        timeout ≥ 0 || throw(ArgumentError("timer cannot have negative timeout of $timeout seconds"))
        interval ≥ 0 || throw(ArgumentError("timer cannot have negative repeat interval of $interval seconds"))

        this = new(Libc.malloc(_sizeof_uv_timer), Condition(), true)
        err = ccall(:uv_timer_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), eventloop(), this)
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_timer_init", err))
        end

        associate_julia_struct(this.handle, this)
        finalizer(uvfinalize, this)

        ccall(:uv_update_time, Cvoid, (Ptr{Cvoid},), eventloop())
        ccall(:uv_timer_start,  Cint,  (Ptr{Cvoid}, Ptr{Cvoid}, UInt64, UInt64),
              this, uv_jl_timercb::Ptr{Cvoid},
              UInt64(round(timeout * 1000)) + 1, UInt64(round(interval * 1000)))
        return this
    end
end

unsafe_convert(::Type{Ptr{Cvoid}}, t::Timer) = t.handle
unsafe_convert(::Type{Ptr{Cvoid}}, async::AsyncCondition) = async.handle

function wait(t::Union{Timer, AsyncCondition})
    isopen(t) || throw(EOFError())
    stream_wait(t, t.cond)
end

isopen(t::Union{Timer, AsyncCondition}) = t.isopen

function close(t::Union{Timer, AsyncCondition})
    if t.handle != C_NULL && isopen(t)
        t.isopen = false
        isa(t, Timer) && ccall(:uv_timer_stop, Cint, (Ptr{Cvoid},), t)
        ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), t)
    end
    nothing
end

function uvfinalize(t::Union{Timer, AsyncCondition})
    if t.handle != C_NULL
        disassociate_julia_struct(t.handle) # not going to call the usual close hooks
        close(t)
        t.handle = C_NULL
    end
    t.isopen = false
    nothing
end

function _uv_hook_close(t::Union{Timer, AsyncCondition})
    uvfinalize(t)
    notify_error(t.cond, EOFError())
    nothing
end

function uv_asynccb(handle::Ptr{Cvoid})
    async = @handle_as handle AsyncCondition
    notify(async.cond)
    nothing
end

function uv_timercb(handle::Ptr{Cvoid})
    t = @handle_as handle Timer
    if ccall(:uv_timer_get_repeat, UInt64, (Ptr{Cvoid},), t) == 0
        # timer is stopped now
        close(t)
    end
    notify(t.cond)
    nothing
end

"""
    sleep(seconds)

Block the current task for a specified number of seconds. The minimum sleep time is 1
millisecond or input of `0.001`.
"""
function sleep(sec::Real)
    sec ≥ 0 || throw(ArgumentError("cannot sleep for $sec seconds"))
    wait(Timer(sec))
    nothing
end

# timer with repeated callback
"""
    Timer(callback::Function, delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object) and
calls the function `callback`.

Waiting tasks are woken and the function `callback` is called after an initial delay of `delay` seconds,
and then repeating with the given `interval` in seconds. If `interval` is equal to `0`, the timer
is only triggered once. The function `callback` is called with a single argument, the timer itself.
When the timer is closed (by [`close`](@ref) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.

# Examples

Here the first number is printed after a delay of two seconds, then the following numbers are printed quickly.

```julia-repl
julia> begin
           i = 0
           cb(timer) = (global i += 1; println(i))
           t = Timer(cb, 2, interval = 0.2)
           wait(t)
           sleep(0.5)
           close(t)
       end
1
2
3
```
"""
function Timer(cb::Function, timeout::Real; interval::Real = 0.0)
    t = Timer(timeout, interval = interval)
    waiter = Task(function()
        while isopen(t)
            success = try
                wait(t)
                true
            catch exc # ignore possible exception on close()
                isa(exc, EOFError) || rethrow()
                false
            end
            success && cb(t)
        end
    end)
    # must start the task right away so that it can wait for the Timer before
    # we re-enter the event loop. this avoids a race condition. see issue #12719
    yield(waiter)
    return t
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Source <: StaticInjection
        name::String
        available::Bool
        bus::Bus
        active_power::Float64
        reactive_power::Float64
        R_th::Float64
        X_th::Float64
        internal_voltage::Float64
        internal_angle::Float64
        dynamic_injector::Union{Nothing, DynamicInjection}
        services::Vector{Service}
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

This struct acts as an infinity bus.

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `active_power::Float64`
- `reactive_power::Float64`
- `R_th::Float64`: Source Thevenin resistance, validation range: `(0, nothing)`
- `X_th::Float64`: Source Thevenin reactance, validation range: `(0, nothing)`
- `internal_voltage::Float64`: Internal Voltage, validation range: `(0, nothing)`
- `internal_angle::Float64`: Internal Angle
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Source <: StaticInjection
    name::String
    available::Bool
    bus::Bus
    active_power::Float64
    reactive_power::Float64
    "Source Thevenin resistance"
    R_th::Float64
    "Source Thevenin reactance"
    X_th::Float64
    "Internal Voltage"
    internal_voltage::Float64
    "Internal Angle"
    internal_angle::Float64
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Source(name, available, bus, active_power, reactive_power, R_th, X_th, internal_voltage=1.0, internal_angle=0.0, dynamic_injector=nothing, services=Device[], ext=Dict{String, Any}(), )
    Source(name, available, bus, active_power, reactive_power, R_th, X_th, internal_voltage, internal_angle, dynamic_injector, services, ext, InfrastructureSystemsInternal(), )
end

function Source(; name, available, bus, active_power, reactive_power, R_th, X_th, internal_voltage=1.0, internal_angle=0.0, dynamic_injector=nothing, services=Device[], ext=Dict{String, Any}(), internal=InfrastructureSystemsInternal(), )
    Source(name, available, bus, active_power, reactive_power, R_th, X_th, internal_voltage, internal_angle, dynamic_injector, services, ext, internal, )
end

# Constructor for demo purposes; non-functional.
function Source(::Nothing)
    Source(;
        name="init",
        available=false,
        bus=Bus(nothing),
        active_power=0.0,
        reactive_power=0.0,
        R_th=0,
        X_th=0,
        internal_voltage=0,
        internal_angle=0,
        dynamic_injector=nothing,
        services=Device[],
        ext=Dict{String, Any}(),
    )
end

"""Get [`Source`](@ref) `name`."""
get_name(value::Source) = value.name
"""Get [`Source`](@ref) `available`."""
get_available(value::Source) = value.available
"""Get [`Source`](@ref) `bus`."""
get_bus(value::Source) = value.bus
"""Get [`Source`](@ref) `active_power`."""
get_active_power(value::Source) = value.active_power
"""Get [`Source`](@ref) `reactive_power`."""
get_reactive_power(value::Source) = value.reactive_power
"""Get [`Source`](@ref) `R_th`."""
get_R_th(value::Source) = value.R_th
"""Get [`Source`](@ref) `X_th`."""
get_X_th(value::Source) = value.X_th
"""Get [`Source`](@ref) `internal_voltage`."""
get_internal_voltage(value::Source) = value.internal_voltage
"""Get [`Source`](@ref) `internal_angle`."""
get_internal_angle(value::Source) = value.internal_angle
"""Get [`Source`](@ref) `dynamic_injector`."""
get_dynamic_injector(value::Source) = value.dynamic_injector
"""Get [`Source`](@ref) `services`."""
get_services(value::Source) = value.services
"""Get [`Source`](@ref) `ext`."""
get_ext(value::Source) = value.ext
"""Get [`Source`](@ref) `internal`."""
get_internal(value::Source) = value.internal

"""Set [`Source`](@ref) `name`."""
set_name!(value::Source, val) = value.name = val
"""Set [`Source`](@ref) `available`."""
set_available!(value::Source, val) = value.available = val
"""Set [`Source`](@ref) `bus`."""
set_bus!(value::Source, val) = value.bus = val
"""Set [`Source`](@ref) `active_power`."""
set_active_power!(value::Source, val) = value.active_power = val
"""Set [`Source`](@ref) `reactive_power`."""
set_reactive_power!(value::Source, val) = value.reactive_power = val
"""Set [`Source`](@ref) `R_th`."""
set_R_th!(value::Source, val) = value.R_th = val
"""Set [`Source`](@ref) `X_th`."""
set_X_th!(value::Source, val) = value.X_th = val
"""Set [`Source`](@ref) `internal_voltage`."""
set_internal_voltage!(value::Source, val) = value.internal_voltage = val
"""Set [`Source`](@ref) `internal_angle`."""
set_internal_angle!(value::Source, val) = value.internal_angle = val
"""Set [`Source`](@ref) `services`."""
set_services!(value::Source, val) = value.services = val
"""Set [`Source`](@ref) `ext`."""
set_ext!(value::Source, val) = value.ext = val


text:
"""
    _A_gain_method(parameters::Parameters1PL, latent::Latent1D)
"""
function _A_gain_method(parameters::Parameters1PL, latent::Latent1D)
    return _item_expected_information(parameters, latent)
end

"""
    _A_gain_method(parameters::Parameters2PL, latent::Latent1D)
"""
function _A_gain_method(parameters::Parameters2PL, latent::Latent1D)
    old_exp_info = copy(parameters.expected_information)
    return LinearAlgebra.tr(old_exp_info + _item_expected_information(parameters, latent)) -
           LinearAlgebra.tr(old_exp_info)
end

"""
    A_gain_method(item::AbstractItem, examinee::AbstractExaminee)

# Description
Computes the gain in the trace of the expected information matrix for an item.

# Arguments
- **`item::AbstractItem`**: The  item.
- **`examinee::AbstractExaminee`**: The examinee at which computing the information.

# Output
It returns a `Float64` scalar.
"""
function A_gain_method(item::AbstractItem, examinee::AbstractExaminee)
    return _A_gain_method(item.parameters, examinee.latent)
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

x_coords = 0.0:1.0:6.0
#[0.0, 6.0, 6.0, 12.0, 12.0, 14.0],
y_coords = repeat(0.0:0.0, outer=7)
#0.0, -4.0, 0.0, -5.0, 0.0],
x_coords = vcat(x_coords, 6.0)
y_coords = vcat(y_coords, -4.0)
x_coords = vcat(x_coords, collect(7.0:1.0:12.0))
y_coords = vcat(y_coords, repeat(0.0:0.0, outer=6))
x_coords = vcat(x_coords, 12.0)
y_coords = vcat(y_coords, -5.0)
x_coords = vcat(x_coords, 13.0:1.0:14.0)
y_coords = vcat(y_coords, repeat(0.0:0.0, outer=2))

xy_coords = hcat(x_coords, y_coords)

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(17, 17, 2, 1, 1, Line(2, 3)),
  :properties => [
      5.0e9 6.0e4;
      1.0e9  2.0e4
    ],
  :etype => vcat(repeat(1:1, outer=14), repeat(2:2, outer=3)),
  :x_coords => x_coords,
  :y_coords => y_coords,
  :g_num => [
      1 2 3 4 5 6 7  9 10 11 12 13 14 16 8  8 15;
      2 3 4 5 6 7 9 10 11 12 13 14 16 17 7 14 14
    ],
  :support => [
      ( 1, [0 0 1]),
      ( 8, [0 0 0]),
      (15, [0 0 0])
    ],
  :loaded_nodes => [
      ( 1, [0.0 -10.0 -20//12]),
      ( 2, [0.0 -20.0     0.0]),
      ( 3, [0.0 -20.0     0.0]),
      ( 4, [0.0 -20.0     0.0]),
      ( 5, [0.0 -20.0     0.0]),
      ( 6, [0.0 -20.0     0.0]),
      ( 7, [0.0 -20.0     0.0]),
      ( 9, [0.0 -20.0     0.0]),
      (10, [0.0 -80.0     0.0]),
      (11, [0.0 -20.0     0.0]),
      (12, [0.0 -80.0     0.0]),
      (13, [0.0 -20.0     0.0]),
      (14, [0.0 -20.0     0.0]),
      (16, [0.0 -20.0     0.0]),
      (17, [0.0 -10.0  20//12])
    ],
  :penalty => 1e19,
  :eq_nodal_forces_and_moments => [
    ( 1, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 2, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 3, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 4, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 5, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 6, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 7, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    ( 9, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (10, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (11, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (12, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (13, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (14, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (16, [0.0 -10.0 -20//12 0.0 -10.0 20//12]),
    (17, [0.0 -10.0 -20//12 0.0 -10.0 20//12])
  ]
)

data |> display
println()

@time fem, dis_df, fm_df = p44(data)
println()

display(dis_df)
println()
display(fm_df)

if VERSION.minor < 6
  using Plots
  gr(size=(400,600))

  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
  titles = ["p4.4.1 translations", "p4.4.1 rotations",
    "p4.4.1 y shear force", "p4.4.1 z moment"]
  indxs = [1,2,3,4,5,6,7,9,10,11,12,13,14,16,17]
  moms = vcat(
    convert(Array, fm_df[:, :z1_Moment]),
    convert(Array, fm_df[:, :z2_Moment])[end]
  )
  fors = vcat(
    convert(Array, fm_df[:, :y1_Force]),
    convert(Array, fm_df[:, :y2_Force])[end]
  )
  x_coords = data[:x_coords]

  indxs = [1,2,3,4,5,6,7,9,10,11,12,13,14,16,17]
  p[1] = plot(fem.displacements[2,indxs], ylim=(-0.015, 0.015),
    xlabel="node", ylabel="y translation [m]", color=:black,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[1], leg=false)
  p[2] = plot(fem.displacements[3,indxs], ylim=(-0.008, 0.008),
    xlabel="node", ylabel="rotation [radians]", color=:red,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[2], leg=false)
  p[3] = plot(fors, lab="y Shear force", ylim=(-250.0, 250.0),
    xlabel="node", ylabel="shear force [N]", color=:blue,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[3], leg=false)
  p[4] = plot(moms, lab="z Moment", ylim=(-150.0, 150.0),
    xlabel="node", ylabel="z moment [Nm]", color=:green,
    line=(:dash,1), marker=(:circle,4,0.8,stroke(1,:black)),
    title=titles[4], leg=false)

  plot(p..., layout=(4, 1))
  savefig(ProjDir*"/Ex44.1a.png")

end


text:
abstract type AbstractResponseModel end

"""
    random_betas(gr::GroupedFeatures, αs)

Suppose `gr` consists of ``K`` groups with ``p_1, \\dotsc, p_K`` features each.
Then this returns a random vector of βs of length ``\\sum p_g``,
where for `j` in the `g`-th group
we draw (independent) ``β_j \\sim N(0, α_g^2/p_g)``.
``α_g`` is the `g`-th element of the vectpr `αs`.
"""
function random_betas(gr::GroupedFeatures, αs)
    ps = gr.ps
    βs = zeros(eltype(αs), gr.p)
    for i = 1:gr.num_groups
        βs[group_idx(gr, i)] .= randn(ps[i]) .* sqrt(αs[i]^2 / ps[i])
    end
    βs
end

Base.@kwdef struct RandomLinearResponseModel <: AbstractResponseModel
    αs::Vector{Float64}
    grp::GroupedFeatures
    iid_measure = Normal()
end

function (resp::RandomLinearResponseModel)(X)
    β = random_betas(resp.grp, resp.αs) #todo, allow other noise dbn.
    Xβ = X * β
    Xβ, β
end

Base.@kwdef struct GroupRidgeSimulationSettings{C,R,D}
    groups::GroupedFeatures
    Σ::C
    response_model::R
    response_noise::D = Normal()
    ntest::Int = 10000
    ntrain::Int
    iid_measure = Normal()
end

Base.@kwdef struct GroupRidgeSimulation
    groups::GroupedFeatures
    X::Matrix{Float64}
    Y::Vector{Float64}
    resampling_idx = nothing
    β = nothing
end

function simulate(group_simulation::GroupRidgeSimulationSettings)
    ntrain = group_simulation.ntrain
    ntest = group_simulation.ntest
    ntotal = ntrain + ntest
    @unpack response_model, response_noise = group_simulation

    X = simulate_rotated_design(
        group_simulation.Σ,
        ntotal;
        rotated_measure = group_simulation.iid_measure,
    )

    Xβ, β = response_model(X)
    Y = Xβ .+ rand(response_noise, ntotal)

    resampling_idx = [(1:ntrain, (ntrain+1):ntotal)]

    GroupRidgeSimulation(;
        groups = group_simulation.groups,
        X=X,
        Y=Y,
        resampling_idx = resampling_idx,
        β = β,
    )
end


text:
using CSV
CSV.File(joinpath(@__DIR__,"sample.csv"))

max_stars_repo_name: jondea/juliab

max_issues_repo_name: jondea/juliab

max_forks_repo_name: jondea/juliab


text:
module Optimisers

    # For now we just use the Knet optimisers and export them for uniform API
    export Adam, SGD, Momentum, Nesterov, Adagrad, Rmsprop, Adadelta

    using Knet: Adam, SGD, Momentum, Nesterov, Adagrad, Rmsprop, Adadelta

end


text:
__precompile__()

module SearchMethods
	# package code goes here
	export bfs, dfs, dls, ids, ucs
	export greedy, a_star

	using DataStructures

	include("uninformed.jl")
	include("informed.jl")

end # module


text:
# Standardization for data arrays containing NAs. The NAs are simply ignored
# when analyzing the population.

function na_mean_and_stddev{T<:Real}(X::AbstractArray{T}, sample::Bool)
  sum = 0.
  n = 0
  for x in X
    if !isna(x)
      n += 1
      sum += x
    end
  end

  # TODO handle n=0
  µ = sum / n

  sum = 0.
  for x in X
    if !isna(x)
      sum += (x - µ)^2
    end
  end

  # Use n - 1 to calculate sample stddev
  if sample
    n -= 1
  end

  # TODO handle n=0
  # TODO could multiply by cached 1/n if available
  σ = √(sum / n)

  (µ, σ)
end

function nazscore{T<:Real}(X::AbstractArray{T}, µ::Real, σ::Real)
  Z = DataArray(Float64, size(X))
  for i = 1:length(X)
    x = X[i]
    if !isna(x)
      @inbounds Z[i] = (x - µ) / σ
    else
      @inbounds Z[i] = x
    end
  end
  Z
end
nazscore{T<:Real}(X::AbstractArray{T}) = ((µ, σ) = na_mean_and_stddev(X, true); nazscore(X, µ, σ))
nazscore_pop{T<:Real}(X::AbstractArray{T}) = ((µ, σ) = na_mean_and_stddev(X, false); nazscore(X, µ, σ))

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


text:
using Test, Dates, Temporal

@testset "Conversion" begin
    data = (time=[DateTime(2018, 11, 21, 12, 0), DateTime(2018, 11, 21, 13, 0)], col1=[10.2, 11.2], col2=[20.2, 21.2])
    result = TS(data)
    for i in 1:size(result,1)
        @test result.index[i] == data[:time][i]
        for j in 1:size(result,2)
            @test result.values[i,j] == data[j+1][i]
        end
    end
end


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

name = "DASKR"
version = v"1.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.netlib.org/ode/daskr.tgz", "5472e98c954062ea69d19bd16d1cb8a8835497d6b0df18fe5753bae61d151774")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd DASKR
install_license LICENSE
mkdir -p "${libdir}"
gfortran -shared -fPIC -o $libdir/libdaskr.${dlext} solver/d*.f
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    unreal(::T) where T <: Real

For real numbers, return an empty array of type `T`.
"""
function unreal(::T) where T <: Real
    T[]
end

"""
    iszerodivisor(a::Real) -> Bool

For real numbers, return `true` if `a` is zero. Equivalent to `iszero(a)`.
"""
function iszerodivisor(a::Real)
    iszero(a)
end

"""
    cloak(a::Real)

For real numbers, equivalent to the minus operation.
"""
function cloak(a::Real)
    -(a)
end

"""
    dagger(a::Real)

For real numbers, equivalent to the identity operation.
"""
function dagger(a::Real)
    a
end

"""
    star(a::Real)

For real numbers, the identity operation.
"""
function star(a::Real)
    a
end

"""
    quadrance(a::Real)

For real numbers, the squared absolute value.
"""
function quadrance(a::Real)
    abs2(a)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Event is a report of an event somewhere in the cluster. It generally denotes some state change in the system.

    IoK8sApiEventsV1beta1Event(;
        action=nothing,
        apiVersion=nothing,
        deprecatedCount=nothing,
        deprecatedFirstTimestamp=nothing,
        deprecatedLastTimestamp=nothing,
        deprecatedSource=nothing,
        eventTime=nothing,
        kind=nothing,
        metadata=nothing,
        note=nothing,
        reason=nothing,
        regarding=nothing,
        related=nothing,
        reportingController=nothing,
        reportingInstance=nothing,
        series=nothing,
        type=nothing,
    )

    - action::String : What action was taken/failed regarding to the regarding object.
    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - deprecatedCount::Int32 : Deprecated field assuring backward compatibility with core.v1 Event type
    - deprecatedFirstTimestamp::IoK8sApimachineryPkgApisMetaV1Time : Deprecated field assuring backward compatibility with core.v1 Event type
    - deprecatedLastTimestamp::IoK8sApimachineryPkgApisMetaV1Time : Deprecated field assuring backward compatibility with core.v1 Event type
    - deprecatedSource::IoK8sApiCoreV1EventSource : Deprecated field assuring backward compatibility with core.v1 Event type
    - eventTime::IoK8sApimachineryPkgApisMetaV1MicroTime : Required. Time when this Event was first observed.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta
    - note::String : Optional. A human-readable description of the status of this operation. Maximal length of the note is 1kB, but libraries should be prepared to handle values up to 64kB.
    - reason::String : Why the action was taken.
    - regarding::IoK8sApiCoreV1ObjectReference : The object this Event is about. In most cases it's an Object reporting controller implements. E.g. ReplicaSetController implements ReplicaSets and this event is emitted because it acts on some changes in a ReplicaSet object.
    - related::IoK8sApiCoreV1ObjectReference : Optional secondary object for more complex actions. E.g. when regarding object triggers a creation or deletion of related object.
    - reportingController::String : Name of the controller that emitted this Event, e.g. `kubernetes.io/kubelet`.
    - reportingInstance::String : ID of the controller instance, e.g. `kubelet-xyzf`.
    - series::IoK8sApiEventsV1beta1EventSeries : Data about the Event series this event represents or nil if it's a singleton Event.
    - type::String : Type of this event (Normal, Warning), new types could be added in the future.
"""
mutable struct IoK8sApiEventsV1beta1Event <: SwaggerModel
    action::Any # spec type: Union{ Nothing, String } # spec name: action
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    deprecatedCount::Any # spec type: Union{ Nothing, Int32 } # spec name: deprecatedCount
    deprecatedFirstTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: deprecatedFirstTimestamp
    deprecatedLastTimestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: deprecatedLastTimestamp
    deprecatedSource::Any # spec type: Union{ Nothing, IoK8sApiCoreV1EventSource } # spec name: deprecatedSource
    eventTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1MicroTime } # spec name: eventTime
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    note::Any # spec type: Union{ Nothing, String } # spec name: note
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    regarding::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: regarding
    related::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: related
    reportingController::Any # spec type: Union{ Nothing, String } # spec name: reportingController
    reportingInstance::Any # spec type: Union{ Nothing, String } # spec name: reportingInstance
    series::Any # spec type: Union{ Nothing, IoK8sApiEventsV1beta1EventSeries } # spec name: series
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiEventsV1beta1Event(;action=nothing, apiVersion=nothing, deprecatedCount=nothing, deprecatedFirstTimestamp=nothing, deprecatedLastTimestamp=nothing, deprecatedSource=nothing, eventTime=nothing, kind=nothing, metadata=nothing, note=nothing, reason=nothing, regarding=nothing, related=nothing, reportingController=nothing, reportingInstance=nothing, series=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("action"), action)
        setfield!(o, Symbol("action"), action)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("deprecatedCount"), deprecatedCount)
        setfield!(o, Symbol("deprecatedCount"), deprecatedCount)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("deprecatedFirstTimestamp"), deprecatedFirstTimestamp)
        setfield!(o, Symbol("deprecatedFirstTimestamp"), deprecatedFirstTimestamp)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("deprecatedLastTimestamp"), deprecatedLastTimestamp)
        setfield!(o, Symbol("deprecatedLastTimestamp"), deprecatedLastTimestamp)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("deprecatedSource"), deprecatedSource)
        setfield!(o, Symbol("deprecatedSource"), deprecatedSource)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("eventTime"), eventTime)
        setfield!(o, Symbol("eventTime"), eventTime)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("note"), note)
        setfield!(o, Symbol("note"), note)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("regarding"), regarding)
        setfield!(o, Symbol("regarding"), regarding)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("related"), related)
        setfield!(o, Symbol("related"), related)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("reportingController"), reportingController)
        setfield!(o, Symbol("reportingController"), reportingController)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("reportingInstance"), reportingInstance)
        setfield!(o, Symbol("reportingInstance"), reportingInstance)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("series"), series)
        setfield!(o, Symbol("series"), series)
        validate_property(IoK8sApiEventsV1beta1Event, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiEventsV1beta1Event

const _property_map_IoK8sApiEventsV1beta1Event = Dict{Symbol,Symbol}(Symbol("action")=>Symbol("action"), Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("deprecatedCount")=>Symbol("deprecatedCount"), Symbol("deprecatedFirstTimestamp")=>Symbol("deprecatedFirstTimestamp"), Symbol("deprecatedLastTimestamp")=>Symbol("deprecatedLastTimestamp"), Symbol("deprecatedSource")=>Symbol("deprecatedSource"), Symbol("eventTime")=>Symbol("eventTime"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("note")=>Symbol("note"), Symbol("reason")=>Symbol("reason"), Symbol("regarding")=>Symbol("regarding"), Symbol("related")=>Symbol("related"), Symbol("reportingController")=>Symbol("reportingController"), Symbol("reportingInstance")=>Symbol("reportingInstance"), Symbol("series")=>Symbol("series"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiEventsV1beta1Event = Dict{Symbol,String}(Symbol("action")=>"String", Symbol("apiVersion")=>"String", Symbol("deprecatedCount")=>"Int32", Symbol("deprecatedFirstTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("deprecatedLastTimestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("deprecatedSource")=>"IoK8sApiCoreV1EventSource", Symbol("eventTime")=>"IoK8sApimachineryPkgApisMetaV1MicroTime", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("note")=>"String", Symbol("reason")=>"String", Symbol("regarding")=>"IoK8sApiCoreV1ObjectReference", Symbol("related")=>"IoK8sApiCoreV1ObjectReference", Symbol("reportingController")=>"String", Symbol("reportingInstance")=>"String", Symbol("series")=>"IoK8sApiEventsV1beta1EventSeries", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiEventsV1beta1Event }) = collect(keys(_property_map_IoK8sApiEventsV1beta1Event))
Swagger.property_type(::Type{ IoK8sApiEventsV1beta1Event }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiEventsV1beta1Event[name]))}
Swagger.field_name(::Type{ IoK8sApiEventsV1beta1Event }, property_name::Symbol) =  _property_map_IoK8sApiEventsV1beta1Event[property_name]

function check_required(o::IoK8sApiEventsV1beta1Event)
    (getproperty(o, Symbol("eventTime")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiEventsV1beta1Event }, name::Symbol, val)
end


text:
alpha = 0.85
accuracy = 0.00001
edgelist = "test/fixtures/dangling.adj"

rank = Rank{Float64,Int64}(alpha, accuracy, edgelist)
r = stationary_distribution(rank)
count = length(r)

# expected rank is  4, 3, 2, 5, 6, 1,
expected = [0.037375, 0.208173, 0.245344, 0.245344, 0.176404, 0.087356]
precision = 5

@assert isequal(floor(expected, precision), floor(r, precision))


text:
# Testing continuous univariate distributions

using Distributions
using Base.Test

### load reference data
#
#   Note
#   -------
#   To generate the reference data:
#   (1) make sure that python, numpy, and scipy are installed in your system
#   (2) enter the sub-directory test
#   (3) run: python discrete_ref.py > discrete_ref.csv
#
#   For most cases, you don't have. You only need to run this when you
#   implement a new distribution and want to add new test cases, then
#   you should add the new test cases to discrete_ref.py and run this
#   procedure to update the reference data.
#

immutable ContinuousRefEntry
    distr::ContinuousUnivariateDistribution
    mean::Float64
    var::Float64
    entropy::Float64
    x25::Float64
    x50::Float64
    x75::Float64
    lp25::Float64
    lp50::Float64
    lp75::Float64
end

function ContinuousRefEntry(row::Vector)
    @assert length(row) == 10
    d = eval(parse(row[1]))
    return ContinuousRefEntry(d, row[2:10]...)
end

csvpath = joinpath(dirname(@__FILE__), "continuous_ref.csv")
table = readcsv(csvpath)

R = [ContinuousRefEntry(vec(table[i,:])) for i = 2:size(table,1)]

### check with references

function verify(e::ContinuousRefEntry)
    d = e.distr

    if isfinite(e.mean)
        @test_approx_eq_eps mean(d)    e.mean    1.0e-12
    end

    if isfinite(e.var)
        @test_approx_eq_eps var(d)     e.var     1.0e-12
    end

    if isfinite(e.entropy) && applicable(entropy, d)
        @test_approx_eq_eps entropy(d) e.entropy 1.0e-6 * (abs(e.entropy) + 1.0)
    end

    @test_approx_eq_eps quantile(d, 0.25) e.x25 5.0e-9
    @test_approx_eq_eps quantile(d, 0.50) e.x50 5.0e-9
    @test_approx_eq_eps quantile(d, 0.75) e.x75 5.0e-9

    @test_approx_eq_eps logpdf(d, e.x25) e.lp25 1.0e-12
    @test_approx_eq_eps logpdf(d, e.x50) e.lp50 1.0e-12
    @test_approx_eq_eps logpdf(d, e.x75) e.lp75 1.0e-12
end

n_tsamples = 10^6

for rentry in R
    println("    testing $(rentry.distr)")
    verify(rentry)
    test_distr(rentry.distr, n_tsamples)
end

# additional distributions that have no direct counterparts in scipy
println("    -----")

for distr in [
    Frechet(0.5, 1.0),
    Frechet(3.0, 1.0),
    Frechet(20.0, 1.0),
    Frechet(120.0, 1.0),
    Frechet(0.5, 2.0),
    Frechet(3.0, 2.0),
    InverseGaussian(1.0, 1.0),
    InverseGaussian(2.0, 7.0),
    Levy(0.0, 1.0),
    Levy(2.0, 8.0),
    Levy(3.0, 3.0),
    LogNormal(0.0, 1.0),
    LogNormal(0.0, 2.0),
    LogNormal(3.0, 0.5),
    LogNormal(3.0, 1.0),
    LogNormal(3.0, 2.0),
    NoncentralBeta(2,2,0),
    NoncentralBeta(2,6,5),
    NoncentralChisq(2,2),
    NoncentralChisq(2,5),
    NoncentralF(2,2,2),
    NoncentralF(8,10,5),
    NoncentralT(2,2),
    NoncentralT(10,2) ]

    println("    testing $(distr)")
    test_distr(distr, n_tsamples)
end

for distr in [
    VonMises(0.0, 1.0),
    VonMises(0.5, 1.0),
    VonMises(0.5, 2.0) ]

    println("    testing $(distr)")
    test_samples(distr, n_tsamples)
end


text:
module TestIntegration

import JSON
using BenchmarkCI
using BenchmarkCI: mktempdir, versionof
using Dates: DateTime
using Test

function check_workspace(workspace = BenchmarkCI.DEFAULT_WORKSPACE)
    @test isfile(joinpath(workspace, "Project.toml"))
    @test isfile(joinpath(workspace, "Manifest.toml"))
    metadata = JSON.parsefile(joinpath(workspace, "metadata.json"))
    run_info = JSON.parsefile(joinpath(workspace, "run_info.json"))
    @test metadata["BenchmarkCI"]["versions"]["BenchmarkCI"] ==
          string(versionof(BenchmarkCI))
    @test metadata["BenchmarkCI"]["format_version"]::Int < 0
    @test VersionNumber(metadata["target_julia_info"]["VERSION"]) == VERSION
    @test metadata["target_julia_info"]["Sys"]["WORD_SIZE"] == Sys.WORD_SIZE
    @test metadata["cpu_info"]["cpubrand"] isa AbstractString
    @testset for side in ["target", "baseline"]
        d = run_info[side]
        @test DateTime(d["start_date"]) isa Any
        @test DateTime(d["end_date"]) isa Any
        @test d["elapsed_time"]::Real > 0
        @test length(d["git_tree_sha1"]::AbstractString) == 40
        @test length(d["git_commit_sha1"]::AbstractString) == 40
        @test length(d["git_tree_sha1_benchmark"]::AbstractString) == 40
        @test length(d["git_tree_sha1_src"]::AbstractString) == 40
    end
end

@testset "BenchmarkCI.jl" begin
    function flushall()
        flush(stderr)
        flush(stdout)
    end

    function printlns(n)
        flushall()
        for _ in 1:n
            println()
        end
        flushall()
    end

    mktempdir(prefix = "BenchmarkCI_jl_test_") do dir
        cd(dir) do
            run(`git clone https://github.com/tkf/BenchmarkCIExample.jl`)
            cd("BenchmarkCIExample.jl")

            # Run a test without $GITHUB_TOKEN
            function runtests(target)
                printlns(2)
                @info "Testing with target = $target"
                flushall()

                @testset "$target" begin
                    run(`git checkout $target`)
                    run(`git clean --force -xd`)
                    withenv("CI" => "true", "GITHUB_EVENT_PATH" => nothing) do
                        @testset "default" begin
                            BenchmarkCI.runall()
                            check_workspace()
                        end
                        @testset "project = benchmark/Project.toml" begin
                            BenchmarkCI.runall(project = "benchmark/Project.toml")
                            check_workspace()
                        end
                    end
                end
            end

            runtests("testcase/0000-with-manifest")
            runtests("testcase/0001-without-manifest")
            printlns(2)

            err = nothing
            @test try
                BenchmarkCI.judge(script = joinpath(dir, "nonexisting", "script.jl"))
                false
            catch err
                true
            end
            @test occursin("One of the following files must exist:", sprint(showerror, err))

            ciresult = BenchmarkCI._loadciresult()

            io = IOBuffer()
            BenchmarkCI.printcommentjson(io, ciresult)
            seekstart(io)
            dict = JSON.parse(io)
            @test dict["body"] isa String

            @testset "kwargs pass-through" begin
                @info "Testing with `BenchmarkCI.judge(retune = true)`"
                @test BenchmarkCI.judge(retune = true) isa Any
                printlns(2)
            end
        end
    end
end

end  # module


text:
include("gumbel_softmax.jl")

using ReverseDiff

const OP = [+, -, *, /]

X = collect(1:5)
gt1(x) = x + 1 #(1,1,2)
gt2(x) = x - 2  #(2,1,3)
gt3(x) = x * 3  #(3,1,4)
gt4(x) = x / 3  #(4,1,4)

gt = gt1

function f1(ps1, ps2, ps3, x)
    T = 0.3
    logits = log.(ps1)
    w = gumbel_softmax(logits, T, true)
    return w' * [op.(f2(ps2,x), f3(ps3,x)) for op in OP]
end
function f2(ps2, x)
    T = 0.3
    logits = log.(ps2)
    w = gumbel_softmax(logits, T, true)
    return w' * [x, 1, 2, 3]
end
function f3(ps3, x)
    T = 0.3
    logits = log.(ps3)
    w = gumbel_softmax(logits, T, true)
    return w' * [x, 1, 2, 3]
end

function objective(ps1, ps2, ps3, x, y_target)
    abs2.(y_target - f1(ps1, ps2, ps3, x))
end

using DataFrames, PGFPlots, TikzPictures
TikzPictures.standaloneWorkaround(true)
function learn(α=0.0002, N=5000; track=false,
               D=track ? DataFrame([Int, Int, Float64, Float64, Float64, Float64], [:iter, :node, :p1, :p2, :p3, :p4], 0) : 0)
    ps1 = ones(length(OP)) ./ length(OP)
    ps2 = ones(4) ./ 4
    ps3 = ones(4) ./ 4
    Y = [gt(x) for x in X]
    for i = 1:N
        @show i
        @show ps1
        @show ps2
        @show ps3
        if track
            push!(D, vcat(i, 1, ps1))
            push!(D, vcat(i, 2, ps2))
            push!(D, vcat(i, 3, ps3))
        end
        for j = 1:length(Y)
            g1, g2, g3, _, _ = ReverseDiff.gradient(objective, (ps1, ps2, ps3, [X[j]], [Y[j]]))
            #@show g
            if !any(isnan.(g1))
                ps1 -= α * g1
            end
            if !any(isnan.(g2))
                ps2 -= α * g2
            end
            if !any(isnan.(g3))
                ps3 -= α * g3
            end
            ps1 = max.(min.(ps1, 1.0), 1e-12)
            normalize!(ps1,1)
            ps2 = max.(min.(ps2, 1.0), 1e-12)
            normalize!(ps2,1)
            ps3 = max.(min.(ps3, 1.0), 1e-12)
            normalize!(ps3,1)
        end
        if isapprox(maximum(ps1),1) && isapprox(maximum(ps2),1) && isapprox(maximum(ps3),1)
            break
        end
    end
    return D, indmax(ps1), indmax(ps2), indmax(ps3)
end

function plot_me(D::DataFrame)
    g = GroupPlot(3, 1, groupStyle="horizontal sep=1.75cm, vertical sep=1.5cm")
    for i = 1:3
        ax = Axis([
                   Plots.Linear(D[D[:node].==i,:][:p1], legendentry="node$i p1"),
                   Plots.Linear(D[D[:node].==i,:][:p2], legendentry="node$i p2"),
                   Plots.Linear(D[D[:node].==i,:][:p3], legendentry="node$i p3"),
                   Plots.Linear(D[D[:node].==i,:][:p4], legendentry="node$i p4"),
                  ])
        push!(g, ax)
    end
    save("main2.pdf", g)
    save("main2.tex", g)
end


text:
using Pkg

Pkg.add("Yao")
Pkg.add("PyCall")

max_stars_repo_name: QuantumBFS/yao-python

max_issues_repo_name: QuantumBFS/yao-python

max_forks_repo_name: QuantumBFS/yao-python


text:
struct MySocket <: IO
    socket
end

import Base.read
read(socket::MySocket,x::Type{UInt8}) = read(socket.socket,x)

import Base.readavailable
readavailable(socket::MySocket) = readavailable(socket.socket)

#import Base.eof
#eof(socket::MySocket) = eof(socket.socket)

import Base.write
write(socket::MySocket,x::UInt8) = write(socket.socket,x)
write(socket::MySocket,x::String) = write(socket.socket,x)
write(socket::MySocket,x::Vector{UInt8}) = write(socket.socket,x)

# serialize seems to be easally replacable with write. Issue however seems to be how can I read what I send to the socket. Seems I need to implement my own serialize and deserialize methods.

using Sockets

server = listen(2000)
@async global socketA = MySocket(accept(server))
socketB = MySocket(connect(2000))


text:
##############################
## determine function names ##
##############################

const VineCppPath = joinpath(Pkg.dir("Copulas"),
                             "src/VineCopulaCPP/libVineCopulaCPP.so.1.0")

pdfName = :_Z13PairCopulaPDFiPKdPdS1_S1_j
cdfName = :_Z13PairCopulaCDFiPKdPdS1_S1_j
hfunName = :_Z14PairCopulaHfuniPKdPdS1_S1_j
vfunName = :_Z14PairCopulaVfuniPKdPdS1_S1_j
hvinName = :_Z17PairCopulaInvHfuniPKdPdS1_S1_j
vinvName = :_Z17PairCopulaInvVfuniPKdPdS1_S1_j

#######################
## VineCpp interface ##
#######################

## copula ID lookup tables
##########################

## list of copula families - order determines copula IDs!
const vineCpp_cops = [:Indep,
                      :AMH,
                      :AsymFGM,
                      :BB1,
                      :BB6,
                      :BB7,
                      :BB8,
                      :Clayton,
                      :FGM,
                      :Frank,
                      :Gaussian,
                      :Gumbel,
                      :IteratedFGM,
                      :Joe,
                      :PartialFrank,
                      :Plackett,
                      :Tawn1,
                      :Tawn2,
                      :Tawn,
                      :t]

## transform to dictionary with families as keys
const vineCpp_IDs = [vineCpp_cops[ii] => ii-1
                     for ii=1:length(vineCpp_cops)]

function getVineCppId(nam::Symbol)
    return vineCpp_IDs[nam]
end

#######################################
## create parametric VineCpp copulas ##
#######################################

copTypeNames_Cpp =
    [symbol(string(vineCpp_cops[ii], "PC_Cpp"))
     for ii=1:length(vineCpp_cops)]

## define copula types through macro
##----------------------------------

macro defineParamCopCpp(nam)
    esc(quote
        type $(nam) <: Copulas.ParamPC_Cpp
            params::FloatVec
        end
    end)
end

for cop in copTypeNames_Cpp
    eval(macroexpand(:(@defineParamCopCpp $cop)))
end

###########################
## VineCpp ID interfaces ##
###########################

## get VineCpp Id for copula object
##---------------------------------

## dict with copula type names as keys
const type2idDictCpp = [copTypeNames_Cpp[ii] => (ii-1)::Int for ii=1:length(copTypeNames_Cpp)]

function getVineCppId(cop::ParamPC_Cpp)
    return type2idDictCpp[symbol(string(typeof(cop)))]
end

## get family name for ID
##-----------------------

function getCopNam(ii::Int)
    return vineCpp_cops[ii+1]
end

## get copula type for ID
##-----------------------

function getCopType_Cpp(ii::Int)
    return copTypeNames_Cpp[ii+1]
end

#############
## display ##
#############

function displayShortString(cop::ParamPC_Cpp)
    copNam = string(typeof(cop))
    par = params(cop)
    nParams = length(par)
    paramsString = string(par[1])
    for ii=2:nParams
        paramsString = join([paramsString, ", ", string(par[ii])])
    end
    return join([copNam, ": ", paramsString])
end

function displayShortString(cop::IndepPC_Cpp)
    return "IndepPC_Cpp"
end

import Base.Multimedia.display
function display(cop::ParamPC_Cpp)
    msg = displayShortString(cop)
    println("$msg")
end

##############
## equality ##
##############

import Base.==
function ==(cop1::ParamPC_Cpp, cop2::ParamPC_Cpp)
    isequal(typeof(cop1), typeof(cop2)) || return false
    return isequal(params(cop1), params(cop2))
end

#######################
## VineCpp functions ##
#######################

## pdf function
##-------------

function pdf(cop::ParamPC_Cpp, u1::FloatVec, u2::FloatVec)
    fam, params = getFamAndParams(cop)

    nObs = length(u1)
    retVals = zeros(Float64, nObs)

    ccall((:_Z13PairCopulaPDFiPKdPdS1_S1_j, VineCppPath),
          Void,
          (Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int),
          fam, params, u1, u2, retVals, nObs)

    return retVals
end

function cdf(cop::ParamPC_Cpp, u1::FloatVec, u2::FloatVec)
    fam, params = getFamAndParams(cop)

    nObs = length(u1)
    retVals = zeros(Float64, nObs)

    ccall((:_Z13PairCopulaCDFiPKdPdS1_S1_j, VineCppPath),
          Void,
          (Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int),
          fam, params, u1, u2, retVals, nObs)

    return retVals
end

function hfun(cop::ParamPC_Cpp, u1::FloatVec, u2::FloatVec)
    fam, params = getFamAndParams(cop)

    nObs = length(u1)
    retVals = zeros(Float64, nObs)

    ccall((:_Z14PairCopulaHfuniPKdPdS1_S1_j, VineCppPath),
          Void,
          (Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int),
          fam, params, u1, u2, retVals, nObs)

    return retVals
end

function vfun(cop::ParamPC_Cpp, u2::FloatVec, u1::FloatVec)
    # u2|u1
    fam, params = getFamAndParams(cop)

    nObs = length(u1)
    retVals = zeros(Float64, nObs)

    ccall((:_Z14PairCopulaVfuniPKdPdS1_S1_j, VineCppPath),
          Void,
          (Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int),
          fam, params, u1, u2, retVals, nObs)

    return retVals
end

function hinv(cop::ParamPC_Cpp, u1::FloatVec, u2::FloatVec)
    fam, params = getFamAndParams(cop)

    nObs = length(u1)
    retVals = zeros(Float64, nObs)

    ccall((:_Z17PairCopulaInvHfuniPKdPdS1_S1_j, VineCppPath),
          Void,
          (Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int),
          fam, params, u1, u2, retVals, nObs)

    return retVals
end

function vinv(cop::ParamPC_Cpp, u2::FloatVec, u1::FloatVec)
    # u2|u1
    fam, params = getFamAndParams(cop)

    nObs = length(u1)
    retVals = zeros(Float64, nObs)

    ccall((:_Z17PairCopulaInvVfuniPKdPdS1_S1_j, VineCppPath),
          Void,
          (Int, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int),
          fam, params, u1, u2, retVals, nObs)

    return retVals
end


text:
using RegressionTables
using FixedEffectModels, GLM, RDatasets, Test

df = dataset("datasets", "iris")
df[!, :SpeciesDummy] = categorical(df[!,:Species])
df[!, :isSmall] = categorical(df[!, :SepalWidth] .< 2.9)
df[!, :isWide] = categorical(df[!, :SepalWidth] .> 2.5)

# FixedEffectModels.jl
rr1 = reg(df, @formula(SepalLength ~ SepalWidth))
rr2 = reg(df, @formula(SepalLength ~ SepalWidth + PetalLength + fe(SpeciesDummy)))
rr3 = reg(df, @formula(SepalLength ~ SepalWidth + PetalLength + PetalWidth + fe(SpeciesDummy) + fe(isSmall)))
rr4 = reg(df, @formula(SepalWidth ~ SepalLength + PetalLength + PetalWidth + fe(SpeciesDummy)))
rr5 = reg(df, @formula(SepalWidth ~ SepalLength + (PetalLength ~ PetalWidth) + fe(SpeciesDummy)))
rr6 = reg(df, @formula(SepalLength ~ SepalWidth + fe(SpeciesDummy)&fe(isWide) + fe(isSmall)))

# GLM.jl
dobson = DataFrame(Counts = [18.,17,15,20,10,20,25,13,12],
    Outcome = categorical(repeat(["A", "B", "C"], outer = 3)),
    Treatment = categorical(repeat(["a","b", "c"], inner = 3)))

lm1 = fit(LinearModel, @formula(SepalLength ~ SepalWidth), df)
lm2 = fit(LinearModel, @formula(SepalLength ~ SepalWidth + PetalWidth), df)
lm3 = fit(LinearModel, @formula(SepalLength ~ SepalWidth * PetalWidth), df) # testing interactions
gm1 = fit(GeneralizedLinearModel, @formula(Counts ~ 1 + Outcome), dobson,
              Poisson())

function checkfilesarethesame(file1::String, file2::String)

    f1 = open(file1, "r")
    f2 = open(file2, "r")

    s1 = read(f1, String)
    s2 = read(f2, String)

    close(f1)
    close(f2)

    # Character-by-character comparison
    for i=1:length(s1)
        if s1[i]!=s2[i]
            println("Character $(i) different: $(s1[i]) $(s2[i])")
        end
    end

    if s1 == s2
        return true
    else
        return false
        println("File 1:")
        @show s1
        println("File 2:")
        @show s2
    end
end

# ASCII TABLES

# default
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput())
#
# # display of statistics below estimates
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), below_statistic = :blank)
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), below_decoration = s -> "[$(s)]")
#
# # ordering of regressors, leaving out regressors
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), regressors = ["SepalLength";"PetalWidth";"SepalWidth"])
#
# # format of the estimates
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), estimformat = "%02.5f")
#
# # replace some variable names by other strings
# regtable(rr1,rr2,rr3; renderSettings = asciiOutput(), labels = Dict("SepalLength" => "My dependent variable: SepalLength", "PetalLength" => "Length of Petal", "PetalWidth" => "Width of Petal", "(Intercept)" => "Const." , "isSmall" => "isSmall Dummies", "SpeciesDummy" => "Species Dummies"))
#
# # do not print the FE block
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), print_fe_section = false)
#
# # re-order fixed effects
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), fixedeffects = ["isSmall", "SpeciesDummy"])
#
# # change the yes/no labels in the fixed effect section, and statistics labels
# regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), labels = Dict("__LABEL_FE_YES__" => "Mhm.", "__LABEL_FE_NO__" => "Nope.", "__LABEL_STATISTIC_N__" => "Number of observations", "__LABEL_STATISTIC_R2__" => "R Squared"))
#
# # full set of available statistics
# regtable(rr1,rr2,rr3,rr5; renderSettings = asciiOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])

# if you want to test locally...
# include("src/RegressionTables.jl")
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.asciiOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(rr1,rr2,rr3,rr5,rr6; renderSettings = RegressionTables.asciiOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(), regression_statistics = [:nobs, :r2], standardize_coef = true)
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.latexOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.latexOutput(), regression_statistics = [:nobs, :r2])
# RegressionTables.regtable(lm1, lm2, lm3, gm1; renderSettings = RegressionTables.latexOutput(), regression_statistics = [:nobs, :r2], transform_labels = RegressionTables.escape_ampersand)
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.htmlOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.htmlOutput(), regression_statistics = [:nobs, :r2])
# --------------
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.latexOutput("test2.tex"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.latexOutput("test4.tex"), regression_statistics = [:nobs, :r2])
# RegressionTables.regtable(lm1, lm2, lm3, gm1; renderSettings = RegressionTables.latexOutput("test6.tex"), regression_statistics = [:nobs, :r2], transform_labels = RegressionTables.escape_ampersand)
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.asciiOutput("test1.txt"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(rr1,rr2,rr3,rr5,rr6; renderSettings = RegressionTables.asciiOutput("test7.txt"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("test3.txt"), regression_statistics = [:nobs, :r2])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("test5.txt"), regression_statistics = [:nobs, :r2], standardize_coef = true)
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.latexOutput("test2.txt"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.latexOutput("test4.txt"), regression_statistics = [:nobs, :r2])
# RegressionTables.regtable(lm1, lm2, lm3, gm1; renderSettings = RegressionTables.latexOutput("test6.txt"), regression_statistics = [:nobs, :r2], transform_labels = RegressionTables.escape_ampersand)
# RegressionTables.regtable(rr1,rr2,rr3,rr5; renderSettings = RegressionTables.htmlOutput("test1.html"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# RegressionTables.regtable(lm1, lm2, gm1; renderSettings = RegressionTables.htmlOutput("test2.html"), regression_statistics = [:nobs, :r2])

# # new tests: all features
# RegressionTables.regtable(rr4,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("ftest1.txt"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# # regressors and labels
# RegressionTables.regtable(rr4,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("ftest2.txt"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof], regressors = ["SepalLength", "PetalWidth"])
# # fixedeffects, estimformat, statisticformat, number_regressions_decoration
# RegressionTables.regtable(rr3,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("ftest3.txt"), fixedeffects = ["SpeciesDummy"], estimformat = "%0.4f", statisticformat = "%0.4f", number_regressions_decoration = i -> "[$i]")
# # estim_decoration, below_statistic, below_decoration, number_regressions
# function dec(s::String, pval::Float64)
#     if pval<0.0
#         error("p value needs to be nonnegative.")
#     end
#     if (pval > 0.05)
#         return "$s"
#     elseif (pval <= 0.05)
#         return "$s<-OMG!"
#     end
# end
# RegressionTables.regtable(rr3,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("ftest4.txt"), estim_decoration = dec, below_statistic = :tstat, below_decoration = s -> "[$s]", number_regressions = false)
# # print_fe_section, print_estimator_section
# RegressionTables.regtable(rr3,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput("ftest5.txt"), print_fe_section = false, print_estimator_section = false)
# # transform_labels and custom labels
# RegressionTables.regtable(rr5,rr6,lm1, lm2, lm3; renderSettings = RegressionTables.asciiOutput("ftest6.txt"), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof], transform_labels = RegressionTables.escape_ampersand,
# labels = Dict("SepalLength" => "My dependent variable: SepalLength", "PetalLength" => "Length of Petal", "PetalWidth" => "Width of Petal", "(Intercept)" => "Const." , "isSmall" => "isSmall Dummies", "SpeciesDummy" => "Species Dummies"))

RegressionTables.regtable(rr4,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "ftest1.txt")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "ftest1.txt"), joinpath(dirname(@__FILE__), "tables", "ftest1_reference.txt"))
# regressors and labels
RegressionTables.regtable(rr4,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "ftest2.txt")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof], regressors = ["SepalLength", "PetalWidth"])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "ftest2.txt"), joinpath(dirname(@__FILE__), "tables", "ftest2_reference.txt"))
# fixedeffects, estimformat, statisticformat, number_regressions_decoration
RegressionTables.regtable(rr3,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "ftest3.txt")), fixedeffects = ["SpeciesDummy"], estimformat = "%0.4f", statisticformat = "%0.4f", number_regressions_decoration = i -> "[$i]")
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "ftest3.txt"), joinpath(dirname(@__FILE__), "tables", "ftest3_reference.txt"))
# estim_decoration, below_statistic, below_decoration, number_regressions
function dec(s::String, pval::Float64)
    if pval<0.0
        error("p value needs to be nonnegative.")
    end
    if (pval > 0.05)
        return "$s"
    elseif (pval <= 0.05)
        return "$s<-OMG!"
    end
end
RegressionTables.regtable(rr3,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "ftest4.txt")), estim_decoration = dec, below_statistic = :tstat, below_decoration = s -> "[$s]", number_regressions = false)
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "ftest4.txt"), joinpath(dirname(@__FILE__), "tables", "ftest4_reference.txt"))
# print_fe_section, print_estimator_section
RegressionTables.regtable(rr3,rr5,lm1, lm2, gm1; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "ftest5.txt")), print_fe_section = false, print_estimator_section = false)
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "ftest5.txt"), joinpath(dirname(@__FILE__), "tables", "ftest5_reference.txt"))
# transform_labels and custom labels
RegressionTables.regtable(rr5,rr6,lm1, lm2, lm3; renderSettings = RegressionTables.asciiOutput(joinpath(dirname(@__FILE__), "tables", "ftest6.txt")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof], transform_labels = RegressionTables.escape_ampersand,
labels = Dict("SepalLength" => "My dependent variable: SepalLength", "PetalLength" => "Length of Petal", "PetalWidth" => "Width of Petal", "(Intercept)" => "Const." , "isSmall" => "isSmall Dummies", "SpeciesDummy" => "Species Dummies"))
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "ftest6.txt"), joinpath(dirname(@__FILE__), "tables", "ftest6_reference.txt"))

regtable(rr1,rr2,rr3,rr5; renderSettings = asciiOutput(joinpath(dirname(@__FILE__), "tables", "test1.txt")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test1.txt"), joinpath(dirname(@__FILE__), "tables", "test1_reference.txt"))

regtable(rr1,rr2,rr3,rr5,rr6; renderSettings = asciiOutput(joinpath(dirname(@__FILE__), "tables", "test7.txt")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test7.txt"), joinpath(dirname(@__FILE__), "tables", "test7_reference.txt"))

regtable(lm1, lm2, gm1; renderSettings = asciiOutput(joinpath(dirname(@__FILE__), "tables", "test3.txt")), regression_statistics = [:nobs, :r2])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test3.txt"), joinpath(dirname(@__FILE__), "tables", "test3_reference.txt"))

regtable(lm1, lm2, gm1; renderSettings = asciiOutput(joinpath(dirname(@__FILE__), "tables", "test5.txt")), regression_statistics = [:nobs, :r2], standardize_coef = true)
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test5.txt"), joinpath(dirname(@__FILE__), "tables", "test5_reference.txt"))

# LATEX TABLES

# # default
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput())
#
# # display of statistics below estimates
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), below_statistic = :blank)
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), below_decoration = s -> "[$(s)]")
#
# # ordering of regressors, leaving out regressors
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), regressors = ["SepalLength";"PetalWidth";"SepalWidth"])
#
# # format of the estimates
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), estimformat = "%02.5f")
#
# # replace some variable names by other strings
# regtable(rr1,rr2,rr3; renderSettings = latexOutput(), labels = Dict("SepalLength" => "My dependent variable: SepalLength", "PetalLength" => "Length of Petal", "PetalWidth" => "Width of Petal", "(Intercept)" => "Const." , "isSmall" => "isSmall Dummies", "SpeciesDummy" => "Species Dummies"))
#
# # do not print the FE block
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), print_fe_section = false)
#
# # re-order fixed effects
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), fixedeffects = ["isSmall", "SpeciesDummy"])
#
# # change the yes/no labels in the fixed effect section, and statistics labels
# regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput(), labels = Dict("__LABEL_FE_YES__" => "Mhm.", "__LABEL_FE_NO__" => "Nope.", "__LABEL_STATISTIC_N__" => "Number of observations", "__LABEL_STATISTIC_R2__" => "R Squared"))
#
# # full set of available statistics
# regtable(rr1,rr2,rr3,rr5; renderSettings = latexOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
#

regtable(rr1,rr2,rr3,rr5; renderSettings = latexOutput(joinpath(dirname(@__FILE__), "tables", "test2.tex")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test2.tex"), joinpath(dirname(@__FILE__), "tables", "test2_reference.tex"))

regtable(lm1, lm2, gm1; renderSettings = latexOutput(joinpath(dirname(@__FILE__), "tables", "test4.tex")), regression_statistics = [:nobs, :r2])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test4.tex"), joinpath(dirname(@__FILE__), "tables", "test4_reference.tex"))

regtable(lm1, lm2, lm3, gm1; renderSettings = latexOutput(joinpath(dirname(@__FILE__), "tables", "test6.tex")), regression_statistics = [:nobs, :r2], transform_labels = escape_ampersand)
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test6.tex"), joinpath(dirname(@__FILE__), "tables", "test6_reference.tex"))

# HTML Tables
regtable(rr1,rr2,rr3,rr5; renderSettings = htmlOutput(joinpath(dirname(@__FILE__), "tables", "test1.html")), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test1.html"), joinpath(dirname(@__FILE__), "tables", "test1_reference.html"))

regtable(lm1, lm2, gm1; renderSettings = htmlOutput(joinpath(dirname(@__FILE__), "tables", "test2.html")), regression_statistics = [:nobs, :r2])
@test checkfilesarethesame(joinpath(dirname(@__FILE__), "tables", "test2.html"), joinpath(dirname(@__FILE__), "tables", "test2_reference.html"))

# clean up
rm(joinpath(dirname(@__FILE__), "tables", "test1.txt"))
rm(joinpath(dirname(@__FILE__), "tables", "test2.tex"))
rm(joinpath(dirname(@__FILE__), "tables", "test3.txt"))
rm(joinpath(dirname(@__FILE__), "tables", "test4.tex"))
rm(joinpath(dirname(@__FILE__), "tables", "test5.txt"))
rm(joinpath(dirname(@__FILE__), "tables", "test1.html"))
rm(joinpath(dirname(@__FILE__), "tables", "test2.html"))


text:
## Re-exports. ##

export

    AbstractEntry,
    Docs,
    Entry,
    Manual,
    Metadata,
    Page

## Module. ##

export documented, isdocumented, metadata

"""
Storage for deprecated ``Metadata`` documentation.

!!set(state:deprecated)
"""
const DOCUMENTED = ObjectIdDict()

"""
Returns the modules that are currently documented by Docile.
"""
documented() = Cache.loadedmodules()

"""
Is the module ``mod`` documented by Docile?
"""
isdocumented(mod::Module) = mod in documented()

"""
Get the ``Metadata`` object associated with a module ``mod``.

!!set(status:deprecated)
"""
function metadata(mod::Module; rebuild = false)
    if isdocumented(mod) && !haskey(DOCUMENTED, mod)
        DOCUMENTED[mod] = Metadata(mod)
    end
    DOCUMENTED[mod]
end

## Metadata. ##

export modulename, manual, entries, root, files, isloaded, isexported

"""
The ``Module`` that a ``Metadata`` object documents.

!!set(status:deprecated)
"""
modulename(meta::Metadata) = meta.modname

"""
The manual files for a ``Metadata`` object ``meta``.

!!set(status:deprecated)
"""
manual(meta::Metadata) = meta.data[:manual]

"""
``ObjectIdDict`` containing documented objects and their associated ``Entry``s.

!!set(status:deprecated)
"""
entries(meta::Metadata) = meta.entries

"""
The rootfile of the module documented by a ``Metadata`` object ``meta``.

!!set(status:deprecated)
"""
root(meta::Metadata) = meta.root

"""
List of all ``include``d files in a module documented by ``Metadata`` object ``meta``.

!!set(status:deprecated)
"""
files(meta::Metadata) = meta.files

"""
Have the docstrings contained in a module been collected yet?

!!set(status:deprecated)
"""
isloaded(meta::Metadata) = meta.loaded

"""
The ``Dict{Symbol, Any}`` containing arbitrary additional data about a ``Metadata`` object.

!!set(status:deprecated)
"""
metadata(meta::Metadata) = meta.data

"""
Is the documented object ``object`` been exported from the given module ``modname``?
"""
function isexported(modname::Module, object)
    meta = metadata(modname)
    data = metadata(meta)

    ent = entries(meta)[object]
    cat = category(ent)

    sym = cat == :macro ? macroname(ent) : name(object)

    sym in data[:exports]
end

isexported(modname::Module, comment::Comment) = false

macroname(ent) = symbol(string("@", metadata(ent)[:signature].args[1]))

name(f::Function) = methods(f).name
name(m::Method)   = Utilities.lsdfield(m, :name)
name(m::Module)   = module_name(m)
name(t::DataType) = t.name.name
name(s::Symbol)   = s
name(s::AbstractString) = s

"""
Get the ``Symbol`` representing an object such as ``Function`` or ``Method``.
"""
name

## Entry. ##

export category, modulename, metadata, docs

"""
What category does an ``Entry`` object belong to?

!!set(status:deprecated)
"""
category{C}(::Entry{C}) = C

"""
Which module does the ``Entry`` object come from?

!!set(status:deprecated)
"""
modulename(e::Entry) = e.modname

"""
Arbitrary additional metadata associated with a particular ``Entry`` ``e``.

!!set(status:deprecated)
"""
metadata(e::Entry) = e.data

"""
The ``Docs`` object for an ``Entry`` object ``e``.

!!set(status:deprecated)
"""
docs(e::Entry) = e.docs

## Docs. ##

export data, format, parsed, parsedocs

"""
Raw docstring associated with a ``Docs`` object ``d``.

!!set(status:deprecated)
"""
data(d::Docs) = d.data

"""
The format that a docstring is written in.

!!set(status:deprecated)
"""
format{F}(d::Docs{F}) = F

"""
Get the parsed docstring for a ``Docs`` object ``d``.

!!set(status:deprecated)
"""
parsed(d::Docs) = isdefined(d, :obj) ? d.obj : (d.obj = parsedocs(d);)

parsedocs{ext}(d::Docs{ext}) = error("Unknown documentation format: $(ext)")

parsedocs(d::Docs{:txt}) = data(d)

"""
Parsing hook for specifying how to parse raw docstrings into formatted text.
"""
parsedocs

## Deprecated. ##

export documentation

"""
The ``Metadata`` object associated with a module ``mod``.

!!set(status:deprecated)
"""
documentation(mod::Module) = metadata(mod)


text:
using GLVisualize, GLAbstraction, GeometryTypes, Reactive

if !isdefined(:runtests)
	window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end

function mgrid(dim1, dim2)
    X = [i for i in dim1, j in dim2]
    Y = [j for i in dim1, j in dim2]
    return X,Y
end

rotation_angle = const_lift(*, timesignal, 2f0*pi)
rotation = map(rotationmatrix_z, rotation_angle)

dphi, dtheta = pi/200.0f0, pi/200.0f0

phi,theta = mgrid(0f0:dphi:(pi+dphi*1.5f0), 0f0:dtheta:(2f0*pi+dtheta*1.5f0));
m0 = 4f0; m1 = 3f0; m2 = 2f0; m3 = 3f0; m4 = 6f0; m5 = 2f0; m6 = 6f0; m7 = 4f0;
a = sin(m0*phi).^m1;
b = cos(m2*phi).^m3;
c = sin(m4*theta).^m5;
d = cos(m6*theta).^m7;
r = a + b + c + d;
x = r.*sin(phi).*cos(theta);
y = r.*cos(phi);
z = r.*sin(phi).*sin(theta);

surface = visualize((x,y,z), :surface, model=rotation)
_view(surface, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# ################################################################
#
# CAVEAT: THIS EXAMPLE DOES NOT WORK YET
#
# REASON: the use of StaticArrays causes some problems:
#         https://github.com/JuliaArrays/StaticArrays.jl/issues/971
#
# NOTE:   Rosenbrock_static is simply an alias of Rosenbrock and not a
#         static version of Rosenbrock as it should be.
#
# -> TO FIX
#
# ################################################################

# An example to count alloc static vs dynamic arrays
#
using Pkg, BenchmarkTools

Pkg.activate("../")

using NLS_Solver

include("../test/TestProblems/TestProblems.jl")

θ = [1.0,1.0]

obj = Rosenbrock()
obj_static = Rosenbrock_Static()

conf = LevenbergMarquardt_Conf()

println("Dynamic : ", @btime solve($obj, $θ, $conf))
println("Static  : ", @btime solve($obj_static, $θ, $conf))

# Sanity check
#
result = solve(obj, θ, conf)
result_static = solve(obj_static, θ, conf)

@assert converged(result)
@assert converged(result_static)
@assert solution(result) ≈ solution(result_static)

println("Sanity check ok")

# Constrained problem

bc = BoundConstraints(2)

conf = LevenbergMarquardt_BC_Conf()

println("Dynamic : ", @btime solve($obj, $θ, $bc, $conf))
println("Static  : ", @btime solve($obj_static, $θ, $bc, $conf))

#
# Debug runtime error for StaticArrays
#
# using Debugger
# breakpoint(joinpath(@__DIR__,"src/Levenberg-Marquardt/lm_bc.jl"),152)
# @enter solve(obj_static, θ, bc, conf)


text:
__precompile__()

using Ensemblearn

export
    DecisionTree,
    fit!,
    predict,
    predict_probs,
    print_tree

# Based on https://github.com/eriklindernoren/ML-From-Scratch/blob/master/mlfromscratch/supervised_learning/decision_tree.py

mutable struct DecisionNode
    feature_i::Union{Int64, Nothing}          # Index for the feature that is tested
    threshold::Union{Float64, Nothing}          # Threshold value for feature
    value#::Union{Float64, Nothing}                  # Value if the node is a leaf in the tree
    true_branch::Union{DecisionNode,Nothing}      # 'Left' subtree
    false_branch::Union{DecisionNode, Nothing}    # 'Right' subtree

    function DecisionNode(;feature_i=nothing::Union{Int64, Nothing}, threshold=nothing::Union{Float64, Nothing},
                          value=nothing, true_branch=nothing::Union{DecisionNode,Nothing}, false_branch=nothing::Union{DecisionNode, Nothing})
        new(feature_i, threshold, value, true_branch, false_branch)
    end
end

# Super class of ClassificationTree
mutable struct DecisionTree
    # Minimum n of samples to justify split
    min_samples_split::Int64
    # The minimum impurity to justify split
    min_impurity::Float64
    # The maximum depth to grow the tree to
    max_depth::Int64
    # Impurity method to use on information gain calculation
    impurity::String
    # Function to determine prediction of y at leaf
    #_leaf_value_calculation = None
    root  # Root node in dec. tree
    classes

    function DecisionTree(min_samples_split=2::Int64, min_impurity=1e-7::Float64, max_depth=10::Int64,
                          impurity="gini"::String, classes=nothing)
        # TODO: here assert impurity ?
        new(min_samples_split, min_impurity, max_depth, impurity, nothing, classes)
    end
end

impurity_methods = Dict("entropy" => entropy, "gini" => gini_index, "cerr" => classification_error)

function calculate_information_gain(y::Array{T,2}, y1::Array{T2,1}, y2::Array{T2,1}, impurity="gini"::String) where{T <: Real, T2 <: Real}  # TODO: y::Array{T,1}

    @assert haskey(impurity_methods, impurity)
    i = impurity_methods[impurity]

    # Calculate information gain
    p = length(y1) / length(y)
    info_gain = i(y) - p * i(y1) - (1 - p) * i(y2)

    return info_gain
end

function fill_dict!(d, ks, v=0)
    for k in ks
        if !haskey(d, k)
            d[k] = v;
        end
    end
end

# TODO: partial fit, by fingerprinting samples to add counts (it is necessary to fingerprint?)

function fit!(tree::DecisionTree, x::Array{T,2}, y::Array{T2,1}; current_depth=0::Int64) where {T <: Real, T2 <: Real}
    """ Recursive method which builds out the decision tree and splits x and respective y
    on the feature of x which (based on impurity) best separates the data"""

    largest_impurity = 0
    best_criteria = nothing    # Feature index and threshold
    best_sets = nothing        # Subsets of the data

    # Check if expansion of y is needed
    # TODO: is this necesary?
    if ndims(y) == 1
        y = reshape(y, length(y), 1)
    end

    # Add y as last column of x
    xy = hcat(x, y)

    if tree.classes == nothing
        tree.classes = unique(y)
    end

    n_samples, n_features = size(x)

    if (n_samples >= tree.min_samples_split) && (current_depth <= tree.max_depth)
        # Calculate the impurity for each feature
        for feature_i in 1:n_features
            # All values of feature_i
            feature_values = x[:, feature_i]
            unique_values = unique(feature_values)

            # Iterate through all unique values of feature column i and
            # calculate the impurity
            for threshold in unique_values
                # Divide X and y depending on if the feature value of X at index feature_i
                # meets the threshold

                # TODO: just get indexes to avoid mixing x and y types in a single array
                xy1, xy2 = divide_on_feature(xy, feature_i, threshold)

                if (size(xy1,1) > 0) && (size(xy2,1) > 0)
                    # Select the y-values of the two sets
                    y1 = xy1[:, end]
                    y2 = xy2[:, end]

                    # Calculate information gain
                    # TODO: impurity (general method) instead of information gain (particular)
                    impurity = calculate_information_gain(y, y1, y2, tree.impurity)

                    # If this threshold resulted in a higher information gain than previously
                    # recorded save the threshold value and the feature
                    # index
                    if impurity > largest_impurity
                        largest_impurity = impurity
                        best_criteria = Dict("feature_i"=> feature_i, "threshold"=> threshold)
                        best_sets = Dict(
                            "leftx"=> xy1[:, 1:end-1],   # X of left subtree
                            "lefty"=> y1,   # y of left subtree
                            "rightx"=> xy2[:, 1:end-1],  # X of right subtree
                            "righty"=> y2   # y of right subtree
                        )
                    end
                end
            end
        end
    end

    #tree.one_dim = size(y, 1) == 1
    #tree.loss = nothing

    if largest_impurity > tree.min_impurity
        # Build subtrees for the right and left branches
        true_branch = DecisionTree(tree.min_samples_split, tree.min_impurity, tree.max_depth, tree.impurity, tree.classes)
        false_branch = DecisionTree(tree.min_samples_split, tree.min_impurity, tree.max_depth, tree.impurity, tree.classes)

        fit!(true_branch, best_sets["leftx"], best_sets["lefty"], current_depth=current_depth + 1)
        fit!(false_branch, best_sets["rightx"], best_sets["righty"], current_depth=current_depth + 1)

        tree.root = DecisionNode(feature_i=best_criteria["feature_i"], threshold=best_criteria["threshold"],
                                 true_branch=true_branch.root, false_branch=false_branch.root)
    else
        # We're at leaf => determine value
        leaf_value = count_distinct(y)
        fill_dict!(leaf_value, tree.classes)

        tree.root = DecisionNode(value=leaf_value)
    end
end

function get_ratios_dict(d)
    n = sum(values(d))
    return Dict([(x[1], x[2]/n) for x in d])
end

function likelihood(n, x, delta=1e-4)
    # TODO: test this!
    # From https://www.datascience.com/blog/introduction-to-bayesian-inference-learn-data-science-tutorials
    """
    likelihood function for a binomial distribution

    n: [int] the number of experiments
    x: [int] the number of successes
    theta: [float] the proposed probability of success
    """
    lk = theta ->  (factorial(n) / (factorial(x) * factorial(n - x))) * (theta ^ x) * ((1 - theta) ^ (n - x))

    probs = 0:delta:1

    likelihoods = map(lk, probs)

    return probs[argmax(likelihoods)]
end

function predict_value(model::DecisionTree, x::Array{T,1}, node=nothing::Union{DecisionNode, Nothing};
                       get_probs=false::Bool) where T <: Real
    """ Do a recursive search down the tree and make a prediction of the data sample by the
        value of the leaf that we end up at """

    if node == nothing
        node = model.root
    end

    # If we have a value (i.e we're at a leaf) => return value as the prediction
    if node.value != nothing
        if get_probs == true
            v = get_ratios_dict(node.value)
        else
            v = findmax(node.value)[2]
        end

        return v
    end

    # Choose the feature that we will test)
    feature_value = x[node.feature_i]

    # Determine if we will follow left or right branch
    branch = node.false_branch
    if feature_value isa Number
        if feature_value >= node.threshold
            branch = node.true_branch
        end
    elseif feature_value == node.threshold
        branch = node.true_branch
    end

    # Test subtree
    return predict_value(model, x, branch, get_probs=get_probs)
end

function predict(model::DecisionTree, x::Array{T,2}) where T <: Real
    """ Classify samples one by one and return the set of labels """
    pred = []

    for i in 1:size(x, 1)
        y  = predict_value(model, x[i,:], get_probs=false)
        push!(pred, y)
    end

    pred
end

function predict_probs(model::DecisionTree, x::Array{T,2}) where T <: Real
    """ Classify samples one by one and return the set of labels """
    pred = Dict{Any,Array{Float64,1}}()

    for i in 1:size(x, 1)
        y  = predict_value(model, x[i,:], get_probs=true)

        for kv in y
            if !haskey(pred, kv[1])
                pred[kv[1]] = []
            end
            push!(pred[kv[1]], kv[2])
        end
    end

    pred
end

function print_tree(tree::DecisionTree, node=nothing::Union{DecisionNode,Nothing}; indent=" "::String)
    """ Recursively print the decision tree """
    if node == nothing
        node = tree.root
    end

    # If we're at leaf => print the label
    if node.value != nothing
        println(node.value)
    # Go deeper down the tree
    else
        # Print test
        # TODO: print column name when trees are fitted with DFs
        println("Feature_$(node.feature_i):$(node.threshold)?")
        # Print the true scenario
        println("$(indent)T->")
        print_tree(tree, node.true_branch, indent=indent*indent)
        # Print the false scenario
        println("$(indent)F->")
        print_tree(tree, node.false_branch, indent=indent*indent)
    end
end