text:
using Test
using ClimaCore:
    Domains, Meshes, Topologies, Geometry, Operators, Spaces, Fields
using ClimaCore.Operators: local_weights, LinearRemap, remap, remap!
using ClimaCore.Topologies: Topology2D
using ClimaCore.Spaces: AbstractSpace
using ClimaCore.DataLayouts: IJFH
using IntervalSets, LinearAlgebra, SparseArrays

FT = Float64

"""
Checks if a linear remapping operator is conservative.

A linear operator `R` is conservative the local weight of each element in the
source mesh is distributed conservatively amongst the target mesh elements.
"""
function conservative(R::LinearRemap)
    Jt = local_weights(R.target)
    Js = local_weights(R.source)
    return (R.map)' * Jt ≈ Js
end

"""
Checks if a linear remapping operator is consistent.

A linear operator `R` is consistent if all its rows sum to one.
"""
function consistent(R::LinearRemap)
    row_sums = sum(R.map, dims = 2)
    v = ones(size(row_sums))
    return row_sums ≈ v
end

"""
Checks if a linear remapping operator is monotone.

A linear operator `R` is monotone if it is consistent and each element
`R_{ij} ≥ 0`.
"""
function monotone(R::LinearRemap)
    consistent(R) && all(i -> i >= 0, R.map)
end

@testset "Linear Operator Properties" begin
    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(-3) .. Geometry.XPoint{FT}(5),
        Geometry.YPoint{FT}(-2) .. Geometry.YPoint{FT}(8),
        x1boundary = (:bottom, :top),
        x2boundary = (:left, :right),
    )
    mesh = Meshes.RectilinearMesh(domain, 1, 1)
    topology = Topologies.Topology2D(mesh)
    space = Spaces.SpectralElementSpace2D(topology, Spaces.Quadratures.GL{1}())

    op = [0.5 0.25 0.0 0.25; 0.0 1.0 0.0 0.0; 0.25 0.25 0.25 0.25]
    R = LinearRemap(space, space, op)
    @test consistent(R)
    @test monotone(R)

    op = [0.5 0.25 0.0 0.25; 1.0 1.0 0.0 -1.0; 0.25 0.25 0.25 0.25]
    R = LinearRemap(space, space, op)
    @test consistent(R)
    @test !monotone(R)

    op = [0.5 0.25 0.0 0.25; 0.0 2.0 0.0 0.0; 0.25 0.25 0.25 0.25]
    R = LinearRemap(space, space, op)
    @test !consistent(R)
    @test !monotone(R)
end

@testset "Finite Volume Remapping" begin
    quad = Spaces.Quadratures.GL{1}() # FV specification

    @testset "1D Domains" begin
        @testset "Aligned Intervals Different Resolutions" begin
            domain = Domains.IntervalDomain(
                Geometry.XPoint(0.0) .. Geometry.XPoint(1.0),
                boundary_tags = (:left, :right),
            )

            n1 = 4
            mesh1 = Meshes.IntervalMesh(domain; nelems = n1)
            source_topo = Topologies.IntervalTopology(mesh1)
            source = Spaces.SpectralElementSpace1D(source_topo, quad)

            n2 = 2
            mesh2 = Meshes.IntervalMesh(domain; nelems = n2)
            target_topo = Topologies.IntervalTopology(mesh2)
            target = Spaces.SpectralElementSpace1D(target_topo, quad)

            @test local_weights(source) ≈ 0.25 * ones(4, 1)
            @test local_weights(target) ≈ 0.5 * ones(2, 1)

            # create remapping operator from 4-elem to 2-elem interval
            R = LinearRemap(target, source)
            R_true = [
                0.5 0.5 0.0 0.0
                0.0 0.0 0.5 0.5
            ]

            @test R.map ≈ R_true
            @test nnz(R.map) == 4 # check R is sparse

            @test conservative(R)
            @test consistent(R)
            @test monotone(R)

            @testset "Scalar Remap Operator Application" begin
                n = length(source.local_geometry)
                source_field =
                    Fields.Field(IJFH{FT, 1}(ones(1, 1, n, 1)), source)

                # test consistent remap
                target_field = remap(R, source_field)
                @test vec(parent(target_field)) ≈
                      ones(length(target.local_geometry))

                # test simple remap
                vec(parent(source_field)) .= [1.0; 2.0; 3.0; 4.0]
                remap!(target_field, R, source_field)
                @test vec(parent(target_field)) ≈ [1.5; 3.5]
            end
        end

        @testset "Unaligned Intervals Same Resolution" begin
            domain1 = Domains.IntervalDomain(
                Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
                boundary_tags = (:left, :right),
            )

            n1 = 3
            mesh1 = Meshes.IntervalMesh(domain1; nelems = n1)
            source_topo = Topologies.IntervalTopology(mesh1)
            source = Spaces.SpectralElementSpace1D(source_topo, quad)

            domain2 = Domains.IntervalDomain(
                Geometry.XPoint(0.0) .. Geometry.XPoint(2.0),
                boundary_tags = (:left, :right),
            )
            n2 = 3
            mesh2 = Meshes.IntervalMesh(domain2; nelems = n2)
            target_topo = Topologies.IntervalTopology(mesh2)
            target = Spaces.SpectralElementSpace1D(target_topo, quad)

            R = LinearRemap(target, source)
            R_true = spzeros(FT, 3, 3)
            R_true[1, 2] = 0.5
            R_true[1, 3] = 0.5
            R_true[2, 3] = 0.5

            @test R.map ≈ R_true
            @test nnz(R.map) == 3

            @test !conservative(R)
            @test !consistent(R)
            @test !monotone(R)
        end

        @testset "Concentric Domains of Different Length" begin
            domain1 = Domains.IntervalDomain(
                Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
                boundary_tags = (:left, :right),
            )

            n1 = 4
            mesh1 = Meshes.IntervalMesh(domain1; nelems = n1)
            source_topo = Topologies.IntervalTopology(mesh1)
            source = Spaces.SpectralElementSpace1D(source_topo, quad)

            domain2 = Domains.IntervalDomain(
                Geometry.XPoint(-0.5) .. Geometry.XPoint(0.5),
                boundary_tags = (:left, :right),
            )
            n2 = 4
            mesh2 = Meshes.IntervalMesh(domain2; nelems = n2)
            target_topo = Topologies.IntervalTopology(mesh2)
            target = Spaces.SpectralElementSpace1D(target_topo, quad)

            R = LinearRemap(target, source)
            R_true = spzeros(4, 4)
            R_true[1, 2] = 1.0
            R_true[2, 2] = 1.0
            R_true[3, 3] = 1.0
            R_true[4, 3] = 1.0

            @test R.map ≈ R_true
            @test nnz(R.map) == 4

            @test !conservative(R)
            @test consistent(R)
            @test monotone(R)

            @testset "Scalar Remap Operator Application" begin
                n = length(source.local_geometry)
                source_field =
                    Fields.Field(IJFH{FT, 1}(ones(1, 1, n, 1)), source)

                # test consistent remap
                target_field = remap(R, source_field)
                @test vec(parent(target_field)) ≈
                      ones(length(target.local_geometry))

                # test simple remap
                vec(parent(source_field)) .= [1.0; 2.0; 3.0; 4.0]
                remap!(target_field, R, source_field)
                @test vec(parent(target_field)) ≈ [2.0; 2.0; 3.0; 3.0]
            end
        end
    end

    @testset "2D Domains" begin
        @testset "Aligned Grids Different Resolutions" begin
            domain = Domains.RectangleDomain(
                Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
                Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0),
                x1boundary = (:bottom, :top),
                x2boundary = (:left, :right),
            )

            m1, n1 = 2, 2
            mesh1 = Meshes.RectilinearMesh(domain, m1, n1)
            source_topo = Topologies.Topology2D(mesh1)
            source = Spaces.SpectralElementSpace2D(source_topo, quad)

            m2, n2 = 3, 3
            mesh2 = Meshes.RectilinearMesh(domain, m2, n2)
            target_topo = Topologies.Topology2D(mesh2)
            target = Spaces.SpectralElementSpace2D(target_topo, quad)

            @test local_weights(source) ≈ ones(4, 1)
            @test local_weights(target) ≈ (2 / 3)^2 * ones(9, 1)

            # create remapping operator from 2x2 to 3x3 grid
            R = LinearRemap(target, source)
            R_true = [
                1.0 0.0 0.0 0.0
                0.5 0.5 0.0 0.0
                0.0 1.0 0.0 0.0
                0.5 0.0 0.5 0.0
                0.25 0.25 0.25 0.25
                0.0 0.5 0.0 0.5
                0.0 0.0 1.0 0.0
                0.0 0.0 0.5 0.5
                0.0 0.0 0.0 1.0
            ]

            @test R.map ≈ R_true
            @test nnz(R.map) == 16 # check R is sparse

            @test conservative(R)
            @test consistent(R)
            @test monotone(R)

            @testset "Scalar Remap Operator Application" begin
                n = length(source.local_geometry)
                source_field =
                    Fields.Field(IJFH{FT, 1}(ones(1, 1, n, 1)), source)

                # test consistent remap
                target_field = remap(R, source_field)
                @test vec(parent(target_field)) ≈
                      ones(length(target.local_geometry))

                # test simple remap
                vec(parent(source_field)) .= [1.0; 2.0; 3.0; 4.0]
                remap!(target_field, R, source_field)
                @test vec(parent(target_field)) ≈
                      [1.0; 1.5; 2.0; 2.0; 2.5; 3.0; 3.0; 3.5; 4.0]
            end
        end

        @testset "Unaligned Grids Same Resolution" begin
            domain1 = Domains.RectangleDomain(
                Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
                Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0),
                x1boundary = (:bottom, :top),
                x2boundary = (:left, :right),
            )

            m1, n1 = 2, 2
            mesh1 = Meshes.RectilinearMesh(domain1, m1, n1)
            source_topo = Topologies.Topology2D(mesh1)
            source = Spaces.SpectralElementSpace2D(source_topo, quad)

            domain2 = Domains.RectangleDomain(
                Geometry.XPoint(0.0) .. Geometry.XPoint(2.0),
                Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0),
                x1boundary = (:bottom, :top),
                x2boundary = (:left, :right),
            )
            m2, n2 = 2, 2
            mesh2 = Meshes.RectilinearMesh(domain2, m2, n2)
            target_topo = Topologies.Topology2D(mesh2)
            target = Spaces.SpectralElementSpace2D(target_topo, quad)

            R = LinearRemap(target, source)
            R_true = spzeros(FT, 4, 4)
            R_true[1, 2] = 1.0
            R_true[3, 4] = 1.0

            @test R.map ≈ R_true
            @test nnz(R.map) == 2

            @test !conservative(R)
            @test !consistent(R)
            @test !monotone(R)
        end

        @testset "Concentric Domains of Different Area" begin
            domain1 = Domains.RectangleDomain(
                Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
                Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0),
                x1boundary = (:bottom, :top),
                x2boundary = (:left, :right),
            )

            m1, n1 = 2, 2
            mesh1 = Meshes.RectilinearMesh(domain1, m1, n1)
            source_topo = Topologies.Topology2D(mesh1)
            source = Spaces.SpectralElementSpace2D(source_topo, quad)

            domain2 = Domains.RectangleDomain(
                Geometry.XPoint(-0.5) .. Geometry.XPoint(0.5),
                Geometry.YPoint(-0.5) .. Geometry.YPoint(0.5),
                x1boundary = (:bottom, :top),
                x2boundary = (:left, :right),
            )
            m2, n2 = 2, 2
            mesh2 = Meshes.RectilinearMesh(domain2, m2, n2)
            target_topo = Topologies.Topology2D(mesh2)
            target = Spaces.SpectralElementSpace2D(target_topo, quad)

            R = LinearRemap(target, source)
            R_true = I

            @test R.map ≈ R_true
            @test nnz(R.map) == 4

            @test !conservative(R)
            @test consistent(R)
            @test monotone(R)

            @testset "Scalar Remap Operator Application" begin
                n = length(source.local_geometry)
                source_field =
                    Fields.Field(IJFH{FT, 1}(ones(1, 1, n, 1)), source)

                # test consistent remap
                target_field = remap(R, source_field)
                @test vec(parent(target_field)) ≈
                      ones(length(target.local_geometry))

                # test simple remap
                vec(parent(source_field)) .= [1.0; 2.0; 3.0; 4.0]
                remap!(target_field, R, source_field)
                @test vec(parent(target_field)) ≈ [1.0; 2.0; 3.0; 4.0]
            end
        end
    end
end

@testset "Finite Volume <-> Spectral Elements Remapping" begin

    @testset "1D Domains" begin
        domain = Domains.IntervalDomain(
            Geometry.XPoint(0.0) .. Geometry.XPoint(1.0),
            boundary_tags = (:left, :right),
        )

        @testset "Single aligned elements" begin
            quad1 = Spaces.Quadratures.GL{1}() # FV specification
            n1 = 1
            mesh1 = Meshes.IntervalMesh(domain; nelems = n1)
            topo1 = Topologies.IntervalTopology(mesh1)
            fv_space = Spaces.SpectralElementSpace1D(topo1, quad1)

            @test local_weights(fv_space) ≈ (1 / n1) * ones(n1, 1)

            for nq2 in [2, 5, 10]

                quad2 = Spaces.Quadratures.GLL{nq2}() # Spectral
                n2 = 1
                mesh2 = Meshes.IntervalMesh(domain; nelems = n2)
                topo2 = Topologies.IntervalTopology(mesh2)
                se_space = Spaces.SpectralElementSpace1D(topo2, quad2)

                _, w = Spaces.Quadratures.quadrature_points(FT, quad2)
                @test local_weights(se_space) ≈ w ./ 2

                # FV -> SE
                R = LinearRemap(se_space, fv_space)
                R_true = ones(nq2, 1)
                @test R.map ≈ R_true
                @test nnz(R.map) == nq2

                @test conservative(R)
                @test consistent(R)
                @test monotone(R)

                # SE -> FV
                R = LinearRemap(fv_space, se_space)
                R_true = local_weights(se_space)'
                @test R.map ≈ R_true
                @test nnz(R.map) == nq2

                @test conservative(R)
                @test consistent(R)
                @test monotone(R)
            end
        end

        @testset "Multiple elements" begin
            quad1 = Spaces.Quadratures.GL{1}() # FV specification
            n1 = 2
            mesh1 = Meshes.IntervalMesh(domain; nelems = n1)
            topo1 = Topologies.IntervalTopology(mesh1)
            fv_space = Spaces.SpectralElementSpace1D(topo1, quad1)

            nq2 = 2
            quad2 = Spaces.Quadratures.GLL{nq2}() # Spectral
            n2 = 3
            mesh2 = Meshes.IntervalMesh(domain; nelems = n2)
            topo2 = Topologies.IntervalTopology(mesh2)
            se_space = Spaces.SpectralElementSpace1D(topo2, quad2)

            # FV -> SE
            R = LinearRemap(se_space, fv_space)
            @test nnz(R.map) == 8
            @test conservative(R)
            @test consistent(R)
            @test monotone(R)

            # SE -> FV
            R = LinearRemap(fv_space, se_space)
            @test nnz(R.map) == 8
            @test conservative(R)
            @test consistent(R)
            @test monotone(R)

            nq2 = 3
            quad2 = Spaces.Quadratures.GLL{nq2}() # Spectral
            n2 = 3
            mesh2 = Meshes.IntervalMesh(domain; nelems = n2)
            topo2 = Topologies.IntervalTopology(mesh2)
            se_space = Spaces.SpectralElementSpace1D(topo2, quad2)

            # FV -> SE
            R = LinearRemap(se_space, fv_space)
            @test nnz(R.map) == 12
            @test conservative(R)
            @test consistent(R)
            @test !monotone(R)

            # SE -> FV
            R = LinearRemap(fv_space, se_space)
            @test nnz(R.map) == 12
            @test conservative(R)
            @test consistent(R)
            @test !monotone(R)
        end
    end
end


text:
@testset "Fizz & Buzz" begin
    # @test fizz(1) == 1
    @test 1 |> fizz == 1
    @test 3 |> fizz == "Fizz"
    @test 1 |> buzz == 1
    @test 5 |> buzz == "Buzz"
    # @test buzz(fizz(1)) == 1
    @test 1 |> fizz |> buzz == 1
    @test 3 |> fizz |> buzz == "Fizz"
    @test 5 |> fizz |> buzz == "Buzz"
    @test 15 |> fizz |> buzz == "FizzBuzz"
    @test 15 |> buzz |> fizz == "BuzzFizz"
end

ext: jl


text:
using Distances

dist = Euclidean()

tic()
tic()
# file lines iterator object
input = eachline(STDIN)

state = start(input)
# read data size N
(n, state) = next(input, state)
n = parse(Int32, n)

# 2xN node's coordinates matrix
# Data[:, k] = [x_k, y_k]
data = Array{Float64}(2, n)
# NxN node's distances matric
d = Array{Float64}(n, n)

# read file data to Data matrix
while !done(input, state)
    (v, state) = next(input, state)
    (i_v, x_v, y_v) = split(v)
    ii = parse(Int32, i_v)
    xx = parse(Float64, x_v)
    yy = parse(Float64, y_v)
    data[:,ii] = [xx, yy]
end
# calculate nodes distances as NxN matrix
pairwise!(d, dist, data)

# show(d)

####### SETTINGS #######
const itermax = Integer(floor((10_000_000*sqrt(n))/n^2))
println(string("iterations: ", itermax))
const penalty = 10 # rounds tabbed
const checkpercent = 1.0 # * sqrt(neighborscount)
const noimprovpercent = 0.1 # * itermax
########################
tabulist = fill(-penalty, n)

currentsol = Array{Int32}(n+1)
currentsol[1] = currentsol[n+1] = Int32(1)
bestsol = Array{Int32}(n+1)
bestsol[1] = bestsol[n+1] = 1

function tab(i::Integer, j::Integer, iteration::Integer)
    tabulist[currentsol[i]] = iteration
    tabulist[currentsol[j]] = iteration
end

function istabbed(i::Integer, j::Integer, iteration::Integer)
    return (iteration - tabulist[currentsol[i]]) < penalty || (iteration - tabulist[currentsol[j]]) < penalty
end

function distance(solution::Array{Int32})
    cost = 0
    for i in 1:length(solution)-1
        cost += d[solution[i],solution[i+1]]
    end
    return cost
end

function distancediff(tour::Array{Int32}, j::Integer, k::Integer)
    return sum([
            -d[tour[j-1], tour[j]],
            -d[tour[j], tour[j+1]],
             d[tour[j-1], tour[k]],
             d[tour[k], tour[j+1]],
            -d[tour[k-1], tour[k]],
            -d[tour[k], tour[k+1]],
             d[tour[k-1], tour[j]],
             d[tour[j], tour[k+1]]
        ])
end

#
# for i in 2:n
#     currentsol[i] = Int32(i)
# end

for i in 1:n
    d[i,i] = Inf
end

d1 = copy(d)
k = 1
for i in 2:length(currentsol)-1
    currentsol[i] = findfirst(d1[k,:].==min(d1[k,:]...))
    k = currentsol[i]
    d1[k,1] = Inf
    for visited_node in 1:i
        d1[k, currentsol[visited_node]] = Inf
    end
end

currentcost = distance(currentsol)

bestsol = copy(currentsol)
bestcost = currentcost

bestneighborsol = copy(currentsol)
bestneighborcost = currentcost

neighborsol = copy(bestneighborsol)
neighborcost = bestneighborcost

noimprov = 0
checked = 0

neighborscount = 0.5 * (length(currentsol)-3) * (length(currentsol)-2)

println(string("initial cost: ", bestcost))
# println(string("initial route: ", bestsol))
println(string("initialization time: ", toq(), "s"))

for iter in 1:itermax
    neighborsol = copy(bestneighborsol)
    neighborcost = bestneighborcost

    for j in shuffle(2:length(neighborsol)-2), k in shuffle(j+1:length(neighborsol)-1)
    # for j in 2:length(neighborsol)-2, k in j+1:length(neighborsol)-1
        if k-j>1 && !istabbed(j,k,iter)
            currentsol = copy(neighborsol)
            currentcost = neighborcost

            currentcost += distancediff(currentsol, j, k)

            currentsol[k], currentsol[j] = currentsol[j], currentsol[k]

            # TODO: tab every neighbor, not only ones checked
            tab(j, k, iter)

            if currentcost <= bestneighborcost
                bestneighborsol = copy(currentsol)
                bestneighborcost = currentcost
            end

            checked+=1
        end
        if checked >= sqrt(neighborscount) * checkpercent
            checked = 0
            break
        end
    end

    for j in 2:length(neighborsol)-2, k in j+1:length(neighborsol)-1
        if k-j>1
            tab(j, k, iter)
        end
    end

    if bestneighborcost < bestcost
        bestsol = copy(bestneighborsol)
        bestcost = bestneighborcost
    end
end

# println(STDOUT, bestcost)
# for n in bestsol
#     print(STDERR, n)
#     print(STDERR, " ")
# end
# println(STDERR)

println(STDOUT, string("best cost: ", bestcost))
println(string("total time: ", toq(), "s"))

max_stars_repo_name: AfroMetal/algorithms-tsp-tabusearch

max_issues_repo_name: AfroMetal/algorithms-tsp-tabusearch

max_forks_repo_name: AfroMetal/algorithms-tsp-tabusearch


text:
module SaveFigs

export savefigs

using Plots

function savefigs(fig, fname::AbstractString, formats...; sizes=nothing)
    isdir(fname) || mkdir(fname)
    if isnothing(sizes)
        dir = fname
        _saveformats(fig, dir, fname, formats)
    else
        dir = pwd()
        _saveformats_sizes(fig, dir, fname, formats, sizes)
    end
end

function savefigs(fname::AbstractString, formats...; sizes=nothing)
    savemultifigs(Plots.current(), fname, formats...; sizes=sizes)
end

function _saveformats(fig, dir, fname, formats)
    for f = formats
        savefig(fig, joinpath(dir, "$(fname).$(f)"))
    end
end

function _saveformats_sizes(fig, dir, fname, formats, sizes)
    # Save original size
    _saveformats(fig, fname, fname, formats)
    # Resize and save
    for (w,h) = sizes
        plot!(fig, size=(w,h))
        _saveformats(fig, fname, "$(fname)_$(w)x$(h)", formats)
    end
end

end


text:
# %% imports

using Statistics
using PyCall
gym = pyimport("gym")

include("../../Modules/Networks/net_arch.jl")
include("../../Modules/Poisson/PoissonStateSpike.jl")

# %% functions

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl >= sr) ? action = 0 : action = 1
    if sl==sr
        action = rand(0:1)
    end
    return action, sl, sr
end

function norm_obs(obs::Array{Float64})
    pos = obs[1]
    theta = mod(obs[3]+pi, 2pi) - pi
    vel = [obs[2], obs[4]]
    pos_scaled = (pos + 2.4)/4.8
    cone = 12*pi/180
    theta_scaled = (theta + pi)/(2pi)
    vel_scaled = (vel .+ 7.5)/15
    obs_scaled = [pos_scaled, vel_scaled[1], theta_scaled, vel_scaled[2]]
    return obs_scaled
end

function norm_cartpole(obs::Array{Float64})
    pos = obs[1]
    theta = mod(obs[3]+pi,2pi)-pi
    vel = [obs[2], obs[4]]
    pos_scaled = (pos + 2.4)/4.8
    cone = 12*pi/180
    theta_scaled = (theta + cone)/(2*cone)
    vel_scaled = (vel .+ 7.5)/15
    obs_scaled = [pos_scaled, vel_scaled[1], theta_scaled, vel_scaled[2]]
    return obs_scaled
end

# %% sim steps

function play_decision(net::Network, state_alpha::Array{Float64}, train::Bool)
    win_size = 40 #ms
    spikes = int_spikes(state_alpha)
    net.v .= net.param["c"]
    count_tot = zeros(net.arch[end])

    @inbounds for msec in 1:win_size
        spiked = step_LIF!(net, spikes, train)
        if msec >= length(net.arch)-1
            count_tot .+= count_spikes(net.arch, spiked, length(net.arch))
            #push!(spike_hist,spiked)
        end
    end

    n_left = div(net.arch[end],2)
    n_right = net.arch[end] - n_left
    action, sl, sr = left_right(count_tot, n_left, n_right)

    return action, sl, sr
end

function reward_net(state::Array{Float64}, action::Int64)
    env = gym.make("CartPole-v1")
    env.reset()
    env.env.state = state
    new_state, r, done, _ = env.step(action)
    env.close()
    env = nothing
    return r,new_state,done
end

function dopa_inc(r_stud::Real, r_teach::Real)
    if r_stud == r_teach
        da = 1.0
    elseif r_stud != r_teach
        da = - 200.0
    end
    return da
end

function play_episode(n_stud::Network, n_teach::Network, matalpha::Array{Float64,2}, n_steps::Int64 = 200)
    # initial point
    env = gym.make("CartPole-v1")
    state = env.reset()
    env.close()
    env = nothing
    Base.GC.gc()

    state_alpha = matalpha * norm_cartpole(state)

    #record da evo
    da_evo = []
    r_tot = 0
    sls = []
    srs = []
    # sim loop
    for j in 1:n_steps
        state_alpha = matalpha * norm_cartpole(state)

        action_stud, sl_stud, sr_stud = play_decision(n_stud, state_alpha, true)
        @show sl_stud-sr_stud
        push!(sls, sl_stud)
        push!(srs, sr_stud)
        action_teach, sl_teach, sr_teach = play_decision(n_teach, state_alpha, false)
        #@show sl_teach,sr_teach
        r_stud, new_state, done= reward_net(state, action_stud)
        r_teach, _,_ = reward_net(state,action_teach)
        # r_teach, new_state, done= reward_net(state, action_teach)
        # r_stud, _,_ = reward_net(state,action_stud)

        # Compare the rewards
        da_inc = dopa_inc(action_stud, rand(0:1))
        n_stud.da += da_inc * n_stud.param["da_inc"]
        push!(da_evo,n_stud.da)
        r_tot += r_stud

        if done
            break
        end
        state = new_state
    end
    ss = [mean(sls), mean(srs)]
    return r_tot,da_evo, ss
end

# %% tests

include("load_net_2.jl")

filepath = "Julia/Actors/cartpoleElites/sNESelites_cartpole_1.jld"

n_teach, n_stud, matalpha = init_teacher_student(filepath, 1)
weight_viz(n_stud)
rewards = []
sls = []
srs = []

for epoch in 1:1000
    @show epoch
    @time r_ep, da_evo, ss= play_episode(n_stud, n_teach, matalpha,200)
    @show r_ep
    push!(sls, ss[1])
    push!(srs, ss[2])
    #weight_viz(n_stud, "heatmap")
    #display(plot(da_evo))
    push!(rewards, r_ep)
end

# %% viz rewards

plot(rewards, label = "Student network", legend = true, xlabel = "Episode", ylabel = "Reward" )
plot!([200], seriestype = "hline", color = :green, label = "Teacher Network")

# %%
plot(sls, label= "Mean action score left", title = "Student action scores", xlabel = "Episode", ylabel = "Score")
plot!(srs, label = "Mean action score right")

# %%
plot(sls-srs)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
module BytePairEncoding

using WordTokenizers
using InternedStrings

export normalize, UtfNormalizer
export BPELearner, learn!, add!, emit
export Bpe, process_line, segment, segment_token

export GenericBPE, ByteLevelBPE

include("./utfnorm.jl")
include("./codemap.jl")
include("./glossary.jl")
include("./mstring.jl")
include("./bpe.jl")
include("./learn.jl")
include("./stats.jl")
include("./defaults.jl")
include("./api.jl")
include("./old_api.jl")

end # module


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Refine prescribed policy for a set of sample paths for all positions
################################################################################

immutable PathPolicy

    ## Members
    policy::Array{Int64, 3}
    pathnum::Int64 # number of paths
    tnum::Int64 # number of time points
    pnum::Int64 # number of positions
    anum::Int64 # number of actions
    snum::Int64 # state size

    function PathPolicy(css::CSS, bellman::Bellman, path::Path)

        ## Get sizes
        pathnum = path.pathnum
        pnum = bellman.pnum
        tnum = bellman.tnum
        anum = bellman.anum
        snum = bellman.snum

        ## Construct tree structure for neighbours
        tree = NearestNeighbors.KDTree(transpose(css.grid))

        ## Container
        policy = Array{Int64}(pathnum, pnum, tnum - 1)
        path_state = Array{Float64}(pathnum, snum)
        container = zero(Array{Float64}(anum, snum))
        neighbours = Array{Int64}(pathnum)

        ## Obtain the prescribed policy via recursion
        t = tnum - 1 # recall that last time is scrap
        if bellman.fullcontrol # positions fully controllable
            while (t > 0)
                path_state = path.sample[t, :, :]
                neighbours = knn(tree, transpose(path_state), 1, true)[1]
                for i = 1:pathnum
                    for p = 1:pnum
                        container = zero(container)
                        for a = 1:anum
                            container[a, :] =
                                bellman.evalue[neighbours[i], :, css.control[p, a], t]
                            container[a, :] += css.Reward(t, path_state[i, :], p, a)
                        end
                        policy[i, p, t] = indmax(container * path_state[i, :])
                    end
                end
                t = t - 1
            end
        else # positions not fully controllable
             while (t > 0)
                path_state = path.sample[t, :, :]
                neighbours = knn(tree, transpose(path_state), 1, true)[1]
                for i = 1:pathnum
                    for p = 1:pnum
                        container = zero(container)
                        for a = 1:anum
                            for pp = 1:pnum
                                container[a, :] += css.control[p, a, pp] *
                                    bellman.evalue[neighbours[i], :, pp, t][:]
                            end
                            container[a, :] += css.Reward(t, path_state[i, :], p, a)
                        end
                        policy[i, p, t] = indmax(container * path_state[i, :])
                    end
                end
                t = t - 1
            end
        end

        return new(policy, pathnum, tnum, pnum, anum, snum)

    end

end


text:
using MCMCChains

@testset "Test Chains Package." begin
    chn = read("../chain-file-10_1.jls", Chains)
end


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

mutable struct TopologyAssociation <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    associationType::Any # spec type: Union{ Nothing, String } # spec name: associationType

    function TopologyAssociation(;name=nothing, resourceId=nothing, associationType=nothing)
        o = new()
        validate_property(TopologyAssociation, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(TopologyAssociation, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(TopologyAssociation, Symbol("associationType"), associationType)
        setfield!(o, Symbol("associationType"), associationType)
        o
    end
end # type TopologyAssociation

const _property_map_TopologyAssociation = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("resourceId")=>Symbol("resourceId"), Symbol("associationType")=>Symbol("associationType"))
const _property_types_TopologyAssociation = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("resourceId")=>"String", Symbol("associationType")=>"String")
Base.propertynames(::Type{ TopologyAssociation }) = collect(keys(_property_map_TopologyAssociation))
Swagger.property_type(::Type{ TopologyAssociation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TopologyAssociation[name]))}
Swagger.field_name(::Type{ TopologyAssociation }, property_name::Symbol) =  _property_map_TopologyAssociation[property_name]

const _allowed_TopologyAssociation_associationType = ["Associated", "Contains"]

function check_required(o::TopologyAssociation)
    true
end

function validate_property(::Type{ TopologyAssociation }, name::Symbol, val)
    if name === Symbol("associationType")
        Swagger.validate_param(name, "TopologyAssociation", :enum, val, _allowed_TopologyAssociation_associationType)
    end
end


text:
###########################################################################
# Least-squares
###########################################################################
function bearing_ls(angles, gains, ref_angles, ref_gains)
	@assert length(angles) == length(gains)
	num_gains = length(gains)
	best_err = Inf
	best_shift = 0
	for shift = 0:359
		err = 0.0
		for j = 1:num_gains
			if validgain(gains[j])
				ref_gain = interp_linear(angles[j], shift, ref_angles, ref_gains)
				err += (gains[j] - ref_gain)^2
			end
		end
		if err < best_err
			best_err = err
			best_shift = shift
		end
	end
	return 360 - best_shift
end

# Uses least-square method
# samples is how many samples you want to use (maybe you dont want whole rotation)
# sense is direction of rotation: -1 for negative sense, 1 for positive sense
function bearing_ls(dir_file::AbstractString, omni_file::AbstractString, ref_file::AbstractString; samples::Int64=0, sense::Int64=1)
	test_gp = GainPattern(dir_file) - GainPattern(omni_file) + 5
	ref_gp = GainPattern(ref_file)

	num_samples = length(test_gp.angles)
	num_exp = length(test_gp.samples[1])

	# Allowing user to specify count
	if (samples > 0) && (samples < num_samples)
		num_samples = samples
	end

	aoa_array = zeros(num_exp)
	angles = zeros(num_samples)
	gains = zeros(num_samples)
	for i = 1:num_exp
		if sense > 0
			for j = 1:num_samples
				# Loop through all samples in this experiment...
				angles[j] = test_gp.angles[j]
				gains[j] = test_gp.samples[j][i]
			end
		elseif sense == -1
			# Backwards loop, starting from last measurement
			for j = 1:num_samples
				angles[j] = test_gp.angles[end - j + 1]
				gains[j] = test_gp.samples[end - j + 1][i]
			end
		elseif sense == -2
			# Backwards loop, starting from origin
			for j = 2:num_samples
				angles[j] = test_gp.angles[end - j + 2]
				gains[j] = test_gp.samples[end - j + 2][i]
			end
			angles[1] = test_gp.angles[1]
			gains[1] = test_gp.samples[1][i]
		end
		aoa_array[i] = bearing_ls(angles, gains, ref_gp.angles, ref_gp.meangains)
	end
	return aoa_array
end


text:
using Statistics
using LinearAlgebra
using Random

using Distributions

include("da.jl")
import .DA

include("ens_forecast.jl")
import .ens_forecast

include("models.jl")
import .Models

include("integrators.jl")
import .Integrators

using BandedMatrices

Random.seed!(1)

D1 = 220
D2 = 20
models = [Models.lorenz96_twoscale_err, Models.lorenz96_half_true]
model_errs_prescribed = [nothing, nothing, nothing, nothing]
model_true = Models.lorenz96_twoscale_true
n_models = length(models)
dd = zeros(D1)
dd[1:11:D1] .= 1
I1to2 = diagm(dd)[Vector{Bool}(dd), :]
I2to1 = I1to2'
obs_ops = [I(D1), I2to1]
mappings = Matrix{AbstractArray}(undef, n_models, n_models)
mappings[1, 1] = I(D1)
mappings[1, 2] = I1to2
mappings[2, 1] = I2to1
mappings[2, 2] = I(D2)

H = I(D1)
ens_sizes = [20, 20]
model_sizes = [D1, D2]
integrators = [Integrators.rk4, Integrators.rk4, Integrators.rk4, Integrators.rk4]
integrator_true = Integrators.rk4
da_method = DA.ensrf
localization = diagm(ones(D1))

indices = reshape(1:220, 11, :)
first_layer_indices = indices[1, :]
second_layer_indices = indices[2:end, :]

c = 4
for (ii, i) in enumerate(first_layer_indices)
    for (ij, j) in enumerate(first_layer_indices)
        r = min(mod(ii - ij, 0:20), mod(ij - ii, 0:20))/c
        localization[i, j] = DA.gaspari_cohn(r)
        localization[j, i] = DA.gaspari_cohn(r)
    end
end

for i=1:20
    layer_indices = indices[:, i]
    for j=layer_indices
        for k=layer_indices
            localization[j, k] = localization[k, j] = 1
        end
    end
end

x0 = rand(D1)
t0 = 0.0
Δt = 0.005
outfreq = 1
transient = 2000
x = integrators[1](models[1], x0, t0, transient*outfreq*Δt, Δt, inplace=false)
R = Symmetric(diagm(var(x, dims=1)[:]*0.1))
ens_errs = [R, R[first_layer_indices, first_layer_indices]]
gen_ensembles = false
assimilate_obs = false
all_orders = false
save_Q_hist = false

leads = 5
ref_model = 1
x0 = x[end, :]

n_cycles = 500*leads
ρ = 1e-3

window = 10

ensembles = [mappings[ref_model, 1]*x0 .+ rand(MvNormal(ens_errs[1]), sum(ens_sizes))]
info_a = ens_forecast.da_cycles(x0=x0, ensembles=ensembles, models=[models[1]],
                                model_true=model_true, obs_ops=[obs_ops[1]], H_true=I,
                                model_errs=[ens_errs[1]],
                                model_errs_prescribed=[model_errs_prescribed[1]],
                                integrators=[integrators[1]],
                                integrator_true=integrator_true, da_method=da_method,
                                localization=mappings[ref_model, 1]*localization*mappings[ref_model, 1]',
                                ens_sizes=[sum(ens_sizes)], Δt=Δt, window=window,
                                n_cycles=n_cycles, outfreq=outfreq,
                                model_sizes=[model_sizes[1]], R=R,
                                ens_errs=[ens_errs[1]], ρ=ρ,
                                gen_ensembles=gen_ensembles,
                                assimilate_obs=true, save_analyses=true, save_trues=true,
                                leads=1, save_Q_hist=save_Q_hist,
                                mappings=mappings[1:1, 1:1])

#x0 = integrators[1](model_true, x0, t0, window*outfreq*Δt, Δt)

infos = Vector(undef, n_models)
for model=1:n_models
    model_errs = [ens_errs[model]]

    ens_size = ens_sizes[model]
    ensembles = [mappings[ref_model, model]*x0 .+ rand(MvNormal(ens_errs[model]), ens_size)]

    if model == 2
        analyses = permutedims(cat([mappings[1, 2]*info_a.analyses[c, :, :] for c=1:n_cycles]..., dims=3), [3, 1, 2])
    else
        analyses = info_a.analyses
    end
    info = ens_forecast.da_cycles(x0=x0, ensembles=ensembles, models=[models[model]],
                                  model_true=model_true, obs_ops=[obs_ops[model]], H_true=I,
                                  model_errs=model_errs,
                                  model_errs_prescribed=[model_errs_prescribed[model]],
                                  integrators=[integrators[model]],
                                  integrator_true=integrator_true, da_method=da_method,
                                  localization=mappings[ref_model, model]*localization*mappings[ref_model, model]',
                                  ens_sizes=[ens_size], Δt=Δt, window=window,
                                  n_cycles=n_cycles, outfreq=outfreq,
                                  model_sizes=[model_sizes[model]], R=R,
                                  ens_errs=[ens_errs[model]], ρ=ρ,
                                  gen_ensembles=gen_ensembles,
                                  assimilate_obs=assimilate_obs, save_analyses=false,
                                  leads=leads, save_Q_hist=save_Q_hist,
                                  mappings=mappings[model:model, model:model],
                                  prev_analyses=analyses)
    infos[model] = info
end

ensembles = [mappings[ref_model, model]*x0 .+ rand(MvNormal(ens_errs[model]), ens_sizes[model]) for model=1:n_models]

model_errs = ens_errs

info_mm = ens_forecast.da_cycles(x0=x0, ensembles=ensembles, models=models,
                                 model_true=model_true, obs_ops=obs_ops,
                                 model_errs=model_errs,
                                 model_errs_prescribed=model_errs_prescribed,
                                 integrators=integrators, integrator_true=integrator_true,
                                 da_method=da_method, localization=localization,
                                 ens_sizes=ens_sizes, Δt=Δt, window=window,
                                 n_cycles=n_cycles, outfreq=outfreq,
                                 model_sizes=model_sizes, R=R, ens_errs=ens_errs, ρ=ρ,
                                 all_orders=all_orders, combine_forecasts=true,
                                 gen_ensembles=gen_ensembles, assimilate_obs=assimilate_obs,
                                 leads=leads, save_Q_hist=save_Q_hist, ref_model=ref_model,
                                 mappings=mappings, prev_analyses=info_a.analyses)

max_stars_repo_name: eviatarbach/mmda

max_issues_repo_name: eviatarbach/mmda

max_forks_repo_name: eviatarbach/mmda


text:
"""
    MostPopular()

Non-personalized baseline model which predicts n-most popular items in the corpus.
"""
mutable struct MostPopular <: AbstractRecommender
    df_popular::Union{Nothing,DataFrame}
    col_item::Union{Nothing,Symbol}
    user_histories::Any
    MostPopular() = new(nothing, nothing, nothing)
end

"""
    fit!(model::MostPopular, table; col_user = :userid, col_item = :itemid)

Fit most popular model.
"""
function fit!(model::MostPopular, table; col_user = :userid, col_item = :itemid, kwargs...)
    df = DataFrame(table)
    model.col_item = col_item
    model.df_popular = sort(combine(groupby(df, col_item), nrow), [:nrow], rev = true)

    users, items = make_u2i_dataset(table, col_user = col_user, col_item = col_item)
    model.user_histories = Dict()
    for (user, item) in zip(users, items)
        model.user_histories[user] = item
    end
    return model
end

"""
    predict_u2i(model::MostPopular, userid::Union{AbstractString,Int}, n::Int64; drop_history::Bool = false)

Make `n` prediction to user by most popular model.
"""
function predict_u2i(
    model::MostPopular,
    userid::Union{AbstractString,Int},
    n::Int64;
    drop_history::Bool = false,
    kwargs...,
)
    pred = model.df_popular[!, model.col_item]
    if drop_history && (userid in keys(model.user_histories))
        pred = filter(p -> !(p in model.user_histories[userid]), pred)
    end
    n = min(n, length(pred))
    return pred[1:n]
end


text:
# sample from a double well in 1D using HMC
using Printf
using Random
using LinearAlgebra
using ForwardDiff
using QuadGK
using BasicMD
using Statistics

include("potentials.jl")

β = 5.0;
x₀ = [-1.0];
M = 1.0;
seed = 100;
Δt = 2e-1;
n_iters = 10^4; # number of samples
nΔt = 10^1; # number of Verlet steps per HMC iteration
n_save_iters=10;

V = x->DoubleWell(x);
cfg = ForwardDiff.GradientConfig(V, x₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

# second moment and energy obsrevables
f₁ = x-> x[1]^2;
f₂ = x-> V(x[1]);
observables = (f₁, f₂);

sampler = HMC(V, gradV!, β, M, Δt, nΔt);

Random.seed!(100);
Xvals,_ = sample_trajectory(x₀, sampler, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
f₁_estimate = mean(f₁.(Xvals));
f₂_estimate = mean(f₂.(Xvals));
@printf("f₁ estimate with %d samples: %g\n",length(Xvals), f₁_estimate);
@printf("f₂ estimate with %d samples: %g\n", length(Xvals), f₂_estimate);
Random.seed!(100);
observable_samples = sample_observables(x₀, sampler, observables, options=MDOptions(n_iters=n_iters,n_save_iters=n_save_iters));
@printf("f₁ estimate with %d samples: %g\n",length(observable_samples), mean(observable_samples[1,:]));
@printf("f₂ estimate with %d samples: %g\n", length(observable_samples), mean(observable_samples[2,:]));

max_stars_repo_name: gideonsimpson/JuBasicMD


text:
depwarn2(old::AbstractString, new::AbstractString, funcsym::Symbol) =
  Base.depwarn("$old is deprecated, use $new instead.", funcsym)


max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# every model interface must implement a `fit` method of the form
# `fit(model, verb::Integer, training_args...) -> fitresult, cache, report`
function fit end

# fallback for static transformations
fit(::Static, ::Integer, a...) = (nothing, nothing, nothing)

# fallbacks for supervised models that don't support sample weights:
fit(m::Supervised, verb::Integer, X, y, w) = fit(m, verb, X, y)

# this operation can be optionally overloaded to provide access to
# fitted parameters (eg, coeficients of linear model):
fitted_params(::Model, fitres) = (fitresult=fitres,)

# each model interface may overload the following refitting method:
update(m::Model, verb::Integer, fitres, cache, a...) = fit(m, verb, a...)

# stub for online learning method update method
function update_data end

# supervised methods must implement the predict operation; additionally,
# probabilistic supervised models may overload one or more of
# `predict_mode`, `predict_median` and `predict_mean`
function predict end
function predict_mean end
function predict_mode end
function predict_median end

# unsupervised methods must implement the transform operation and
# may implement the inverse_transform operation
function transform end
function inverse_transform end

# models can optionally overload these for enable serialization in a
# custom format:
function save end
function restore end

# operations implemented by some meta-models:
function evaluate end


text:
using SpeedMapping, Test, LinearAlgebra, ForwardDiff

function f(x) # Rosenbrock objective
    f_out = 0.0
    for i ∈ 1:Int(length(x)/2)
        f_out += 100(x[2i - 1]^2 - x[2i])^2 + (x[2i - 1] - 1)^2
    end
    return f_out
end

function g!(∇,x) # Rosenbrock gradient
    for i ∈ 1:Int(length(x)/2)
        ∇[2i] = -200 * (x[2i - 1]^2 - x[2i])
        ∇[2i - 1] =  400 * (x[2i - 1]^2 - x[2i]) * x[2i - 1] + 2(x[2i - 1] - 1)
    end
    return nothing
end

function g_matrix!(∇,x) # Rosenbrock gradient for N × 2 input
    for i ∈ 1:size(x, 1)
        ∇[i, 2] = -200 * (x[i, 1]^2 - x[i, 2])
        ∇[i, 1] =  400 * (x[i, 1]^2 - x[i, 2]) * x[i, 1] + 2(x[i, 1] - 1)
    end
    return nothing
end

# Power method
C = [1 2 3; 4 5 6; 7 8 9]
A = C + C'

function m!(x_out, x_in) # map for the power method
    mul!(x_out, A, x_in)
    x_out ./= norm(x_out, Inf)
end

function m_horizontal!(x_out, x_in) # map for the power method
    x_out .= (A * x_in')'
    x_out ./= norm(x_out, Inf)
end

# Testing
function exception(expr, error)
    goodexception = false
    try
        eval(expr)
    catch e
        goodexception = isa(e, error)
    end

    return goodexception
end

@testset "SpeedMapping.jl" begin
    # Testing the complete algorithm with f, g!, ForwardDiff, m!, lower, upper
    @test speedmapping([-2.0,5.0]; g!).minimizer ≈ [1,1]
    @test speedmapping([-2.0,5.0]; g!, store_info = true).info.x[1][1] == -2
    @test speedmapping(zeros(2); f, g!).minimizer ≈ [1,1]
    @test speedmapping(zeros(2); f, g!, check_obj = true).minimizer ≈ [1,1]
    @test speedmapping(zeros(4); f).minimizer ≈ [1,1,1,1]
    @test speedmapping([5.0,5.0]; f, g!, lower = [1.5,-Inf]).minimizer ≈ [1.5, 2.25]
    @test speedmapping([0.0,0.0]; f, g!, upper = [Inf,0.25]).minimizer ≈ [0.5048795423935251, 0.25]
    @test speedmapping(ones(3); m!).minimizer' * A[:,3] ≈ 32.916472867168096
    @test speedmapping(ones(3); m!, stabilize = true).minimizer' * A[:,3] ≈ 32.91647286145264
    @test speedmapping(ones(3); m!, Lp = Inf).minimizer' * A[:,3] ≈ 32.916472867168096

    # Exceptions
    # Starting point outside boundary
    @test exception(:(speedmapping([0.0, 0.0]; f, g!, lower = [1,1])), DomainError)
    # Can't provide both g! and m!
    @test exception(:(speedmapping([0.0, 0.0]; g!, m!)), ArgumentError)
    # eltype(x_in) is Int
    #using parametric types,it fails before check_arguments
    @test_broken exception(:(speedmapping([0, 0]; g!)), ArgumentError)
    # check_obj without providing f
    @test exception(:(speedmapping([0.0, 0.0]; g!, check_obj = true)), ArgumentError)
    # The gradient is zero
    @test exception(:(speedmapping([1.0, 1.0]; g!)), DomainError)
end

@testset "other number types" begin
    #support for bigfloats
    @test speedmapping(zeros(BigFloat,2); g!).minimizer isa Vector{BigFloat}
    #support for ForwardDiff.Dual
    @test ForwardDiff.jacobian(x -> speedmapping(x; g!).minimizer
    ,[0.0,0.0]) isa Matrix{Float64}
    #support for float32 and matrices
    @test speedmapping(Float32.([-2.0 5.0]); g! = g_matrix!, tol = 1e-4).minimizer ≈ [1 1]
    @test (speedmapping(Float32.(ones(3)'); m! = m_horizontal!).minimizer * A[:, 3])[1] ≈ 32.916473f0
end


text:
module Upwind

function generate_solver(f₀, f, c)
   function upwind⁺!(g, f, i, i⁻)
        g[i] = f[i] - c*(f[i]-f[i⁻])
    end

    function solve⁺!()
        upwind⁺!(f, f₀, 1, length(f))
        @inbounds @fastmath @simd for i = 2:length(f)
            upwind⁺!(f, f₀, i, i-1)
        end
    end

    function upwind⁻!(g, f, i, i⁺)
        g[i] = f[i] - c*(f[i⁺]-f[i])
    end

    function solve⁻!()
        @inbounds @fastmath @simd for i = 1:length(f)-1
            upwind⁻!(f, f₀, i, i+1)
        end
        upwind⁻!(f, f₀, length(f), 1)
    end

    if c > 0
        return solve⁺!
    else
        return solve⁻!
    end
end

function generate_solver(f₀, f)
   function upwind⁺!(g, f, i, i⁻, c)
        g[i] = f[i] - c*(f[i]-f[i⁻])
    end

    function upwind⁻!(g, f, i, i⁺, c)
        g[i] = f[i] - c*(f[i⁺]-f[i])
    end

    function solve⁺!(c)
        upwind⁺!(f, f₀, 1, length(f), c)
        @inbounds @fastmath @simd for i = 2:length(f)
            upwind⁺!(f, f₀, i, i-1, c)
        end
    end

    function solve⁻!(c)
        @inbounds @fastmath @simd for i = 1:length(f)-1
            upwind⁻!(f, f₀, i, i+1, c)
        end
        upwind⁻!(f, f₀, length(f), 1, c)
    end

    function solve!(c)
        if c > 0
            solve⁺!(c)
        else
            solve⁻!(c)
        end
    end

    return solve!
end

end # module


text:
const OBJTYPE_TO_SIZENAME = Dict(
    MJCore.mjOBJ_BODY => :nbody,
    MJCore.mjOBJ_JOINT => :njnt,
    MJCore.mjOBJ_GEOM => :ngeom,
    MJCore.mjOBJ_SITE => :nsite,
    MJCore.mjOBJ_CAMERA => :ncam,
    MJCore.mjOBJ_LIGHT => :nlight,
    MJCore.mjOBJ_MESH => :nmesh,
    MJCore.mjOBJ_SKIN => :nskin,
    MJCore.mjOBJ_HFIELD => :nhfield,
    MJCore.mjOBJ_TEXTURE => :ntex,
    MJCore.mjOBJ_MATERIAL => :nmat,
    MJCore.mjOBJ_PAIR => :npair,
    MJCore.mjOBJ_EXCLUDE => :nexclude,
    MJCore.mjOBJ_EQUALITY => :neq,
    MJCore.mjOBJ_TENDON => :ntendon,
    MJCore.mjOBJ_ACTUATOR => :nu,
    MJCore.mjOBJ_SENSOR => :nsensor,
    MJCore.mjOBJ_NUMERIC => :nnumeric,
    MJCore.mjOBJ_TEXT => :ntext,
    MJCore.mjOBJ_TUPLE => :ntuple,
    MJCore.mjOBJ_KEY => :nkey,
)

function expandnq(m::jlModel, root, id)
    jtype = m.jnt_type[id]
    if jtype == MJCore.mjJNT_FREE
        map(suffix -> Symbol(root, :_, suffix), (:x, :y, :z, :qw, :qx, :qy, :qz))
    elseif jtype == MJCore.mjJNT_BALL
        map(suffix -> Symbol(root, :_, suffix), (:qw, :qx, :qy, :qz))
    else
        (root,)
    end
end

function expandnv(m::jlModel, root, id)
    jtype = m.jnt_type[id]
    if jtype == MJCore.mjJNT_FREE
        map(suffix -> Symbol(root, :_, suffix), (:x, :y, :z, :rx, :ry, :rz))
    elseif jtype == MJCore.mjJNT_BALL
        map(suffix -> Symbol(root, :_, suffix), (:rx, :ry, :rz))
    else
        (root,)
    end
end

function ragged2linearaxis(expandfn, m::jlModel, donor_labels, ragged_idxs)
    linearlabels = Symbol[]
    for (id, root) in enumerate(donor_labels)
        append!(linearlabels, expandfn(m, root, id)::Tuple)
    end
    linearlabels
end

const RAGGED_SIZENAME_TO_ADREXPANDER = Dict(
    :nq => (:jnt_qposadr, expandnq),
    :nv => (:jnt_dofadr, expandnv),
    #:nsensordata=>:sensor_adr, # can label this instead of _
    #:nnumericdata=>:numeric_adr
)

jl_id2name(m::jlModel, objtype::MJCore.mjtObj, id::Integer) =
    MJCore.mj_id2name(m, objtype, Int32(id - 1))
jl_type2szname(objtype::MJCore.mjtObj) =
    haskey(OBJTYPE_TO_SIZENAME, objtype) ? OBJTYPE_TO_SIZENAME[objtype] : nothing

makeunique(x::Vector{Symbol}; kwargs...) = (makeunique!(copy(x); kwargs...))
makeunique(x::Tuple{Vararg{Symbol}}; kwargs...) = Tuple(makeunique!(collect(x); kwargs...))

function makeunique!(x::AbstractVector{Symbol}; sep = :_)
    uniqs = unique(x)
    for uniq in uniqs
        idxs = findall(isequal(uniq), x)
        if length(idxs) > 1
            for (count, idx) in enumerate(idxs)
                x[idx] = Symbol(uniq, sep, count)
            end
        end
    end
    x
end

abstract type AbstractAxis end

struct LinearAxis <: AbstractAxis
    name::Symbol
    labels::Vector{Symbol}
end

struct UnnamedAxis <: AbstractAxis
    name::Symbol
    len::Int
end

Base.length(ax::LinearAxis) = length(ax.labels)
Base.length(ax::UnnamedAxis) = ax.len

# e.g. `:nsensor` -> MJAxis([:mysensor1, :mysensor2], [1:3, 4:4])
function sizenames2axes(m::jlModel)
    szname2axes = Dict{Symbol,LinearAxis}()
    sinfo = STRUCT_INFO[:mjModel]

    # Get the labels for each kind of object.
    for objtype in instances(MJCore.mjtObj)

        # get the szname corresponding to `objtype`
        szname = jl_type2szname(objtype)
        isnothing(szname) && continue # not all types have sizes (e.g. mjOBJ_UNKNOWN)
        sz = getproperty(m, szname) # get the integer value of `szname`

        # lookup the names for `objtype`
        elnames = map(1:sz) do id
            elname = jl_id2name(m, objtype, id)
            if isnothing(elname)
                # For unnamed objects, default to `typename_id` (e.g. :geom_3)
                Symbol(MJCore.mju_type2Str(objtype), :_, id)
            else
                Symbol(elname)
            end
        end

        # Now every dimension of size `szname` can be labeled with `elnames`.
        szname2axes[szname] = LinearAxis(szname, elnames)
    end

    # "Ragged" axes may inherit their element names from other "non-ragged" axes.
    # For example, the element names for "nv" axis come from "njnt".
    # This is because for system containing 2 joints, 1 slide and 1 free, njnt=2,
    # but nv = 7.
    for (ragged_szname, (addr_fieldname, expandfn)) in pairs(RAGGED_SIZENAME_TO_ADREXPANDER)
        # e.g. `jnt_qposadr` has size `njnt`
        donor_szname = first(sinfo[addr_fieldname].mjsize::Tuple{Symbol}) # for ragged_szname=nv, donor_szname=njnt
        ragged_idxs = getproperty(m, addr_fieldname) .+ 1 # convert to one-indxed
        donor_labels = szname2axes[donor_szname].labels

        linearlabels = ragged2linearaxis(expandfn, m, donor_labels, ragged_idxs)

        szname2axes[ragged_szname] = LinearAxis(ragged_szname, linearlabels)
    end

    # Arrays with dimension `na` correspond to stateful actuators. MuJoCo's
    # compiler requires that these are always defined after stateless actuators,
    # so we only need the final `na` elements in the list of all actuator names.
    if m.na != 0
        act_names = szname2axes[:nu].labels
        szname2axes[:na] = LinearAxis(:na, act_names[end-m.na+1:end])
    end

    # And a sanity check
    for (szname, axis) in szname2axes
        sz = getproperty(m, szname)
        @assert axis.name === szname
        @assert length(axis) == sz
    end

    # Lastly, add the labels for our custom axes
    for jlsize in JLSizes
        szname2axes[nameof(jlsize)] = LinearAxis(nameof(jlsize), collect(jlsize.labels))
    end

    szname2axes
end

# Associate axes in `sznames2axes` with each field of x
function associate_axes(m::jlModel, x::Union{jlModel,jlData}, sznames2axes)
    # Used to create the namedtuple.
    named_axes = Dict{Symbol,Tuple{Vararg{<:AbstractAxis}}}()
    sinfo = x isa jlModel ? STRUCT_INFO[:mjModel] : STRUCT_INFO[:mjData]

    for name in propertynames(x)
        # we only care about arrays that we wrap
        field = getproperty(x, name)
        if field isa Array && haskey(sinfo, name)
            finfo = sinfo[name]

            sznames = reverse(finfo.jlsize)
            sizes = map(sznames) do szname
                if szname isa Symbol
                    Int(hasproperty(m, szname) ? getproperty(m, szname) :
                        getproperty(x, szname))
                elseif szname isa JLSize
                    length(szname)
                else
                    szname
                end
            end

            sznames = map(enumerate(sznames)) do (i, dim)
                # replace Integer sizes with :dim_i
                if dim isa Symbol
                    dim
                elseif dim isa JLSize
                    nameof(dim)
                else
                    Symbol(:dim, i)
                end
            end

            uniquenames = makeunique(sznames)

            # `field`'s axes
            axes = map(sizes, sznames, uniquenames) do sz, szname, uniqname
                if haskey(sznames2axes, szname)
                    LinearAxis(uniqname, sznames2axes[szname].labels)
                else
                    UnnamedAxis(uniqname, sz)
                end
            end
            named_axes[name] = Tuple(axes)
        end
    end
    named_axes
end

function namify_array(A, newaxes)
    axisarray_axes = map(newaxes) do newax
        axname = newax.name
        if newax isa UnnamedAxis
            Axis{axname}(1:length(newax))
        else
            Axis{axname}(SVector{length(newax),Symbol}(newax.labels))
        end
    end
    newsize = length.(axisarray_axes)
    AxisArray(reshape(A, newsize), axisarray_axes)
end

function namify_one(x::Union{jlModel,jlData}, axismapping)
    newfields = []
    fieldnames = Symbol[]
    for fieldname in propertynames(x)
        field = getproperty(x, fieldname)
        if field isa Array && haskey(axismapping, fieldname)
            newaxes = axismapping[fieldname]
            push!(newfields, namify_array(field, newaxes))
            push!(fieldnames, fieldname)
        end
    end
    NamedTuple{Tuple(fieldnames)}(Tuple(newfields))
end

function namify(m::jlModel, d::jlData)
    sz2ax = sizenames2axes(m)
    d_axmappings = associate_axes(m, d, sz2ax)
    m_axmappings = associate_axes(m, m, sz2ax)
    namify_one(m, m_axmappings), namify_one(d, d_axmappings)
end


text:
# ================================
# Abstract types
# ================================

"""
    AbstractFlowpipe

Abstract type representing a flowpipe.

### Notes

A flowpipe is the set union of an array of reach-sets.
"""
abstract type AbstractFlowpipe end

"""
    basetype(T::Type{<:AbstractFlowpipe})

Return the base type of the given flowpipe type (i.e., without type parameters).

### Input

- `T` -- flowpipe type, used for dispatch

### Output

The base type of `T`.
"""
basetype(T::Type{<:AbstractFlowpipe}) = Base.typename(T).wrapper
basetype(fp::AbstractFlowpipe) = basetype(typeof(fp))

# LazySets interface: fallback behaves like UnionSetArray

"""
    LazySets.ρ(d::AbstractVector, fp::AbstractFlowpipe)

### Input

- `d`  -- direction
- `fp` -- flowpipe

### Output

The support function of the flowpipe along the given direction `d`.

### Notes

In this fallback implementation, the flowpipe behaves like the union of the
reach-sets, i.e. the implementation is analogue to that of a `LazySet.UnionSetArray`.
"""
function LazySets.ρ(d::AbstractVector, fp::AbstractFlowpipe)
    return map(Ri -> ρ(d, set(Ri)), array(fp)) |> maximum
end

"""
    LazySets.σ(d::AbstractVector, fp::AbstractFlowpipe)

### Input

- `d`  -- direction
- `fp` -- flowpipe

### Output

The support vector of the flowpipe along the given direction `d`.

### Notes

In this fallback implementation, the flowpipe behaves like the union of the
reach-sets, i.e. the implementation is analogue to that of a `LazySet.UnionSetArray`.
"""
function LazySets.σ(d::AbstractVector, fp::AbstractFlowpipe)
    return _σ_vec(d, array(fp))
end

"""
    LazySets.dim(fp::AbstractFlowpipe)

### Input

- `fp` -- flowpipe

### Output

An integer representing the ambient dimension of the flowpipe.
"""
function LazySets.dim(fp::AbstractFlowpipe)
    length(fp) > 0 || throw(ArgumentError("the dimension is not defined because this flowpipe is empty"))
    return dim(first(fp)) # assumes that the first set is representative
end

# iteration interface
@inline Base.iterate(fp::AbstractFlowpipe) = iterate(array(fp))
@inline Base.iterate(fp::AbstractFlowpipe, state) = iterate(array(fp), state)
@inline Base.length(fp::AbstractFlowpipe) = length(array(fp))
#@inline Base.size(fp::AbstractFlowpipe) = (length(array(fp)),)
@inline Base.first(fp::AbstractFlowpipe) = getindex(fp, 1)
@inline Base.last(fp::AbstractFlowpipe) = getindex(fp, lastindex(fp))
@inline Base.firstindex(fp::AbstractFlowpipe) = 1
@inline Base.lastindex(fp::AbstractFlowpipe) = length(array(fp))
@inline Base.eachindex(fp::AbstractFlowpipe) = eachindex(array(fp))

# fallback implementations of set getter functions

"""
    set(fp::AbstractFlowpipe)

Return the geometric set represented by this flowpipe as the union of reach-sets.

### Input

- `fp`  -- flowpipe

### Output

The set union of the array of reach-sets of the flowpipe.

## Notes

To retrieve the array of sets stored in the flowpipe use `array(fp)`. To get
a set at a particular index, use `set(F[ind])` or `set(F, ind)`.
"""
set(fp::AbstractFlowpipe) = UnionSetArray([set(R) for R in array(fp)])

"""
    set(fp::AbstractFlowpipe, ind::Integer)

Return the geometric set represented by this flowpipe at the given index.

## Input

- `fp`  -- flowpipe
- `ind` -- index (from `1` to `length(flowpipe)`)

## Output

The set wrapped by the flowpipe at the given index.
"""
set(fp::AbstractFlowpipe, ind::Integer) = set(getindex(array(fp), ind))

"""
    set(fp::AbstractFlowpipe, ind::AbstractVector)

Return the union of set represented by this flowpipe at the given indices.

## Input

- `fp`  -- flowpipe
- `ind` -- vector of indices

## Output

The set union stored in the flowpipe at the given indices.
"""
set(fp::AbstractFlowpipe, ind::AbstractVector) = UnionSetArray([set(fp, i) for i in ind])

# time domain interface

"""
    tstart(fp::AbstractFlowpipe)

Return the initial time of this flowpipe.

### Input

- `fp` -- flowpipe

### Output

A float representing the initial time of the given flowpipe. The fallback is
computed by taking the initial time of the first reach-set.
"""
@inline tstart(fp::AbstractFlowpipe) = tstart(first(fp))

"""
    tend(fp::AbstractFlowpipe)

Return the final time of this flowpipe.

### Input

- `fp` -- flowpipe

### Output

A float representing the initial time of the given flowpipe. The fallback is
computed by taking the final time of the last reach-set.
"""
@inline tend(fp::AbstractFlowpipe) = tend(last(fp))

"""
    tspan(fp::AbstractFlowpipe)

Return time span of this flowpipe.

### Input

- `fp` -- flowpipe

### Output

The interval representing the time span of the given flowpipe. The fallback
is computed as `(tstart(fp), tend(fp))`, see `tstart(::AbstractFlowpipe)` and
`tend(::AbstractFlowpipe)` for details.
"""
@inline tspan(fp::AbstractFlowpipe) = TimeInterval(tstart(fp), tend(fp))

"""
    vars(fp::AbstractFlowpipe)

Return the tuple of variable indices of the flowpipe.

### Input

- `fp` -- flowpipe

### Output

Tuple of integers with the variable indices of the flowpipe, typically ``1, 2, …, n``
where ``n`` is the dimension of the flowpipe.

### Notes

The fallback implementation assumes first reach-set is representative.
"""
vars(fp::AbstractFlowpipe) = vars(first(fp))

# indexing with ranges or with vectors of integers
Base.getindex(fp::AbstractFlowpipe, i::Int) = getindex(array(fp), i)
Base.getindex(fp::AbstractFlowpipe, i::Number) = getindex(array(fp), convert(Int, i))
Base.getindex(fp::AbstractFlowpipe, I::AbstractVector) = getindex(array(fp), I)

# ------------------------------
# Methods to check disjointness
# ------------------------------

# interface
isdisjoint(F1::AbstractFlowpipe, F2::AbstractFlowpipe, method::AbstractDisjointnessMethod=FallbackDisjointness()) = all(X -> isdisjoint(F1, X, method), array(F2))
@commutative isdisjoint(F::AbstractFlowpipe, Y::SetOrReachSet, method::AbstractDisjointnessMethod=FallbackDisjointness()) = all(X -> _is_intersection_empty(X, Y, method), array(F))

# flowpipe and vector of reach-sets TODO

# internal
_is_intersection_empty(F::AbstractFlowpipe, Y::SetOrReachSet, method::AbstractDisjointnessMethod) = all(X -> _is_intersection_empty(X, Y, method), array(F))

# ------------------------------
# Methods to check inclusion
# ------------------------------

Base.:⊆(F::AbstractFlowpipe, X::LazySet) = all(R ⊆ X for R in F)
Base.:⊆(F::AbstractFlowpipe, Y::AbstractLazyReachSet) = all(R ⊆ set(Y) for R in F)


text:
using EngEconomics, Roots, Plots

# Given


text:
using AutoTest

Package = "BlackBoxOptim"

using BlackBoxOptim

function run(packagename, srcdir = "src", testdir = "test";
    testfileregexp = r"^test_.*\.jl$",
    srcfileregexp = r"^.*\.jl$")

    testfiles = AutoTest.findfiles(testdir, testfileregexp; recursive = true) # in AutoTest this is false
    srcfiles = AutoTest.findfiles(srcdir, srcfileregexp; recursive = true)

    ts = AutoTest.TestSuite(testfiles, srcfiles, "$packagename test suite")

    AutoTest.runtestsuite(ts)

end

if length(ARGS) > 0 && ARGS[1] == "continuous"
    AutoTest.autorun(Package, "src", "test/autotests")
else
    run(Package, "src", "test/autotests")
end


text:
"""
    _build_node

Function to construct an internal node of a probabilistic circuit.
"""
function _build_node(::Type{N}, params::T, ns::AbstractNode...) where {T<:AbstractVector{<:Real},N<:NodeType}
    children = collect(ns)
    s_ = reduce(union, scope.(children))
    scope_ = length(s_) == 1 ? first(s_) : s_
    prop_ = NodeProperties[]

    check_smooth(children) && push!(prop_, Smoothness())
    check_decomposable(children) && push!(prop_, Decomposability())
    check_determinism(children) && push!(prop_, Determinism())
    supports_ = N == SumNode ? Dict(reduce(joinsupport, children)) : Dict(reduce(vcat, support.(children)))

    support_ = scope_ isa Int ? supports_[scope_] : [supports_[s] for s in scope_]

    return Node{typeof(scope_),typeof(support_),T,typeof((prop_...,)),N}(scope_, support_, children, params, (prop_...,))
end

# --
# Default Internal Nodes
# --

# product node
Prod(t::Type{<:AbstractNode}, partition::Vector) = Prod([t(s) for s in partition]...)
Prod(f::Function, partition::Vector) = Prod([f(s) for s in partition]...)
Prod(ns::AbstractNode...) = Prod(Float32, ns...)
function Prod(::Type{T}, ns::AbstractNode...) where {T<:Real}
    nodes = mapreduce(n -> isproduct(n) ? children(n) : n, vcat, ns)
    return _build_node(ProductNode, T[], nodes...)
end

isproduct(n::Node{V,S,T,P,ProductNode}) where {V,T,S,P} = true
isproduct(n::Node{V,S,T,P,N}) where {V,T,S,P,N} = false
isproduct(n::T) where {T<:AbstractLeaf} = false

logpdf(n::Node{T,V,S,P,ProductNode}, x) where {T,V,S,P} = sum(logpdf.(children(n),Ref(x)))

# sum node
Sum(t::Type{<:AbstractNode}, scope, K::Int) = Sum([t(scope) for k in 1:K]...)
Sum(f::Function, K::Int) = Sum([f(k) for k in 1:K]...)

Sum(ns::AbstractNode...) = Sum(Float32, ns...)
function Sum(::Type{T}, ns::AbstractNode...; init = (y) -> log.(rand(length(y)))) where {T<:Real}
    return _build_node(SumNode, T.(init(ns)), ns...)
end

"""
    _fast_logpdf(n, x)

More efficient evaluation of a sum node using log(sum(exp(x + log(w)))) based on:
[Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)
"""
function _fast_logpdf(n::Node{T,V,S,P,SumNode}, x) where {T,V,S,P}
    cs = children(n)
    xmax_r = mapreduce(i -> logpdf(cs[i], x) .+ n.params[i], _logsumexp_onepass_op, 1:length(cs))
    return first(xmax_r) .+ log1p.(last(xmax_r))
end
function logpdf(n::Node{T,V,S,P,SumNode}, x) where {T,V,S,P}
    return logsumexp(n.params[i]+logpdf(c,x) for (i,c) in enumerate(n.children))
end
logpdf(n::Node{T,V,S,P,SumNode}, x::AbstractMatrix) where {T,V,S,P} = _fast_logpdf(n,x)

# variational sum node
VISum(t::Type{<:AbstractNode}, scope, K::Int) = VISum([t(scope) for k in 1:K]...)
VISum(f::Function, K::Int) = VISum([f(k) for k in 1:K]...)

VISum(ns::AbstractNode...) = VISum(Float32, ns...)
function VISum(::Type{T}, ns::AbstractNode...; init = (y) -> log.(rand(length(y)))) where {T<:Real}
    return _build_node(VISumNode, T.(init(ns)), ns...)
end

"""
    _fast_logpdf(n, x)

More efficient evaluation of a sum node using log(sum(exp(x + log(w)))) based on:
[Sebastian Nowozin: Streaming Log-sum-exp Computation.](http://www.nowozin.net/sebastian/blog/streaming-log-sum-exp-computation.html)
"""
function _fast_logpdf(n::Node{T,V,S,P,VISumNode}, x) where {T,V,S,P}
    cs = children(n)
    logbeta = n.params
    logweights = digamma.(exp.(logbeta)) .- digamma.(exp(logsumexp(logbeta)))
    xmax_r = mapreduce(i -> logpdf(cs[i], x) .+ logweights[i], _logsumexp_onepass_op, 1:length(cs))
    return first(xmax_r) .+ log1p.(last(xmax_r))
end
function logpdf(n::Node{T,V,S,P,VISumNode}, x) where {T,V,S,P}
    cs = children(n)
    logbeta = n.params
    logweights = digamma.(exp.(logbeta)) .- digamma.(exp(logsumexp(logbeta)))
    return logsumexp(logweights[i]+logpdf(c,x) for (i,c) in enumerate(n.children))
end
logpdf(n::Node{T,V,S,P,VISumNode}, x::AbstractMatrix) where {T,V,S,P} = _fast_logpdf(n,x)

# --
# Default leaf nodes
# --

logpdfnormal(μ, σ, x) = -log(σ) - (log2π + ((x-μ) / σ)^2)/2

@leaf Normal(μ=0.0, σ=1.0) RealInterval(-Inf, Inf)
logpdf(n::Normal{P}, x::Real) where {P<:NamedTuple{(:μ, :σ)}} = logpdfnormal(n.params.μ, n.params.σ, x)

@leaf Indicator(v = 1)
logpdf(n::Indicator{P}, x::Real) where {P<:Real} = n.params == x ? 0 : -Inf
support(n::Indicator{P}) where {P<:Real} = (n.scope => [n.params])

@leaf TruncatedNormal(μ = 0.0, σ = 1.0, min = -Inf, max = Inf)
function logpdf(n::TruncatedNormal{P}, x::Real) where {P<:NamedTuple{(:μ, :σ, :min, :max)}}
    return (x >= n.params.min) && (x < n.params.max) ? logpdfnormal(n.params.μ, n.params.σ, x) : -Inf
end
function support(n::TruncatedNormal{P}) where {P<:NamedTuple{(:μ, :σ, :min, :max)}}
    return (n.scope => RealInterval(n.params.min, n.params.max))
end

# --
# Helper functions
# --

function descendants(n::Node{T,V,S,P,N}, N2::Type{<:NodeType}; level = Inf) where {T,V,S,P,N}
    l = N == N2 ? level-1 : level
    nodes = l > 0 ? mapreduce(c -> descendants(c, N2), vcat, n.children) : []
    N == N2 && push!(nodes, n)
    return unique(nodes)
end
descendants(::AbstractLeaf, ::Type{<:NodeType}; level = Inf) = []

leaves(n::Node) = unique(mapreduce(c -> leaves(c), vcat, n.children))
leaves(n::AbstractLeaf) = [n]


text:
# Significance tests
using Test
using CircStats

@testset "Tests" begin
    # Watson's U²n test for goodness-of-fit to a distribution
    # Test 96 of Kanji (2006) 100 Statistical Tests
    let θ = [20, 135, 145, 165, 170, 200, 300, 325, 335, 350, 350, 350, 355]
        fit_is_significant, U², U²crit = watson_U2n(θ, x->x/2π, 0.05, true)
        @test !fit_is_significant
        @test U² ≈ 0.1361 atol=0.0001
        @test U²crit ≈ 0.184 atol=0.002
    end

    # V-test (modified Rayleigh) to test if a distribution has a preferred direction
    # Test 95 of Kanji (2006)
    let θ = [250, 275, 285, 285, 290, 290, 295, 300, 305, 310, 315, 320, 330, 330, 5],
            θ₀ = 265, α = 0.01, V = 3.884, Vcrit = 2.302, is_not_random = true
        is_not_random′, V′, Vcrit′ = V_test(θ, θ₀, α, true)
        @test is_not_random′ == is_not_random
        @test V′ ≈ V atol=0.1
        @test Vcrit′ ≈ Vcrit atol=0.01
    end

    # Watson's U² test for the difference of two samples
    # Data and values from worked example from
    #    http://webspace.ship.edu/pgmarr/Geo441/Examples/Watsons%20U2%20Test.pdf
    let θ = [38,45,46,52,53,54,56,57,60,64],
            ϕ = [36,40,44,45,51,51,52,54,54,55,55,56,67,78,89,314], α = 0.05,
            reject = false, U² = 0.0427, U²crit = 0.1856
        reject′, U²′, U²crit′ = watson_U2(θ, ϕ, α, true)
        @test reject′ == reject
        @test U²′ ≈ U² atol=0.0001
        @test U²crit′ ≈ U²crit atol=0.02
    end

    # Data of Example 8.2 of Mardia and Jupp (2000); test values of Example 8.4
    # Currently fails because Mardia and Jupp inexplicably use the value of
    # U²crit for n = m = ∞, and their application of expression (8.3.8) does not
    # seem to give the same answer as as my application of (8.3.7).
    # let θ = [50, 290, 300, 300, 305, 320, 330, 330, 335, 340, 340, 355],
    #         ϕ = [70, 155, 190, 195, 215, 235, 235, 240, 255, 260, 290, 300, 300, 300],
    #         α = 0.01, U² = 0.320, U²crit = 0.268, reject = true
    #     reject′, U²′, U²crit′ = watson_U2(θ, ϕ, α, true)
    #     @test reject′ == reject
    #     println(U²′)
    #     @test U²′ ≈ U² atol=0.001
    #     println(U²crit′)
    #     @test U²crit′ ≈ U²crit atol=0.001
    # end
end


text:
# Constructors
# ============
#
# Constructor methods for Biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

function (::Type{BioSequence{A}}){A<:Alphabet}(len::Integer)
    return BioSequence{A}(Vector{UInt64}(seq_data_len(A, len)), 1:len, false)
end

BioSequence(::Type{DNA}) = DNASequence()
BioSequence(::Type{RNA}) = RNASequence()
BioSequence(::Type{AminoAcid}) = AminoAcidSequence()
BioSequence(::Type{Char}) = CharSequence()

function BioSequence()
    return BioSequence{VoidAlphabet}(Vector{UInt64}(0), 0:-1, false)
end

function (::Type{BioSequence{A}}){A<:Alphabet}(
        src::Union{AbstractString,AbstractVector},
        startpos::Integer=1,
        stoppos::Integer=length(src))
    len = stoppos - startpos + 1
    seq = BioSequence{A}(len)
    return encode_copy!(seq, 1, src, startpos, len)
end

# create a subsequence
function BioSequence{A,T<:Integer}(other::BioSequence{A}, part::UnitRange{T})
    checkbounds(other, part)
    start = other.part.start + part.start - 1
    stop = start + length(part) - 1
    subseq = BioSequence{A}(other.data, start:stop, true)
    other.shared = true
    return subseq
end

function (::Type{BioSequence{A}}){A}(other::BioSequence{A}, part::UnitRange)
    return BioSequence(other, part)
end

# concatenate chunks
function (::Type{BioSequence{A}}){A}(chunks::BioSequence{A}...)
    len = 0
    for chunk in chunks
        len += length(chunk)
    end
    seq = BioSequence{A}(len)
    offset = 1
    for chunk in chunks
        copy!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

function Base.repeat{A}(chunk::BioSequence{A}, n::Integer)
    seq = BioSequence{A}(length(chunk) * n)
    offset = 1
    for _ in 1:n
        copy!(seq, offset, chunk, 1)
        offset += length(chunk)
    end
    return seq
end

# operators for concat and repeat
Base.:*{A}(chunk::BioSequence{A}, chunks::BioSequence{A}...) =
    BioSequence{A}(chunk, chunks...)
Base.:^(chunk::BioSequence, n::Integer) = repeat(chunk, n)

function Base.similar{A}(seq::BioSequence{A}, len::Integer=length(seq))
    return BioSequence{A}(len)
end


text:
# This file was generated, do not modify it. # hide
knnc.K = 3
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
accuracy(ŷ, y[test]) |> r3


text:
using NonconvexSearch, Test

include("problems.jl")

@testset "test function sanity checks" begin
    for F in TEST_FUNCTIONS
        @test F(minimum_location(F, 10)) ≈ 0
    end
end

d_tol = Dict(SHIFTED_QUADRATIC=>1e-6, GRIEWANK=>1e-6, LEVY_MONTALVO2=>1e-6, RASTRIGIN=>1e-6, ROSENBROCK=>1e-1)
tol(F) = (d_tol[F])

test_dim = 2

@testset "MTS" begin
    # Temporary disable ROSENBROCK
    println("Testing MTS... ")
    for F in setdiff(TEST_FUNCTIONS, (ROSENBROCK, ))
        println("Testing nonconvex function: ", F)
        m = Model(x -> F(x))
        lb = [lu_bounds(F)[1] for _ in 1:test_dim]
        ub = [lu_bounds(F)[2] for _ in 1:test_dim]
        addvar!(m, lb, ub)
        alg = MTSAlg()
        r = optimize(m, alg, options=MTSOptions())
        println(r.minimizer)
        println(r.minimum)
        @test abs(r.minimum) < tol(F)
    end
end

@testset "Localsearch1" begin
    println("Testing Localsearch1... ")
    for F in setdiff(TEST_FUNCTIONS, (ROSENBROCK, ))
        println("Testing nonconvex function: ", F)
        m = Model(x -> F(x))
        lb = [lu_bounds(F)[1] for _ in 1:test_dim]
        ub = [lu_bounds(F)[2] for _ in 1:test_dim]
        addvar!(m, lb, ub)
        alg = LS1Alg()
        r = optimize(m, alg, options=LS1Options())
        println(r.minimizer)
        println(r.minimum)
        @test abs(r.minimum) < tol(F)
    end
end


text:
#!/usr/bin/env julia --project
#=
# This experiment file establishes the initial conditions, boundary conditions,
# source terms and simulation parameters (domain size + resolution) for
# a dry neutrally stratified Ekman layer.
#
# The initial conditions are given by constant horizontal velocity of 1 m/s,
# and a constant potential temperature profile. The bottom boundary condition
# results in momentum drag, and there is no exchange of heat from the surface
# since the fluxes are zero and the temperature is homogeneous.
#
=#

using ArgParse
using Distributions
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra
using Printf
using UnPack

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.ODESolvers
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.TurbulenceConvection
using ClimateMachine.VariableTemplates
using ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws: source, prognostic_vars

using CLIMAParameters
using CLIMAParameters.Planet: cp_d, cv_d, grav, T_surf_ref
using CLIMAParameters.Atmos.SubgridScale: C_smag, C_drag
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()
import CLIMAParameters

using ClimateMachine.Atmos: altitude, recover_thermo_state, density

"""
  EkmanLayer Geostrophic Forcing (Source)
"""
struct EkmanLayerGeostrophic{FT} <: TendencyDef{Source}
    "Coriolis parameter [s⁻¹]"
    f_coriolis::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
prognostic_vars(::EkmanLayerGeostrophic) = (Momentum(),)

function source(::Momentum, s::EkmanLayerGeostrophic, m, args)
    @unpack state, aux = args
    @unpack f_coriolis, u_geostrophic, u_slope, v_geostrophic = s

    z = altitude(m, aux)
    # Note z dependence of eastward geostrophic velocity
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    fkvector = f_coriolis * ẑ
    # Accumulate sources
    return -fkvector × (state.ρu .- state.ρ * u_geo)
end

"""
  EkmanLayer Sponge (Source)
"""
struct EkmanLayerSponge{FT} <: TendencyDef{Source}
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Sponge exponent"
    γ::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
prognostic_vars(::EkmanLayerSponge) = (Momentum(),)

function source(::Momentum, s::EkmanLayerSponge, m, args)
    @unpack state, aux = args

    @unpack z_max, z_sponge, α_max, γ = s
    @unpack u_geostrophic, u_slope, v_geostrophic = s

    z = altitude(m, aux)
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    # Accumulate sources
    if z_sponge <= z
        r = (z - z_sponge) / (z_max - z_sponge)
        β_sponge = α_max * sinpi(r / 2)^s.γ
        return -β_sponge * (state.ρu .- state.ρ * u_geo)
    else
        FT = eltype(state)
        return SVector{3, FT}(0, 0, 0)
    end
end

add_perturbations!(state, localgeo) = nothing

"""
  Initial Condition for EkmanLayer simulation
"""
function init_problem!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord
    # Problem floating point precision
    param_set = parameter_set(bl)
    FT = eltype(state)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v
    # Initialise speeds [u = Eastward, v = Northward, w = Vertical]
    u::FT = 1
    v::FT = 0
    w::FT = 0
    # Assign constant θ profile and equal to surface temperature
    θ::FT = T_surf_ref(param_set)

    p = aux.ref_state.p
    TS = PhaseDry_pθ(param_set, p, θ)

    compress = compressibility_model(bl) isa Compressible
    ρ = compress ? air_density(TS) : aux.ref_state.ρ
    # Compute momentum contributions
    ρu = ρ * u
    ρv = ρ * v
    ρw = ρ * w

    # Compute energy contributions
    e_kin = FT(1 // 2) * (u^2 + v^2 + w^2)
    e_pot = _grav * z
    ρe_tot = ρ * total_energy(e_kin, e_pot, TS)

    # Assign initial conditions for prognostic state variables
    state.ρ = ρ
    state.ρu = SVector(ρu, ρv, ρw)
    state.energy.ρe = ρe_tot
    add_perturbations!(state, localgeo)
    init_state_prognostic!(turbconv_model(bl), bl, state, aux, localgeo, t)
end

function ekman_layer_model(
    ::Type{FT},
    config_type,
    zmax,
    surface_flux;
    turbulence = ConstantKinematicViscosity(FT(0.1)),
    turbconv = NoTurbConv(),
    compressibility = Compressible(),
    ref_state = HydrostaticState(DryAdiabaticProfile{FT}(param_set),),
) where {FT}

    ics = init_problem!     # Initial conditions

    C_drag_::FT = C_drag(param_set) # FT(0.001)    # Momentum exchange coefficient
    u_star = FT(0.30)
    z_0 = FT(0.1)          # Roughness height

    z_sponge = FT(300)     # Start of sponge layer
    α_max = FT(0.75)       # Strength of sponge layer (timescale)
    γ = 2                  # Strength of sponge layer (exponent)

    u_geostrophic = FT(1)        # Eastward relaxation speed
    u_slope = FT(0)              # Slope of altitude-dependent relaxation speed
    v_geostrophic = FT(0)        # Northward relaxation speed
    f_coriolis = FT(1.39e-4) # Coriolis parameter at 73N

    q_sfc = FT(0)
    θ_sfc = T_surf_ref(param_set)
    g = compressibility isa Compressible ? (Gravity(),) : ()

    # Assemble source components
    source_default = (
        g...,
        EkmanLayerSponge{FT}(
            zmax,
            z_sponge,
            α_max,
            γ,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        EkmanLayerGeostrophic(
            f_coriolis,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        turbconv_sources(turbconv)...,
    )
    source = source_default

    # Set up problem initial and boundary conditions
    if surface_flux == "prescribed"
        energy_bc = PrescribedEnergyFlux((state, aux, t) -> FT(0))
    elseif surface_flux == "bulk"
        energy_bc = BulkFormulaEnergy(
            (bl, state, aux, t, normPu_int) -> C_drag_,
            (bl, state, aux, t) -> (θ_sfc, q_sfc),
        )
    elseif surface_flux == "custom_sbl"
        energy_bc = PrescribedTemperature((state, aux, t) -> θ_sfc)
    elseif surface_flux == "Nishizawa2018"
        energy_bc = NishizawaEnergyFlux(
            (bl, state, aux, t, normPu_int) -> z_0,
            (bl, state, aux, t) -> (θ_sfc, q_sfc),
        )
    else
        @warn @sprintf(
            """
%s: unrecognized surface flux; using 'prescribed'""",
            surface_flux,
        )
    end

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = turbulence,
        moisture = DryModel(),
        turbconv = turbconv,
        compressibility = compressibility,
    )

    moisture_bcs = ()
    boundary_conditions = (
        AtmosBC(
            physics;
            momentum = Impenetrable(DragLaw(
                # normPu_int is the internal horizontal speed
                # P represents the projection onto the horizontal
                (state, aux, t, normPu_int) -> (u_star / normPu_int)^2,
            )),
            energy = energy_bc,
            moisture_bcs...,
            turbconv = turbconv_bcs(turbconv)[1],
        ),
        AtmosBC(physics; turbconv = turbconv_bcs(turbconv)[2]),
    )

    problem = AtmosProblem(
        init_state_prognostic = ics,
        boundaryconditions = boundary_conditions,
    )

    # Assemble model components
    model =
        AtmosModel{FT}(config_type, physics; problem = problem, source = source)

    return model
end

function config_diagnostics(driver_config)
    default_dgngrp = setup_atmos_default_diagnostics(
        AtmosLESConfigType(),
        "60ssecs",
        driver_config.name,
    )
    core_dgngrp = setup_atmos_core_diagnostics(
        AtmosLESConfigType(),
        "60ssecs",
        driver_config.name,
    )
    return ClimateMachine.DiagnosticsConfiguration([
        default_dgngrp,
        core_dgngrp,
    ])
end


text:
# This script plots the objective
# function with different views

# Define the floating point type used across the script
F = Float64

# Define the integer type used across the script
I = UInt64

# Define default values for optional arguments
POSTFIX = ""

"Check if the value of the option is the last argument"
function check_last(i)
    if i + 1 == length(ARGS)
        println("The last argument should be the input directory.")
        exit(1)
    end
end

# Parse the options
for i in eachindex(ARGS)
    # A postfix for the names of output files
    if ARGS[i] == "--postfix"
        check_last(i)
        try
            global POSTFIX = " ($(ARGS[i+1]))"
        catch
            println("Couldn't parse the value of the `--postfix` argument.")
            exit(1)
        end
    end
end

# Check for required arguments
if length(ARGS) == 0
    println("""
        Usage:
        julia --project=. scripts/plot.jl [--postfix <POSTFIX>] <INPUT> """
    )
    exit(1)
end

# Define the input directory
INPUT_DIR = ARGS[end]

println('\n', " "^4, "> Loading the packages...")

using LaTeXStrings
using Plots

# Use the GR backend for plots
gr()

# Change some of the default parameters for plots
default(fontfamily = "Computer Modern", dpi = 300, legend = :topright, size = (400, 400))

# Define the paths to output directories
CURRENT_DIR = @__DIR__
ROOT_DIR = basename(CURRENT_DIR) == "scripts" ? dirname(CURRENT_DIR) : CURRENT_DIR
PLOTS_DIR = joinpath(ROOT_DIR, "plots")

# Make sure the needed directories exist
mkpath(PLOTS_DIR)

# Define the paths to the binary files
maximum_path = joinpath(INPUT_DIR, "maximum.bin")
point_path = joinpath(INPUT_DIR, "point.bin")
θ_path = joinpath(INPUT_DIR, "theta.bin")
φ_path = joinpath(INPUT_DIR, "phi.bin")
f_path = joinpath(INPUT_DIR, "obj.bin")

"Read binary files in the `bincode` format"
function read_bincode(path::AbstractString)::Tuple{I, Vector{F}}
    open(path, "r") do io
        n = read(io, I)
        a = reinterpret(F, read(io, sizeof(F) * n))
        n, a
    end
end

# Plot if the corresponding data files exist
if isfile(maximum_path) &&
    isfile(point_path) &&
    isfile(θ_path) &&
    isfile(φ_path) &&
    isfile(f_path)

    # Read the data
    _, maximum = read_bincode(maximum_path)
    _, point = read_bincode(point_path)
    n, θ = read_bincode(θ_path)
    _, φ = read_bincode(φ_path)
    _, f = read_bincode(f_path)

    # Normalize the vectors
    point ./= 2π
    θ ./= 2π
    φ ./= 2π

    println(" "^4, "> Plotting the heatmap...")

    # Plot the heatmap
    p = heatmap(
        θ,
        φ,
        reshape(f, (n, n));
        label = "",
        xlabel = L"\theta \; \mathrm{[2 \pi]}",
        ylabel = L"\varphi \; \mathrm{[2 \pi]}",
        zlabel = L"f",
    )

    # Point out the maximum
    scatter!([point[1]], [point[2]]; label = "maximum = $(maximum[1])", legend = :outertop, legendfontsize = 6)

    # Save the figure
    savefig(p, joinpath(PLOTS_DIR, "Heatmap$(POSTFIX).png"))

    println(" "^4, "> Plotting the surface plot...")

    # Plot the surface plot
    p = surface(
        θ,
        φ,
        reshape(f, (n, n));
        label = "",
        xlabel = L"\theta \; \mathrm{[2 \pi]}",
        ylabel = L"\varphi \; \mathrm{[2 \pi]}",
        zlabel = L"f",
    )

    # Point out the maximum
    scatter!([point[1]], [point[2]], maximum; label = "maximum = $(maximum[1])", legend = :outertop, legendfontsize = 6)

    # Save the figure
    savefig(p, joinpath(PLOTS_DIR, "Surface$(POSTFIX).png"))

    println(" "^4, "> Plotting the θ profile...")

    # Plot the surface plot
    p = surface(
        θ,
        φ,
        reshape(f, (n, n));
        camera = (0, 0),
        label = "",
        xlabel = L"\hspace{4.5} \theta \; \mathrm{[2 \pi]}",
        yaxis = false,
        yticks = false,
        zlabel = L"f",
        colorbar = false,
    )

    # Point out the maximum
    scatter!([point[1]], [point[2]], maximum; label = "maximum = $(maximum[1])", legend = :outertop, legendfontsize = 6)

    # Save the figure
    savefig(p, joinpath(PLOTS_DIR, "θ Profile$(POSTFIX).png"))

    println(" "^4, "> Plotting the φ profile...")

    # Plot the surface plot
    p = surface(
        θ,
        φ,
        reshape(f, (n, n));
        camera = (90, 0),
        label = "",
        xaxis = false,
        xticks = false,
        ylabel = L"\varphi \; \mathrm{[2 \pi]} \hspace{4.5}",
        zlabel = L"f",
        colorbar = false,
    )

    # Point out the maximum
    scatter!([point[1]], [point[2]], maximum; label = "maximum = $(maximum[1])", legend = :outertop, legendfontsize = 6)

    # Save the figure
    savefig(p, joinpath(PLOTS_DIR, "φ Profile$(POSTFIX).png"))
end

println()


text:
function add_dual_cone_constraint(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    ci::CI{F,S},
) where {F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}
    vi, con_index = MOI.add_constrained_variable(
        dual_model,
        _dual_set(get_set(primal_model, ci)),
    )
    return [vi], con_index
end

function add_dual_cone_constraint(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    ci::CI{F,MOI.EqualTo{T}},
) where {T,F<:MOI.AbstractScalarFunction}
    vi = MOI.add_variable(dual_model)
    return [vi], nothing
end

function add_dual_cone_constraint(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    ci::CI{F,S},
) where {F<:MOI.AbstractVectorFunction,S<:MOI.AbstractVectorSet}
    return MOI.add_constrained_variables(
        dual_model,
        _dual_set(get_set(primal_model, ci)),
    )
end

function add_dual_cone_constraint(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    ci::CI{F,MOI.Zeros},
) where {F<:MOI.AbstractVectorFunction}
    # Add as many variables as the dimension of the constraint
    vis = MOI.add_variables(dual_model, get_ci_row_dimension(primal_model, ci))
    return vis, nothing
end


text:
export BrickletIndustrialQuadRelayMonoflop
struct BrickletIndustrialQuadRelayMonoflop
    value::Integer
    time::Integer
    time_remaining::Integer
end

export BrickletIndustrialQuadRelayIdentity
struct BrickletIndustrialQuadRelayIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialQuadRelay
"""
4 galvanically isolated solid state relays
"""
mutable struct BrickletIndustrialQuadRelay <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialQuadRelay(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_quad_relay")
        deviceInternal = package.BrickletIndustrialQuadRelay(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_value
"""
    $(SIGNATURES)

Sets the output value with a bitmask (16bit). A 1 in the bitmask means relay
closed and a 0 means relay open.

For example: The value 3 or 0b0011 will close the relay of pins 0-1 and open
the other pins.

If no groups are used (see :func:`Set Group`), the pins correspond to the
markings on the Industrial Quad Relay Bricklet.

If groups are used, the pins correspond to the element in the group.
Element 1 in the group will get pins 0-3, element 2 pins 4-7, element 3
pins 8-11 and element 4 pins 12-15.

All running monoflop timers will be aborted if this function is called.
"""
function set_value(device::BrickletIndustrialQuadRelay, value_mask)
    device.deviceInternal.set_value(value_mask)
end

export get_value
"""
    $(SIGNATURES)

Returns the bitmask as set by :func:`Set Value`.
"""
function get_value(device::BrickletIndustrialQuadRelay)
    return device.deviceInternal.get_value()
end

export set_monoflop
"""
    $(SIGNATURES)

Configures a monoflop of the pins specified by the first parameter
bitmask.

The second parameter is a bitmask with the desired value of the specified
pins. A 1 in the bitmask means relay closed and a 0 means relay open.

The third parameter indicates the time that the pins should hold
the value.

If this function is called with the parameters (9, 1, 1500) or
(0b1001, 0b0001, 1500): Pin 0 will close and pin 3 will open. In 1.5s pin 0
will open and pin 3 will close again.

A monoflop can be used as a fail-safe mechanism. For example: Lets assume you
have a RS485 bus and a Quad Relay Bricklet connected to one of the slave
stacks. You can now call this function every second, with a time parameter
of two seconds and pin 0 closed. Pin 0 will be closed all the time. If now
the RS485 connection is lost, then pin 0 will be opened in at most two seconds.
"""
function set_monoflop(device::BrickletIndustrialQuadRelay, selection_mask, value_mask, time)
    device.deviceInternal.set_monoflop(selection_mask, value_mask, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns (for the given pin) the current value and the time as set by
:func:`Set Monoflop` as well as the remaining time until the value flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletIndustrialQuadRelay, pin)
    return device.deviceInternal.get_monoflop(pin)
end

export set_group
"""
    $(SIGNATURES)

Sets a group of Quad Relay Bricklets that should work together. You can
find Bricklets that can be grouped together with :func:`Get Available For Group`.

The group consists of 4 elements. Element 1 in the group will get pins 0-3,
element 2 pins 4-7, element 3 pins 8-11 and element 4 pins 12-15.

Each element can either be one of the ports ('a' to 'd') or 'n' if it should
not be used.

For example: If you have two Quad Relay Bricklets connected to port A and
port B respectively, you could call with ``['a', 'b', 'n', 'n']``.

Now the pins on the Quad Relay on port A are assigned to 0-3 and the
pins on the Quad Relay on port B are assigned to 4-7. It is now possible
to call :func:`Set Value` and control two Bricklets at the same time.
"""
function set_group(device::BrickletIndustrialQuadRelay, group)
    device.deviceInternal.set_group(group)
end

export get_group
"""
    $(SIGNATURES)

Returns the group as set by :func:`Set Group`
"""
function get_group(device::BrickletIndustrialQuadRelay)
    return device.deviceInternal.get_group()
end

export get_available_for_group
"""
    $(SIGNATURES)

Returns a bitmask of ports that are available for grouping. For example the
value 5 or 0b0101 means: Port A and port C are connected to Bricklets that
can be grouped together.
"""
function get_available_for_group(device::BrickletIndustrialQuadRelay)
    return device.deviceInternal.get_available_for_group()
end

export set_selected_values
"""
    $(SIGNATURES)

Sets the output value with a bitmask, according to the selection mask.
The bitmask is 16 bit long, *true* refers to a closed relay and
*false* refers to an open relay.

For example: The values (3, 1) or (0b0011, 0b0001) will close the relay of
pin 0, open the relay of pin 1 and leave the others untouched.

If no groups are used (see :func:`Set Group`), the pins correspond to the
markings on the Industrial Quad Relay Bricklet.

If groups are used, the pins correspond to the element in the group.
Element 1 in the group will get pins 0-3, element 2 pins 4-7, element 3
pins 8-11 and element 4 pins 12-15.

Running monoflop timers for the selected relays will be aborted if this function
is called.
"""
function set_selected_values(device::BrickletIndustrialQuadRelay, selection_mask, value_mask)
    device.deviceInternal.set_selected_values(selection_mask, value_mask)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialQuadRelay)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialQuadRelay, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
include("../imports.jl")

@testset "`Checkpointer`" begin
    checkpointer = Checkpointer(mktempdir())
    learner = testlearner(Recorder(), Metrics(), checkpointer)
    @test length(readdir(checkpointer.folder)) == 0
    @test_nowarn fit!(learner, 5)
    @test length(readdir(checkpointer.folder)) == 5
end


text:
# Pkg.add("DataFrames")
# Pkg.add("Dates")
Pkg.add("CSV")

using PlotlyJS, DataFrames, CSV, Dates
function dumbell()
    # reference: https://plot.ly/r/dumbbell-plots/
    # read Data into dataframe
    nm = tempname()
    url = "https://raw.githubusercontent.com/plotly/datasets/master/school_earnings.csv"
    download(url, nm)
    df = CSV.read(nm, DataFrame)
    rm(nm)

    df = sort(df, :Men, rev=false)

    men = PlotlyJS.scatter(;y=df[:School], x=df[:Men], mode="markers", name="Men",
                   marker=attr(color="blue", size=12))
    women = PlotlyJS.scatter(;y=df[:School], x=df[:Women], mode="markers", name="Women",
                     marker=attr(color="pink", size=12))

    lines = map(eachrow(df)) do r
        PlotlyJS.scatter(y=fill(r[:School], 2), x=[r[:Women], r[:Men]], mode="lines",
                name=r[:School], showlegend=false, line_color="gray")
    end

    data = Base.typed_vcat(GenericTrace, men, women, lines)
    layout = Layout(width=650, height=650, margin_l=100, yaxis_title="School",
                    xaxis_title="Annual Salary (thousands)",
                    title="Gender earnings disparity")

    plot(data, layout)
end
dumbell()

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations