text:
# Composition
# ===========
#
# Sequence composition counter.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

type Composition{A<:Alphabet}
    counts::Vector{Int}
    function Composition()
        new(zeros(Int, length(alphabet(A))))
    end
end

typealias DNANucleotideCounts{n} Composition{DNAAlphabet{n}}
typealias RNANucleotideCounts{n} Composition{RNAAlphabet{n}}
typealias AminoAcidCounts Composition{AminoAcidAlphabet}

function Composition{A}(seq::BioSequence{A})
    comp = Composition{A}()
    @inbounds for x in seq
        comp.counts[encode(A, x)+1] += 1
    end
    return comp
end

function Composition(seq::ReferenceSequence)
    comp = Composition{DNAAlphabet{4}}()
    @inbounds for x in seq
        comp.counts[encode(DNAAlphabet{4}, x)+1] += 1
    end
    return comp
end

function Composition{K}(kmer::DNAKmer{K})
    A = DNAAlphabet{2}
    comp = Composition{A}()
    x = UInt64(kmer)
    @inbounds begin
        comp.counts[encode(A, DNA_A)+1] = count_a(x) - (32 - K)
        comp.counts[encode(A, DNA_C)+1] = count_c(x)
        comp.counts[encode(A, DNA_G)+1] = count_g(x)
        comp.counts[encode(A, DNA_T)+1] = count_t(x)
    end
    return comp
end

function Composition{K}(kmer::RNAKmer{K})
    A = RNAAlphabet{2}
    comp = Composition{A}()
    x = UInt64(kmer)
    @inbounds begin
        comp.counts[encode(A, RNA_A)+1] = count_a(x) - (32 - K)
        comp.counts[encode(A, RNA_C)+1] = count_c(x)
        comp.counts[encode(A, RNA_G)+1] = count_g(x)
        comp.counts[encode(A, RNA_U)+1] = count_t(x)
    end
    return comp
end

"""
    composition(seq)

Calculate composition of biological symbols in `seq`.
"""
function composition(seq::Sequence)
    return Composition(seq)
end

function Base.getindex{A}(comp::Composition{A}, x)
    try
        i = encode(A, x) + 1
        return Int(comp.counts[i])
    catch ex
        if isa(ex, EncodeError{A})
            return 0
        end
        rethrow()
    end
end

Base.summary(::DNANucleotideCounts) = "DNA Nucleotide Counts"
Base.summary(::RNANucleotideCounts) = "RNA Nucleotide Counts"
Base.summary(::AminoAcidCounts)     = "Amino Acid Counts"

function Base.show{A}(io::IO, comp::Composition{A})
    print(io, summary(comp), ':')
    for x in alphabet(A)
        print(io, "\n  ", x, " => ", comp[x])
    end
end


text:
using Pkg
pkg"activate ."

using Genie, Genie.Router, Genie.Renderer

markup = """<button class="nes-btn" id="submit-button">Go!</button>"""
markup *= """<label for="whatever"></label>"""

route("/") do
  html(markup)
end

Genie.AppServer.startup(async = false)

ext: jl


text:
module TestIteration

  using GBIF
  using Test

  set = occurrences()

  @test iterate(set) == (set[1], 2)
  @test iterate(set, 2) == (set[2], 3)

  plotor = taxon("Procyon lotor")
  plotor_occ = occurrences(plotor)
  occurrences!(plotor_occ)
  for o in plotor_occ
      @test typeof(o) <: GBIFRecord
      @test o.taxon.species.second == plotor.species.second
  end

end


text:
const Menu = Ptr{Cvoid}
const Control = Ptr{Cvoid}
const Label = Ptr{Cvoid}
const Button = Ptr{Cvoid}

@enum MenuOrientation begin
  Vertical = 0
  Horizontal = 1
  Central = 2
end

# Menu
function createMenu()
  ccall((:createMenu, @fullLibraryPath), Menu, ())
end

function setMenuOrientation(menu, orientation)
  ccall((:setMenuOrientation, @fullLibraryPath), Cvoid, (Menu, MenuOrientation), menu, orientation)
end

function setMenuPosition(menu, position)
  ccall((:setMenuPosition, @fullLibraryPath), Cvoid, (Menu, Vector2), menu, position)
end

function setMenuSize(menu, size)
  ccall((:setMenuSize, @fullLibraryPath), Cvoid, (Menu, Vector2), menu, size)
end

function setMenuColor(menu, color)
  ccall((:setMenuColor, @fullLibraryPath), Cvoid, (Menu, Vector4), menu, color)
end

function addControlToMenu(menu, control)
  ccall((:addControlToMenu, @fullLibraryPath), Cvoid, (Menu, Control), menu, control)
end

function updateMenu(menu)
  ccall((:updateMenu, @fullLibraryPath), Cvoid, (Menu,), menu)
end

# Label
function createLabel()
  ccall((:createLabel, @fullLibraryPath), Label, ())
end

function setLabelZ(label, z)
  ccall((:setLabelZ, @fullLibraryPath), Cvoid, (Label, Cint), label, z)
end

function setLabelPosition(label, position)
  ccall((:setLabelPosition, @fullLibraryPath), Cvoid, (Label, Vector2), label, position)
end

function setLabelSize(label, size)
  ccall((:setLabelSize, @fullLibraryPath), Cvoid, (Label, Vector2), label, size)
end

function setLabelColor(label, color)
  ccall((:setLabelColor, @fullLibraryPath), Cvoid, (Label, Vector4), label, color)
end

function setLabelTextHeight(label, textHeight)
  ccall((:setLabelTextHeight, @fullLibraryPath), Cvoid, (Label, Cint), label, textHeight)
end

function setLabelText(label, text)
  ccall((:setLabelText, @fullLibraryPath), Cvoid, (Label, Cstring), label, text)
end

# Button
function createButton()
  ccall((:createButton, @fullLibraryPath), Button, ())
end

function setButtonZ(button, z)
  ccall((:setButtonZ, @fullLibraryPath), Cvoid, (Button, Cint), button, z)
end

function setButtonPosition(button, position)
  ccall((:setButtonPosition, @fullLibraryPath), Cvoid, (Button, Vector2), button, position)
end

function setButtonSize(button, size)
  ccall((:setButtonSize, @fullLibraryPath), Cvoid, (Button, Vector2), button, size)
end

function setButtonTextHeight(button, textHeight)
  ccall((:setButtonTextHeight, @fullLibraryPath), Cvoid, (Button, Cint), button, textHeight)
end

function setButtonText(button, text)
  ccall((:setButtonText, @fullLibraryPath), Cvoid, (Button, Cstring), button, text)
end

function setButtonAction(button, action)
  ccall((:setButtonAction, @fullLibraryPath), Cvoid, (Button, Ptr{Cvoid}), button, @cfunction($action, Cvoid, (Button,)))
end

function setButtonTextColor(button, textColor)
  ccall((:setButtonTextColor, @fullLibraryPath), Cvoid, (Button, Vector4), button, textColor)
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
# bootstrapped class of curve methods

# Forward curves

"""
    ForwardStarting(curve,forwardstart)

Rebase a `curve` so that `discount`/`accumulation`/etc. are re-based so that time zero from the new curves perspective is the given `forwardstart` time.

# Examples

```julia-repl
julia> zero = [5.0, 5.8, 6.4, 6.8] ./ 100
julia> maturity = [0.5, 1.0, 1.5, 2.0]
julia> curve = Yields.Zero(zero, maturity)
julia> fwd = Yields.ForwardStarting(curve, 1.0)

julia> discount(curve,1,2)
0.9275624570410582

julia> discount(fwd,1) # `curve` has effectively been reindexed to `1.0`
0.9275624570410582
```

# Extended Help

While `ForwardStarting` could be nested so that, e.g. the third period's curve is the one-period forward of the second period's curve, it will be more efficient to reuse the initial curve from a runtime and compiler perspective.

`ForwardStarting` is not used to construct a curve based on forward rates. See  [`Forward`](@ref) instead.
"""
struct ForwardStarting{T,U} <: AbstractYield
    curve::U
    forwardstart::T
end

function discount(c::ForwardStarting, to)
    discount(c.curve, c.forwardstart, to + c.forwardstart)
end

function Base.zero(c::ForwardStarting, to,cf::C) where {C<:CompoundingFrequency}
    z = forward(c.curve,c.forwardstart,to+c.forwardstart)
    return convert(cf,z)
end

"""
    zero(curve,time)
    zero(curve,time,CompoundingFrequency)

Return the zero rate for the curve at the given time. If not specified, will use `Periodic(1)` compounding.
"""
Base.zero(c::YC, time) where {YC<:AbstractYield} = zero(c, time, Periodic(1))
function Base.zero(c::YC, time, cf::C) where {YC<:AbstractYield,C<:CompoundingFrequency}
    return convert(cf, Continuous(c.zero(time))) # c.zero is a curve of continuous rates represented as floats. explicitly wrap in continuous before converting
end

"""
    Constant(rate::Real, cf::CompoundingFrequency=Periodic(1))
    Constant(r::Rate)

Construct a yield object where the spot rate is constant for all maturities. If `rate` is not a `Rate` type, will assume `Periodic(1)` for the compounding frequency

# Examples

```julia-repl
julia> y = Yields.Constant(0.05)
julia> discount(y,2)
0.9070294784580498     # 1 / (1.05) ^ 2
```
"""
struct Constant{T} <: AbstractYield
    rate::T
end

function Constant(rate::T, cf::C = Periodic(1)) where {T<:Real,C<:CompoundingFrequency}
    return Constant(Rate(rate, cf))
end

Base.zero(c::Constant, time) = c.rate
Base.zero(c::Constant, time, cf::CompoundingFrequency) = convert(cf, c.rate)
rate(c::Constant) = c.rate
rate(c::Constant, time) = c.rate
discount(r::Constant, time) = discount(r.rate, time)
discount(r::Constant, from, to) = discount(r.rate, to - from)
accumulation(r::Constant, time) = accumulation(r.rate, time)
accumulation(r::Constant, from, to) = accumulation(r.rate, to - from)

"""
    Step(rates,times)

Create a yield curve object where the applicable rate is the effective rate of interest applicable until corresponding time. If `rates` is not a `Vector{Rate}`, will assume `Periodic(1)` type.

# Examples

```julia-repl
julia>y = Yields.Step([0.02,0.05], [1,2])

julia>rate(y,0.5)
0.02

julia>rate(y,1.5)
0.05

julia>rate(y,2.5)
0.05
```
"""
struct Step <: AbstractYield
    rates
    times
end

Step(rates) = Step(rates, collect(1:length(rates)))

function rate(y::Step, time)
    i = findfirst(t -> time <= t, y.times)
    if isnothing(i)
        return y.rates[end]
    else
        return y.rates[i]
    end
end

function discount(y::Step, time)
    v = 1 / (1 + y.rates[1])^min(y.times[1], time)

    if y.times[1] >= time
        return v
    end

    for i = 2:length(y.times)

        if y.times[i] >= time
            # take partial discount and break
            v /= (1 + y.rates[i])^(time - y.times[i-1])
            break
        else
            # take full discount and continue
            v /= (1 + y.rates[i])^(y.times[i] - y.times[i-1])
        end

    end

    return v
end

"""
    Zero(rates, maturities; interpolation=CubicSpline())

Construct a yield curve with given zero-coupon spot `rates` at the given `maturities`. If `rates` is not a `Vector{Rate}`, will assume `Periodic(1)` type.

See [`bootstrap`](@ref) for more on the `interpolation` parameter, which is set to `CubicSpline()` by default.
"""
function Zero(rates::Vector{<:Rate}, maturities; interpolation=CubicSpline())
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))
    return _zero_inner(rates,maturities,interpolation)
end

# zero is different than the other boostrapped curves in that it doesn't actually need to bootstrap
# because the rate are already zero rates. Instead, we just cut straight to the
# appropriate interpolation function based on the type dispatch.
function _zero_inner(rates::Vector{<:Rate}, maturities, interp::CubicSpline)
    continuous_zeros = rate.(convert.(Continuous(), rates))
    return YieldCurve(
        rates,
        maturities,
        cubic_interp([0.0; maturities],[first(continuous_zeros); continuous_zeros])
    )
end

function _zero_inner(rates::Vector{<:Rate}, maturities, interp::LinearSpline)
    continuous_zeros = rate.(convert.(Continuous(), rates))
    return YieldCurve(
        rates,
        maturities,
        linear_interp([0.0; maturities],[first(continuous_zeros); continuous_zeros])
    )
end

# fallback for user provied interpolation function
function _zero_inner(rates::Vector{<:Rate}, maturities, interp)
    continuous_zeros = rate.(convert.(Continuous(), rates))
    return YieldCurve(
        rates,
        maturities,
        interp([0.0; maturities],[first(continuous_zeros); continuous_zeros])
    )
end

#fallback if `rates` aren't `Rate`s. Assume `Periodic(1)` per Zero docstring
function Zero(rates, maturities; interpolation=CubicSpline())
    Zero(Periodic.(rates, 1), maturities; interpolation)
end

function Zero(rates; interpolation=CubicSpline())
    # bump to a constant yield if only given one rate
    maturities = collect(1:length(rates))
    return Zero(rates, maturities; interpolation)
end

"""
    Par(rates, maturities; interpolation=CubicSpline())

Construct a curve given a set of bond equivalent yields and the corresponding maturities. Assumes that maturities <= 1 year do not pay coupons and that after one year, pays coupons with frequency equal to the CompoundingFrequency of the corresponding rate.

See [`bootstrap`](@ref) for more on the `interpolation` parameter, which is set to `CubicSpline()` by default.

# Examples

```julia-repl

julia> par = [6.,8.,9.5,10.5,11.0,11.25,11.38,11.44,11.48,11.5] ./ 100
julia> maturities = [t for t in 1:10]
julia> curve = Par(par,maturities);
julia> zero(curve,1)
Rate(0.06000000000000005, Periodic(1))

```
"""
function Par(rates::Vector{<:Rate}, maturities; interpolation=CubicSpline())
    # bump to a constant yield if only given one rate
    if length(rates) == 1
        return Constant(rate[1])
    end
    return YieldCurve(
        rates,
        maturities,
        # assume that maturities less than or equal to 12 months are settled once, otherwise semi-annual
        # per Hull 4.7
        bootstrap(rates, maturities, [m <= 1 ? nothing : 1 / r.compounding.frequency for (r, m) in zip(rates, maturities)], interpolation)
    )
end

function Par(rates::Vector{T}, maturities; interpolation=CubicSpline()) where {T<:Real}
    return Par(Rate.(rates), maturities; interpolation)
end

"""
    Forward(rate_vector,maturities)

Takes a vector of 1-period forward rates and constructs a discount curve. If rate_vector is not a vector of `Rates` (ie is just a vector of `Float64` values), then
the assumption is that each value is `Periodic` rate compounded once per period.

# Examples

```julia-repl
julia> Yields.Forward( [0.01,0.02,0.03] );

julia> Yields.Forward( Yields.Continuous.([0.01,0.02,0.03]) );

```
"""
function Forward(rates::Vector{<:Rate}, maturities)
    # convert to zeros and pass to Zero
    disc_v = Vector{Float64}(undef, length(rates))

    v = 1.0

    for i = 1:length(rates)
        Δt = maturities[i] - (i == 1 ? 0 : maturities[i-1])
        v *= discount(rates[i], Δt)
        disc_v[i] = v
    end

    z = (1.0 ./ disc_v) .^ (1 ./ maturities) .- 1 # convert disc_v to zero
    return Zero(z, maturities)
end

# if rates isn't a vector of Rates, then we assume periodic rates compounded once per period.
function Forward(rates, maturities)
    return Forward(Yields.Periodic.(rates, 1), maturities)
end

Forward(rates) = Forward(rates, collect(1:length(rates)))

"""
    Yields.CMT(rates, maturities; interpolation=CubicSpline())

Takes constant maturity (treasury) yields (bond equivalent), and assumes that instruments <= one year maturity pay no coupons and that the rest pay semi-annual.

See [`bootstrap`](@ref) for more on the `interpolation` parameter, which is set to `CubicSpline()` by default.

# Examples

```
# 2021-03-31 rates from Treasury.gov
rates =[0.01, 0.01, 0.03, 0.05, 0.07, 0.16, 0.35, 0.92, 1.40, 1.74, 2.31, 2.41] ./ 100
mats = [1/12, 2/12, 3/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30]

Yields.CMT(rates,mats)
```
"""
function CMT(rates::Vector{T}, maturities; interpolation=CubicSpline()) where {T<:Real}
    rs = map(zip(rates, maturities)) do (r, m)
        if m <= 1
            Rate(r, Periodic(1 / m))
        else
            Rate(r, Periodic(2))
        end
    end

    CMT(rs, maturities;interpolation)
end

function CMT(rates::Vector{<:Rate}, maturities; interpolation=CubicSpline())
    return YieldCurve(
        rates,
        maturities,
        # assume that maturities less than or equal to 12 months are settled once, otherwise semi-annual
        # per Hull 4.7
        bootstrap(rates, maturities, [m <= 1 ? nothing : 0.5 for m in maturities], interpolation)
    )
end

"""
    OIS(rates,maturities)
Takes Overnight Index Swap rates, and assumes that instruments <= one year maturity are settled once and other agreements are settled quarterly with a corresponding CompoundingFrequency.

See [`bootstrap`](@ref) for more on the `interpolation` parameter, which is set to `CubicSpline()` by default.

"""
function OIS(rates::Vector{T}, maturities; interpolation=CubicSpline()) where {T<:Real}
    rs = map(zip(rates, maturities)) do (r, m)
        if m <= 1
            Rate(r, Periodic(1 / m))
        else
            Rate(r, Periodic(4))
        end
    end

    return OIS(rs, maturities; interpolation)
end
function OIS(rates::Vector{<:Rate}, maturities ; interpolation=CubicSpline())
    return YieldCurve(
        rates,
        maturities,
        # assume that maturities less than or equal to 12 months are settled once, otherwise quarterly
        # per Hull 4.7
        bootstrap(rates, maturities, [m <= 1 ? nothing : 1 / 4 for m in maturities], interpolation)
    )
end


text:
using JuMP, MOModelingExtensions
using Test, GLPK

@testset "Modeling extensions" begin
    @testset "Maximum" begin
        @testset "No bounds" begin
            m = Model()
            @variable(m, x <= 5)
            @variable(m, y <= 5)

            @test_throws ErrorException @constraint(m, max(x, y) <= 4)
            @test_throws ErrorException @constraint(m, maximum([x, y]) <= 4)
        end

        @testset "max" begin
            @testset "Both variables nonnegative" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, 0 <= y <= 5)
                @constraint(m, max(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 4.0
                @test value(y) ≈ 4.0
            end

            @testset "First variable nonpositive" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, -5 <= x <= 0)
                @variable(m, 0 <= y <= 5)
                @constraint(m, max(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 0.0
                @test value(y) ≈ 4.0
            end

            @testset "Second variable nonpositive" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, -5 <= y <= 0)
                @constraint(m, max(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 4.0
                @test value(y) ≈ 0.0
            end

            @testset "Both variables mixed sign" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, -5 <= x <= 5)
                @variable(m, -5 <= y <= 5)
                @constraint(m, max(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 4.0
                @test value(y) ≈ 4.0
            end
        end

        @testset "maximum" begin
            @testset "maximum([x, y])" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, 0 <= y <= 5)
                @constraint(m, maximum([x, y]) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 4.0
                @test value(y) ≈ 4.0
            end

            @testset "maximum([x, y, z])" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, 0 <= y <= 5)
                @variable(m, 4 <= z <= 5)
                @constraint(m, maximum([x, y, z]) <= 4)
                @objective(m, Max, x + y + z)
                optimize!(m)

                @test value(x) ≈ 4.0
                @test value(y) ≈ 4.0
                @test value(y) ≈ 4.0
            end
        end
    end

    @testset "Minimum" begin
        @testset "No bounds" begin
            m = Model()
            @variable(m, x <= 5)
            @variable(m, y <= 5)

            @test_throws ErrorException @constraint(m, min(x, y) <= 4)
            @test_throws ErrorException @constraint(m, minimum([x, y]) <= 4)
        end

        @testset "min" begin
            @testset "Both variables nonnegative" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, 0 <= y <= 5)
                @constraint(m, min(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                if value(x) ≈ 4.0
                    @test value(x) ≈ 4.0
                    @test value(y) ≈ 5.0
                else
                    @test value(x) ≈ 5.0
                    @test value(y) ≈ 4.0
                end
            end

            @testset "First variable nonpositive" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, -5 <= x <= 0)
                @variable(m, 0 <= y <= 5)
                @constraint(m, min(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 0.0
                @test value(y) ≈ 5.0
            end

            @testset "Second variable nonpositive" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, -5 <= y <= 0)
                @constraint(m, min(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                @test value(x) ≈ 5.0
                @test value(y) ≈ 0.0
            end

            @testset "Both variables mixed sign" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, -5 <= x <= 5)
                @variable(m, -5 <= y <= 5)
                @constraint(m, min(x, y) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                if value(x) ≈ 4.0
                    @test value(x) ≈ 4.0
                    @test value(y) ≈ 5.0
                else
                    @test value(x) ≈ 5.0
                    @test value(y) ≈ 4.0
                end
            end
        end

        @testset "minimum" begin
            @testset "minimum([x, y])" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, 0 <= y <= 5)
                @constraint(m, minimum([x, y]) <= 4)
                @objective(m, Max, x + y)
                optimize!(m)

                if value(x) ≈ 4.0
                    @test value(x) ≈ 4.0
                    @test value(y) ≈ 5.0
                else
                    @test value(x) ≈ 5.0
                    @test value(y) ≈ 4.0
                end
            end

            @testset "minimum([x, y, z])" begin
                m = Model(with_optimizer(GLPK.Optimizer))
                @variable(m, 0 <= x <= 5)
                @variable(m, 0 <= y <= 5)
                @variable(m, 4 <= z <= 5)
                @constraint(m, minimum([x, y, z]) <= 4)
                @objective(m, Max, x + y + z)
                optimize!(m)

                if value(x) ≈ 4.0
                    @test value(x) ≈ 4.0
                    @test value(y) ≈ 5.0
                    @test value(z) ≈ 5.0
                elseif value(y) ≈ 4.0
                    @test value(x) ≈ 5.0
                    @test value(y) ≈ 4.0
                    @test value(z) ≈ 5.0
                else
                    @test value(x) ≈ 5.0
                    @test value(y) ≈ 5.0
                    @test value(z) ≈ 4.0
                end
            end
        end
    end

    @testset "Boolean operators" begin
        @testset "&" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, x & y == true)
            @objective(m, Max, x + y)
            optimize!(m)

            @test value(x) ≈ 1.0
            @test value(y) ≈ 1.0

            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, x & y == true)
            @constraint(m, x + y <= 1)
            @objective(m, Max, x + y)
            optimize!(m)

            @test termination_status(m) == MOI.INFEASIBLE
        end

        @testset "all([x, y])" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, all([x, y]) == true)
            @objective(m, Max, x + y)
            optimize!(m)

            @test value(x) ≈ 1.0
            @test value(y) ≈ 1.0

            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, all([x, y]) == true)
            @constraint(m, x + y <= 1)
            @objective(m, Max, x + y)
            optimize!(m)

            @test termination_status(m) == MOI.INFEASIBLE
        end

        @testset "all([x, y, z])" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @variable(m, z, Bin)
            @constraint(m, all([x, y, z]) == true)
            @objective(m, Max, x + y)
            optimize!(m)

            @test value(x) ≈ 1.0
            @test value(y) ≈ 1.0
            @test value(z) ≈ 1.0

            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @variable(m, z, Bin)
            @constraint(m, all([x, y, z]) == true)
            @constraint(m, x + y <= 1)
            @objective(m, Max, x + y + z)
            optimize!(m)

            @test termination_status(m) == MOI.INFEASIBLE
        end

        @testset "|" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, x | y == true)
            @objective(m, Min, x + y)
            optimize!(m)

            if value(x) ≈ 1.0
                @test value(x) ≈ 1.0
                @test value(y) ≈ 0.0
            else
                @test value(x) ≈ 0.0
                @test value(y) ≈ 1.0
            end

            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, x | y == true)
            @constraint(m, x + y <= 1)
            @objective(m, Max, x + y)
            optimize!(m)

            if value(x) ≈ 1.0
                @test value(x) ≈ 1.0
                @test value(y) ≈ 0.0
            else
                @test value(x) ≈ 0.0
                @test value(y) ≈ 1.0
            end
        end

        @testset "any([x, y])" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, any([x, y]) == true)
            @objective(m, Min, x + y)
            optimize!(m)

            if value(x) ≈ 1.0
                @test value(x) ≈ 1.0
                @test value(y) ≈ 0.0
            else
                @test value(x) ≈ 0.0
                @test value(y) ≈ 1.0
            end

            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, any([x, y]) == true)
            @constraint(m, x + y <= 1)
            @objective(m, Max, x + y)
            optimize!(m)

            if value(x) ≈ 1.0
                @test value(x) ≈ 1.0
                @test value(y) ≈ 0.0
            else
                @test value(x) ≈ 0.0
                @test value(y) ≈ 1.0
            end
        end

        @testset "any([x, y, z])" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @variable(m, z, Bin)
            @constraint(m, any([x, y, z]) == true)
            @objective(m, Min, x + y + z)
            optimize!(m)

            @test value(x) + value(y) + value(z) ≈ 1.0
        end

        @testset "⊻" begin
            m = Model(with_optimizer(GLPK.Optimizer))
            @variable(m, x, Bin)
            @variable(m, y, Bin)
            @constraint(m, x ⊻ y == true)
            @objective(m, Min, x + y)
            optimize!(m)

            if value(x) ≈ 1.0
                @test value(x) ≈ 1.0
                @test value(y) ≈ 0.0
            else
                @test value(x) ≈ 0.0
                @test value(y) ≈ 1.0
            end
        end
    end
end


text:
using LabelledArrays

export rossler_rossler_bidir

@inline @inbounds function eom_rossler_rossler_bidir(u, p, t)
    ω₁, ω₂, c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃ = (p...,)
    x1, x2, x3, y1, y2, y3 = (u...,)

    dx1 = -ω₁*(x2 + x3) + c_yx*(y1 - x1)
    dx2 = ω₁*x1 + a₁*x2
    dx3 = a₂ + x3*(x1 - a₃)

    dy1 = -ω₂*(y2 + y3) + c_xy*(x1 - y1)
    dy2 = ω₂*y1 + b₁*y2
    dy3 = b₂ + y3*(y1 - b₃)

    return SVector{6}(dx1, dx2, dx3, dy1, dy2, dy3)
end

"""
    rossler_rossler_bidir(; u0 = rand(6),
        ω₁ = 1.015, ω₂ = 0.985,
        c_xy = 0.1, c_yx = 0.1,
        a₁ = 0.15, a₂ = 0.2, a₃ = 10,
        b₁ = 0.15, b₂ = 0.2, b₃ = 10)

Initialise a system of two bidirectionally coupled 3D Rössler systems.
This system has been modified from [1] to allow other parameterisations,
but default parameters are as in [1].

The ``X`` and ``Y`` subsystems are mostly synchronized for
`c_xy > 0.1` or `c_yx > 0.1`.

## Equations of motion

The dynamics is generated by the following vector field

```math
\\begin{aligned}
\\dot{x_1} &= -\\omega_1(x_2 + x_3) + c_{yx}(y_1 - x_1) \\\\
\\dot{x_2} &= \\omega_1 x_1 + a_1 x_2 \\\\
\\dot{x_3} &= a_2 + x_3 (x_1 - a_3) \\\\
\\dot{y_1} &= -\\omega_2 (y_2 + y_3) + c_{xy}(x_1 - y_1) \\\\
\\dot{y_2} &= \\omega_2 y_1 + b_1 y_2 \\\\
\\dot{y_3} &= b_2 + y_3 (y_1 - b_3)
\\end{aligned}
```

## References

1. Amigó, José M., and Yoshito Hirata. "Detecting directional couplings from
    multivariate flows by the joint distance distribution." Chaos: An
    Interdisciplinary Journal of Nonlinear Science 28.7 (2018): 075302.
"""
function rossler_rossler_bidir(; u0 = rand(6),
        ω₁ = 1.015, ω₂ = 0.985,
        c_xy = 0.1, c_yx = 0.1,
        a₁ = 0.15, a₂ = 0.2, a₃ = 10,
        b₁ = 0.15, b₂ = 0.2, b₃ = 10)

    p = @LArray [ω₁, ω₂, c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃] (:ω₁, :ω₂, :c_xy, :c_yx, :a₁, :a₂, :a₃, :b₁, :b₂, :b₃)
    ContinuousDynamicalSystem(eom_rossler_rossler_bidir, u0, p)
end

function rossler_rossler_bidir_trajectory(npts, sample_dt; n_transient = 2000, dt = 0.2,
        u0 = rand(6), ω₁ = 1.015, ω₂ = 0.985,
        c_xy = 0.2, c_yx = 0.2,
        a₁ = 0.15, a₂ = 0.2, a₃ = 10,
        b₁ = 0.15, b₂ = 0.2, b₃ = 10)

    s = rossler_rossler_bidir(u0 = u0,
        c_xy = c_xy, c_yx = c_yx,
        ω₁ = ω₁, ω₂ = ω₂,
        a₁ = a₁, a₂ = a₂, a₃ = a₃,
        b₁ = b₁, b₂ = b₂, b₃ = b₃)

    # the system is recorded at times t0:dt:T
    T = npts*dt*sample_dt
    o = trajectory(s, T, dt = dt, Ttr = n_transient*dt, alg = SimpleDiffEq.SimpleATsit5())[1:sample_dt:end-1, :] #alg = SimpleDiffEq.SimpleATsit5()
end

function good_rossler_rossler_bidir_trajectory(npts; sample_dt = 1,
        Da₁ = Uniform(0.12, 0.17),
        Da₂ = Uniform(0.18, 0.22),
        Da₃ = Uniform(9.0, 11.0),
        Db₁ = Uniform(0.10, 0.20),
        Db₂ = Uniform(0.18, 0.22),
        Db₃ = Uniform(9.0, 11.0),
        Dω₁ = Uniform(0.95, 0.999), #Uniform(0.97, 1.03)
        Dω₂ = Uniform(1.001, 1.05), #Uniform(0.97, 1.03)
        a₁ = nothing,
        a₂ = nothing,
        a₃ = nothing,
        b₁ = nothing,
        b₂ = nothing,
        b₃ = nothing,
        ω₁ = nothing,
        ω₂ = nothing,
        c_xy = 0.2,  c_yx = 0.2,
        Ttr = 5000, dt = 0.2,
        u0 = rand(6),
        n_maxtries = 300)

    n_tries = 0

    while n_tries <= n_maxtries
        ω₁ == nothing ? ω₁ = rand(Dω₁) : ω₁ = ω₁
        ω₂ == nothing ? ω₂ = rand(Dω₂) : ω₂ = ω₂
        a₁ == nothing ? a₁ = rand(Da₁) : a₁ = a₁
        a₂ == nothing ? a₂ = rand(Da₂) : a₂ = a₂
        a₃ == nothing ? a₃ = rand(Da₃) : a₃ = a₃
        b₁ == nothing ? b₁ = rand(Db₁) : b₁ = b₁
        b₂ == nothing ? b₂ = rand(Db₂) : b₂ = b₂
        b₃ == nothing ? b₃ = rand(Db₃) : b₃ = b₃

        pts = rossler_rossler_bidir_trajectory(npts, sample_dt, dt = dt, n_transient = Ttr,
            c_xy = c_xy,  c_yx = c_yx, ω₁ = ω₁, ω₂ = ω₂, a₁ = a₁, a₂ = a₂,a₃ = a₃,b₁ = b₁,  b₂ = b₂, b₃ = b₃)

        if all(Matrix(pts) .< 1e10) && length(unique(pts)) > npts/2
            return pts
        end
        println("no attractor found. trying with new initial condition and parameters")
        n_tries += 1
    end
end


text:
# Constants and types go here

const portaudio = "libportaudio"
const sndfile = "libsndfile"

type DeviceInfo
    structVersion::Int32
    name::Ptr{Uint8}
    hostApi::Int32
    maxInputChannels::Int32
    maxOutputChannels::Int32
    defaultLowInputLatency::Float64
    defaultLowOutputLatency::Float64
    defaultHighInputLatency::Float64
    defaultHighOutputLatency::Float64
    defaultSampleRate::Float64
end

type StreamParameters
    device::Int32
    channelCount::Int32
    sampleFormat::Uint
    suggestedLatency::Float64
    hostApiSpecificStreamInfo::Ptr{Void}
end

typealias PaStream Ptr{Void}

type SF_INFO
    frames::Int64
    samplerate::Int32
    channels::Int32
    format::Int32
    sections::Int32
    seekable::Int32
end


text:
module PotentialWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AstronomicalAmplitudeReader = @jimport org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader
const BodiesElements = @jimport org.orekit.data.BodiesElements
const CachedNormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.CachedNormalizedSphericalHarmonicsProvider
const Class = @jimport java.lang.Class
const ConstantSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.ConstantSphericalHarmonics
const DataLoader = @jimport org.orekit.data.DataLoader
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const EGMFormatReader = @jimport org.orekit.forces.gravity.potential.EGMFormatReader
const Enum = @jimport java.lang.Enum
const FESCHatEpsilonReader = @jimport org.orekit.forces.gravity.potential.FESCHatEpsilonReader
const FESCnmSnmReader = @jimport org.orekit.forces.gravity.potential.FESCnmSnmReader
const GRGSFormatReader = @jimport org.orekit.forces.gravity.potential.GRGSFormatReader
const GravityFieldFactory = @jimport org.orekit.forces.gravity.potential.GravityFieldFactory
const GravityFields = @jimport org.orekit.forces.gravity.potential.GravityFields
const ICGEMFormatReader = @jimport org.orekit.forces.gravity.potential.ICGEMFormatReader
const InputStream = @jimport java.io.InputStream
const JString = @jimport java.lang.String
const LazyLoadedGravityFields = @jimport org.orekit.forces.gravity.potential.LazyLoadedGravityFields
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const NormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider
const NormalizedSphericalHarmonicsProvider_NormalizedSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider$NormalizedSphericalHarmonics
const Normalizer = @jimport org.orekit.forces.gravity.potential.Normalizer
const Object = @jimport java.lang.Object
const OceanLoadDeformationCoefficients = @jimport org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients
const OceanTidesReader = @jimport org.orekit.forces.gravity.potential.OceanTidesReader
const OceanTidesWave = @jimport org.orekit.forces.gravity.potential.OceanTidesWave
const Optional = @jimport java.util.Optional
const PotentialCoefficientsReader = @jimport org.orekit.forces.gravity.potential.PotentialCoefficientsReader
const PulsatingSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.PulsatingSphericalHarmonics
const RawSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.RawSphericalHarmonicsProvider
const RawSphericalHarmonicsProvider_RawSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.RawSphericalHarmonicsProvider$RawSphericalHarmonics
const SHMFormatReader = @jimport org.orekit.forces.gravity.potential.SHMFormatReader
const SecularTrendSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.SecularTrendSphericalHarmonics
const SphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.SphericalHarmonicsProvider
const TideSystem = @jimport org.orekit.forces.gravity.potential.TideSystem
const TideSystemProvider = @jimport org.orekit.forces.gravity.potential.TideSystemProvider
const TimeScale = @jimport org.orekit.time.TimeScale
const UnnormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider
const UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider$UnnormalizedSphericalHarmonics
const Unnormalizer = @jimport org.orekit.forces.gravity.potential.Unnormalizer
const WrappingNormalizedProvider = @jimport org.orekit.forces.gravity.potential.WrappingNormalizedProvider
const WrappingUnnormalizedProvider = @jimport org.orekit.forces.gravity.potential.WrappingUnnormalizedProvider

export AstronomicalAmplitudeReader
export CachedNormalizedSphericalHarmonicsProvider
export ConstantSphericalHarmonics
export EGMFormatReader
export FESCHatEpsilonReader
export FESCnmSnmReader
export GRGSFormatReader
export GravityFieldFactory
export GravityFields
export ICGEMFormatReader
export LazyLoadedGravityFields
export NormalizedSphericalHarmonicsProvider
export Normalizer
export OceanLoadDeformationCoefficients
export OceanTidesReader
export OceanTidesWave
export PotentialCoefficientsReader
export PulsatingSphericalHarmonics
export RawSphericalHarmonicsProvider
export SHMFormatReader
export SecularTrendSphericalHarmonics
export SphericalHarmonicsProvider
export TideSystem
export TideSystemProvider
export UnnormalizedSphericalHarmonicsProvider
export Unnormalizer
export WrappingNormalizedProvider
export WrappingUnnormalizedProvider
export add_contribution
export add_default_ocean_tides_readers
export add_default_potential_coefficients_readers
export add_ocean_tides_reader
export add_potential_coefficients_reader
export can_add
export clear_ocean_tides_readers
export clear_potential_coefficients_readers
export compare_to
export configure_ocean_load_deformation_coefficients
export describe_constable
export equals
export get_ae
export get_astronomical_amplitudes_map
export get_class
export get_coefficients
export get_constant_normalized_provider
export get_constant_unnormalized_provider
export get_declaring_class
export get_doodson
export get_gravity_fields
export get_max_available_degree
export get_max_available_order
export get_max_degree
export get_max_order
export get_max_parse_degree
export get_max_parse_order
export get_mu
export get_normalized_provider
export get_ocean_load_deformation_coefficients
export get_ocean_tides_waves
export get_offset
export get_provider
export get_reference_date
export get_supported_names
export get_tide_system
export get_unnormalization_factors
export get_unnormalized_provider
export get_waves
export hash_code
export load_data
export missing_coefficients_allowed
export name
export notify
export notify_all
export on_date
export ordinal
export read_gravity_field
export set_max_parse_degree
export set_max_parse_order
export still_accepts_data
export to_string
export value_of
export values
export wait

include("astronomical_amplitude_reader.jl")
include("cached_normalized_spherical_harmonics_provider.jl")
include("constant_spherical_harmonics.jl")
include("egm_format_reader.jl")
include("fes_cnm_snm_reader.jl")
include("fesc_hat_epsilon_reader.jl")
include("gravity_field_factory.jl")
include("gravity_fields.jl")
include("grgs_format_reader.jl")
include("icgem_format_reader.jl")
include("lazy_loaded_gravity_fields.jl")
include("normalized_spherical_harmonics_provider.jl")
include("normalizer.jl")
include("ocean_load_deformation_coefficients.jl")
include("ocean_tides_reader.jl")
include("ocean_tides_wave.jl")
include("potential_coefficients_reader.jl")
include("pulsating_spherical_harmonics.jl")
include("raw_spherical_harmonics_provider.jl")
include("secular_trend_spherical_harmonics.jl")
include("shm_format_reader.jl")
include("spherical_harmonics_provider.jl")
include("tide_system.jl")
include("tide_system_provider.jl")
include("unnormalized_spherical_harmonics_provider.jl")
include("unnormalizer.jl")
include("wrapping_normalized_provider.jl")
include("wrapping_unnormalized_provider.jl")

end


text:
function make_instance(client_matrix, station_matrix; p_dist=2, threashold=1e-3)
    A, Λ = init_alpha_matrix(client_matrix, station_matrix)

    for (j, (x, y, v, d_csr, d_cli)) in enumerate(eachrow(client_matrix))
        p_csr, p_cli = make_probability(x, y, d_csr, d_cli, p_dist)
        for (i, (sx, sy)) in enumerate(eachrow(station_matrix))
            A[i, j] = p_cli(sx, sy) * v
            Λ[i] += p_csr(sx, sy) * v
        end
    end
    A[A .<= threashold] .= 0
    A, Λ
end

function make_probability(i, j, d_csr, d_client, p_dist)
    d_csr /= log(p_dist)
    d_client /= log(p_dist)
    p_csr = (x, y) -> exp(-dist(i, j, x, y) / d_csr)
    p_client = (x, y) -> exp(-dist(i, j, x, y) / d_client)
    p_csr, p_client
end

function init_alpha_matrix(cm, sm)
    x = size(cm, 1)
    y = size(sm, 1)
    zeros(y, x), zeros(y)
end

dist = (x₁, y₁, x₂, y₂) -> √((x₁ - x₂)^2 + (y₁ - y₂)^2)


text:
"""
    Greyscale

    Greyscale(min=nothing, max=nothing)

A greeyscale scheme ith better performance than using a
Colorschemes.jl scheme as there is not array access or interpolation.

`min` and `max` are values between `0.0` and `1.0` that define the range of greys used.
"""
struct Greyscale{M1,M2}
    min::M1
    max::M2
end
Greyscale(; min=nothing, max=nothing) = Greyscale(min, max)

Base.get(scheme::Greyscale, x::Real) = scale(x, scheme.min, scheme.max)

const Grayscale = Greyscale

"""
    ObjectScheme

    ObjectScheme()

Default colorscheme. Similar to `GreyScale` for `Number`.

Other grid objects can define a custom method to return colors from composite objects:

```julia
DynamicGrids.to_rgb(::ObjectScheme, obj::MyObjectType) = ...
```

Which must return an `ARGB32` value.
"""
struct ObjectScheme end

to_rgb(scheme::ObjectScheme, x::Real) = to_rgb(x)


text:
"""
    find_global_maximum(graph::Dict{Int64,Array{Int64,1}})

Find the total global maximum based on a comparsion between the intial vertex
(`global_maximum = 0`), the current dictionary key, and the accesible vertexes from
the array-list (value).

# Arguments
- `graph::Dict{Int64,Array{Int64,1}}`: Graph of the connected nodes
"""
function find_global_maximum(graph::Dict{Int64,Array{Int64,1}})
    global_maximum = 0
    # Going through the dictionary
    for (key, value) in graph
        # Check the keys as new global maximum
        if key > global_maximum
            global_maximum = key
        end

        # Check the values as new local maximum
        local_maximum = maximum(value)

        # Check if new local maximum becomes global maximum
        if local_maximum > global_maximum
            global_maximum = local_maximum
        end
    end
    return global_maximum
end

"""
    find_parent_in_cycle(parent::Array{Int64,1}, i::Int64, ref_key::Int64)

The find parent algorithm is a recursive function to find the subeset of an item i for
graph test of being cycling.

# Arguments
- `parent::Array{Int64,1}`: Array of the subset of the items
- `i::Int64`: Index of the parents-item
- `ref_key::Int64`: Reference Key, which can be for example the total size of the subset
"""
function find_parent_in_cycle(parent::Array{Int64,1}, i::Int64, ref_key::Int64)
    if parent[i] == ref_key
        return i
    end
    return find_parent_in_cycle(parent, parent[i], ref_key)
end

"""
    graph_cycle_check(graph::Dict{Int64,Array{Int64,1}})

The disjoint-set data structure principle is used to check if a direct or undirect
graph contains a cycle. For this reason, the algorithm keeps the first track of a set
of items partitioned into several disjoint (non-overlapping) subsets to find which
subset a particular item is kept. This procedure is essential to figure out if two
items are in the same subgroup. Next, the two subsets have to be merged into a single
subset. For more information see: [https://en.wikipedia.org/wiki/Disjoint-set_data_structure](https://en.wikipedia.org/wiki/Disjoint-set_data_structure)

# Arguments
- `graph::Dict{Int64,Array{Int64,1}}`: Graph of the connected nodes

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> graph_cycle_true = Dict(1 => [2], 2 => [3], 3 => [1, 4])
julia> ClassicAlgorithmsCollections.graph_cycle_check(graph_cycle_true)
true
julia> graph_cycle_false = Dict(1 => [2], 2 => [5], 3 => [1, 4])
julia> ClassicAlgorithmsCollections.graph_cycle_check(graph_cycle_false)
false
```
"""
function graph_cycle_check(graph::Dict{Int64,Array{Int64,1}})
    # Get the maximum length of the vertex
    vertex_size = find_global_maximum(graph)

    # Create array-list for finding the parents and intialize all subsets to vertex_size, because
    # negative indixing is forbidden in Julia.
    parent = fill(vertex_size, vertex_size)

    # Iterate through all edges of graph for finding the subset of both vertices of every edge.
    # If the two subsets are equal, then graph has a cycle inside.
    for (i, values) in graph
        for j in values
            x = find_parent_in_cycle(parent, i, vertex_size)
            y = find_parent_in_cycle(parent, j, vertex_size)

            # Check if equal, then graph has a cycle
            if x == y
                return true
            end

            x_set = find_parent_in_cycle(parent, x, vertex_size)
            y_set = find_parent_in_cycle(parent, y, vertex_size)
            # Update the parent-set
            parent[x_set] = y_set
        end
    end
    return false

end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# write your own tests here
N = 30
ph = discrete_phs(N, 0, 1);
eigval, eigvec = eigen(ph);
num_freq = frequencies(eigval)
exact_freq = pi/2*collect(1:2:2*N);
err = exact_freq - num_freq;

@test abs(err[2]) < 1e-13

@test abs(minimum(damp(ph)[:,1])-pi/2)<1e-13

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
module InterPol

export interpolateAt, interpolateDerivativeAt

# returns the value of the interpolation polynomial at the point x0
function interpolateAt{T<:Number,N}(x::Vector{T}, y::AbstractArray{T,N}, x0::T)

    if length(x)!=size(y,1)
        error("x and y have to be of the same size.")
    end

    n = length(x)
    p = zero(T)

    for i=1:n
        Li =one(T)
        dLi=zero(T)
        for j=1:n
            if j==i
                continue
            else
                Li*=(x0-x[j])/(x[i]-x[j])
            end
        end
        p.+=Li.*y[i, [1:size(y,n) for n=2:N]...]
    end
    return p
end

# returns the value of the derivative of the interpolation polynomial
# at the point x0
function interpolateDerivativeAt{T<:Number,N}(x::Vector{T}, y::AbstractArray{T,N}, x0::T)

    if length(x)!=size(y,1)
        error("x and y have to be of the same size.")
    end

    n = length(x)
    p = zero(T)

    for i=1:n
        dLi=zero(T)
        for k=1:n
            if k==i
                continue
            else
                dLi1=one(T)
                for j=1:n
                    if j==k || j==i
                        continue
                    else
                        dLi1*=(x0-x[j])/(x[i]-x[j])
                    end
                end
                dLi+=dLi1/(x[i]-x[k])
            end
        end
        p.+=dLi.*y[i, [1:size(y,n) for n=2:N]...]
        # p+=dLi*y[i]
    end
    return p
end

end # module


text:
using Base: tail

using Flux: applychain

"""
    @toNd f(x, y, z...; a=a, b=b, c=c...) n

macro for calling 2-d array function on N-d array by reshape input with reshape(x, size(x, 1), :)
and reshape back with reshape(out, :, input[n][2:end]...) where n is the n-th input(default=1).

"""
macro toNd(ex, outref::Int=1)
    fname = esc(ex.args[1])
    fkw = ex.args[2] isa Expr && ex.args[2].head == :parameters ? ex.args[2] : nothing
    _targs = Tuple(ex.args)
    fargs = esc.(fkw === nothing ? tail(_targs) : tail(tail(_targs)))
    fsize = map(x->Expr(:call, :size, x), fargs)
    rfargs = map((x, s) -> Expr(:call, :reshape, x, Expr(:ref, s, 1), :(:)), fargs, fsize)
    func = fkw === nothing ? Expr(:call, fname, rfargs...) : Expr(:call, fname, fkw, rfargs...)
    rsize = Expr(:call, :tail, fsize[outref])
    ret = Expr(:call, :reshape, func, :(:), Expr(:..., rsize))
    sT = gensym(:T)
    Expr(:(::), ret, Expr(:where,
                          Expr(:curly, :AbstractArray,
                               sT,
                               Expr(:call, :ndims, fargs[outref])),
                          sT)
         )
end

"""
    Positionwise(layers)

just like Flux.Chain, but reshape input to 2d and reshape back when output.
"""
struct Positionwise{T<:Tuple}
    models::T
    Positionwise(xs...) = new{typeof(xs)}(xs)
end

Flux.children(pw::Positionwise) = pw.models
Flux.mapchildren(f, pw::Positionwise) = Positionwise(f.(pw.models)...)

function (pw::Positionwise)(x)
    insize = size(x)
    y = applychain(pw.models, reshape(x, insize[1], :))
    reshape(y, :, Base.tail(insize)...)
end


text:
using Base.Test
# inputlayerstate.jl
import WaveletScattering: InputLayerState, kthrange, setup
# inputlayer.jl
import WaveletScattering: InputLayer

backend = Mocha.CPUBackend()
inputlayer = InputLayer(
        tops = [:signal],
        symbols = [:time, :chunk],
        data = rand(Float32, 256, 2))
inputlayerstate_constructor =
    InputLayerState(backend, inputlayer)
@test inputlayerstate_constructor.layer == inputlayer
inputlayerstate_setup = Mocha.setup(backend,
    inputlayer, Mocha.Blob[], Mocha.Blob[])
@test inputlayerstate_setup.layer == inputlayer


text:
## reader

function Base.parse(::Type{Module}, ir::String, ctx::Context=GlobalContext())
    data = unsafe_wrap(Vector{UInt8}, ir)
    membuf = MemoryBuffer(data, "", false)

    out_ref = Ref{API.LLVMModuleRef}()
    out_error = Ref{Cstring}()
    status = convert(Core.Bool, API.LLVMParseIRInContext(ref(ctx), ref(membuf), out_ref, out_error))

    if status
        error = unsafe_string(out_error[])
        API.LLVMDisposeMessage(out_error[])
        throw(LLVMException(error))
    end

    Module(out_ref[])
end

## writer

Base.convert(::Type{String}, mod::Module) = unsafe_string(API.LLVMPrintModuleToString(ref(mod)))


text:
using DifferentialEquations
using SimpleDiffEq
using DiffEqSensitivity
using Random
using Distributions
using DiffEqParamEstim
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R,C) = u
    (β,c,γ) = p
    N = S+I+R
    infection = β*c*I/N*S
    recovery = γ*I
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection
    end
    nothing
end;

δt = 1.0
tmax = 40.0
tspan = (0.0,tmax)
obstimes = 1.0:1.0:tmax;

u0 = [990.0,10.0,0.0,0.0]; # S,I.R,Y

p = [0.05,10.0,0.25]; # β,c,γ

prob_ode = ODEProblem(sir_ode!,u0,tspan,p)
sol_ode = solve(prob_ode,Tsit5(),saveat=δt);

out = Array(sol_ode)
C = out[4,:];

X = C[2:end] .- C[1:(end-1)];

Random.seed!(1234);

Y = rand.(Poisson.(X));

using Optim

function ss1(β)
    prob = remake(prob_ode,u0=[990.0,10.0,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),saveat=δt)
    out = Array(sol)
    C = out[4,:]
    X = C[2:end] .- C[1:(end-1)]
    nonpos = sum(X .<= 0)
    if nonpos > 0
        return Inf
    end
    return(sum((X .- Y) .^2))
end;

function nll1(β)
    prob = remake(prob_ode,u0=[990.0,10.0,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),saveat=δt)
    out = Array(sol)
    C = out[4,:]
    X = C[2:end] .- C[1:(end-1)]
    nonpos = sum(X .<= 0)
    if nonpos > 0
        return Inf
    end
    -sum(logpdf.(Poisson.(X),Y))
end;

lower1 = 0.0
upper1 = 1.0
initial_x1 = 0.1;

opt1_ss = Optim.optimize(ss1,lower1,upper1)

opt1_nll = Optim.optimize(nll1,lower1,upper1)

function ss2(x)
    (i0,β) = x
    I = i0*1000.0
    prob = remake(prob_ode,u0=[1000.0-I,I,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),saveat=δt)
    out = Array(sol)
    C = out[4,:]
    X = C[2:end] .- C[1:(end-1)]
    nonpos = sum(X .<= 0)
    if nonpos > 0
        return Inf
    end
    return(sum((X .- Y) .^2))
end;

function nll2(x)
    (i0,β) = x
    I = i0*1000.0
    prob = remake(prob_ode,u0=[1000.0-I,I,0.0,0.0],p=[β,10.0,0.25])
    sol = solve(prob,Tsit5(),saveat=δt)
    out = Array(sol)
    C = out[4,:]
    X = C[2:end] .- C[1:(end-1)]
    nonpos = sum(X .<= 0)
    if nonpos > 0
        return Inf
    end
    -sum(logpdf.(Poisson.(X),Y))
end;

lower2 = [0.0,0.0]
upper2 = [1.0,1.0]
initial_x2 = [0.01,0.1];

opt2_ss = Optim.optimize(ss2,lower2,upper2,initial_x2)

opt2_nll = Optim.optimize(nll2,lower2,upper2,initial_x2)

function loss_function(sol)
    out = Array(sol)
    C = out[4,:]
    X = C[2:end] .- C[1:(end-1)]
    nonpos = sum(X .<= 0)
    if nonpos > 0
        return Inf
    end
    -sum(logpdf.(Poisson.(X),Y))
end;

prob_generator = (prob,q) -> remake(prob,
    u0=[1000.0-(q[1]*1000),q[1]*1000,0.0,0.0],
    p=[q[2],10.0,0.25]);

cost_function = build_loss_objective(prob_ode,
    Tsit5(),
    loss_function,
    saveat=δt,
    prob_generator = prob_generator,
    maxiters=100,
    verbose=false);

opt_pe1 = Optim.optimize(cost_function,lower2,upper2,initial_x2)

using NLopt
opt = Opt(:LD_MMA, 2)
opt.lower_bounds = lower2
opt.upper_bounds = upper2
opt.min_objective = cost_function
opt.maxeval = 10000
(minf,minx,ret) = NLopt.optimize(opt,initial_x2)

using BlackBoxOptim
bound1 = Tuple{Float64, Float64}[(0.0,1.0),(0.0, 1.0)]
result = bboptimize(cost_function;SearchRange = bound1, MaxSteps = 1e4)


text:
using Statistics
using Elefridge

i = [0,0,0,0];

round_better(A::Array{Float32}) = abs(mean(A-round(A,7))) <= abs(mean(A-groom(A,7)))

for j in 1:100
    A = 4*rand(Float32,10000000)
    i[1] += round_better(A) ? 1 : 0
    A = A .+ 2
    i[2] += round_better(A) ? 1 : 0
    A = randn(Float32,10000000)
    i[3] += round_better(A) ? 1 : 0
    A = A .+ 2
    i[4] += round_better(A) ? 1 : 0
 end


text:
#for now just the main functions, no wrapper
#remove all unnecessary allocations
mutable struct PowerBuffer #one for each *power calculation* (can be > 1 for each frequency/incident wave)
    iₚ::Integer #which problem (β) is this??
    Ez::Vector{Complex{Float64}}
    Hx::Vector{Complex{Float64}}
    Hy::Vector{Complex{Float64}}
    Ez_inc::Vector{Complex{Float64}}
    Hx_inc::Vector{Complex{Float64}}
    Hy_inc::Vector{Complex{Float64}}
    pow::Float64 #power
    HEz::Vector{Vector{Complex{Float64}}}
    HHx::Vector{Vector{Complex{Float64}}}
    HHy::Vector{Vector{Complex{Float64}}}
    nhat::Array{Float64,2}
    S::Array{Complex{Float64},2}
    ∂pow::Vector{Complex{Float64}} # ∂power/∂βᵀ
    l::Float64 #arc length
    PowerBuffer(Ns::Integer, P::Integer, nhat, iₚ) =
            new(iₚ,
                Array{Complex{Float64}}(undef, size(nhat, 1)),
                Array{Complex{Float64}}(undef, size(nhat, 1)),
                Array{Complex{Float64}}(undef, size(nhat, 1)),
                Array{Complex{Float64}}(undef, size(nhat, 1)),
                Array{Complex{Float64}}(undef, size(nhat, 1)),
                Array{Complex{Float64}}(undef, size(nhat, 1)),
                0.0,
                [Array{Complex{Float64}}(undef, Ns*(2P+1)) for i=1:size(nhat, 1)],
                [Array{Complex{Float64}}(undef, Ns*(2P+1)) for i=1:size(nhat, 1)],
                [Array{Complex{Float64}}(undef, Ns*(2P+1)) for i=1:size(nhat, 1)],
                nhat,
                Array{Complex{Float64},2}(undef, size(nhat)),
                Array{Complex{Float64}}(undef, Ns*(2P+1)),
                0.0)
end

#TODO: take k0,ui from opb
function optMatrixPwr(points, centers, Ns, P, k0, ui, iₚ, nhat, l)
    sv = PowerBuffer(Ns, P, nhat, iₚ)
    cf = 1/(1im*k0*eta0)
    pt = Array{Float64}(undef, 2)
    for ip = 1:size(points,1)
        for ic = 1:Ns
            pt[1] = points[ip,1] - centers[ic,1]
            pt[2] = points[ip,2] - centers[ic,2]
            θ = atan(pt[2], pt[1])
            R = hypot(pt[1], pt[2])

            ind = (ic-1)*(2*P+1) .+ P + 1
            Hₚ₋₂ = besselh(-1, 1, k0*R)
            Hₚ₋₁ = besselh(0, 1, k0*R)

            sv.HEz[ip][ind] = Hₚ₋₁
            ∂H∂R = k0*Hₚ₋₂
            sv.HHx[ip][ind] = cf*∂H∂R*pt[2]/R
            sv.HHy[ip][ind] = -cf*∂H∂R*pt[1]/R
            for p = 1:P
                Hₚ = (2*(p-1)/(k0*R))*Hₚ₋₁ - Hₚ₋₂
                ∂H∂R = k0*(Hₚ₋₁ - (p/k0/R)*Hₚ)
                sv.HEz[ip][ind - p] = (-1)^p*Hₚ*exp(-1im*p*θ)
                sv.HHy[ip][ind - p] = -cf*exp(-1im*p*θ)*(-1)^p*(∂H∂R*pt[1] - Hₚ*1im*p*(-pt[2])/R)/R
                sv.HHx[ip][ind - p] = cf*exp(-1im*p*θ)*(-1)^p*(∂H∂R*pt[2] - Hₚ*1im*p*(pt[1])/R)/R

                sv.HEz[ip][ind + p] = Hₚ*exp(1im*p*θ)
                sv.HHy[ip][ind + p] = -cf*exp(1im*p*θ)*(∂H∂R*pt[1] + Hₚ*1im*p*(-pt[2])/R)/R
                sv.HHx[ip][ind + p] = cf*exp(1im*p*θ)*(∂H∂R*pt[2] + Hₚ*1im*p*(pt[1])/R)/R
                Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ
            end
        end
    end
    sv.l = l
    sv.Ez_inc = uinc(k0, points, ui)
    sv.Hx_inc = hxinc(k0, points, ui)
    sv.Hy_inc = hyinc(k0, points, ui)
    sv
end

mutable struct OptimProblemBuffer#one for each problem, e.g. 3 if same incident wave but 3 different frequencies, or 2 if same frequency but different incident wave
    k0::Float64
    kin::Float64
    α::Vector{Complex{Float64}}
    ui::Einc
    β::Vector{Complex{Float64}}
    rhs_grad::Vector{Complex{Float64}}
    λadj::Vector{Complex{Float64}}
end

function OptimProblemBuffer(k0::Float64, kin::Float64, centers::Array{Float64,2}, ui::Einc, P::Integer)
    Ns = size(centers, 1)
    α = ParticleScattering.u2α(k0, ui, centers, P)
    OptimProblemBuffer(k0, kin, α, ui,
        Array{Complex{Float64}}(undef, Ns*(2*P+1)),
        Array{Complex{Float64}}(undef, Ns*(2*P+1)),
        Array{Complex{Float64}}(undef, Ns*(2*P+1)))
end

function optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM)
    if rs != last_rs
        copyto!(last_rs, rs)
        #do whatever common calculations and save to shared_var
        for i in eachindex(opb)
            #TODO: if multiple problems have same k0, these calcs are duplicated
            for id in unique(ids)
                ParticleScattering.updateCircleScatteringDerivative!(scatteringMatrices[i][id], dS_S[i][id], opb[i].k0, opb[i].kin, rs[id], P)
            end
            for ic = 1:Ns
                rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
                fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*opb[i].α[rng]
            end

            MVP = LinearMap{eltype(fmmbuf.rhs)}(
                        (output_, x_) -> ParticleScattering.FMM_mainMVP_pre!(output_, x_,
                                            scatteringMatrices[i], φs, ids, P, mFMM[i],
                                            fmmbuf.pre_agg, fmmbuf.trans),
                        Ns*(2*P+1), Ns*(2*P+1), ismutating = true)

            opb[i].β[:] .= 0
            opb[i].β,ch = gmres!(opb[i].β, MVP, fmmbuf.rhs,
                                    restart = Ns*(2*P+1), tol = fmmopts.tol,
                                    log = true, initially_zero = true) #no restart
            !ch.isconverged && error("""FMM process did not converge, normalized
                residual: $(norm(MVP*opb[i].β - fmmbuf.rhs)/norm(fmmbuf.rhs))""")
        end
        #calculate power for each power calculation arc
        calc_multi_pwr!(power_buffer, opb)
    end
end

function optimize_pwr_rs_f(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM, fobj_pwr)
    optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM)

    fobj_pwr(power_buffer)
end

function optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM, gobj_pwr!)
    optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM)

    fill!(grad_stor, 0) #gradient is the sum of both adjoint gradients
    #compute all ∂P/∂βᵀ
    dPdβ_pwr!.(power_buffer)
    #build rhs of all adjoint problems (-∂f/∂βᵀ)
    gobj_pwr!(power_buffer, opb)

    #for each problem, solve adjoint problem and add to total gradient
    for i in eachindex(opb)
        MVP = LinearMap{eltype(fmmbuf.rhs)}(
                (output_, x_) -> ParticleScattering.FMM_mainMVP_transpose!(output_, x_,
                                        scatteringMatrices[i], φs, ids, P, mFMM[i],
                                        fmmbuf.pre_agg, fmmbuf.trans),
                Ns*(2*P+1), Ns*(2*P+1), ismutating = true)

        #solve adjoint problem
        opb[i].λadj[:] .= 0
        opb[i].λadj, ch = gmres!(opb[i].λadj, MVP, opb[i].rhs_grad,
                            restart = Ns*(2*P+1), tol = fmmopts.tol, log = true,
                            initially_zero = true) #no restart
        for n = 1:length(rs)
            #compute n-th gradient - here we must pay the price for symmetry
            #as more than one β is affected. Overwrites rhs_grad (ok because
            #there is seperate one for each i). TODO:
            #rhs_grad is still usually sparse - utilize this to reduce complexity
            #here O(N^2) -> O(N)
            for ic = 1:Ns
                rng = (ic-1)*(2*P+1) .+ (1:2*P+1)
                if ids[ic] == n
                    opb[i].rhs_grad[rng] = dS_S[i][n]*opb[i].β[rng]
                else
                    opb[i].rhs_grad[rng] .= 0.0
                end
            end
            grad_stor[n] += -2*real(transpose(opb[i].λadj)*opb[i].rhs_grad)
        end
    end
end

function calc_multi_pwr!(power_buffer, opb)
    for sv in power_buffer
        len = length(sv.HEz)
        Sn = Array{Complex{Float64}}(undef, len) # S⋅n
        βi = opb[sv.iₚ].β
        for ip = 1:len
            sv.Ez[ip] = transpose(sv.HEz[ip])*βi + sv.Ez_inc[ip]
            sv.Hx[ip] = transpose(sv.HHx[ip])*βi + sv.Hx_inc[ip]
            sv.Hy[ip] = transpose(sv.HHy[ip])*βi + sv.Hy_inc[ip]
            Sn[ip] = -0.5*real(sv.Ez[ip]*conj(sv.Hy[ip]))*sv.nhat[ip,1]
            Sn[ip] += 0.5*real(sv.Ez[ip]*conj(sv.Hx[ip]))*sv.nhat[ip,2]
        end
        sv.pow = (sum(Sn) - 0.5*Sn[1] - 0.5*Sn[end])*sv.l/(len-1) #trapezoidal rule, correct up to arc length for simple integrals
    end
end

function dPdβ_pwr!(sv::PowerBuffer)
    len = length(sv.Ez) #num of points
    fill!(sv.∂pow, 0.0)
    #this is a sum of real and imaginary parts, hence the additional 1/2 factor
    for ip = 1:len
        cf = ifelse(ip == 1 || ip == len, 0.5, 1.0)*sv.l/(len-1) #trapezoidal rule constant
        sv.∂pow .+= (-0.25*cf*sv.nhat[ip,1])*(conj(sv.Hy[ip])*sv.HEz[ip] +
                                    conj(sv.Ez[ip])*sv.HHy[ip])
        sv.∂pow .+= (0.25*cf*sv.nhat[ip,2])*(conj(sv.Hx[ip])*sv.HEz[ip] +
                                    conj(sv.Ez[ip])*sv.HHx[ip])
    end
end


text:
using Plots, LaTeXStrings, JLD, Statistics
TotalTime=1000
probability=[0.25,0.5,0.75,0.95]
function RandomWalk(TotalTime, probability)
    arr=zeros((TotalTime,length(probability)))
    for j in 1:length(probability)
        for i in 1:TotalTime
            if probability[j]>rand()
                arr[i,j]=1
            else
                arr[i,j]=-1
            end
        end
    end
    return cumsum(arr,dims=1)
end
#TheMatrix=RandomWalk(TotalTime, probability)
run_num=10000
# a list cinsisting of lists to the number of probabilities,
#each list having lists of each position in each run number per time:
RandomWalkers=[[[] for j in 1:TotalTime] for i in 1:length(probability)]
for n in 1:run_num
    TheMatrix=RandomWalk(TotalTime, probability)
    for col in 1:length(probability)
        for i in 1:TotalTime
            push!(RandomWalkers[col][i],TheMatrix[i,col])
        end
    end
end
SavedData_mean=[[] for i in 1:length(probability)]
SavedData_var=[[] for i in 1:length(probability)]
for i in 1:length(RandomWalkers)
    for j in RandomWalkers[i]
        push!(SavedData_var[i],var(j))
        push!(SavedData_mean[i],mean(j))
    end
end
save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.2\\chapter5_4.2_MeanAndStdData.jld",
 "MeanData", SavedData_mean, "VarData", SavedData_var)
 #plots:
load("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.2\\chapter5_4.2_MeanAndStdData.jld")
scatter(dpi=400)
for i in 1:length(probability)
    p=probability[i]
    scatter!(1:TotalTime, SavedData_mean[i],label=L"%$p",markersize=1,markerstrokewidth=0, alpha=0.7)
end
scatter!(title="<x> per time for different probabilities", xlabel="t", ylabel="<x>")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.2\\Figs\\MeanPlot.png")
scatter(dpi=400)
for i in 1:length(probability)
    p=probability[i]
    scatter!(1:TotalTime, SavedData_var[i],label=L"%$p",markersize=1,markerstrokewidth=0, alpha=0.7)
end
scatter!(title="Var[x] per time  for different probabilities", xlabel="t", ylabel="Var[x]")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.2\\Figs\\VarPlot.png")


text:
using SoleViz
using Documenter

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

makedocs(;
    modules=[SoleViz],
    authors="Eduard I. STAN, Giovanni PAGLIARINI, Federico MANZELLA",
    repo="https://github.com/aclai-lab/SoleViz.jl/blob/{commit}{path}#{line}",
    sitename="SoleViz.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://aclai-lab.github.io/SoleViz.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/aclai-lab/SoleViz.jl",
    devbranch="main",
)

max_stars_repo_name: aclai-lab/SoleVisualizations

max_issues_repo_name: aclai-lab/SoleVisualizations

max_forks_repo_name: aclai-lab/SoleVisualizations


text:
using DataFrames

function df2mm(df::DataFrame)
    n = size(df, 1)
    mm_raw = [fill(1.0, n, 1)]
    mm_name = ["const"]
    for (name, value) in eachcol(df)
        if eltype(value) <: Real
            push!(mm_raw, hcat(Float64.(value)))
            push!(mm_name, string(name))
        else
            uvalue = unique(value)
            length(uvalue) == 1 && continue
            dvalue = Dict(v=>i for (i, v) in enumerate(uvalue))
            mvalue = zeros(n, length(uvalue))
            for i in 1:n
                mvalue[i, dvalue[value[i]]] = 1.0
            end
            push!(mm_raw, mvalue[:, 2:end])
            append!(mm_name, string.(name, "_", uvalue[2:end]))
        end
    end
    (data=hcat(mm_raw...), names=mm_name)
end

function lm(df, y, xs)
    yv = Float64.(df[y])
    xv, xn = df2mm(df[[xs;]])
    params = (transpose(xv)*xv)\(transpose(xv)*yv)
    DataFrame(name = xn, estimate=params)
end

max_stars_repo_path: Chapter04/09. Estimating linear regression/lm.jl


text:
__precompile__()

module BooleanMatrixFactorization

using Parameters
using Distributions
using StatsFuns

export BeliefPropagationOptions, ProblemSpecification, makeData, makeMask,
    inference, experiment

"""
    BeliefPropagationOptions

Pararmeters define the operation of belief propagation.
"""
@with_kw immutable BeliefPropagationOptions
    learningRate = 0.5     # damping parameter
    maxIter = 200          # max number of message updates
    verbose = true
    tol = 1e-6         # convergence criterion on message diffs
    minSum = true          # true for minSum, false for sumProd
    @assert tol > 0.0
    @assert learningRate > 0.0
    @assert maxIter > 0
end

"""
    ProblemSpecification

Structure representing te parameters defining a factorization problem
"""
@with_kw immutable ProblemSpecification
    xor = false                                # is problem xor or and
    K = 1                                      # desired rank
    PX1 = 0.5                                  # prior probability of X[i,k] = 1
    PY1 = 0.5                                  # prior probability of Y[j,k] = 1
    PO1GivenZ1 = 0.99                          # noise model, P(O[i,j]=1|Z[i,j]=1)
    PO0GivenZ0 = 0.99                          # noise model, P(O[i,j]=0|Z[i,j]=0)
    Obs = Array{Bool}(0,0)                     # observation matrix
    mask = Array{Bool}(0,0)                    # only use observation O[i,j] where mask[i,j]==true
    numObs = count(mask)                       # number of observations in O
    posObs = Obs[mask]                         # positive observations
    negObs = Vector{Bool}(.!(posObs))          # negative observations
    logOddsX = log(PX1) - log(1-PX1)           # log odds for X[i,j]=1
    logOddsY = log(PY1) - log(1-PY1)           # log odds for X[i,j]=1
    co0 = log(1-PO1GivenZ1) - log(PO0GivenZ0)  # log(P(O[i,j]=0|Z[i,j]=1)/P(O[i,j]=0|Z[i,j]=0))
    co1 = log(PO1GivenZ1) - log(1-PO0GivenZ0)  # log(P(O[i,j]=1|Z[i,j]=1)/P(O[i,j]=1|Z[i,j]=0))
    @assert 0 < K
    @assert 0.0 ≤ PX1 ≤ 1.0
    @assert 0.0 ≤ PY1 ≤ 1.0
    @assert 0.0 ≤ PO1GivenZ1 ≤ 1.0
    @assert 0.0 ≤ PO0GivenZ0 ≤ 1.0
    @assert size(Obs) == size(mask)
end

"""
    setData(probSPec::ProblemSpec,O,m=trues(O))

Set the data and mask in probSpec to O and m.
"""
function setData(probSpec::ProblemSpecification, O, m=trues(O), K=probSpec.K)
    size(O) ≠ size(m) && error("Observation and mask matrices must be the same size")
    (any(K .> size(O)) || (K<0)) && error("require 0 < K ≤ size(O)")
    posObs = O[m]
    ProblemSpecification(probSpec; Obs=O, mask=m, K=K, posObs=posObs, negObs=Vector{Bool}(.!(posObs)))
end

"""
    X,Y,Z,O = makeData(M,N,probSpec)

Generate test data X,Y,Z,O such that O is a noisy version of Z = X*Y'. Z is M×N
"""
function makeData(M, N, probSpec::ProblemSpecification)
    X = Matrix{Bool}(rand(M, probSpec.K) .< probSpec.PX1)
    Y = Matrix{Bool}(rand(M, probSpec.K) .< probSpec.PY1)
    Z = Matrix{Bool}(probSpec.xor ? ((X * Y' % 2) .== 1) :  (X * Y' .> 0))
    # apply the noise model
    B0 = Bernoulli(1.0-probSpec.PO0GivenZ0) # Prob(O=1|Z=0)
    B1 = Bernoulli(probSpec.PO1GivenZ1)     # Prob(O=1|Z=1)
    O = similar(Z)
    for (i,z) ∈ enumerate(Z)
        O[i] = z ? rand(B1) : rand(B0)
    end
    X, Y, Z, O
end

"""
    makeMask(M,N, pObserved=1.0, permMatrix=false)

Generate a random mask of observations within an M×N observation matrix
"""
makeMask(M, N, pObserved=1.0) = Matrix{Bool}(rand(M,N) .< pObserved)
makeMask(O::Matrix{Bool}, pObserved=1.0) = makeMask(size(O,1),size(O,2),pObserved)

phi(x::Float64) = max(0, x)
psi(x::Float64) = log1pexp.(x)

function getAdjList(A, bip_connections = false)# A is m x n, nbM: list of nonzeros in each row
    (M, N) = findn(A);  m,n = size(A)
    nbM = [Int64[] for i=1:m];  nbN = [Int64[] for i=1:n]
    for z=1:length(N)
        push!(nbN[N[z]], z);  push!(nbM[M[z]], z)
    end
    if bip_connections
        nbMN = [Int64[] for i=1:m];  nbNM = [Int64[] for i=1:n]
        for z=1:length(N)
            push!(nbNM[N[z]], M[z]);  push!(nbMN[M[z]], N[z])
        end
        nbM, nbN, nbMN, nbNM
    else
        nbM, nbN
    end
end

"""
    updateMinSum!(outX,outY,inZ,outZ,newInX,newInY,ps::ProblemSpecification)
"""
function updateMinSum!(outX::Array{Float64, 2}, outY::Array{Float64, 2}, inZ::Array{Float64, 2},
                       outZ::Array{Float64, 2}, newInX::Array{Float64, 2}, newInY::Array{Float64, 2},
                       probSpec::ProblemSpecification)
    inZ[:,:] = min.(outX + outY, outX, outY)
    inZpos = phi.(deepcopy(inZ))
    inZ_max, inZ_maxind = findmax(inZ, 2)   # max over k
    inZ[inZ_maxind[:]] = -Inf
    inZ_max_sec = max.(-maximum(inZ, 2), 0) # new
    inZ_max = phi.(-inZ_max)
    sumval = sum(inZpos, 2)
    sumval[probSpec.posObs,1] += probSpec.co1
    sumval[probSpec.negObs,1] += probSpec.co0
    inZ_maxind_k = map(x -> ind2sub(size(outX), x)[2], inZ_maxind)
    tmp_inZ_max = deepcopy(inZ_max)
    inZpos = -inZpos .+ sumval
    for k = 1:size(outZ,2)
        self_maxind = (inZ_maxind_k .== k)
        tmp_inZ_max[self_maxind] = inZ_max_sec[self_maxind]
        outZ[:, k] = min.( tmp_inZ_max[:], inZpos[:,k]) #sumval_edge[:] - inZpos[posObs,K])
        tmp_inZ_max[self_maxind] = inZ_max[self_maxind]
    end
    newInX[:,:] = phi.(outZ + outY) - phi.(outY)
    newInY[:,:] = phi.(outZ + outX) - phi.(outX)
end

"""
    updateMinsumNoiseless(outX, outY, inZ, outZ, newInX, newInY, probSpec)
"""
function updateMinSumNoiseLess!(outX::Array{Float64, 2}, outY::Array{Float64, 2}, inZ::Array{Float64, 2},
                                outZ::Array{Float64, 2}, newInX::Array{Float64, 2}, newInY::Array{Float64, 2},
                                prodSpec::ProblemSpecification)
    inZ[:,:] = min(outX[probSpec.posObs,:]+outY[probSpec.posObs,:], outX[probSpec.posObs,:], outY[probSpec.posObs,:])
    inZ_max, inZ_maxind = findmax(inZ, 2)
    inZ[inZ_maxind[:]] = -Inf
    inZ_max_sec = phi(-maximum(inZ, 2))
    inZ_max = phi(-inZ_max)
    inZ_maxind_k = map(x -> ind2sub(size(outZ), x)[2], inZ_maxind)
    tmp_inZ_max = deepcopy(inZ_max)
    for k = 1:size(outZ,2)
        self_maxind = (inZ_maxind_k .== k)
        tmp_inZ_max[self_maxind] = inZ_max_sec[self_maxind]
        outZ[:, k] = tmp_inZ_max[:]
        tmp_inZ_max[self_maxind] = inZ_max[self_maxind]
    end
    newInY[:,:] = -phi(outX)
    newInX[:,:] = -phi(outY)
    newInX[probSpec.posObs,:] += phi(outZ + outY[probSpec.posObs,:])
    newInY[probSpec.posObs,:] += phi(outZ + outX[probSpec.posObs,:])
end

"""
    updateSumProd!(outX, outY, inZ, outZ, newInX, newInY, probSpec::ProblemSpecification)
"""
function updateSumProd!(outX::Array{Float64, 2}, outY::Array{Float64, 2}, inZ::Array{Float64, 2},
                        outZ::Array{Float64, 2}, newInX::Array{Float64, 2}, newInY::Array{Float64, 2},
                        probSpec::ProblemSpecification)
    inZ[:,:] = psi(outX + outY - log(1 + exp(outX) + exp(outY)))
    sumvalLog = sum(inZ, 2)
    inZ[:,:] = -inZ .+ sumvalLog
    outZ[:,:] = exp(inZ) - 1
    inZ[probSpec.posObs,:] += probSpec.co1
    inZ[probSpec.negObs,:] += probSpec.co0
    outZ[probSpec.posObs,:] *= exp(probSpec.co1)
    outZ[probSpec.negObs,:] *= exp(probSpec.co0)
    outZ[:,:] = inZ - log(1 + outZ)
    newInX[:,:]= psi(outZ + outY) - psi(outY)
    newInY[:,:] = psi(outZ + outX) - psi(outX)
end

"""
    updateMargs!()
"""
function updateMargs!(margX::Array{Float64, 2}, margY::Array{Float64, 2}, outX::Array{Float64, 2},
                      outY::Array{Float64, 2}, inX::Array{Float64, 2}, inY::Array{Float64, 2},
                      nbM::Vector{Vector{Int64}}, nbN::Vector{Vector{Int64}}, probSpec::ProblemSpecification)
    for M=1:length(nbM)
        margX[M,:] = sum(inX[nbM[M],:], 1) .+ probSpec.logOddsX
        outX[nbM[M], :] = -inX[nbM[M],:] .+ margX[M,:]'
    end
    for N=1:length(nbN)
        margY[N,:] = sum(inY[nbN[N],:], 1) .+ probSpec.logOddsY
        outY[nbN[N], :] = -inY[nbN[N],:] .+ margY[N,:]'
    end
end

"""
    inference(O,mask,)
"""
function inference(probSpec::ProblemSpecification, bpOpts::BeliefPropagationOptions)
    O, mask, K = probSpec.Obs, probSpec.mask, probSpec.K
    M, N = size(O);  nObs = probSpec.numObs
    nObs = countnz(mask)   # number of observations in O
    nbM, nbN = getAdjList(mask)
    posObs = O[mask];  negObs = Vector{Bool}(.!(posObs))
    margX = zeros(M, K);  margY = zeros(N, K)
    localX = probSpec.logOddsX + log.(rand(M,K))
    localY = probSpec.logOddsY + log.(rand(N,K))
    inX = zeros(nObs, K); outX = log.(rand(nObs, K))
    newInX = zeros(nObs, K);  newInY = zeros(nObs, K)
    inY = zeros(nObs, K);  outY = log.(rand(nObs, K))
    inZ = zeros(nObs, K);  outZ = zeros(nObs, K)
    iters = 1;  diffMsg = Inf
    while (diffMsg > bpOpts.tol) && iters ≤ bpOpts.maxIter
        if bpOpts.minSum
            updateMinSum!(outX, outY, inZ, outZ, newInX, newInY, probSpec)
        else
            updateSumProd!(outX, outY, inZ, outZ, newInX, newInY, probSpec)
        end
        diffMsg = maximum(abs,newInX - inX)
        inX .= (1. - bpOpts.learningRate)*inX .+ bpOpts.learningRate*newInX
        inY .= (1. - bpOpts.learningRate)*inY .+ bpOpts.learningRate*newInY
        updateMargs!(margX, margY, outX, outY, inX, inY, nbM, nbN, probSpec)
        bpOpts.verbose && println("iter $iters, diff:$diffMsg")
        iters += 1
    end
    X = Matrix{Bool}(margX .> 0);  Y = Matrix{Bool}(margY .> 0)
    X, Y, Matrix{Bool}(probSpec.xor ? (X * Y' % 2) .== 1 : X * Y' .> 0)
end

"""
    compare(A,B,mask)

Compare two Boolean matrices A and B by reporting the fraction of elements
differing both across all elements, and across those flagged as true within mask.
"""
function compare(A, B, mask=trues(A))
    any(size(A)≠size(B)) && error("Matrices must be the same size")
    d1, d2, obs = countnz(A)/length(A), countnz(B)/length(B), countnz(mask)/length(mask)
    println("densities: $d1,$d2    observations: $obs")
    dist = countnz(A .!== B) / length(A)
    distObs = countnz(A[mask] .!== B[mask]) / sum(mask)
    println("diff (all elements): $dist,  diff (on observations): $distObs")
    dist, distObs
end

"""
    experiment(M,N)

Run an experiment generating random test data, running BP and comparing inference
to  test data.
"""
function experiment(M,N,k=1)
    srand(1234)
    probSpec = ProblemSpecification(K=k)
    X, Y, Z, O = makeData(M, N, probSpec)
    pObserved = 0.99
    mask = makeMask(O, pObserved)
    probSpec = setData(probSpec, O, mask, k)
    bpOpts = BeliefPropagationOptions()
    Xh,Yh,Zh = inference(probSpec, bpOpts)
    println("comparing Zh and O");  compare(Zh, O, mask)
    println("comparing X and Xh");  compare(X, Xh)
    print("X: ");  show(STDOUT, "text/plain", X);  println("")
    print("Xh: ");  show(STDOUT, "text/plain", Xh);  println("")
    println("comparing Y and Yh");  compare(Y, Yh)
    print("Y: ");  show(STDOUT, "text/plain", Y);  println("")
    print("Yh: ");  show(STDOUT, "text/plain", Yh);  println("")
    println("comparing Z and Zh");  compare(Z, Zh, mask)
    print("Z: ");  show(STDOUT, "text/plain", Z);  println("")
    print("Zh: ");  show(STDOUT, "text/plain", Zh);  println("")
    O, Z, Zh, mask
end

end # module BooleanMatrixFactorization


text:
const AllowedDateTimeField = Union{BLPDateTime, DateTime, Date}

"""
    bdh_intraday_ticks(session, security, event_types, date_start, date_end;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        )

Runs a `IntradayTickRequest` request.

# Example

```julia
# use a recent date
d0 = DateTime(2022, 3, 22, 13)
d1 = DateTime(2022, 3, 22, 13, 5)
res = BLPData.bdh_intraday_ticks(session, "PETR4 BS Equity", ["TRADE", "BID", "ASK"], d0, d1)
df = DataFrame(res)
show(df)
```
"""
function bdh_intraday_ticks(session::Session, security::AbstractString, event_types::Vector{T}, date_start::T1, date_end::T2;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        ) where {T<:AbstractString, T1<:AllowedDateTimeField, T2<:AllowedDateTimeField}

    @assert !isempty(event_types) "Event types vector is empty."

    queue, corr_id = send_request(session, "//blp/refdata", "IntradayTickRequest") do req
        req["security"] = security
        req["startDateTime"] = date_start
        req["endDateTime"] = date_end
        append!(req["eventTypes"], event_types)

        if options != nothing
            for (k, v) in options
                req[k] = v
            end
        end
    end

    result = Vector()

    for_each_response_message_element(queue, corr_id, timeout_milliseconds=timeout_milliseconds, verbose=verbose) do element
        @assert has_name(element, "IntradayTickResponse")
        response_element = get_choice(element)

        if has_name(response_element, "responseError")
            error("Got responseError. \n$response_element")
        end

        @assert has_name(response_element, "tickData")
        tick_data_array = response_element["tickData"]
        @assert isa(tick_data_array, Element{true, BLPAPI_DATATYPE_SEQUENCE})
        push_named_tuples!(result, tick_data_array)
    end

    return result
end

function bdh_intraday_ticks(session::Session, securities::Vector{T1}, event_types::Vector{T2}, date_start::DT1, date_end::DT2;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        ) where {T1<:AbstractString, T2<:AbstractString, DT1<:AllowedDateTimeField, DT2<:AllowedDateTimeField}

    result = Dict()

    @sync for security in securities
        @async result[security] = bdh_intraday_ticks($session, $security, $event_types, $date_start, $date_end, options=$options, verbose=$verbose, timeout_milliseconds=$timeout_milliseconds)
    end

    return result
end

function bdh_intraday_ticks(session::Session, security::AbstractString, event_type::AbstractString, date_start::T1, date_end::T2;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        ) where {T1<:AllowedDateTimeField, T2<:AllowedDateTimeField}

    bdh_intraday_ticks(session, security, [event_type], date_start, date_end, options=options, verbose=verbose, timeout_milliseconds=timeout_milliseconds)
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
using DiffEqSensitivity, Flux, Optimization, OptimizationFlux, OptimizationOptimJL, OrdinaryDiffEq, Test

x = Float32[0.8; 0.8]
tspan = (0.0f0,10.0f0)

ann = Chain(Dense(2,10,tanh), Dense(10,1))
p = Float32[-2.0,1.1]
p2,re = Flux.destructure(ann)
_p = [p;p2]
θ = [x;_p]

function dudt2_(u,p,t)
    x, y = u
    [(re(p[3:end])(u)[1]),p[1]*y + p[2]*x]
end

prob = ODEProblem(dudt2_,x,tspan,_p)
solve(prob,Tsit5())

function predict_rd(θ)
  Array(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:end],abstol=1e-7,reltol=1e-5,sensealg=TrackerAdjoint()))
end
loss_rd(p) = sum(abs2,x-1 for x in predict_rd(p))
l = loss_rd(θ)

cb = function (θ,l)
  @show l
  # display(plot(solve(remake(prob,u0=Flux.data(_x),p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
cb(θ,l)

loss1 = loss_rd(θ)
optfunc = Optimization.OptimizationFunction((x, p) -> loss_rd(x), Optimization.AutoZygote())
optprob = Optimization.OptimizationProblem(optfunc, θ)
res = Optimization.solve(optprob, BFGS(initial_stepnorm = 0.01), callback = cb)
loss2 = res.minimum
@test 3loss2 < loss1

## Partial Neural Adjoint

u0 = Float32[0.8; 0.8]
tspan = (0.0f0,25.0f0)

ann = Chain(Dense(2,10,tanh), Dense(10,1))

p1,re = Flux.destructure(ann)
p2 = Float32[-2.0,1.1]
p3 = [p1;p2]
θ = [u0;p3]

function dudt_(du,u,p,t)
    x, y = u
    du[1] = re(p[1:41])(u)[1]
    du[2] = p[end-1]*y + p[end]*x
end
prob = ODEProblem(dudt_,u0,tspan,p3)
solve(prob,Tsit5(),abstol=1e-8,reltol=1e-6)

function predict_adjoint(θ)
  Array(solve(prob,Tsit5(),u0=θ[1:2],p=θ[3:end],saveat=0.0:1:25.0))
end
loss_adjoint(θ) = sum(abs2,x-1 for x in predict_adjoint(θ))
l = loss_adjoint(θ)

cb = function (θ,l)
  @show l
  # display(plot(solve(remake(prob,p=Flux.data(p3),u0=Flux.data(u0)),Tsit5(),saveat=0.1),ylim=(0,6)))
  false
end

# Display the ODE with the current parameter values.
cb(θ,l)

loss1 = loss_adjoint(θ)
optfunc = Optimization.OptimizationFunction((x, p) -> loss_adjoint(x), Optimization.AutoZygote())
optprob = Optimization.OptimizationProblem(optfunc, θ)
res1 = Optimization.solve(optprob, ADAM(0.01), callback = cb, maxiters = 100)

optprob = Optimization.OptimizationProblem(optfunc, res1.minimizer)
res = Optimization.solve(optprob, BFGS(initial_stepnorm = 0.01), callback = cb)
loss2 = res.minimum
@test 3loss2 < loss1


text: # Common functions go here


text:
import SortingAlgorithms:uint_mapping
using SortingAlgorithms, StrFs

#export uint_mapping, radixsort

SortingAlgorithms.uint_mapping(::Base.Order.ReverseOrdering, s::StrF{S}) where S = begin
    res = s |> Ref |> pointer_from_objref |> Ptr{UInt128} |> unsafe_load

    # this might contain rubbish so need to clear them out
    ~(((res << (8*(16-sizeof(s)))) >> (8*(16-sizeof(s)))) |> ntoh)
end

SortingAlgorithms.uint_mapping(::Base.Order.ForwardOrdering, s::StrF{S}) where S  = begin
    res = s |> Ref |> pointer_from_objref |> Ptr{UInt128} |> unsafe_load
    ((res << (8*(16-sizeof(s)))) >> (8*(16-sizeof(s)))) |> ntoh
end

radixsort(v::AbstractVector{StrF{S}}; rev = false) where S = begin
    if S <=16
        return sort(v, alg = RadixSort, rev = rev)
    else
        throw(ErrorException("radix sort of Vector{StrF(S)} is only possible for S <= 16 for now"))
    end
end

radixsort!(v::AbstractVector{StrF{S}}; rev = false) where S = begin
    if S <=16
        return sort!(v, alg = RadixSort, rev = rev)
    else
        throw(ErrorException("radix sort of Vector{StrF(S)} is only possible for S <= 16 for now"))
    end
end


text:
to_discrete_position(pos, pathfinder) =
    floor.(Int, pos ./ pathfinder.dims .* size(pathfinder.walkmap)) .+ 1
to_continuous_position(pos, pathfinder) =
    pos ./ size(pathfinder.walkmap) .* pathfinder.dims .-
    pathfinder.dims ./ size(pathfinder.walkmap) ./ 2.
sqr_distance(from, to, pathfinder::AStar{D,true}) where {D} =
    sum(min.(abs.(from .- to), pathfinder.dims .- abs.(from .- to)) .^ 2)
sqr_distance(from, to, pathfinder::AStar{D,false}) where {D} =
    sum((from .- to) .^ 2)

"""
    find_continuous_path(pathfinder, from, to)
Functions like `find_path`, but uses the output of `find_path` and converts it to the coordinate
space used by the corresponding `ContinuousSpace`. Performs checks on the last two waypoints
in the discrete path to ensure continuous path is optimal.
"""
function find_continuous_path(
    pathfinder::AStar{D},
    from::NTuple{D,Float64},
    to::NTuple{D,Float64},
) where {D}
    discrete_from = to_discrete_position(from, pathfinder)
    discrete_to = to_discrete_position(to, pathfinder)
    discrete_path = find_path(pathfinder, discrete_from, discrete_to)
    # find_path returns nothing if no path exists
    isnothing(discrete_path) && return
    # if discrete_path is empty, `from` and `to` are in the same grid cell,
    # so `to` is the only waypoint
    isempty(discrete_path) && return Path{D,Float64}(to)

    cts_path = Path{D,Float64}()
    for pos in discrete_path
        push!(cts_path, to_continuous_position(pos, pathfinder))
    end

    # Handles an edge case to prevent backtracking for a fraction of a grid cell
    # Consider case where each grid cell is (1., 1.), and the path is to be calculated from
    # (0.5, 0.5) to (0.5, 1.2). Without this, the resultant path would look like
    # [(0.5, 0.5), (0.5, 1.5), (0.5, 1.2)], causing the agent to go to the last waypoint and
    # then backtrack to the target
    last_pos = last(cts_path)
    pop!(cts_path)
    # It's possible there's only one waypoint in the path, in which case the second last
    # position is the starting position
    second_last_pos = isempty(cts_path) ? from : last(cts_path)
    last_to_end = sqr_distance(last_pos, to, pathfinder)
    second_last_to_end = sqr_distance(second_last_pos, to, pathfinder)
    if last_to_end < second_last_to_end
        push!(cts_path, last_pos)
    end
    # If `to` is already at the center of a grid cell, there's no need
    # to push it to the path
    last_to_end ≈ 0. || push!(cts_path, to)
    return cts_path
end

function Agents.plan_route!(
    agent::A,
    dest::NTuple{D,Float64},
    pathfinder::AStar{D,P,M,Float64},
) where {A<:AbstractAgent,D,P,M}
    path = find_continuous_path(pathfinder, agent.pos, dest)
    isnothing(path) && return
    pathfinder.agent_paths[agent.id] = path
end

function Agents.plan_best_route!(
    agent::A,
    dests,
    pathfinder::AStar{D,P,M,Float64};
    condition::Symbol = :shortest,
) where {A<:AbstractAgent,D,P,M}
    @assert condition ∈ (:shortest, :longest)
    compare = condition == :shortest ? (a, b) -> a < b : (a, b) -> a > b
    best_path = Path{D,Float64}()
    best_target = nothing
    for target in dests
        path = find_continuous_path(pathfinder, agent.pos, target)
        isnothing(path) && continue
        if isempty(best_path) || compare(length(path), length(best_path))
            best_path = path
            best_target = target
        end
    end

    isnothing(best_target) && return
    pathfinder.agent_paths[agent.id] = best_path
    return best_target
end

"""
    move_along_route!(agent, model::ABM{<:ContinuousSpace{D}}, pathfinder::AStar{D}, speed, dt = 1.0)
Move `agent` for one step along the route toward its target set by
[`plan_route!`](@ref) at the given `speed` and timestep `dt`.

For pathfinding in models with [`ContinuousSpace`](@ref)

If the agent does not have a precalculated path or the path is empty, it remains stationary.
"""
function Agents.move_along_route!(
    agent::A,
    model::ABM{<:ContinuousSpace{D},A},
    pathfinder::AStar{D},
    speed::Float64,
    dt::Real = 1.0,
) where {D,A<:AbstractAgent}
    isempty(agent.id, pathfinder) && return
    from = agent.pos
    next_pos = agent.pos
    while true
        next_waypoint = first(pathfinder.agent_paths[agent.id])
        dir = get_direction(from, next_waypoint, model)
        dist_to_target = norm(dir)
        # edge case
        if dist_to_target ≈ 0.
            from = next_waypoint
            popfirst!(pathfinder.agent_paths[agent.id])
            if isempty(agent.id, pathfinder)
                next_pos = next_waypoint
                break
            end
            continue
        end
        dir = dir ./ dist_to_target
        next_pos = from .+ dir .* (speed * dt)
        # overshooting means we reached the waypoint
        dist_to_next = edistance(from, next_pos, model)
        if dist_to_next > dist_to_target
            # change from and dt so it appears we're coming from the waypoint just skipped, instead
            # of directly where the agent was. E.g:
            # agent.pos = (0, 0)
            # pathfinder.agent_paths[agent.id] = [(1, 0), (1, 1)]
            # speed = 1
            # dt = 1.2
            # without this, agent would end up at (0.85, 0.85) instead of (1, 0.2)
            from = next_waypoint
            dt -= dist_to_target / speed
            popfirst!(pathfinder.agent_paths[agent.id])
            # if the path is now empty, go directly to the end
            if isempty(agent.id, pathfinder)
                next_pos = next_waypoint
                break
            end
        else
            break
        end
    end

    move_agent!(agent, next_pos, model)
end

function random_walkable(model::ABM{<:ContinuousSpace{D}}, pathfinder::AStar{D}) where {D}
    discrete_pos = Tuple(rand(
        model.rng,
        filter(x -> pathfinder.walkmap[x], CartesianIndices(pathfinder.walkmap))
    ))
    half_cell_size = model.space.extent ./ size(pathfinder.walkmap) ./ 2.
    return to_continuous_position(discrete_pos, pathfinder) .+
        Tuple(rand(model.rng, D) .- 0.5) .* half_cell_size
end

walkable_cells_in_radius(pos, r, pathfinder::AStar{D,false}) where {D} =
    Iterators.filter(
        x -> all(1 .<= x .<= size(pathfinder.walkmap)) &&
            pathfinder.walkmap[x...] &&
            sum(((x .- pos) ./ r) .^ 2) <= 1,
        Iterators.product([(pos[i]-r[i]):(pos[i]+r[i]) for i in 1:D]...)
    )

walkable_cells_in_radius(pos, r, pathfinder::AStar{D,true}) where {D} =
    Iterators.map(
        x -> mod1.(x, size(pathfinder.walkmap)),
            Iterators.filter(
                x -> pathfinder.walkmap[mod1.(x, size(pathfinder.walkmap))...] && sum(((x .- pos) ./ r) .^ 2) <= 1,
                Iterators.product([(pos[i]-r[i]):(pos[i]+r[i]) for i in 1:D]...)
            )
    )

"""
    Pathfinding.random_walkable(pos, model::ABM{<:ContinuousSpace{D}}, pathfinder::AStar{D}, r = 1.0)
Return a random position within radius `r` of `pos` which is walkable, as specified by `pathfinder`.
Return `pos` if no such position exists.
"""
function random_walkable(
    pos,
    model::ABM{<:ContinuousSpace{D}},
    pathfinder::AStar{D},
    r = 1.0,
) where {D}
    discrete_r = to_discrete_position(r, pathfinder) .- 1
    discrete_pos = to_discrete_position(pos, pathfinder)
    options = collect(walkable_cells_in_radius(discrete_pos, discrete_r, pathfinder))
    isempty(options) && return pos
    discrete_rand = rand(
        model.rng,
        options
    )
    half_cell_size = model.space.extent ./ size(pathfinder.walkmap) ./ 2.
    cts_rand = to_continuous_position(discrete_rand, pathfinder) .+
        Tuple(rand(model.rng, D) .- 0.5) .* half_cell_size
    dist = edistance(pos, cts_rand, model)
    dist > r && (cts_rand = mod1.(
        pos .+ get_direction(pos, cts_rand, model) ./ dist .* r,
        model.space.extent
    ))
    return cts_rand
end


text:
#=
Given a array of numbers representing the stock prices of a company in chronological order, write a function that calculates the maximum profit you could have made from buying and selling that stock once. You must buy before you can sell it.

For example, given [9, 11, 8, 5, 7, 10], you should return 5, since you could buy the stock at 5 dollars and sell it at 10 dollars.
=#

using Test
include("Solutions/problem47_find_maximum_profit.jl")

@test maximumprofit([9, 11, 8, 5, 7, 10]) == 5

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module MLJModelInterface

# ------------------------------------------------------------------------
# Dependencies (ScientificTypes itself does not have dependencies)
using ScientificTypes
using Random

# ------------------------------------------------------------------------
# exports

# mode
export LightInterface, FullInterface

# MLJ model hierarchy
export MLJType, Model, Supervised, Unsupervised,
       Probabilistic, JointProbabilistic, Deterministic, Interval, Static,
       UnivariateFinite

# parameter_inspection:
export params

# model constructor + metadata
export @mlj_model, metadata_pkg, metadata_model

# model api
export fit, update, update_data, transform, inverse_transform,
    fitted_params, predict, predict_mode, predict_mean, predict_median,
    predict_joint, evaluate, clean!

# model traits
export input_scitype, output_scitype, target_scitype,
       is_pure_julia, package_name, package_license,
       load_path, package_uuid, package_url,
       is_wrapper, supports_weights, supports_online,
       docstring, name, is_supervised,
       prediction_type, implemented_methods, hyperparameters,
       hyperparameter_types, hyperparameter_ranges

# data operations
export matrix, int, classes, decoder, table,
       nrows, selectrows, selectcols, select

# equality
export is_same_except, isrepresented

# re-exports from ScientificTypes
export Scientific, Found, Unknown, Known, Finite, Infinite,
       OrderedFactor, Multiclass, Count, Continuous, Textual,
       Binary, ColorImage, GrayImage, Image, Table
export scitype, scitype_union, elscitype, nonmissing, trait, info

# ------------------------------------------------------------------------
# To be extended

import Base.==
import Base: in, isequal
#
# ------------------------------------------------------------------------
# Mode trick

abstract type Mode end
struct LightInterface <: Mode end
struct FullInterface  <: Mode end

const INTERFACE_MODE = Ref{Mode}(LightInterface())

set_interface_mode(m::Mode) = (INTERFACE_MODE[] = m)

get_interface_mode() = INTERFACE_MODE[]

struct InterfaceError <: Exception
    m::String
end

# ------------------------------------------------------------------------
# Model types

abstract type MLJType end

abstract type Model   <: MLJType end

abstract type   Supervised <: Model end
abstract type Unsupervised <: Model end

abstract type Probabilistic <: Supervised end
abstract type Deterministic <: Supervised end
abstract type      Interval <: Supervised end

abstract type Static <: Unsupervised end

abstract type JointProbabilistic <: Probabilistic end

# ------------------------------------------------------------------------
# includes

include("utils.jl")
include("parameter_inspection.jl")
include("data_utils.jl")
include("metadata_utils.jl")
include("model_traits.jl")
include("model_def.jl")
include("model_api.jl")
include("equality.jl")

end # module


text:
using WeatherStation

basedir = get(ENV,"WEATHERSTATION_DIR","/var/lib/WeatherStation.jl/")
WeatherStation.server(basedir = basedir)