text:
function (control::PlotlinesControl)(plot_context::PlotContext, interval_context::NestedIntervalContext)
    plot_control = plot_context.control
    plot_model = plot_context.model
    plot_display_properties = plot_context.display_properties

    interval_control = interval_context.control
    interval_model = interval_context.model
    interval_display_properties = interval_context.display_properties

    data = get_data(plot_model)
    i₀ = Int(get_start(interval_model))
    i₁ = Int(get_stop(interval_model))
    plot_model₂ = PlotlinesModel(data[i₀:i₁])

    #CImGui.Begin("EDA",C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus)
    isenabled(plot_control) ? plot_control(plot_model₂, plot_display_properties) : nothing
    CImGui.Dummy(ImVec2(0, 20))
    isenabled(interval_control) ? interval_control(interval_model, interval_display_properties) : nothing
    #CImGui.End()
end


text:
export molteno_dim, molteno_boxing

################################################################################
# Molteno histogram based dimension by boxing values
################################################################################
"""
    molteno_dim(data::Dataset; k0 = 10, q = 1.0, base = ℯ)
Calculate the generalized dimension using the algorithm for box division defined
by Molteno[^Molteno1993].

## Description
Divide the data into boxes with each new box having half the side length of the
former box using [`molteno_boxing`](@ref). Break if the number of points over
the number of filled boxes falls below `k0`. Then the generalized dimension can
be calculated by using [`genentropy`](@ref) to calculate the sum over the
logarithm also considering possible approximations and fitting this to the
logarithm of one over the boxsize using [`linear_region`](@ref).

This algorithm is faster than the traditional approach of using [`probabilities`](@ref),
but it is only suited for low dimensional data since it divides each
box into `2^D` new boxes if `D` is the dimension. For large `D` this leads to low numbers of
box divisions before the threshold is passed and the divison stops. This results
to a low number of data points to fit the dimension to and thereby a poor estimate.

[^Molteno1993]: Molteno, T. C. A., [Fast O(N) box-counting algorithm for estimating dimensions. Phys. Rev. E 48, R3263(R) (1993)](https://doi.org/10.1103/PhysRevE.48.R3263)
"""
function molteno_dim(data; k0 = 10, α = nothing, q=1.0, base = ℯ)
    if α ≠ nothing
        @warn "Keyword `α` is deprecated in favor of `q`."
        q = α
    end
    boxes, εs = molteno_boxing(data; k0)
    dd = genentropy.(boxes; q, base)
    return linear_region(-log.(base, εs), dd)[2]
end

"""
    molteno_boxing(data::Dataset; k0 = 10) → (probs, εs)
Distribute the `data` into boxes whose size is halved in each step. Stop if the
average number of points per filled box falls below the threshold `k0`.

Return `probs`, a vector of `Propabilities` for different box sizes and the
corresponding box sizes `εs`.

## Description
Project the `data` onto the whole interval of numbers that is covered by
`UInt64`. This projected data is distributed into boxes whose size
decreases by factor 2 in each step. For each box that contains more than one
point `2^D` new boxes are created where `D` is the dimension of the data.

The process of dividing the data into new boxes stops when the number of points
over the number of filled boxes falls below `k0`. The box sizes `εs` are
calculated and returned together with the `probs`.

See[^Molteno1993] for more.
"""
function molteno_boxing(data::AbstractDataset; k0 = 10)
    integers, ε0 = real_to_uint64(data)
    boxes = _molteno_boxing(integers; k0)
    εs = ε0 ./ (2 .^ (1:length(boxes)))
    return boxes, εs
end

"""
    real_to_uint64(data::Dataset{D,T}) where {D, T}
Calculate maximum and minimum value of `data` to then project the values onto
``[0, M - ε(M)]`` where ``\\epsilon`` is the
precision of the used Type and ``M`` is the maximum value of the UInt64 type.
"""
function real_to_uint64(data::AbstractDataset{D,T}) where {D,T<:Real}
    N = length(data)
    # The maximum value needs to be smaller than the absolute typemax due to
    # the 12 bits used by the float type for sign and exponent.
    max_val = typemax(UInt64) - 2(UInt64 |> typemax |> T |> eps)
    mins, maxs = minmaxima(data)
    spans = maxs .- mins
    ε0 = maximum(spans)
    # Let f:[min,max] -> [0,typemax(UInt64)] linearly, then f(x) = m*x + b.
    m = max_val ./ ε0
    b = - mins .* m
    res = Vector{SVector{D,UInt64}}()
    sizehint!(res, N)
    for x in data
        int_val = floor.(UInt64, m .* x .+ b)
        push!(res, int_val)
    end
    Dataset(res), ε0
end

function _molteno_boxing(data::Dataset{D,T}; k0 = 10) where {D,T<:UInt}
    N = length(data)
    box_probs = Vector{Float64}[]
    iteration = 1
    boxes = [[1:N;]]
    while N / length(boxes) > k0
        l = length(boxes)
        for t in 1:l
            # Take the first box.
            box = popfirst!(boxes)
            # Check if only one element is contained.
            if length(box) == 1
                push!(boxes, box)
                continue
            end
            # Append a new partitioned box.
            append!(boxes, molteno_subboxes(box, data, iteration))
        end
        # Calculate probabilities by dividing the number of elements in a box by N.
        push!(box_probs, length.(boxes) ./ N)
        iteration += 1
    end
    Probabilities.(box_probs)
end

"""
    molteno_subboxes(box, data, iteration)
Divide a `box` containing indices into `data` to `2^D` smaller boxes and sort
the points contained in the former box into the new boxes. Implemented according to
Molteno[^Molteno]. Sorts the elements of the former box into the smaller boxes
using cheap bit shifting and `&` operations on the value of `data` at each box
element. `iteration` determines which bit of the array should be shifted to the
last position.
"""
function molteno_subboxes(box, data::AbstractDataset{D,UInt64}, iteration) where {D}
    new_boxes = [UInt64[] for i in 1:2^D]
    index_multipliers = [2^i for i in 0:D-1]
    sorting_number = 64-iteration
    for elem in box
        index = one(UInt64)
        for (i, multi) in enumerate(index_multipliers)
            # index shifting magic
            index += ((data[elem][i] >> sorting_number) & 1) * multi
        end
        push!(new_boxes[index], elem)
    end
    filter!(!isempty, new_boxes)
end

ext: jl


text:
@testset "FlipX" begin
    @test typeof(@inferred(FlipX())) <: FlipX <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test @inferred(FlipX(0.7)) === Either(FlipX(), 0.7)
        @test str_show(FlipX()) == "Augmentor.FlipX()"
        @test str_showconst(FlipX()) == "FlipX()"
        @test str_showcompact(FlipX()) == "Flip the X axis"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(FlipX(), nothing)
        @test Augmentor.supports_eager(FlipX) === true
        res1 = reverse(rect; dims=2)
        imgs = [
            (rect, res1),
            (Augmentor.prepareaffine(rect), res1),
            (OffsetArray(rect, -2, -1), res1),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res1),
        ]
        @testset "single image" begin
            for (img_in, img_out) in imgs
                res = @inferred(Augmentor.applyeager(FlipX(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for (img_in1, img_out1) in imgs, (img_in2, img_out2) in imgs
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(FlipX(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(FlipX) === true
        @test_throws MethodError Augmentor.applyaffine(FlipX(), nothing)
        @test @inferred(Augmentor.toaffinemap(FlipX(), rect)) ≈ AffineMap([1. 0.; 0. -1.], [0.0,4.0])
        @testset "single image" begin
            wv = @inferred Augmentor.applyaffine(FlipX(), Augmentor.prepareaffine(square))
            @test parent(wv).itp.coefs === square
            @test wv == reverse(square; dims=2)
            @test typeof(wv) <: InvWarpedView{eltype(square),2}
        end
        @testset "multiple images" begin
            img_in = Augmentor.prepareaffine.((rgb_rect, square))
            res1, res2 = @inferred(Augmentor.applyaffine(FlipX(), img_in))
            # make sure affine map is specific to image
            @test res1 == reverse(rgb_rect; dims=2)
            @test res2 == reverse(square; dims=2)
            @test typeof(res1) <: InvWarpedView{eltype(rgb_rect),2}
            @test typeof(res2) <: InvWarpedView{eltype(square),2}
        end
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(FlipX) === true
        @test_throws MethodError Augmentor.applyaffineview(FlipX(), nothing)
        wv = @inferred Augmentor.applyaffineview(FlipX(), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == reverse(square; dims=2)
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(FlipX) === true
        @testset "single image" begin
            v = @inferred Augmentor.applylazy(FlipX(), rect)
            @test v === view(rect, 1:1:2, 3:-1:1)
            @test v == reverse(rect; dims=2)
            @test typeof(v) <: SubArray
            wv = @inferred Augmentor.applylazy(FlipX(), Augmentor.prepareaffine(square))
            @test parent(wv).itp.coefs === square
            @test wv == reverse(square; dims=2)
            @test typeof(wv) <: InvWarpedView{eltype(square),2}
        end
        @testset "multiple images" begin
            img_in = (rgb_rect, square)
            res1, res2 = @inferred(Augmentor.applylazy(FlipX(), img_in))
            @test res1 == reverse(rgb_rect; dims=2)
            @test res2 == reverse(square; dims=2)
            @test typeof(res1) <: SubArray{eltype(rgb_rect),2}
            @test typeof(res2) <: SubArray{eltype(square),2}
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(FlipX) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(FlipX) === true
        v = @inferred Augmentor.applylazy(FlipX(), rect)
        @test v === view(rect, 1:1:2, 3:-1:1)
        @test v == reverse(rect; dims=2)
        @test typeof(v) <: SubArray
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(FlipX) === false
    end
end

# --------------------------------------------------------------------

@testset "FlipY" begin
    @test typeof(@inferred(FlipY())) <: FlipY <: Augmentor.AffineOperation
    @testset "constructor" begin
        @test @inferred(FlipY(0.7)) === Either(FlipY(), 0.7)
        @test str_show(FlipY()) == "Augmentor.FlipY()"
        @test str_showconst(FlipY()) == "FlipY()"
        @test str_showcompact(FlipY()) == "Flip the Y axis"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(FlipY(), nothing)
        @test Augmentor.supports_eager(FlipY) === true
        res1 = reverse(rect; dims=1)
        imgs = [
            (rect, res1),
            (Augmentor.prepareaffine(rect), res1),
            (OffsetArray(rect, -2, -1), res1),
            (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res1),
        ]
        @testset "single image" begin
            for (img_in, img_out) in imgs
                res = @inferred(Augmentor.applyeager(FlipY(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
        @testset "multiple images" begin
            for (img_in1, img_out1) in imgs, (img_in2, img_out2) in imgs
                img_in = (img_in1, img_in2)
                img_out = (img_out1, img_out2)
                res = @inferred(Augmentor.applyeager(FlipY(), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(FlipY) === true
        @test_throws MethodError Augmentor.applyaffine(FlipY(), nothing)
        @test @inferred(Augmentor.toaffinemap(FlipY(), rect)) ≈ AffineMap([-1. 0.; 0. 1.], [3.0,0.0])
        @testset "single image" begin
            wv = @inferred Augmentor.applyaffine(FlipY(), Augmentor.prepareaffine(square))
            @test parent(wv).itp.coefs === square
            @test wv == reverse(square; dims=1)
            @test typeof(wv) <: InvWarpedView{eltype(square),2}
        end
        @testset "multiple images" begin
            img_in = Augmentor.prepareaffine.((rgb_rect, square))
            res1, res2 = @inferred(Augmentor.applyaffine(FlipY(), img_in))
            # make sure affine map is specific to image
            @test res1 == reverse(rgb_rect; dims=1)
            @test res2 == reverse(square; dims=1)
            @test typeof(res1) <: InvWarpedView{eltype(rgb_rect),2}
            @test typeof(res2) <: InvWarpedView{eltype(square),2}
        end
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(FlipY) === true
        @test_throws MethodError Augmentor.applyaffineview(FlipY(), nothing)
        wv = @inferred Augmentor.applyaffineview(FlipY(), Augmentor.prepareaffine(square))
        @test typeof(wv) <: SubArray{eltype(square),2}
        @test typeof(parent(wv)) <: InvWarpedView
        @test parent(parent(wv)).itp.coefs === square
        @test wv == reverse(square; dims=1)
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(FlipY) === true
        @testset "single image" begin
            v = @inferred Augmentor.applylazy(FlipY(), rect)
            @test v === view(rect, 2:-1:1, 1:1:3)
            @test v == reverse(rect; dims=1)
            @test typeof(v) <: SubArray
            wv = @inferred Augmentor.applylazy(FlipY(), Augmentor.prepareaffine(square))
            @test parent(wv).itp.coefs === square
            @test wv == reverse(square; dims=1)
            @test typeof(wv) <: InvWarpedView{eltype(square),2}
        end
        @testset "multiple images" begin
            img_in = (rgb_rect, square)
            res1, res2 = @inferred(Augmentor.applylazy(FlipY(), img_in))
            @test res1 == reverse(rgb_rect; dims=1)
            @test res2 == reverse(square; dims=1)
            @test typeof(res1) <: SubArray{eltype(rgb_rect),2}
            @test typeof(res2) <: SubArray{eltype(square),2}
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(FlipY) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(FlipY) === true
        v = @inferred Augmentor.applylazy(FlipY(), rect)
        @test v === view(rect, 2:-1:1, 1:1:3)
        @test v == reverse(rect; dims=1)
        @test typeof(v) <: SubArray
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(FlipY) === false
    end
end


text:
module TravellingSalesmanProblem

using JuMP, GLPK

export solve_TSP

export hacematrixdedistancias

include("solving_TSP.jl")

end # module


text: print("hello, world!\n")

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct Simulation{Td, Tl}
    name::Symbol
    data::Td
    latent::Tl
end

# Base functions
function Base.convert(::Type{DataFrames.DataFrame}, S::AnimalBehavior.Simulation)
    df = hcat(unpack(S.data), unpack(S.latent))
    return df
end

function Base.show(io::IO, mime::MIME"text/plain", S::AnimalBehavior.Simulation)
    table_conf = set_pt_conf(tf = tf_markdown, alignment = :c)
    println(io, "Simulation of one $(S.name) agent")
    println(io)
    pretty_table_with_conf(table_conf,
        collect(values(S.latent[1]))';
        header=collect(keys(S.latent[1])),
        title="Initial latent variables")

    println(io)
    vals = hcat(collect(StructArrays.components(S.data))..., collect(StructArrays.components(S.latent))...)
    header = vcat(["State", "Action", "Feedback", "Hidden"], keys(S.latent[1])...)
    pretty_table_with_conf(table_conf,
        vals;
        header=header,
        title="Simulation of $(length(S.data)) trials")
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
function mass_matrix_avr_entries!(
    mass_matrix,
    avr::T,
    global_index::Dict{Symbol, Int64},
) where {T <: PSY.AVR}
    @debug "Using default mass matrix entries $T)"
end

function mass_matrix_avr_entries!(
    mass_matrix,
    avr::PSY.SEXS,
    global_index::Dict{Symbol, Int64},
)
    mass_matrix[global_index[:Vf], global_index[:Vf]] = PSY.get_Te(avr)
    mass_matrix[global_index[:Vr], global_index[:Vr]] = PSY.get_Tb(avr)
end

function mdl_avr_ode!(
    device_states,
    output_ode,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRFixed, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #TODO Change EMF name for Vf in PowerSystems
    #Update Vf voltage on inner vars
    get_inner_vars(dynamic_device)[Vf_var] = PSY.get_Vf(PSY.get_avr(dynamic_device))

    return
end

function mdl_avr_ode!(
    device_states,
    output_ode,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRSimple, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AVRSimple)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]

    #Define external states for device
    V_th = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )

    #Get Parameters
    Kv = PSY.get_Kv(PSY.get_avr(dynamic_device))

    #Compute ODEs
    output_ode[local_ix[1]] = Kv * (V_ref - V_th)

    #Update inner_vars
    get_inner_vars(dynamic_device)[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states,
    output_ode,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRTypeI, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AVRTypeI)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Vm = internal_states[4]

    #Define external states for device
    V_th = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )
    Vs = get_inner_vars(dynamic_device)[V_pss_var]

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Ka = PSY.get_Ka(avr)
    Ke = PSY.get_Ke(avr)
    Kf = PSY.get_Kf(avr)
    Ta = PSY.get_Ta(avr)
    Te = PSY.get_Te(avr)
    Tf = PSY.get_Tf(avr)
    Tr = PSY.get_Tr(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)

    #Compute auxiliary parameters
    Se_Vf = Ae * exp(Be * abs(Vf)) #16.13
    V_ref = V0_ref + Vs

    #Set anti-windup for Vr1. #TODO in callbacks
    #if Vr1 > Vr_max
    #    Vr1 = Vr_max
    #elseif Vr1 < Vr_min
    #    Vr1 = Vr_min
    #end

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = -(1.0 / Te) * (Vf * (Ke + Se_Vf) - Vr1) #16.12c
    output_ode[local_ix[2]] = (1.0 / Ta) * (Ka * (V_ref - Vm - Vr2 - (Kf / Tf) * Vf) - Vr1) #16.12a
    output_ode[local_ix[3]] = -(1.0 / Tf) * ((Kf / Tf) * Vf + Vr2) #16.12b
    output_ode[local_ix[4]] = (1.0 / Tr) * (V_th - Vm) #16.11

    #Update inner_vars
    get_inner_vars(dynamic_device)[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states,
    output_ode,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.AVRTypeII, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.AVRTypeII)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Vm = internal_states[4]

    #Define external states for device
    V_th = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )
    Vs = get_inner_vars(dynamic_device)[V_pss_var]

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    K0 = PSY.get_K0(avr)
    T1 = PSY.get_T1(avr)
    T2 = PSY.get_T2(avr)
    T3 = PSY.get_T3(avr)
    T4 = PSY.get_T4(avr)
    Te = PSY.get_Te(avr)
    Tr = PSY.get_Tr(avr)
    Ae = PSY.get_Ae(avr)
    Be = PSY.get_Be(avr)

    #Compute auxiliary parameters
    Se_Vf = Ae * exp(Be * abs(Vf)) #16.13
    V_ref = V0_ref + Vs
    Vr = K0 * Vr2 + (T4 / T3) * (Vr1 + K0 * (T2 / T1) * (V_ref - Vm)) #16.21

    #Set anti-windup for Vr1. #TODO in callbacks
    #if Vr > Vr_max
    #    Vr = Vr_max
    #elseif Vr < Vr_min
    #    Vr = Vr_min
    #end

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = -(1.0 / Te) * (Vf * (1.0 + Se_Vf) - Vr) #16.18
    output_ode[local_ix[2]] = (1.0 / T1) * (K0 * (1.0 - (T2 / T1)) * (V_ref - Vm) - Vr1) #16.14
    output_ode[local_ix[3]] =
        (1.0 / (K0 * T3)) *
        ((1.0 - (T4 / T3)) * (Vr1 + K0 * (T2 / T1) * (V_ref - Vm)) - K0 * Vr2)  #16.20
    output_ode[local_ix[4]] = (1.0 / Tr) * (V_th - Vm) #16.11

    #Update inner_vars
    get_inner_vars(dynamic_device)[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states,
    output_ode,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.ESAC1A, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]

    #Obtain avr
    avr = PSY.get_avr(dynamic_device)

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, typeof(avr))

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vm = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Ve = internal_states[4]
    Vr3 = internal_states[5]

    #Define external states for device
    V_th = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )
    Vs = get_inner_vars(dynamic_device)[V_pss_var]
    Xad_Ifd = get_inner_vars(dynamic_device)[Xad_Ifd_var]

    #Get parameters
    inv_Tr = PSY.get_Tr(avr) < eps() ? 1.0 : 1 / PSY.get_Tr(avr)
    inv_Tb = PSY.get_Tb(avr) < eps() ? 1.0 : 1 / PSY.get_Tb(avr)
    Tc_Tb_ratio = PSY.get_Tb(avr) < eps() ? 0.0 : PSY.get_Tc(avr) * inv_Tb
    Ka = PSY.get_Ka(avr)
    inv_Ta = PSY.get_Ta(avr) < eps() ? 1.0 : 1 / PSY.get_Ta(avr)
    #Va_min, Va_max = PSY.get_Va_lim(avr) #Not used without UEL or OEL
    Te = PSY.get_Te(avr) # Te > 0
    Kf = PSY.get_Kf(avr)
    Tf = PSY.get_Tf(avr) # Te > 0
    Kc = PSY.get_Kc(avr)
    Kd = PSY.get_Kd(avr)
    Ke = PSY.get_Ke(avr)
    E1, E2 = PSY.get_E_sat(avr)
    SE1, SE2 = PSY.get_Se(avr)
    Vr_min, Vr_max = PSY.get_Vr_lim(avr)
    #Obtain saturation
    Se = saturation_function(avr, Ve)

    #Compute auxiliary parameters
    I_N = Kc * Xad_Ifd / Ve
    V_FE = Kd * Xad_Ifd + Ke * Ve + Se * Ve
    V_F = Vr3 + (Kf / Tf) * V_FE
    V_in = V_ref - Vm - V_F
    V_out = Vr1 + (Tc_Tb_ratio) * V_in
    Vf = Ve * rectifier_function(I_N)
    V_R = Vr2

    #Set anti-windup for Vr2.
    if Vr2 > Vr_max
        V_R = Vr_max
    elseif Vr2 < Vr_min
        V_R = Vr_min
    end

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = inv_Tr * (V_th - Vm) #dVm/dt
    output_ode[local_ix[2]] = inv_Tb * (V_in * (1 - Tc_Tb_ratio) - Vr1) #dVr1/dt
    output_ode[local_ix[3]] = inv_Ta * (Ka * V_out - Vr2) #dVr2/dt
    output_ode[local_ix[4]] = (1.0 / Te) * (V_R - V_FE) #dVe/dt
    output_ode[local_ix[5]] = (1.0 / Tf) * (-(Kf / Tf) * V_FE - Vr3) #dVr3/dt

    #Update inner_vars
    get_inner_vars(dynamic_device)[Vf_var] = Vf

    return
end

function mdl_avr_ode!(
    device_states,
    output_ode,
    dynamic_device::PSY.DynamicGenerator{M, S, PSY.SEXS, TG, P},
) where {M <: PSY.Machine, S <: PSY.Shaft, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #Obtain references
    V0_ref = PSY.get_ext(dynamic_device)[CONTROL_REFS][V_ref_index]

    #Obtain indices for component w/r to device
    local_ix = get_local_state_ix(dynamic_device, PSY.SEXS)

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr = internal_states[2]

    #Define external states for device
    V_th = sqrt(
        get_inner_vars(dynamic_device)[VR_gen_var]^2 +
        get_inner_vars(dynamic_device)[VI_gen_var]^2,
    )
    Vs = get_inner_vars(dynamic_device)[V_pss_var]

    #Get parameters
    avr = PSY.get_avr(dynamic_device)
    Ta_Tb = PSY.get_Ta_Tb(avr)
    K = PSY.get_K(avr)
    V_min, V_max = PSY.get_V_lim(avr)

    #Compute auxiliary parameters
    V_in = V0_ref + Vs - V_th
    V_LL = Vr + Ta_Tb * V_in

    #Set anti-windup for Vf. #TODO in callbacks
    #if Vf > V_max
    #    Vf = V_max
    #elseif Vf < V_min
    #    Vf = V_min
    #end

    #Compute 2 States AVR ODE:
    output_ode[local_ix[1]] = K * V_LL - Vf
    output_ode[local_ix[2]] = V_in * (1 - Ta_Tb) - Vr

    #Update inner_vars
    get_inner_vars(dynamic_device)[Vf_var] = Vf

    return
end


text:
# ---
# title: 1444. Number of Ways of Cutting a Pizza
# id: problem1444
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-ways-of-cutting-a-pizza/description/>
# hidden: true
# ---
#
# Given a rectangular pizza represented as a `rows x cols` matrix containing the
# following characters: `'A'` (an apple) and `'.'` (empty cell) and given the
# integer `k`. You have to cut the pizza into `k` pieces using `k-1` cuts.
#
# For each cut you choose the direction: vertical or horizontal, then you choose
# a cut position at the cell boundary and cut the pizza into two pieces. If you
# cut the pizza vertically, give the left part of the pizza to a person. If you
# cut the pizza horizontally, give the upper part of the pizza to a person. Give
# the last piece of pizza to the last person.
#
# _Return the number of ways of cutting the pizza such that each piece contains
# **at least** one apple. _Since the answer can be a huge number, return this
# modulo 10^9 + 7.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/04/23/ways_to_cut_apple_1.png)**
#
#
#
#     Input: pizza = ["A..","AAA","..."], k = 3
#     Output: 3
#     Explanation: The figure above shows the three ways to cut the pizza. Note that pieces must contain at least one apple.
#
#
# **Example 2:**
#
#
#
#     Input: pizza = ["A..","AA.","..."], k = 3
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: pizza = ["A..","A..","..."], k = 1
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= rows, cols <= 50`
#   * `rows == pizza.length`
#   * `cols == pizza[i].length`
#   * `1 <= k <= 10`
#   * `pizza` consists of characters `'A'` and `'.'` only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using DataFrameTools
using Test

function targetlist(targetexpr::Regex, dir2walk::AbstractString)
    # targetexpr = r"(\.jl)$";
    fulllist = String[];
    for (root, folders, files) in walkdir(dir2walk)
        for file in files
            push!(fulllist, joinpath(root, file));
        end
    end
    allfiles = basename.(fulllist);
    targetid = occursin.(targetexpr, allfiles);

    targetpaths = fulllist[targetid];
    return targetpaths, allfiles
end

println("Running tests:")
(testfiles, allnames) = targetlist(r"^(?!runtests).*(\.jl)$", "./");
# pwd() here should be "./test/"
for f in testfiles
   # write test in the file of the same name.
   # E.g., `test/mycode.jl` for testing `src/mycode.jl`.
   include(f);
end


text:
citation(year, vv) = @match my_varname(vv) begin
    "DSi" || "DIP" || "DIN" => citation_Nutrients(year)
    "Temp" => citation_Temperature(year)
    "Salt" => citation_Salinity(year)
    "O2" || "O2sat" || "AOU" => citation_Oxygen(year)
    "Dens" || "Cond" => citation(year)
    _ => error("Not sure what you are trying to cite.")
end

citation(year) = @match my_product_year(year) begin
    "13" => "Boyer, T. P., J. I. Antonov, O. K. Baranova, C. Coleman, H. E. Garcia, A. Grodsky, D. R. Johnson, R. A. Locarnini, A. V. Mishonov, T. D. O'Brien, C. R. Paver, J. R. Reagan, D. Seidov, I. V. Smolyar, and M. M. Zweng, 2013: World Ocean Database 2013, NOAA Atlas NESDIS 72, S. Levitus, Ed., A. Mishonov, Technical Ed.; Silver Spring, MD, 209 pp., doi:10.7289/V5NZ85MT"
    "18" => "Boyer, T. P., O. K. Baranova, C. Coleman, H. E. Garcia, A. Grodsky, R. A. Locarnini, A. V. Mishonov, C. R. Paver, J. R. Reagan, D. Seidov, I. V. Smolyar, K. Weathers, M. M. Zweng,(2018): World Ocean Database 2018. A. V. Mishonov, Technical Editor, NOAA Atlas NESDIS 87."
    _ => "I could not find a citation for WOD$(my_product_year(year)). Add it if you know where to find it!"
end

citation_Temperature(year) = @match my_product_year(year) begin
    "09" => "Locarnini, R. A., A. V. Mishonov, J. I. Antonov, T. P. Boyer, H. E. Garcia, O. K. Baranova, M. M. Zweng, and D. R. Johnson, 2010. World Ocean Atlas 2009, Volume 1: Temperature. S. Levitus, Ed. NOAA Atlas NESDIS 68, U.S. Government Printing Office, Washington, D.C., 184 pp."
    "13" => "Locarnini, R. A., A. V. Mishonov, J. I. Antonov, T. P. Boyer, H. E. Garcia, O. K. Baranova, M. M. Zweng, C. R. Paver, J. R. Reagan, D. R. Johnson, M. Hamilton, and D. Seidov, 2013. World Ocean Atlas 2013, Volume 1: Temperature. S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 73, 40 pp."
    "18" => "Locarnini, R. A., A. V. Mishonov, O. K. Baranova, T. P. Boyer, M. M. Zweng, H. E. Garcia, J. R. Reagan, D. Seidov, K. Weathers, C. R. Paver, and I. Smolyar, 2018. World Ocean Atlas 2018, Volume 1: Temperature. A. Mishonov Technical Ed.; in preparation."
    _ => error("No citation for WOA$(my_product_year(year)) Temperature. Add it if it should be there!")
end

citation_Salinity(year) = @match my_product_year(year) begin
    "09" => "Antonov, J. I., D. Seidov, T. P. Boyer, R. A. Locarnini, A. V. Mishonov, H. E. Garcia, O. K. Baranova, M. M. Zweng, and D. R. Johnson, 2010. World Ocean Atlas 2009, Volume 2: Salinity. S. Levitus, Ed. NOAA Atlas NESDIS 69, U.S. Government Printing Office, Washington, D.C., 184 pp."
    "13" => "Zweng, M. M., J. R. Reagan, J. I. Antonov, R. A. Locarnini, A. V. Mishonov, T. P. Boyer, H. E. Garcia, O. K. Baranova, D. R. Johnson, D. Seidov, M. M. Biddle, 2013. World Ocean Atlas 2013, Volume 2: Salinity. S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 74, 39 pp."
    "18" => "Zweng, M. M., J. R. Reagan, D. Seidov, T. P. Boyer, R. A. Locarnini, H. E. Garcia, A. V. Mishonov, O. K. Baranova, K. Weathers, C. R. Paver, and I. Smolyar, 2018. World Ocean Atlas 2018, Volume 2: Salinity. A. Mishonov Technical Ed.; in preparation."
    _ => error("No citation for WOA$(my_product_year(year)) Salinity. Add it if it should be there!")
end

citation_Oxygen(year) = @match my_product_year(year) begin
    "09" => "Garcia, H. E., R. A. Locarnini, T. P. Boyer, J. I. Antonov, O. K. Baranova, M. M. Zweng, and D. R. Johnson, 2010. World Ocean Atlas 2009, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation. S. Levitus, Ed. NOAA Atlas NESDIS 70, U.S. Government Printing Office, Washington, D.C., 344 pp."
    "13" => "Garcia, H. E., R. A. Locarnini, T. P. Boyer, J. I. Antonov, O. K. Baranova, M. M. Zweng, J. R. Reagan, D. R. Johnson, 2014. World Ocean Atlas 2013, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation. S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 75, 27 pp."
    "18" => "Garcia, H. E., K. Weathers, C. R. Paver, I. Smolyar, T. P. Boyer, R. A. Locarnini, M. M. Zweng, A. V. Mishonov, O. K. Baranova, D. Seidov, and J. R. Reagan, 2018. World Ocean Atlas 2018, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation. A. Mishonov Technical Ed.; in preparation."
    _ => error("No citation for WOA$(my_product_year(year)) AOU. Add it if it should be there!")
end

citation_Nutrients(year) = @match my_product_year(year) begin
    "09" => "Garcia, H. E., R. A. Locarnini, T. P. Boyer, J. I. Antonov, M. M. Zweng, O. K. Baranova, and D. R. Johnson, 2010. World Ocean Atlas 2009, Volume 4: Nutrients (phosphate, nitrate, silicate). S. Levitus, Ed. NOAA Atlas NESDIS 71, U.S. Government Printing Office, Washington, D.C., 398 pp."
    "13" => "Garcia, H. E., R. A. Locarnini, T. P. Boyer, J. I. Antonov, O. K. Baranova, M. M. Zweng, J. R. Reagan, D. R. Johnson, 2014. World Ocean Atlas 2013, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate, silicate). S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 76, 25 pp."
    "18" => "Garcia, H. E., K. Weathers, C. R. Paver, I. Smolyar, T. P. Boyer, R. A. Locarnini, M. M. Zweng, A. V. Mishonov, O. K. Baranova, D. Seidov, and J. R. Reagan, 2018. World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate and nitrate+nitrite, silicate). A. Mishonov Technical Ed.; in preparation."
    _ => error("No citation for WOA$(my_product_year(year)) Nutrients. Add it if it should be there!")
end

# TODO Add other citations for other years!


text:
using Gnuplot
@gp 2:3 "notit"  xrange = (0,1) yrange =(0,1)
@gp :- "set object 1 circle at 0.5,0.5 size 0.2"
# 4 more circles
x, y, r = [0.2,0.8,0.8,0.2], [0.2,0.2,0.8,0.8], 0.13
color = ["#4063D8","#9558B2","#389826","#CB3C33"]
for c in 1:4
    @gp :- "set obj $(c+1) circle at $(x[c]),$(y[c]) size $(r)" :-
    @gp :- "set obj $(c+1) fc '$(color[c])' lw 0" :-
end
@gp
save(term ="pngcairo size 600,400", output ="circles.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# MIT license (c) 2019 by Andrew Lyasoff

# Julia 1.1.0.

# The code illustrates some basic operations with large chunks of market data
# (extracted from public sources), such as
# creating histograms from the returns, and other similar operations.

using Dates, JuliaDB

# The historical quotes are downloaded in .csv format from nasdaq.com and are placed
# in the directory "HistoricalQuotes." Before those files can be used here they
# must be modified (slightly) in Excel, LibreOffice Calc, or a text editor (e.g., Emacs):
# there should be no empty line, or a line that contains a time-stamp (all lines except
# the first one must be identically formatted -- this may involve removing the second
# line in the spreadhseet). The data for all stocks can be stored in a single variable
# (note that HistoricalQuotes is the name of the sub-directory that contains all .csv files).

stocksdata = loadndsparse("HistoricalQuotes"; filenamecol = :ticker, indexcols = [:ticker, :date])

# The same database can now be saved in a special binary format that makes reloading
# at a later time very fast.
save(stocksdata, "stocksdata.jdb")
@time reloaded_stocksdata = load("stocksdata.jdb")

# Test the saved database for consistency.
stocksdata == reloaded_stocksdata

# Look up the data associated only with Apple (symbol AAPL):
stocksdata["AAPL",:]
length(stocksdata["AAPL",:])

# Similarly, look up the data linked to Microsoft.
stocksdata["MSFT",:]
length(stocksdata["MSFT",:])

# Look up the price of Google on a specific date in the past.
stocksdata["GOOGL", Date(2009,9,9)]
stocksdata["GOOGL", Date(2009,9,9)].close

# Extract the closing prices of Amazon *only*.
selectvalues(stocksdata,:close)["AMZN",:]
keytype(stocksdata)

# List all dates for which the closing prices for Apple are available.
AAPLdate=columns(stocksdata["AAPL",:])[1]

#Extract the closing prices of Apple on those days.
AAPLclose=columns(stocksdata["AAPL",:])[2]

stocksdata["AAPL", Date(2008,8,11)].close

# Produce some plots.
using Plots
pyplot()
plot(AAPLdate, AAPLclose,label="AAPL closing price")

# Calculate and plot the daily returns for AAPL.
begin
    llc=length(AAPLclose);
    AAPLreturns=(AAPLclose[2:llc]-AAPLclose[1:llc-1])./AAPLclose[1:llc-1];
    scatter(AAPLdate[2:llc],AAPLreturns,label="AAPL daily returns",markersize=2)
end

# Now build the histogram from the returns using a custom made 'histogram' function
# (called 'hstgram' to avoid the confucion with the standard 'histogram').
# It takes as an input a single 1-dimensional array of data. The number of bins in
# the histogram is determined automatically by using the Diaconis-Friedman rule.
# The function returns two arrays: the mid-points of the bins and the (unnormalized)
# heights of the bars.

using StatsBase

function hstgram(data_sample::Array{Float64,1})
    data_sorted=sort(data_sample)
    first=data_sorted[1]
    last=data_sorted[end]
    nmb=length(data_sorted)
    IQR=percentile(data_sorted,75)-percentile(data_sorted,25)
    bin_size_loc = 2*IQR*(nmb^(-1.0/3))
    num_bins=Int(floor((last-first)/bin_size_loc))
    bin_size=(last-first)/(num_bins)
    bin_end_points=[first+(i-1)*bin_size for i=1:(num_bins+1)]
    ahist_val=[length(data_sorted[data_sorted .< u]) for u in bin_end_points]
    hist_val=[ahist_val[i+1]-ahist_val[i] for i=1:num_bins]
    mid_bins=[first-bin_size/2+i*bin_size for i=1:num_bins]
    return mid_bins, hist_val
end

# Normalize the bars so that the area of the histogram equals 1.
begin
    U,V=hstgram(AAPLreturns);
    VV=V/(sum(V)*(U[2]-U[1]));
end

begin
    plot(U,VV,line=(:sticks,0.75),label="")
    xlabel!("returns")
    ylabel!("normalized frequency")
    title!("Histogram from the 10y daily returns from AAPL.")
end

# Test that the area of the histogram is indeed 1.
sum(VV)*(U[2]-U[1])

# Another variation of the same histogram.
begin
    plot(U.+(U[2]-U[1])/2,VV,label="",line=(:steppre,1),linewidth=0.05)
    xlabel!("samples")
    ylabel!("normalized frequency")
end

# Normalize the bars to give the probabilities for hitting the bins.
begin
    VVV=V/sum(V);
    plot(U.+(U[2]-U[1])/2,VVV,label="",line=(:steppre,1),linewidth=0.05)
    xlabel!("samples")
    ylabel!("probability")
end

# Test that the probabilities do sum to 1.
sum(VVV)

max_stars_repo_name: AndrewLyasoff/SMAP

max_issues_repo_name: lhyzh/SMAP

max_forks_repo_name: lhyzh/SMAP


text:
using SampleJuliaProj
using Documenter

DocMeta.setdocmeta!(SampleJuliaProj, :DocTestSetup, :(using SampleJuliaProj); recursive=true)

makedocs(;
    modules=[SampleJuliaProj],
    authors="Nivedita Rethnakar et al.",
    repo="https://github.com/nivupai/SampleJuliaProj.jl/blob/{commit}{path}#{line}",
    sitename=" Documentation for SampleJuliaProj.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nivupai.github.io/SampleJuliaProj.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nivupai/SampleJuliaProj.jl",
    devbranch="main",
)


text:
####################################################
# Script to train a patchy SC network with 1 hidden layer (SGD)
# patch connectivity is controled via "n_patches" & "n_hidden_per_patch" variables
# Algorithm for sparse coding as proposed in:
#   Brito CSN, GerstnerW(2016) PLoS Comput Biol 12(9).
# OR (depending on "Földiak_model" boolian)
#   Földiak 1990 Biol. Cybernetics

###################################################
# Include libraries and core

using Pkg; Pkg.activate("./../core/"); Pkg.instantiate()
using ProgressMeter, LinearAlgebra, Statistics, Distributions
include("./../core/src/autoencoders/autoencoders.jl")

###################################################
# Define task and network

dataset = "MNIST" # dataset to be used: "MNIST" or "CIFAR10"
colorimage = dataset == "CIFAR10"
Földiak_model = false

patch_size = 10 # (linear) patch size of receptive field in first layer p ∈ [1,28] or [1,32] for MNIST/CIFAR10 resp.
n_patches = 100 # number of feature extraction patches
n_hidden_per_patch = 10 # hidden neurons/features per patch. Overall number of hidden neurons is: n_patches * n_hidden_per_patch

iterations_SC = 10^5 # iterations for SC learning
n_inits_c = 3 # number of initializations for classifier
iterations_c = 10^6 # iterations for classifier learning
learningrate_c = 1e-3*[1] # learningrate for classifier

###################################################
# Load and preprocess data

smallimgs, labels, smallimgstest, labelstest = import_data(dataset)
subtractmean!(smallimgs) # substract pixel-wise mean
subtractmean!(smallimgstest)

############################################################################
# Set up network for sparse coding layer

layers = [[size(smallimgs)[1],n_hidden_per_patch]]
nAE = size(layers, 1)[1]
learningrates = [1e-3*[1] for i in 1:nAE] # learningrates for SC learning
lambdas = [[1]*1e-2 for i in 1:nAE] # sparsity parameters for SC learning

nonlinearity, nonlinearity_diff = init_nonlin([lin!],[lin_diff!]) # nonlinearities
nonlinearity = [nonlinearity for i in 1:nAE]
nonlinearity_diff = [nonlinearity_diff for i in 1:nAE]
confSAE = configSAE_sparse("_",[size(smallimgs, 2),size(smallimgstest, 2)],
            nAE, layers, 1, convert(Array{Int64, 1}, [iterations_SC for i in 1:nAE]), learningrates, lambdas, 0. ,
            "none", [0. for i in 1:nAE], nonlinearity, nonlinearity_diff, true)

############################################################################
# Set up classifier

hidden_sizes_c = [0] # shallow classifier (no hidden layer)
confClass = configClass("_", hidden_sizes_c,n_inits_c,iterations_c,
        learningrate_c, 0., [relu!], [relu_diff!], true)
confSAEc = configSAEc_sparse(confSAE,confClass)

############################################################################
# Train network

if colorimage
    im_size = Int(sqrt(size(smallimgs)[1]/3))
else
    im_size = Int(sqrt(size(smallimgs)[1]))
end
sae, errors = trainSAEc_full_patchy!(smallimgs, smallimgstest, labels, labelstest, n_hidden_per_patch, confSAEc;
                                        patchsize = patch_size, n_of_populations = n_patches,
                                        stride = 0, im_size = im_size,
                                        Földiak_model = Földiak_model,
                                        colorimage = colorimage)


text:
__precompile__()

module PortAudio

using SampledSignals
using RingBuffers
#using Suppressor

import Base: eltype, show
import Base: close, isopen
import Base: read, read!, write, flush

# Get binary dependencies loaded from BinDeps
include("../deps/deps.jl")
include("suppressor.jl")
include("pa_shim.jl")
include("libportaudio.jl")

function __init__()
    init_pa_shim()
    global const notifycb_c = cfunction(notifycb, Cint, (Ptr{Void}, ))
    # initialize PortAudio on module load
    @suppress_err Pa_Initialize()
end

export PortAudioStream

# These sizes are all in frames

# the block size is what we request from portaudio if no blocksize is given.
# The ringbuffer and pre-fill will be twice the blocksize
const DEFAULT_BLOCKSIZE=4096

# data is passed to and from the ringbuffer in chunks with this many frames
# it should be at most the ringbuffer size, and must evenly divide into the
# the underlying portaudio buffer size. E.g. if PortAudio is running with a
# 2048-frame buffer period, the chunk size can be 2048, 1024, 512, 256, etc.
const CHUNKSIZE=128

# ringbuffer to receive errors from the audio processing thread
const ERR_BUFSIZE=512

function versioninfo(io::IO=STDOUT)
    println(io, Pa_GetVersionText())
    println(io, "Version: ", Pa_GetVersion())
    println(io, "Shim Source Hash: ", shimhash()[1:10])
end

type PortAudioDevice
    name::String
    hostapi::String
    maxinchans::Int
    maxoutchans::Int
    defaultsamplerate::Float64
    idx::PaDeviceIndex
end

PortAudioDevice(info::PaDeviceInfo, idx) = PortAudioDevice(
        unsafe_string(info.name),
        unsafe_string(Pa_GetHostApiInfo(info.host_api).name),
        info.max_input_channels,
        info.max_output_channels,
        info.default_sample_rate,
        idx)

function devices()
    ndevices = Pa_GetDeviceCount()
    infos = PaDeviceInfo[Pa_GetDeviceInfo(i) for i in 0:(ndevices - 1)]
    PortAudioDevice[PortAudioDevice(info, idx-1) for (idx, info) in enumerate(infos)]
end

# not for external use, used in error message printing
devnames() = join(["\"$(dev.name)\"" for dev in devices()], "\n")

##################
# PortAudioStream
##################

type PortAudioStream{T}
    samplerate::Float64
    blocksize::Int
    stream::PaStream
    sink # untyped because of circular type definition
    source # untyped because of circular type definition
    errbuf::RingBuffer{pa_shim_errmsg_t} # used to send errors from the portaudio callback
    bufinfo::pa_shim_info_t # data used in the portaudio callback

    # this inner constructor is generally called via the top-level outer
    # constructor below
    function PortAudioStream{T}(indev::PortAudioDevice, outdev::PortAudioDevice,
                                inchans, outchans, sr, blocksize, synced) where {T}
        inchans = inchans == -1 ? indev.maxinchans : inchans
        outchans = outchans == -1 ? outdev.maxoutchans : outchans
        inparams = (inchans == 0) ?
            Ptr{Pa_StreamParameters}(0) :
            Ref(Pa_StreamParameters(indev.idx, inchans, type_to_fmt[T], 0.0, C_NULL))
        outparams = (outchans == 0) ?
            Ptr{Pa_StreamParameters}(0) :
            Ref(Pa_StreamParameters(outdev.idx, outchans, type_to_fmt[T], 0.0, C_NULL))
        this = new(sr, blocksize, C_NULL)
        finalizer(this, close)
        this.sink = PortAudioSink{T}(outdev.name, this, outchans, blocksize*2)
        this.source = PortAudioSource{T}(indev.name, this, inchans, blocksize*2)
        this.errbuf = RingBuffer{pa_shim_errmsg_t}(1, ERR_BUFSIZE)
        if synced && inchans > 0 && outchans > 0
            # we've got a synchronized duplex stream. initialize with the output buffer full
            write(this.sink, SampleBuf(zeros(T, blocksize*2, outchans), sr))
        end
        # pass NULL for input/output we're not using
        this.bufinfo = pa_shim_info_t(
                inchans > 0 ? bufpointer(this.source) : C_NULL,
                outchans > 0 ? bufpointer(this.sink) : C_NULL,
                pointer(this.errbuf),
                synced, notifycb_c,
                inchans > 0 ? notifyhandle(this.source) : C_NULL,
                outchans > 0 ? notifyhandle(this.sink) : C_NULL,
                notifyhandle(this.errbuf))
        this.stream = @suppress_err Pa_OpenStream(inparams, outparams,
                                                  float(sr), blocksize,
                                                  paNoFlag, shim_processcb_c,
                                                  this.bufinfo)

        Pa_StartStream(this.stream)
        @async handle_errors(this)

        this
    end
end

"""
    PortAudioStream(inchannels=2, outchannels=2; options...)
    PortAudioStream(duplexdevice, inchannels=2, outchannels=2; options...)
    PortAudioStream(indevice, outdevice, inchannels=2, outchannels=2; options...)

Audio devices can either be `PortAudioDevice` instances as returned
by `PortAudio.devices()`, or strings with the device name as reported by the
operating system. If a single `duplexdevice` is given it will be used for both
input and output. If no devices are given the system default devices will be
used.

Options:

* `eltype`:       Sample type of the audio stream (defaults to Float32)
* `samplerate`:   Sample rate (defaults to device sample rate)
* `blocksize`:    Size of the blocks that are written to and read from the audio
                  device. (Defaults to $DEFAULT_BLOCKSIZE)
* `synced`:       Determines whether the input and output streams are kept in
                  sync. If `true`, you must read and write an equal number of
                  frames, and the round-trip latency is guaranteed constant. If
                  `false`, you are free to read and write separately, but
                  overflow or underflow can affect the round-trip latency.
"""
# this is the top-level outer constructor that all the other outer constructors
# end up calling
function PortAudioStream(indev::PortAudioDevice, outdev::PortAudioDevice,
        inchans=2, outchans=2; eltype=Float32, samplerate=-1, blocksize=DEFAULT_BLOCKSIZE, synced=false)
    if samplerate == -1
        sampleratein = indev.defaultsamplerate
        samplerateout = outdev.defaultsamplerate
        if inchans > 0 && outchans > 0 && sampleratein != samplerateout
            error("""
            Can't open duplex stream with mismatched samplerates (in: $sampleratein, out: $samplerateout).
                   Try changing your sample rate in your driver settings or open separate input and output
                   streams""")
        elseif inchans > 0
            samplerate = sampleratein
        else
            samplerate = samplerateout
        end
    end
    PortAudioStream{eltype}(indev, outdev, inchans, outchans, samplerate, blocksize, synced)
end

# handle device names given as streams
function PortAudioStream(indevname::AbstractString, outdevname::AbstractString, args...; kwargs...)
    indev = nothing
    outdev = nothing
    for d in devices()
        if d.name == indevname
            indev = d
        end
        if d.name == outdevname
            outdev = d
        end
    end
    if indev == nothing
        error("No device matching \"$indevname\" found.\nAvailable Devices:\n$(devnames())")
    end
    if outdev == nothing
        error("No device matching \"$outdevname\" found.\nAvailable Devices:\n$(devnames())")
    end

    PortAudioStream(indev, outdev, args...; kwargs...)
end

# if one device is given, use it for input and output, but set inchans=0 so we
# end up with an output-only stream
function PortAudioStream(device::PortAudioDevice, inchans=2, outchans=2; kwargs...)
    PortAudioStream(device, device, inchans, outchans; kwargs...)
end
function PortAudioStream(device::AbstractString, inchans=2, outchans=2; kwargs...)
    PortAudioStream(device, device, inchans, outchans; kwargs...)
end

# use the default input and output devices
function PortAudioStream(inchans=2, outchans=2; kwargs...)
    inidx = Pa_GetDefaultInputDevice()
    indevice = PortAudioDevice(Pa_GetDeviceInfo(inidx), inidx)
    outidx = Pa_GetDefaultOutputDevice()
    outdevice = PortAudioDevice(Pa_GetDeviceInfo(outidx), outidx)
    PortAudioStream(indevice, outdevice, inchans, outchans; kwargs...)
end

function close(stream::PortAudioStream)
    if stream.stream != C_NULL
        Pa_StopStream(stream.stream)
        Pa_CloseStream(stream.stream)
        close(stream.source)
        close(stream.sink)
        stream.stream = C_NULL
    end

    nothing
end

isopen(stream::PortAudioStream) = stream.stream != C_NULL

SampledSignals.samplerate(stream::PortAudioStream) = stream.samplerate
eltype{T}(stream::PortAudioStream{T}) = T

read(stream::PortAudioStream, args...) = read(stream.source, args...)
read!(stream::PortAudioStream, args...) = read!(stream.source, args...)
write(stream::PortAudioStream, args...) = write(stream.sink, args...)
write(sink::PortAudioStream, source::PortAudioStream, args...) = write(sink.sink, source.source, args...)
flush(stream::PortAudioStream) = flush(stream.sink)

function show(io::IO, stream::PortAudioStream)
    println(io, typeof(stream))
    println(io, "  Samplerate: ", samplerate(stream), "Hz")
    print(io, "  Buffer Size: ", stream.blocksize, " frames")
    if nchannels(stream.sink) > 0
        print(io, "\n  ", nchannels(stream.sink), " channel sink: \"", name(stream.sink), "\"")
    end
    if nchannels(stream.source) > 0
        print(io, "\n  ", nchannels(stream.source), " channel source: \"", name(stream.source), "\"")
    end
end

"""
    handle_errors(stream::PortAudioStream)

Handle errors coming over the error stream from PortAudio. This is run as an
independent task while the stream is active.
"""
function handle_errors(stream::PortAudioStream)
    err = Vector{pa_shim_errmsg_t}(1)
    while true
        nread = read!(stream.errbuf, err)
        nread == 1 || break
        if err[1] == PA_SHIM_ERRMSG_ERR_OVERFLOW
            warn("Error buffer overflowed on stream $(stream.name)")
        elseif err[1] == PA_SHIM_ERRMSG_OVERFLOW
            # warn("Input overflowed from $(name(stream.source))")
        elseif err[1] == PA_SHIM_ERRMSG_UNDERFLOW
            # warn("Output underflowed to $(name(stream.sink))")
        else
            error("""
                Got unrecognized error code $(err[1]) from audio thread for
                stream "$(stream.name)". Please file an issue at
                https://github.com/juliaaudio/portaudio.jl/issues""")
        end
    end
end

##################################
# PortAudioSink & PortAudioSource
##################################

# Define our source and sink types
for (TypeName, Super) in ((:PortAudioSink, :SampleSink),
                          (:PortAudioSource, :SampleSource))
    @eval type $TypeName{T} <: $Super
        name::String
        stream::PortAudioStream{T}
        chunkbuf::Array{T, 2}
        ringbuf::RingBuffer{T}
        nchannels::Int

        function $TypeName{T}(name, stream, channels, ringbufsize) where {T}
            # portaudio data comes in interleaved, so we'll end up transposing
            # it back and forth to julia column-major
            chunkbuf = zeros(T, channels, CHUNKSIZE)
            ringbuf = RingBuffer{T}(channels, ringbufsize)
            new(name, stream, chunkbuf, ringbuf, channels)
        end
    end
end

SampledSignals.nchannels(s::Union{PortAudioSink, PortAudioSource}) = s.nchannels
SampledSignals.samplerate(s::Union{PortAudioSink, PortAudioSource}) = samplerate(s.stream)
SampledSignals.blocksize(s::Union{PortAudioSink, PortAudioSource}) = s.stream.blocksize
eltype(::Union{PortAudioSink{T}, PortAudioSource{T}}) where {T} = T
close(s::Union{PortAudioSink, PortAudioSource}) = close(s.ringbuf)
isopen(s::Union{PortAudioSink, PortAudioSource}) = isopen(s.ringbuf)
RingBuffers.notifyhandle(s::Union{PortAudioSink, PortAudioSource}) = notifyhandle(s.ringbuf)
bufpointer(s::Union{PortAudioSink, PortAudioSource}) = pointer(s.ringbuf)
name(s::Union{PortAudioSink, PortAudioSource}) = s.name

function show(io::IO, stream::T) where {T <: Union{PortAudioSink, PortAudioSource}}
    println(io, T, "(\"", stream.name, "\")")
    print(io, nchannels(stream), " channels")
end

flush(sink::PortAudioSink) = flush(sink.ringbuf)

function SampledSignals.unsafe_write(sink::PortAudioSink, buf::Array, frameoffset, framecount)
    nwritten = 0
    while nwritten < framecount
        towrite = min(framecount-nwritten, CHUNKSIZE)
        # make a buffer of interleaved samples
        transpose!(view(sink.chunkbuf, :, 1:towrite),
                   view(buf, (1:towrite)+nwritten+frameoffset, :))
        n = write(sink.ringbuf, sink.chunkbuf, towrite)
        nwritten += n
        # break early if the stream is closed
        n < towrite && break
    end

    nwritten
end

function SampledSignals.unsafe_read!(source::PortAudioSource, buf::Array, frameoffset, framecount)
    nread = 0
    while nread < framecount
        toread = min(framecount-nread, CHUNKSIZE)
        n = read!(source.ringbuf, source.chunkbuf, toread)
        # de-interleave the samples
        transpose!(view(buf, (1:toread)+nread+frameoffset, :),
                   view(source.chunkbuf, :, 1:toread))

        nread += toread
        # break early if the stream is closed
        n < toread && break
    end

    nread
end

# this is called by the shim process callback to notify that there is new data.
# it's run in the audio context so don't do anything besides wake up the
# AsyncCondition handle associated with that ring buffer
notifycb(handle) = ccall(:uv_async_send, Cint, (Ptr{Void}, ), handle)

end # module PortAudio


text:
mutable struct FeaturePlot
    xmin::Float64
    ymin::Float64
    xscale::Float64
    yscale::Float64
    xc_i::Float64
end

FeaturePlot()=FeaturePlot(0.0,0.0,1.0,1.0,1.0)

mutable struct Buffer
    count::Int64
    ind::Int64
    spikes::Array{Int16,2}
    clus::Array{UInt8,1}
    mask::Array{Bool,1}
    selected_clus::UInt8
    replot::Bool
    selected::Array{Bool,1}
    c_changed::Bool
end

Buffer(wave_points)=Buffer(500,1,zeros(Int16,wave_points,500),zeros(UInt8,500),trues(500),1,false,falses(500),false)

struct Vec2
    x::Float64
    y::Float64
end

mutable struct RubberBand
    pos0::Vec2
    pos1::Vec2
    pos2::Vec2
    polygon::Array{Vec2,1}
    moved::Bool
    minpixels::Int
end

mutable struct SortView
    win::Gtk.GtkWindowLeaf

    c::Gtk.GtkCanvasLeaf

    b1::Gtk.GtkButtonLeaf

    features::Dict{String,Array{Float64,1}}

    pca::PCA{Float64}
    pca_calced::Bool

    n_col::Int64
    n_row::Int64

    popup_axis::Gtk.GtkMenuLeaf

    selected_plot::Int64

    axes::Array{Bool,2}
    axes_name::Array{String,2}

    col_sb::Gtk.GtkSpinButtonLeaf

    plots::Array{FeaturePlot,1}

    h::Float64
    w::Float64

    buf::Buffer
end

mutable struct Thres_Widgets
    sb::Gtk.GtkLabelLeaf
    slider::Gtk.GtkScaleLeaf
    adj::Gtk.GtkAdjustmentLeaf
    all::Gtk.GtkCheckButtonLeaf
    show::Gtk.GtkCheckButtonLeaf
end

mutable struct Gain_Widgets
    gainbox::Gtk.GtkSpinButtonLeaf
    multiply::Gtk.GtkCheckButtonLeaf
    all::Gtk.GtkCheckButtonLeaf
end

mutable struct Single_Channel
    c2::Gtk.GtkCanvasLeaf
    c3::Gtk.GtkCanvasLeaf
    ctx2::Cairo.CairoContext
    ctx2s::Cairo.CairoContext
    rb_active::Bool
    rb::RubberBand
    click_button::Int64
    selected::Array{Bool,1}
    plotted::Array{Bool,1}
    hold::Bool
    pause::Bool
    pause_button::Gtk.GtkToggleButtonLeaf
    rb_buttons::Array{Gtk.GtkRadioButton,1}
    pause_state::Int64
    mi::NTuple{2,Float64} #saved x,y position of mouse input
    show_thres::Bool
    w2::Int64
    h2::Int64
    wave_points::Int64
    s::Float64
    o::Float64
    buf::Buffer
    thres::Float64
    old_thres::Float64
    temp::ClusterTemplate

    total_clus::Int64
    spike::Int64
    sort_cb::Bool
    sort_list::Gtk.GtkListStoreLeaf
    sort_tv::Gtk.GtkTreeViewLeaf
    adj_sort::Gtk.GtkAdjustmentLeaf

    adj_thres::Gtk.GtkAdjustmentLeaf
    thres_slider::Gtk.GtkScaleLeaf
    thres_changed::Bool
    thres_widgets::Thres_Widgets
    gain_widgets::Gain_Widgets
end


text:
using DelayEmbeddings
using TransferEntropy

# Estimate transfer entropy from scratch from a random
# set of points `n_realizations` times. If everything
# works as expected, there should be no negative transfer
# entropy values.
ts_length = 100

@testset "3D #$i" for i in 1:n_realizations
	x = rand(100)
	y = rand(100)
	D = Dataset(x, y)
	pts = customembed(D, Positions(2, 2, 1), Lags(1, 0, 0))
	v = TEVars(Tf = [1], Tpp = [2], Spp = [3])
	binning_scheme = RectangularBinning([2, 4, 5])

	@test transferentropy_visitfreq(pts, binning_scheme, v, b = 2) >= 0
end

@testset "4D #$i" for i in 1:n_realizations
	x = rand(100)
	y = rand(100)
	D = Dataset(x, y)
	pts = customembed(D, Positions(2, 2, 2, 1), Lags(1, 0, -2, 0))
	v = TEVars(Tf = [1], Tpp = [2, 3], Spp = [4])
	binning_scheme = RectangularBinning([2, 4, 5, 5])

	@test transferentropy_visitfreq(pts, binning_scheme, v, b = 2) >= 0
end

@testset "5D #$i" for i in 1:n_realizations
	x = rand(100)
	y = rand(100)
	D = Dataset(x, y)
	pts = customembed(D, Positions(2, 2, 2, 2, 1), Lags(1, 0, -2, -5, 0))
	v = TEVars(Tf = [1], Tpp = [2, 3, 4], Spp = [5])
	binning_scheme = RectangularBinning([2, 4, 5, 5, 3])

	@test transferentropy_visitfreq(pts, binning_scheme, v, b = 2) >= 0
end


text:
using DashBootstrapComponents, DashHtmlComponents

row = html_div([
    dbc_row(
        [
            dbc_col(html_div("One of three columns")),
            dbc_col(html_div("One of three columns")),
            dbc_col(html_div("One of three columns")),
        ],
        align="start",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of three columns")),
            dbc_col(html_div("One of three columns")),
            dbc_col(html_div("One of three columns")),
        ],
        align="center",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of three columns")),
            dbc_col(html_div("One of three columns")),
            dbc_col(html_div("One of three columns")),
        ],
        align="end",
    ),
    dbc_row([
        dbc_col(html_div("One of three columns"), align="start"),
        dbc_col(html_div("One of three columns"), align="center"),
        dbc_col(html_div("One of three columns"), align="end"),
    ]),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
"""
    SweepDetails

Details of a single DMRG sweep.
"""
struct SweepDetails
    "Settings for the sweep."
    set::SweepSettings

    "Total energy of the state."
    energy::Float64
    "Normalized fluctuation of the energy: `<(H - <H>)^2>/<H>^2`."
    dH2::Float64
    "Eigenvalues across the middle bond."
    middle_eigvals::Vector{Float64}

    "Maximum wavefunction bond dimension."
    mps_max_bond_dim::Int

    "Duration of the sweep in seconds."
    duration::Float64
end

"""
    SweepDetails(set::SweepSettings, state::SweepState)

Create a `SweepDetails` for settings `set` using the values in `state`.
"""
SweepDetails(set::SweepSettings, state::SweepState) =
    SweepDetails(set, state.energy, state.dH2, state.middle_eigvals,
                 state.mps_max_bond_dim, state.duration)

"""
    SweepHistory

Information about a sequence DMRG sweeps.
"""
struct SweepHistory
    "Details from all completed sweeps."
    details::Vector{SweepDetails}
    "Whether the state is converged to the desired tolerance."
    converged::Bool
end

Base.getindex(hist::SweepHistory, i::Int) = hist.details[i]
Base.lastindex(hist::SweepHistory) = length(hist.details)
Base.length(hist::SweepHistory) = length(hist.details)


text:
using DataDepsGenerators
using Test

using ReferenceTests

@testset "DataOneV1 test" begin
    @test_reference "references/DataOneV1 Ecology.txt" remove_cite_version(generate(DataOneV1(), "https://datadryad.org/mn/object/http://dx.doi.org/10.5061/dryad.74699"))

    #Checking files with multiple files available for download
    @test_reference "references/DataOneV1 Plasticity.txt"  remove_cite_version(generate(DataOneV1(), "https://datadryad.org/mn/object/http://dx.doi.org/10.5061/dryad.f9s4424"))
    @test_reference "references/DataOneV1 Drought.txt"  remove_cite_version(generate(DataOneV1(), "https://datadryad.org/mn/object/http://dx.doi.org/10.5061/dryad.cc8834s"))
end


text:
import QuantumCircuitOpt as QCOpt
using JuMP
using Gurobi
# using Cbc
# using CPLEX

include("optimizers.jl")
include("2qubit_gates.jl")
include("3qubit_gates.jl")
include("4qubit_gates.jl")
include("5qubit_gates.jl")

decompose_gates = ["decompose_hadamard",
                   "decompose_controlled_Z",
                   "decompose_controlled_V",
                   "decompose_controlled_H",
                   "decompose_controlled_H_with_R",
                   "decompose_magic",
                   "decompose_magic_using_SHCnot",
                   "decompose_S",
                   "decompose_revcnot",
                   "decompose_revcnot_with_U",
                   "decompose_swap",
                   "decompose_W",
                   "decompose_W_using_HCnot",
                   "decompose_GroverDiffusion_using_Clifford",
                   "decompose_GroverDiffusion_using_U3",
                   "decompose_iSwap",
                   "decompose_qft2_using_HT",
                   "decompose_RX_on_q3"]

# decompose_gates = ["decompose_quantum_fulladder"]

#----------------------------------------------#
#      Quantum Circuit Optimization model      #
#----------------------------------------------#
result = Dict{String,Any}()

for gates in decompose_gates

    params = getfield(Main, Symbol(gates))()

    model_options = Dict{Symbol, Any}(:convex_hull_gate_constraints => false,
                                      :unit_magnitude_constraints => false,
                                      :idempotent_gate_constraints => false)

    qcopt_optimizer = get_gurobi()
    global result = QCOpt.run_QCModel(params, qcopt_optimizer; options = model_options)
end


text:
#================================================================

Definition of MinHash, an LSH function for hashing on Jaccard similarity.

================================================================#

using Random: shuffle!

#========================
MinHash struct definition and constructors
========================#

struct MinHash{T, I <: Union{UInt32,UInt64}} <: SymmetricLSHFunction
    # Flag that indicates whether or not we were provided all of the symbols
    # we will possibly see in the MinHash constructor. This is used during
    # hashing to improve efficiency in the case where we were already provided
    # all of the symbols.
    fixed_symbols :: Bool

    # A list of dictionaries mapping symbols to integers. To compute a single
    # hash function, we pass every element of an input set through one of
    # the dictionaries and use the smallest output as our hash.
    mappings :: Vector{Dict{T,I}}
end

"""
    MinHash(n_hashes::Integer = $(DEFAULT_N_HASHES);
            dtype::Type = Any,
            symbols::Union{Vector,Set} = Set())

Construct a locality-sensitive hash function for Jaccard similarity.

# Arguments
- $(N_HASHES_DOCSTR())

# Keyword parameters
- `dtype::Type` (default: `Any`): the type of symbols in the sets you're hashing. This is overriden by the data type contained in `symbols` when `symbols` is non-empty.
- `symbols::Union{Vector,Set}`: a `Vector` or `Set` containing all of the possible elements ("symbols") of the sets that you will be hashing. If left empty, `MinHash` will instead expand its dictionary when it sees new symbols (at small additional computational expense).

# Examples
Construct a hash function to hash sets whose elements are integers between `1` and `50`:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash(40; symbols = Set(1:50));

julia> n_hashes(hashfn) == 40 && similarity(hashfn) == jaccard
true

julia> hashfn(Set([38, 14, 29, 48, 11]));

julia> hashfn([1, 1, 2, 3, 4]); # You can also hash Vectors

julia> hashfn(Set([100]))
ERROR: Symbol 100 not found
```

If you aren't sure ahead of time exactly what kinds of elements will be in the sets you're hashing, you can opt not to specify `symbols`, in which case `MinHash` will lazily update its hash functions as it encounters new symbols:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash();

julia> hashfn(Set([1, 2, 3]));

julia> hashfn(Set(["a", "b", "c"]));

```

If you don't know what elements you'll encounter, but you know that they'll all be of a specific data type, you can specify the `dtype` argument for increased efficiency:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash(10; dtype = String);

julia> hashfn(Set(["a", "b", "c"]));

```

# References
- Broder, A. *On the resemblance and containment of documents*. Compression and Complexity of Sequences: Proceedings, Positano, Amalfitan Coast, Salerno, Italy, June 11-13, 1997. doi:10.1109/SEQUEN.1997.666900.

See also: [`jaccard`](@ref)
"""
function MinHash(args...;
                 dtype::Type = Any,
                 symbols::C = Set{Any}()) where {T, C <: Union{Vector{T},Set{T}}}

    if length(symbols) > 0
        MinHash{T}(args...; symbols=symbols)
    else
        MinHash{dtype}(args...; symbols=Set{dtype}())
    end
end

function MinHash{T}(n_hashes::Integer = DEFAULT_N_HASHES;
                    symbols::C = Set{T}()) where {T, C <: Union{Vector{<:T},Set{<:T}}}

    fixed_symbols = (length(symbols) > 0)

    # If fixed_symbols is true, and the symbol set is sufficiently small, then
    # we make our mappings from symbols -> integers map to UInt32 rather than
    # UInt64.
    # Note that we shouldn't do this if the symbol set is unspecified, because
    # in that case new mappings are determined randomly. With only 32  bits of
    # randomness, the probability that two symbols map to the same integer gets
    # high relatively quickly.
    hash_type =
        (fixed_symbols && length(symbols) ≤ typemax(UInt32)) ?
        UInt32 : UInt64

    mappings = Vector{Dict{T, hash_type}}(undef, n_hashes)

    if !fixed_symbols
        for ii = 1:length(mappings)
            mappings[ii] = Dict{T,hash_type}()
        end
    else
        # Create a new random mapping (i.e. a new Dict) that maps symbols to
        # integers in the range 1:length(symbols).
        mapping_range = convert.(hash_type, 1:length(symbols))

        for ii = 1:length(mappings)
            new_mapping = Dict{T,hash_type}()
            shuffle!(mapping_range)

            for (sym,h) in zip(symbols, mapping_range)
                new_mapping[sym] = h
            end

            mappings[ii] = new_mapping
        end
    end

    MinHash{T, hash_type}(fixed_symbols, mappings)
end

#========================
LSHFunction and SymmetricLSHFunction API compliance
========================#
n_hashes(hashfn::MinHash) = length(hashfn.mappings)
hashtype(::MinHash{T, I}) where {T, I} = I
similarity(::MinHash) = jaccard

single_hash_collision_probability(::MinHash, sim) = sim

### Hash computation

function (hashfn :: MinHash{T,I})(x) where {T,I}
    # The iith hash is the smallest output of hashfn.mappings[ii] for all
    # inputs in the set x.
    if hashfn.fixed_symbols
        # If the symbols are fixed, throw an error whenever we encounter
        # a symbol we don't recognize.
        map(mapping ->
            minimum(get(mapping, xjj) do
                        "Symbol $(xjj) not found" |>
                        ErrorException |>
                        throw
                    end
                    for xjj in x),
            hashfn.mappings)
    else
        # If the symbols are non-fixed, assign new integer labels to
        # new symbols.
        map(mapping ->
            minimum(get!(mapping, xjj) do
                        rand(I)
                    end
                    for xjj in x),
            hashfn.mappings)
    end
end


text:
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("FVadvection2d");

useLog("FVadvection2dlog")

# Configuration setup
domain(2)
solverType(FV)
timeStepper(RK4)

# Mesh
n = 15 # number of elements in each direction
mesh(QUADMESH, elsperdim=n, bids=4)

# Variables and BCs
u = variable("u", SCALAR, CELL)
boundary(u, 1, FLUX, "(abs(y-0.2) < 0.11) ? sin(2*pi*t)^2 : 0") # x=0
boundary(u, 2, NO_BC) # x=1
boundary(u, 3, NO_BC) # y=0
boundary(u, 4, NO_BC) # y=1

# Time interval and initial condition
T = 1;
timeInterval(T)
initial(u, "0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
coefficient("a", ["cos(pi*x/2)","sin(pi*x/2)"], VECTOR) # advection velocity
# The "upwind" function applies upwinding to the term (a.n)*u with flow velocity a.
# The optional third parameter is for tuning. Default upwind = 0, central = 1. Choose something between these.
flux(u, "upwind(a,u)")
# Note that there is no source() for this problem

exportCode("fvad2dcode") # uncomment to export generated code to a file
#importCode("fvad2dcode") # uncomment to import code from a file

solve(u)

finalize_femshop()

##### Uncomment below to plot

xy = Femshop.fv_info.cellCenters

using Plots
pyplot();
display(plot(xy[1,:], xy[2,:], u.values[:], st=:surface))


text:
x = [1, 2, 4, -3]

y = [1.1, π, -4/3, 1e-1]

z = Float64[1, 2, 4, -3]    # Without the `Float64`, Julia would make this an integer array

push!(z, 10)

append!(z, y)

y[3]   # Access the 3rd element in the array y

z[1] = z[end] * z[end-1]      # Set the 1st element in z to the product of the last two elements
z                             # To print the entire array

for i = 1:length(x)
    println("Element ", i, " of the array x has the value ", x[i])
end

for element in x
    println("This element of x has the value ", element)
end

function SieveOfEratosthenes(n)
    # Create a boolean array of length n, initialize all entries to as true
    # After the algorithm finishes, prime[i] will be true is i≥2 is a prime
    prime = trues(n)

    for p = 2:floor(Int, sqrt(n))
        if prime[p]
            for i = p^2:p:n
                prime[i] = false
            end
        end
    end

    # Return an array with all prime numbers
    primes = Int64[]
    for i = 2:n
        if prime[i]
            push!(primes, i)
        end
    end
    primes
end

println(SieveOfEratosthenes(100))            # All prime numbers up to 100

println(SieveOfEratosthenes(1000000)[end])   # The largest prime number less than 1,000,000

x = [1, 3, 5]
y = [0.2, 4.0 , 3.1]
z = x .* y   # Element-by-element multiplication

sqrt.(x)     # Element-by-element square roots

w = x .+ sqrt.(y .^ x)      # More complex mathematical expression, element-by-element

@. x + sqrt(y^x)       # Same as before

push!(y, 7)
x .* y      # Error: arrays must have same number of elements

x .* -3

-3x

function pick_largest(x, y)
    if x > y
        return x
    else
        return y
    end
end

println("z = ", z)
println("w = ", w)
println("elementwise largest = ", pick_largest.(z,w))

println(y)               # Original vector
println(y[1:3])          # First 3 elements
println(y[1:2:4])        # All odd-numbered elements
println(y[end:-1:2])     # From end back to second element in reverse
println(y[4:3])          # Empty subset
println(y[:])            # All elements (same as original vector)
println(y[[4,2,4,3,3]])  # Index by vector - pick elements 4,2,4,3,3

x1 = collect(1:3)

x2 = collect(10:-2.5:-3)

x = [1, -2, 3]
y = x
y[2] = 123
println("y = ", y)
println("x = ", x)

z = copy(x)
z[3] = -53
println("z = ", z)
println("x = ", x)

function modify_scalar(x)
    x = 111
    return nothing
end

function modify_vector!(x)
    x[:] .= 111
    return nothing
end

x = 0
modify_scalar(x)
println("x = ", x)    # Still 0 - function does not modify a scalar

x = zeros(1,5)
modify_vector!(x)
println("x = ", x)    # Function modifies the original vector


text:
using RandomWalks
using BenchmarkTools
using EmpiricalCDFs: EmpiricalCDF, get_data
using Statistics

nsamples = 10^4
nsteps = 10^3
tolerance = 1e-2

ecdf = EmpiricalCDF()
trial_loop = SampleLoopActor(nsamples, ECDFActor(get_x, ecdf))
@time trial!(WalkPlan(ContinuousWalk(WalkB(Point(0.0))), StepLimitActor(nsteps)), trial_loop)

x = get_data(ecdf)
x ./= sqrt(nsteps)

@info "Test standard deviation and mean of an ordinary RW"
@show isapprox(std(ecdf), 1; atol = tolerance)
@show isapprox(mean(ecdf), 0; atol = tolerance)