text:
"""
    FourierDerivativeOperator{T<:Real, Grid, RFFT, BRFFT}

A derivative operator on a periodic grid with scalar type `T` computing the
first derivative using a spectral Fourier expansion via real discrete Fourier
transforms.
"""
struct FourierDerivativeOperator{T<:Real, Grid, RFFT, BRFFT} <: AbstractPeriodicDerivativeOperator{T}
    jac::T
    Δx::T
    grid_compute::Grid   # N-1 nodes, including the left and excluding the right boundary
    grid_evaluate::Grid #  N  nodes, including both boundaries
    tmp::Vector{Complex{T}}
    rfft_plan::RFFT
    brfft_plan::BRFFT

    function FourierDerivativeOperator(jac::T, Δx::T, grid_compute::Grid, grid_evaluate::Grid,
                                        tmp::Vector{Complex{T}}, rfft_plan::RFFT, brfft_plan::BRFFT) where {T<:Real, Grid, RFFT, BRFFT}
        @argcheck length(brfft_plan) == length(tmp) DimensionMismatch
        @argcheck length(brfft_plan) == (length(rfft_plan)÷2)+1 DimensionMismatch
        @argcheck length(grid_compute) == length(rfft_plan) DimensionMismatch
        @argcheck length(grid_compute) == length(grid_evaluate)-1 DimensionMismatch
        @argcheck first(grid_compute) == first(grid_evaluate)
        @argcheck step(grid_compute) ≈ step(grid_evaluate)
        @argcheck last(grid_compute) < last(grid_evaluate)

        new{T, Grid, RFFT, BRFFT}(jac, Δx, grid_compute, grid_evaluate, tmp, rfft_plan, brfft_plan)
    end
end

"""
    FourierDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}

Construct the `FourierDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N÷2+1` complex Fourier modes.
"""
function FourierDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}
    @argcheck N >= 1

    jac = 2*T(π) / (xmax - xmin) / N # / N because of brfft instead of BRFFT
    Δx = (xmax - xmin) / N
    grid_evaluate = range(xmin, stop=xmax, length=N+1) # two boundary nodes
    grid_compute = range(xmin, stop=grid_evaluate[end-1], length=N)
    u = zero.(grid_compute)
    rfft_plan = plan_rfft(u)
    uhat = rfft_plan*u
    brfft_plan = plan_brfft(uhat, N)

    FourierDerivativeOperator(jac, Δx, grid_compute, grid_evaluate, uhat, rfft_plan, brfft_plan)
end

function fourier_derivative_operator(xmin::T, xmax::T, N::Int) where {T<:Real}
    FourierDerivativeOperator(xmin, xmax, N)
end

derivative_order(D::FourierDerivativeOperator) = 1
LinearAlgebra.issymmetric(D::FourierDerivativeOperator) = false

function Base.show(io::IO, D::FourierDerivativeOperator{T}) where {T}
    grid = D.grid_evaluate
    print(io, "Periodic 1st derivative Fourier operator {T=", T, "} \n")
    print(io, "on a grid in [", first(grid), ", ", last(grid),
                "] using ", length(D.rfft_plan), " nodes and ",
                length(D.brfft_plan), " modes. \n")
end

function mul!(dest::AbstractVector{T}, D::FourierDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack jac, tmp, rfft_plan, brfft_plan = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        tmp[j] *= (j-1)*im * jac
    end
    @inbounds tmp[end] = 0
    mul!(dest, brfft_plan, tmp)

    nothing
end

"""
    fourier_derivative_matrix(N, xmin::Real=0.0, xmax::Real=2π)

Compute the Fourier derivative matrix with respect to the corresponding nodal
basis using `N` nodes, see
Kopriva (2009) Implementing Spectral Methods for PDEs, Algorithm 18.
"""
function fourier_derivative_matrix(N, xmin::Real=0.0, xmax::Real=2π)
    T = promote_type(typeof(xmin), typeof(xmax))
    jac_2 = T(π) / (xmax - xmin)
    D = Array{T}(N, N)
    @inbounds for j in 1:N, i in 1:N
        j == i && continue
        D[i,j] = (-1)^(i+j) * cot((i-j)*T(π)/N) * jac_2
        D[i,i] -= D[i,j]
    end
    D
end

"""
    ConstantFilter(D::FourierDerivativeOperator, filter, TmpEltype=T)

Create a modal filter with constant parameters adapted to the Fourier
derivative operator `D` with parameters given by the filter function `filter`.
"""
function ConstantFilter(D::FourierDerivativeOperator{T}, filter) where {T}
    Np1 = length(D.brfft_plan)
    coefficients = Array{T}(undef, Np1)
    set_filter_coefficients!(coefficients, filter)
    tmp = copy(D.tmp)
    modal2nodal = plan_irfft(D.tmp, length(D.rfft_plan))
    nodal2modal = D.rfft_plan
    ConstantFilter(coefficients, nodal2modal, modal2nodal, tmp, filter)
end

abstract type AbstractFourierViscosity{T} <: AbstractDerivativeOperator{T} end

@inline source_of_coefficients(Di::AbstractFourierViscosity) = (Di.source_of_coefficients)

LinearAlgebra.issymmetric(Di::AbstractFourierViscosity) = true
grid(Di::AbstractFourierViscosity) = grid(Di.D)

function mul!(dest::AbstractVector{T}, Di::AbstractFourierViscosity{T}, u::AbstractVector{T}) where {T}
    @unpack coefficients, D = Di
    @unpack tmp, rfft_plan, brfft_plan = D
    N = size(D, 1)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
        @argcheck length(tmp) == length(coefficients)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp))
        tmp[j] *= coefficients[j]
    end
    mul!(dest, brfft_plan, tmp)

    nothing
end

"""
    FourierConstantViscosity

Fourier viscosity operator with constant coefficients for the periodic 1st
derivative Fourier operator.
"""
struct FourierConstantViscosity{T<:Real, Grid, RFFT, BRFFT} <: AbstractFourierViscosity{T}
    coefficients::Vector{T}
    D::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}
    parameters
    source_of_coefficients

    function FourierConstantViscosity(D::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}, parameters, source_of_coefficients) where {T<:Real, Grid, RFFT, BRFFT}
        # precompute coefficients
        N = size(D,1)
        jac = N * D.jac # *N: brfft instead of irfft
        coefficients = Vector{T}(undef, length(D.brfft_plan))
        set_filter_coefficients!(coefficients, jac, N, parameters, source_of_coefficients)
        new{T,Grid,RFFT,BRFFT}(coefficients, D, parameters, source_of_coefficients)
    end
end

function Base.show(io::IO, Di::FourierConstantViscosity{T}) where {T}
    grid = Di.D.grid_evaluate
    print(io, "Fourier viscosity operator with constant coefficients for the periodic 1st\n")
    print(io, "derivative Fourier operator {T=", T, "} on a grid in [",
            first(grid), ", ", last(grid), "]\n")
    print(io, "using ", length(Di.D.rfft_plan), " nodes and ",
            length(Di.D.brfft_plan), " modes with coefficients from\n")
    print(io, Di.source_of_coefficients)
end

function dissipation_operator(source_of_coefficients, D::FourierDerivativeOperator; kwargs...)
    parameters = get_parameters(source_of_coefficients, D; kwargs...)
    FourierConstantViscosity(D, parameters, source_of_coefficients)
end

"""
    Tadmor1989

Coefficients of the Fourier spectral viscosity given in
  Tadmor (1989)
  Convergence of Spectral Methods for Nonlinear Conservation Laws.
  SIAM Journal on Numerical Analysis 26.1, pp. 30-44.
"""
struct Tadmor1989 <: SourceOfCoefficients end

function Base.show(io::IO, ::Tadmor1989)
    print(io,
        "  Tadmor (1989) \n",
        "  Convergence of Spectral Methods for Nonlinear Conservation Laws. \n",
        "  SIAM Journal on Numerical Analysis 26.1, pp. 30-44. \n")
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::Tadmor1989) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac
    end
end

"""
    MadayTadmor1989

Coefficients of the Fourier spectral viscosity given in
  Maday, Tadmor (1989)
  Analysis of the Spectral Vanishing Viscosity Method for Periodic Conservation
    Laws.
  SIAM Journal on Numerical Analysis 26.4, pp. 854-870.
"""
struct MadayTadmor1989 <: SourceOfCoefficients end

function Base.show(io::IO, ::MadayTadmor1989)
    print(io,
        "  Maday, Tadmor (1989) \n",
        "  Analysis of the Spectral Vanishing Viscosity Method for Periodic Conservation\n",
        "    Laws. \n",
        "  SIAM Journal on Numerical Analysis 26.4, pp. 854-870. \n")
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::MadayTadmor1989) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:min(2cutoff,length(coefficients))
        coefficients[j] = -strength * (j-1)^2 * jac * (j-cutoff)/cutoff
    end
    @inbounds @simd for j in 2cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac
    end
end

"""
    TadmorWaagan2012Standard

Coefficients of the Fourier spectral viscosity given in
  Tadmor, Waagan (2012)
  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws.
  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009.
"""
struct TadmorWaagan2012Standard <: SourceOfCoefficients end

function Base.show(io::IO, ::TadmorWaagan2012Standard)
    print(io,
        "  Tadmor, Waagan (2012) \n",
        "  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws. \n",
        "  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009. \n")
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::TadmorWaagan2012Standard) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff+1:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac * exp(-((length(coefficients)-j)/(j-cutoff))^2)
    end
end

"""
    TadmorWaagan2012Convergent

Coefficients of the Fourier spectral viscosity given in
  Tadmor, Waagan (2012)
  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws.
  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009.
See also
  Schochet (1990)
  The Rate of Convergence of Spectral-Viscosity Methods for Periodic Scalar
    Conservation Laws.
  SIAM Journal on Numerical Analysis 27.5, pp. 1142-1159.
"""
struct TadmorWaagan2012Convergent <: SourceOfCoefficients end

function Base.show(io::IO, ::TadmorWaagan2012Convergent)
    print(io,
        "  Tadmor, Waagan (2012) \n",
        "  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws. \n",
        "  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009. \n")
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::TadmorWaagan2012Convergent) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:min(2cutoff,length(coefficients))
        coefficients[j] = -strength * (j-1)^2 * jac * exp(-((2cutoff-j)/(j-cutoff))^2)
    end
    @inbounds @simd for j in 2cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac
    end
end

const FourierSpectralViscosityCoefficients = Union{Tadmor1989,MadayTadmor1989,TadmorWaagan2012Standard,TadmorWaagan2012Convergent}

function get_parameters(source_of_coefficients::FourierSpectralViscosityCoefficients,
                        D::FourierDerivativeOperator;
                        strength::Real=eltype(D)(1)/size(D,2),
                        cutoff::Int=1+round(Int, sqrt(size(D,2))), #+1: 1 based indexing
                        kwargs...)
    @argcheck cutoff >= 1

    Dict(:strength=>strength, :cutoff=>cutoff)
end

"""
    Tadmor1993

Coefficients of the Fourier super spectral viscosity given in
  Tadmor (1993)
  Super Viscosity and Spectral Approximations of Nonlinear Conservation Laws.
  Numerical Methods for Fluid Dynamics IV, pp. 69-82.
"""
struct Tadmor1993 <: SourceOfCoefficients end

function Base.show(io::IO, ::Tadmor1993)
    print(io,
        "  Tadmor (1993) \n",
        "  Super Viscosity and Spectral Approximations of Nonlinear Conservation Laws. \n",
        "  Numerical Methods for Fluid Dynamics IV, pp. 69-82. \n")
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::Tadmor1993) where {T<:Real}
    @unpack order, strength, cutoff = parameters
    @argcheck order >= 1
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2order / N # ^2order: order-th derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2order * jac * (1 - (cutoff/j)^2order)
    end
end

const FourierSuperSpectralViscosityCoefficients = Tadmor1993

function get_parameters(source::FourierSuperSpectralViscosityCoefficients,
                        D::FourierDerivativeOperator;
                        order::Int=1,
                        strength=eltype(D)(1)/size(D,2)^(2order-1),
                        cutoff::Int=1+round(Int, size(D,2)^(1-1/2order)), #+1: 1 based indexing
                        kwargs...)
    @argcheck order >= 1
    @argcheck cutoff >= 1

    Dict(:order=>order, :strength=>strength, :cutoff=>cutoff)
end


text:
function mdl_converter_ode!( device::DynInverter{AvgCnvFixedDC,O,VC,DC,P,F}) where {O <: OuterControl,
                                                   VC<: VSControl,
                                                   DC<: DCSource,
                                                   P <: FrequencyEstimator,
                                                   F <: Filter}

    #Obtain external states inputs for component

    #Obtain inner variables for component
     md = device.inner_vars[md_var]
     mq = device.inner_vars[mq_var]
    VDC = device.inner_vars[Vdc_var]

    #Get parameters

    #Obtain indices for component w/r to device

    #Define internal states for converter

    #Compute states... none?

    #Update inner_vars
    device.inner_vars[Vdcnv_var] = md*VDC
    device.inner_vars[Vqcnv_var] = mq*VDC
    #device.inner_vars[Vdcnv_var] = 1.02
    #device.inner_vars[Vqcnv_var] = 0.01
end

#TODO: Same as above, but:
## VDC is state of DC source
## IDC is algebraic inner_var output from Vcnv real output power

ext: jl


text:
using Test
using MeshCat
using MeshCatMechanisms
using RigidBodyDynamics
using RigidBodyDynamics.OdeIntegrators
using CoordinateTransformations: Translation
using ValkyrieRobot
using NBInclude
using StaticArrays

vis = Visualizer()

@testset "MeshCatMechanisms" begin
    @testset "URDF mechanism" begin
        urdf = joinpath(@__DIR__, "urdf", "Acrobot.urdf")
        robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
        mvis = MechanismVisualizer(robot, URDFVisuals(urdf), vis)
        if !haskey(ENV, "CI")
            open(mvis)
            wait(mvis)
        end
        set_configuration!(mvis, [0, -0.5])
        set_configuration!(mvis, findjoint(robot, "shoulder"), 1.0)
        @test configuration(mvis) == [1.0, -0.5]
        @test configuration(mvis, findjoint(robot, "shoulder")) == [1.0]
        copyto!(mvis, [0.1, -0.6, 0.0, 0.0])
        @test configuration(mvis) == [0.1, -0.6]

        setelement!(mvis, default_frame(bodies(robot)[end]))
        setelement!(mvis, Point3D(default_frame(bodies(robot)[3]), 0.2, 0.2, 0.2))

        @testset "simulation and animation" begin
            state = MechanismState(robot, randn(2), randn(2))
            t, q, v = simulate(state, 5.0);
            animate(mvis, t, q)
            animation = Animation(mvis, t, q)
            setanimation!(mvis, animation)
        end

        @testset "deprecated `setanimation!`" begin
            state = MechanismState(robot, randn(2), randn(2))
            t, q, v = simulate(state, 5.0);
            animate(mvis, t, q)
            setanimation!(mvis, t, q)
        end
    end

    @testset "URDF with fixed joints" begin
        urdf = joinpath(@__DIR__, "urdf", "Acrobot_fixed.urdf")
        robot = parse_urdf(urdf)
        delete!(vis)
        mvis = MechanismVisualizer(robot, URDFVisuals(urdf), vis)
        set_configuration!(mvis, [0.5])

        @testset "simulation and animation" begin
            state = MechanismState(robot, randn(1), randn(1))
            t, q, v = simulate(state, 5.0);
            animate(mvis, t, q)
            setanimation!(mvis, Animation(mvis, t, q))
        end
    end

    @testset "Valkyrie" begin
        val = Valkyrie();
        delete!(vis)
        mvis = MechanismVisualizer(
           val.mechanism,
           URDFVisuals(ValkyrieRobot.urdfpath(), package_path=[ValkyrieRobot.packagepath()]),
           vis[:val_urdf])
    end

    @testset "valkyrie inertias" begin
        val = Valkyrie();
        mvis = MechanismVisualizer(
           val.mechanism,
           Skeleton(),
           vis[:val_inertia])
        settransform!(vis[:val_inertia], Translation(0, 2, 0))
    end

    @testset "visualization during simulation" begin
        urdf = joinpath(@__DIR__, "urdf", "Acrobot.urdf")
        robot = parse_urdf(urdf, remove_fixed_tree_joints=false)
        mvis = MechanismVisualizer(robot, URDFVisuals(urdf), vis[:acrobot][:robot])
        settransform!(vis[:acrobot], Translation(0, -2, 0))
        result = DynamicsResult{Float64}(robot)
        function damped_dynamics!(vd::AbstractArray, sd::AbstractArray, t, state)
            damping = 2.
            τ = -damping * velocity(state)
            dynamics!(result, state, τ)
            copyto!(vd, result.v̇)
            copyto!(sd, result.ṡ)
            nothing
        end
        state = MechanismState(robot, randn(2), randn(2))
        integrator = MuntheKaasIntegrator(state, damped_dynamics!, runge_kutta_4(Float64), MeshCatSink(mvis))
        integrate(integrator, 10., 1e-3, max_realtime_rate = 1.)
    end

    @testset "position bounds resolution" begin
        j = Joint("foo", Revolute(SVector(1., 0, 0)))
        RigidBodyDynamics.position_bounds(j) .= RigidBodyDynamics.Bounds(-1.234, 0.0)
        MeshCatMechanisms.Manipulate.sliders(j)
        MeshCatMechanisms.Manipulate.widget(j)
    end

    @testset "manipulation" begin
        delete!(vis)
        mechanism = rand_chain_mechanism(Float64, [QuaternionFloating{Float64}; [Revolute{Float64} for i = 1:5]]...)
        mvis = MechanismVisualizer(mechanism)
        widget = manipulate!(mvis)
    end

    @testset "manipulation - negative bounds issue" begin
        mechanism = rand_chain_mechanism(Float64, Revolute{Float64})
        joint = first(joints(mechanism))
        joint.position_bounds .= RigidBodyDynamics.Bounds(-2.0, -1.0)
        mvis = MechanismVisualizer(mechanism)
        widget = manipulate!(mvis)
    end

    @testset "examples" begin
        dir = joinpath(@__DIR__, "..", "examples")
        for file in readdir(dir)
            base, ext = splitext(file)
            if ext == ".jl"
                if file == "manipulation_with_blink.jl"
                    # The interaction between MeshCat and Interact
                    # causes weird task deadlocks that I don't want
                    # to deal with right now.
                    @warn "Skipping $file"
                else
                    @info "Running $file"
                    include(joinpath(dir, file))
                end
            elseif ext == ".ipynb"
                @info "Running notebook $file"
                @nbinclude(joinpath(dir, file))
            end
        end
    end
end


text:
### Data structures ###
abstract type InputData end
abstract type TSData <:InputData end
abstract type OptData <: InputData end
abstract type AbstractClustResult end

struct FullInputData <: TSData
 region::String
 years::Array{Int,1}
 N::Int
 data::Dict{String,Array}
end

"""
      ClustData <: TSData

Contains time series data by attribute (e.g. wind, solar, electricity demand) and respective information.

Fields:
- region::String: optional information to specify the region data belongs to
- K::Int: number of periods
- T::Int: time steps per period
- data::Dict{String,Array}: Dictionary with an entry for each attribute `[file name (attribute: e.g technology)]-[column name (node: e.g. location)]`, Each entry of the dictionary is a 2-dimensional `time-steps T x periods K`-Array holding the data
- weights::Array{Float64,2}: 1-dimensional `periods K`-Array with the absolute weight for each period. The weight of a period corresponds to the number of days it representes. E.g. for a year of 365 days, sum(weights)=365
- mean::Dict{String,Array}: Dictionary with a entry for each attribute `[file name (e.g technology)]-[column name (e.g. location)]`, Each entry of the dictionary is a 1-dimensional `periods K`-Array holding the shift of the mean. This is used internally for normalization.
- sdv::Dict{String,Array}:  Dictionary with an entry for each attribute `[file name (e.g technology)]-[column name (e.g. location)]`, Each entry of the dictionary is a 1-dimensional `periods K`-Array holding the standard deviation. This is used internally for normalization.
- delta_t::Array{Float64,2}: 2-dimensional `time-steps T x periods K`-Array with the temporal duration Δt for each timestep. The default is that all timesteps have the same length.
- k_ids::Array{Int}: 1-dimensional `original periods I`-Array with the information, which original period is represented by which period K. E.g. if the data is a year of 365 periods, the array has length 365. If an original period is not represented by any period within this ClustData the entry will be `0`.
"""
struct ClustData <: TSData
 region::String
 years::Array{Int}
 K::Int
 T::Int
 data::Dict{String,Array}
 weights::Array{Float64}
 mean::Dict{String,Array}
 sdv::Dict{String,Array}
 delta_t::Array{Float64,2}
 k_ids::Array{Int}
end

"""
      ClustDataMerged <: TSData

Contains time series data by attribute (e.g. wind, solar, electricity demand) and respective information.

Fields:
- region::String: optional information to specify the region data belongs to
- K::Int: number of periods
- T::Int: time steps per period
- data::Array: Array of the dimension `(time-steps T * length(data_types)  x periods K`. The first T rows are data_type 1, the second T rows are data_type 2, ...
- data_type::Array{String}: The data types (attributes) of the data.
- weights::Array{Float64,2}: 1-dimensional `periods K`-Array with the absolute weight for each period. E.g. for a year of 365 days, sum(weights)=365
- mean::Dict{String,Array}: Dictionary with a entry for each attribute `[file name (e.g technology)]-[column name (e.g. location)]`, Each entry of the dictionary is a 1-dimensional `periods K`-Array holding the shift of the mean
- sdv::Dict{String,Array}:  Dictionary with an entry for each attribute `[file name (e.g technology)]-[column name (e.g. location)]`, Each entry of the dictionary is a 1-dimensional `periods K`-Array holding the standard deviation
- delta_t::Array{Float64,2}: 2-dimensional `time-steps T x periods K`-Array with the temporal duration Δt for each timestep in [h]
- k_ids::Array{Int}: 1-dimensional `original periods I`-Array with the information, which original period is represented by which period K. If an original period is not represented by any period within this ClustData the entry will be `0`.
"""
struct ClustDataMerged <: TSData
 region::String
 years::Array{Int}
 K::Int
 T::Int
 data::Array
 data_type::Array{String}
 weights::Array{Float64}
 mean::Dict{String,Array}
 sdv::Dict{String,Array}
 delta_t::Array{Float64,2}
 k_ids::Array{Int}
end

"""
    ClustResult <: AbstractClustResult
Contains the results from a clustering run: The data, the cost in terms of the clustering algorithm, and a config file describing the clustering method used.

Fields:
- clust_data::ClustData
- cost::Float64: Cost of the clustering algorithm
- config::Dict{String,Any}: Details on the clustering method used
"""
struct ClustResult <: AbstractClustResult
 clust_data::ClustData
 cost::Float64
 config::Dict{String,Any}
end

"""
    ClustResultAll <: AbstractClustResult
Contains the results from a clustering run for all locally converged solutions

Fields:
- clust_data::ClustData: The best centers, weights, clustids in terms of cost of the clustering algorithm
- cost::Float64: Cost of the clustering algorithm
- config::Dict{String,Any}: Details on the clustering method used
- centers_all::Array{Array{Float64},1}
- weights_all::Array{Array{Float64},1}
- clustids_all::Array{Array{Int,1},1}
- cost_all::Array{Float64,1}
- iter_all::Array{Int,1}
"""
struct ClustResultAll <: AbstractClustResult
 clust_data::ClustData
 cost::Float64
 config::Dict{String,Any}
 centers_all::Array{Array{Float64},1}
 weights_all::Array{Array{Float64},1}
 clustids_all::Array{Array{Int,1},1}
 cost_all::Array{Float64,1}
 iter_all::Array{Int,1}
end

"""
    SimpleExtremeValueDescr

Defines a simple extreme day by its characteristics

Fields:

- data_type::String : Choose one of the attributes from the data you have loaded into ClustData
- extremum::String : `min`,`max`
- peak_def::String : `absolute`,`integral`
- consecutive_periods::Int: For a single extreme day, set as 1
"""
struct SimpleExtremeValueDescr
  # TODO: make this one constructor, with consecutive_periods as optional argument
   data_type::String
   extremum::String
   peak_def::String
   consecutive_periods::Int
   "Replace default constructor to only allow certain entries"
   function SimpleExtremeValueDescr(data_type::String,
                                    extremum::String,
                                    peak_def::String,
                                    consecutive_periods::Int)
       # only allow certain entries
       if !(extremum in ["min","max"])
         error("extremum - "*extremum*" - not defined")
       elseif !(peak_def in ["absolute","integral"])
         error("peak_def - "*peak_def*" - not defined")
       end
       new(data_type,extremum,peak_def,consecutive_periods)
   end
end

"""
    SimpleExtremeValueDescr(data_type::String,
                                 extremum::String,
                                 peak_def::String)

Defines a simple extreme day by its characteristics

Input options:
- data_type::String : Choose one of the attributes from the data you have loaded into ClustData
- extremum::String : `min`,`max`
- peak_def::String : `absolute`,`integral`
"""
function SimpleExtremeValueDescr(data_type::String,
                                 extremum::String,
                                 peak_def::String)
   return SimpleExtremeValueDescr(data_type, extremum, peak_def, 1)
end

#### Constructors for data structures###

"""
    ClustData(region::String,
                      years::Array{Int,1},
                      K::Int,
                      T::Int,
                      data::Dict{String,Array},
                      weights::Array{Float64},
                      delta_t::Array{Float64,2},
                      k_ids::Array{Int,1};
                      mean::Dict{String,Array}=Dict{String,Array}(),
                      sdv::Dict{String,Array}=Dict{String,Array}()
                      )
constructor 1 for ClustData: provide data as dict
"""
function ClustData(region::String,
                       years::Array{Int,1},
                       K::Int,
                       T::Int,
                       data::Dict{String,Array},
                       weights::Array{Float64},
                       k_ids::Array{Int,1};
                       delta_t::Array{Float64,2}=ones(T,K),
                       mean::Dict{String,Array}=Dict{String,Array}(),
                       sdv::Dict{String,Array}=Dict{String,Array}()
                       )
 isempty(data) && error("Need to provide at least one input data stream")
 mean_sdv_provided = ( !isempty(mean) && !isempty(sdv))
 if !mean_sdv_provided
   for (k,v) in data
     mean[k]=zeros(T)
     sdv[k]=ones(T)
   end
 end
 # TODO check if right keywords are used
 ClustData(region,years,K,T,data,weights,mean,sdv,delta_t,k_ids)
end

"""
    ClustData(data::ClustDataMerged)
constructor 2: Convert ClustDataMerged to ClustData
"""
function ClustData(data::ClustDataMerged)
 data_dict=Dict{String,Array}()
 i=0
 for k in data.data_type
   i+=1
   data_dict[k] = data.data[(1+data.T*(i-1)):(data.T*i),:]
 end
 ClustData(data.region,data.years,data.K,data.T,data_dict,data.weights,data.mean,data.sdv,data.delta_t,data.k_ids)
end

"""
    ClustData(data::FullInputData,K,T)
constructor 3: Convert FullInputData to ClustData
"""
function ClustData(data::FullInputData,
                                 K::Int,
                                 T::Int)
  data_reshape = Dict{String,Array}()
  for (k,v) in data.data
     data_reshape[k] =  reshape(v,T,K)
  end
  return ClustData(data.region,data.years,K,T,data_reshape,ones(K),collect(1:K))
end

"""
    ClustDataMerged(region::String,
                        years::Array{Int,1},
                        K::Int,
                        T::Int,
                        data::Array,
                        data_type::Array{String},
                        weights::Array{Float64},
                        k_ids::Array{Int,1};
                        delta_t::Array{Float64,2}=ones(T,K),
                        mean::Dict{String,Array}=Dict{String,Array}(),
                        sdv::Dict{String,Array}=Dict{String,Array}()
                        )
constructor 1: construct ClustDataMerged
"""
function ClustDataMerged(region::String,
                       years::Array{Int,1},
                       K::Int,
                       T::Int,
                       data::Array,
                       data_type::Array{String},
                       weights::Array{Float64},
                       k_ids::Array{Int,1};
                       delta_t::Array{Float64,2}=ones(T,K),
                       mean::Dict{String,Array}=Dict{String,Array}(),
                       sdv::Dict{String,Array}=Dict{String,Array}()
                       )
 mean_sdv_provided = ( !isempty(mean) && !isempty(sdv))
 if !mean_sdv_provided
   for dt in data_type
     mean[dt]=zeros(T)
     sdv[dt]=ones(T)
   end
 end
 ClustDataMerged(region,years,K,T,data,data_type,weights,mean,sdv,delta_t,k_ids)
end

"""
    ClustDataMerged(data::ClustData)
constructor 2: convert ClustData into merged format
"""
function ClustDataMerged(data::ClustData)
 n_datasets = length(keys(data.data))
 data_merged= zeros(data.T*n_datasets,data.K)
 data_type=String[]
 i=0
 for (k,v) in data.data
   i+=1
   data_merged[(1+data.T*(i-1)):(data.T*i),:] = v
   push!(data_type,k)
 end
 if maximum(data.delta_t)!=1
   error("You cannot recluster data with different Δt")
 end
 ClustDataMerged(data.region,data.years,data.K,data.T,data_merged,data_type,data.weights,data.mean,data.sdv,data.delta_t,data.k_ids)
end


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

"""
    integrate(data, var; rank=nothing)

Integrate `data` for variable `var` over geometries of given `rank`.
Default rank is the parametric dimension of the underlying domain
where the data is georeferenced.
"""
function integrate(data::Data, vars...; rank=nothing)
  # domain and vertex table
  𝒟 = domain(data)
  𝒯 = values(data, 0)

  valid = Tables.schema(𝒯).names
  @assert vars ⊆ valid "invalid variables for vertex table"

  # vertices and topology
  vert = vertices(𝒟)
  topo = topology(𝒟)

  # retrieve columns
  cols = Tables.columns(𝒯)
  vals = [Tables.getcolumn(cols, var) for var in vars]

  # rank of integration
  R = isnothing(rank) ? paramdim(𝒟) : rank

  # loop over faces
  table = map(faces(topo, R)) do face
    # perform integration of all variables
    ints = _integrate(face, vert, vals)

    # row of table with results
    (; zip(vars, ints)...)
  end

  meshdata(𝒟, Dict(R => table))
end

# The surface integral ∫fdA over a 2D geometry can be
# expressed as ∫ᵤ∫ᵥf(u,v)||rᵤ×rᵥ||dudv where the vector
# r = [x(u,v), y(u,v), z(u,v)] lives on the geometry and
# where rᵤ = ∂r/∂u and rᵥ = ∂r/∂v are partial derivatives
# with respect to parameters u and v.
#
# For triangles, we can approximate functions linearly
# as f(u,v) = θ₀ + θ₁u + θ₂v using the values at the
# three vertices. This is a 3x3 linear system with
# analytical solution hard-coded below:
#
# |1 0 0| |θ₀|   |f₁|
# |1 1 0| |θ₁| = |f₂|
# |1 0 1| |θ₂|   |f₃|
#
# θ₀ = f₁, θ₁ = f₂-f₁, θ₂ = f₃-f₁
#
# f(u,v) = f₁ + (f₂-f₁)u + (f₃-f₁)v
#
# Coordinate functions can be approximated with the same
# system of equations (isometric approximation):
#
# x(u,v) = x₁ + (x₂-x₁)u + (x₃-x₁)v
# y(u,v) = y₁ + (y₂-y₁)u + (y₃-y₁)v
# z(u,v) = z₁ + (z₂-z₁)u + (z₃-z₁)v
#
# Consequently, we have the following constant:
#
# rᵤ = [(x₂-x₁), (y₂-y₁), (z₂-z₁)]
# rᵥ = [(x₃-x₁), (y₃-y₁), (z₃-z₁)]
#
# ||rᵤ×rᵥ|| = ||(p₂-p₁)×(p₃-p₁)|| = c
#
# where p₁, p₂ and p₃ are the three vertices.
#
# Finally, for the limits of integration u ∈ [0,1]
# and v ∈ [0,1-u] we can solve the integrand as:
#
# ∫ᵤ(∫ᵥf(u,v)dv)dv = ∫ᵤ(f₁(1-u) + (f₂-f₁)u(1-u) + (f₃-f₁)(1-u)²/2)du
#                  = ∫ᵤ((f₁+f₃)/2 + (f₂-f₃-f₁)u + (f₁+f₃-2f₂)u²/2)du
#                  = (f₁+f₃)/2 + (f₂-f₃-f₁)/2 + (f₁+f₃-2f₂)/6
#                  = (f₁+f₂+f₃)/6
#
# which leads to:
#
# ∫ᵤ∫ᵥf(u,v)||rᵤ×rᵥ||dudv = c(f₁+f₂+f₃)/6
function _integrate(face::Connectivity{<:Triangle}, vert, vals)
  i, j, k = indices(face)
  pᵢ, pⱼ, pₖ = vert[[i,j,k]]
  c = norm((pⱼ - pᵢ) × (pₖ - pᵢ))
  [c*(f[i] + f[j] + f[k])/6 for f in vals]
end

# fallback method ignores geometry and simply averages
# values of the variables at the vertices
function _integrate(face, vert, vals)
  inds = collect(indices(face))
  [mean(val[inds]) for val in vals]
end


text:
module Debugger

using Highlights
using Crayons
import InteractiveUtils

using Markdown
using Base.Meta: isexpr
using REPL
using REPL.LineEdit
using REPL.REPLCompletions

using CodeTracking
using JuliaInterpreter: JuliaInterpreter, Frame, @lookup, FrameCode, BreakpointRef, debug_command, leaf, root, BreakpointState,
                        finish_and_return!, Compiled

using JuliaInterpreter: pc_expr, moduleof, linenumber, extract_args,
                        root, caller, whereis, get_return, nstatements, getargs

const SEARCH_PATH = []
function __init__()
    append!(SEARCH_PATH,[joinpath(Sys.BINDIR,"../share/julia/base/"),
            joinpath(Sys.BINDIR,"../include/")])
    return nothing
end

export @enter, @run

# reexport some useful things from JuliaInterpreter
using JuliaInterpreter: @bp, @breakpoint, breakpoint, break_on, break_off
export @bp, @breakpoint, breakpoint, break_on, break_off, remove, disable

include("LineNumbers.jl")
using .LineNumbers: SourceFile, compute_line

# We make WATCH_LIST a global since it is likely useful to keep
# watch expressions between invocations of the debugger interface
const WATCH_LIST = []

Base.@kwdef mutable struct DebuggerState
    frame::Union{Nothing, Frame}
    level::Int = 1
    broke_on_error::Bool = false
    watch_list::Vector = WATCH_LIST
    lowered_status::Bool = false
    mode = finish_and_return!
    repl = nothing
    terminal = nothing
    main_mode = nothing
    julia_prompt::Ref{LineEdit.Prompt} = Ref{LineEdit.Prompt}()
    standard_keymap = nothing
    overall_result = nothing
end

function toggle_mode(state)
    state.mode = (state.mode === finish_and_return! ? (state.mode = Compiled()) : (state.mode = finish_and_return!))
end

toggle_lowered(state) = state.lowered_status = !state.lowered_status

function active_frame(state)
    frame = state.frame
    for i in 1:(state.level - 1)
        frame = caller(frame)
    end
    @assert frame !== nothing
    return frame
end

@deprecate break_on_error(v::Bool) (v ? break_on(:error) : break_off(:error))

maybe_quote(x) = (isa(x, Expr) || isa(x, Symbol)) ? QuoteNode(x) : x

include("locationinfo.jl")
include("repl.jl")
include("commands.jl")
include("printing.jl")
include("watch.jl")

function _make_frame(mod, arg)
    args = try
        extract_args(mod, arg)
    catch e
        return :(throw($e))
    end
    quote
        theargs = $(esc(args))
        frame = JuliaInterpreter.enter_call_expr(Expr(:call,theargs...))
        frame === nothing && error("failed to enter the function, perhaps it is set to run in compiled mode")
        frame = JuliaInterpreter.maybe_step_through_kwprep!(frame)
        frame = JuliaInterpreter.maybe_step_through_wrapper!(frame)
        JuliaInterpreter.maybe_next_call!(frame)
        frame
    end
end

_check_is_call(arg) = !(arg isa Expr && arg.head == :call) && throw(ArgumentError("@enter and @run should be applied to a function call"))

macro make_frame(arg)
    _make_frame(__module__, arg)
end

macro enter(arg)
    _check_is_call(arg)
    quote
        let frame = $(_make_frame(__module__, arg))
            RunDebugger(frame)
        end
    end
end

macro run(arg)
    _check_is_call(arg)
    quote
        let frame = $(_make_frame(__module__, arg))
            RunDebugger(frame; initial_continue=true)
        end
    end
end

end # module


text:
#=
3D advection using FV
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVadvection3d");

useLog("FVadvection3dlog")

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

# Mesh
n = 10 # number of elements in each direction
mesh(HEXMESH, elsperdim=n, bids=2)

# Variables and BCs
u = variable("u", location=CELL)
boundary(u, 1, FLUX, "sin(pi*y)*sin(pi*z)*sin(2*pi*t)^2") # x=0
boundary(u, 2, NO_BC) # other

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

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

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

solve(u)

finalize_finch()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using MAT
using Meganet

nImg = vec([32 32])
sK  =  vec([3 3 32 32])
nex = 100;

K1 = getConvGEMMKernel(Float64,nImg,sK)
K2 = getConvFFTKernel(Float64,nImg,sK)

theta = randn(nTheta(K1));

Y = zeros(tuple([nImg;sK[3];nex]...))
Y[2:end-1,2:end-1,:] = randn(tuple([nImg-2;sK[3];nex]...));

t1 = Amv(K1,theta,Y)
@time t1 = Amv(K1,theta,Y)

t2 = Amv(K2,theta,Y)
@time t2 = Amv(K2,theta,Y)

println(norm(t1[:]-t2[:])/norm(t1[:]))


text:
module TestVariableZeros

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_zeros()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Variable.Zeros{Float64}(mock)

    x, cx = MOI.add_constrained_variable(bridged_mock, MOI.GreaterThan(0.0))
    MOI.set(bridged_mock, MOI.VariableName(), x, "x")
    yz, cyz = MOI.add_constrained_variables(bridged_mock, MOI.Zeros(2))
    MOI.set(bridged_mock, MOI.VariableName(), yz, ["y", "z"])
    MOI.set(bridged_mock, MOI.ConstraintName(), cyz, "cyz")
    y, z = yz
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(bridged_mock, MOI.ObjectiveFunction{typeof(x)}(), x)
    @test MOI.get(bridged_mock, MOI.ObjectiveFunction{typeof(x)}()) == x

    # Test before adding affine constraints are affine expressions cannot be
    # unbridged when `Variable.ZerosBridge` is used.
    s = """
    variables: x, y, z
    x >= 0.0
    cyz: [y, z] in VecMathOptInterface.Zeros(2)
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        ["x", "y", "z"],
        ["cyz"],
        [("x", MOI.GreaterThan{Float64}(0.0))],
    )

    c1, c2 = MOI.add_constraints(
        bridged_mock,
        [1.0y + 1.0z, 1.0x + 1.0y + 1.0z],
        [MOI.EqualTo(0.0), MOI.GreaterThan(1.0)],
    )
    MOI.set(bridged_mock, MOI.ConstraintName(), c1, "con1")
    MOI.set(bridged_mock, MOI.ConstraintName(), c2, "con2")
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    obj = 1.0x - 1.0y - 1.0z
    MOI.set(bridged_mock, MOI.ObjectiveFunction{typeof(obj)}(), obj)

    @test MOI.Bridges.Variable.unbridged_map(
        MOI.Bridges.bridge(bridged_mock, y),
        y,
        MOI.Bridges.IndexInVector(1),
    ) === nothing
    @test MOI.Bridges.Variable.unbridged_map(
        MOI.Bridges.bridge(bridged_mock, z),
        z,
        MOI.Bridges.IndexInVector(2),
    ) === nothing

    err = ErrorException(
        "Cannot delete constraint index of bridged constrained variables. Delete" *
        " the scalar variable or the vector of variables instead.",
    )
    @test_throws err MOI.delete(bridged_mock, cyz)

    err = ErrorException(
        "Cannot unbridge function because some variables are bridged by" *
        " variable bridges that do not support reverse mapping, e.g.," *
        " `ZerosBridge`.",
    )
    @test_throws err MOI.get(bridged_mock, MOI.ObjectiveFunction{typeof(obj)}())
    # With `c1`, the function does not contain any variable so it tests that it
    # also throws an error even if it never calls `variable_unbridged_function`.
    @test_throws err MOI.get(bridged_mock, MOI.ConstraintFunction(), c1)
    @test_throws err MOI.get(bridged_mock, MOI.ConstraintFunction(), c2)

    err = ArgumentError(
        "Variable bridge of type `VecMathOptInterface.Bridges.Variable.ZerosBridge{Float64}`" *
        " does not support accessing the attribute `VecMathOptInterface.Test.UnknownVariableAttribute()`.",
    )
    @test_throws err MOI.get(
        bridged_mock,
        MOI.Test.UnknownVariableAttribute(),
        y,
    )

    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                0.0,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                1.0,
        ),
    )
    MOI.optimize!(bridged_mock)
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), x) == 1.0
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), y) == 0.0
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), z) == 0.0

    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cyz) == zeros(2)

    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cx) == 0.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c1) == 0.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c2) == 1.0
    bridge = VecMathOptInterface.Bridges.Variable.ZerosBridge{Float64}
    attr = MOI.ConstraintDual()
    err = ArgumentError(
        "Bridge of type `$(bridge)` does not support accessing " *
        "the attribute `$attr`. If you encountered this error " *
        "unexpectedly, it probably means your model has been " *
        "reformulated using the bridge, and you are attempting to query " *
        "an attribute that we haven't implemented yet for this bridge. " *
        "Please open an issue at https://github.com/jump-dev/MathOptInterface.jl/issues/new " *
        "and provide a reproducible example explaining what you were " *
        "trying to do.",
    )
    @test_throws err MOI.get(bridged_mock, attr, cyz)

    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), cyz).variables == yz
    @test MOI.get(mock, MOI.NumberOfVariables()) == 1
    @test MOI.get(mock, MOI.ListOfVariableIndices()) == [x]
    @test MOI.get(bridged_mock, MOI.NumberOfVariables()) == 3
    @test MOI.get(bridged_mock, MOI.ListOfVariableIndices()) == [x, y, z]
    @test MOI.get(
        mock,
        MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.Zeros}(),
    ) == 0
    @test MOI.get(
        bridged_mock,
        MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.Zeros}(),
    ) == 1
    @test MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Zeros}(),
    ) == [cyz]

    s = """
    variables: x
    x >= 0.0
    con1: 0.0 == 0.0
    con2: x + 0.0 >= 1.0
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        ["x"],
        ["con1", "con2"],
        [("x", MOI.GreaterThan{Float64}(0.0))],
    )

    _test_delete_bridged_variables(
        bridged_mock,
        yz,
        MOI.Zeros,
        3,
        ((MOI.VariableIndex, MOI.GreaterThan{Float64}, 1),),
    )
    @test MOI.is_valid(bridged_mock, x)
    @test !MOI.is_valid(bridged_mock, y)
    @test !MOI.is_valid(bridged_mock, z)
    return
end

end  # module

TestVariableZeros.runtests()


text:
# crossing monads

# Maybe => Either
either(x::None) = left(x)

# if x is a value, return right of x.  Else, return left of y.
right(x::T, y) where {T} = right(MaybeTypeTrait(T), x, y)
right(::IsJust, x, y) = right(x.value)
right(::IsNone, x, y) = left(y)
right(::IsVal,  x, y) = right(x)

# if x is a value, return left of x.  Elwse, return right of y.
left(x::T, y) where {T} = left(MaybeTypeTrait(T), x, y)
left(::IsJust, x, y) = left(x.value)
left(::IsNone, x, y) = right(y)
left(::IsVal,  x, y) = left(x)

# Either => Maybe

just(x::LeftEither) = NONE
just(x::RightEither) = just(right_value(x))


text:
using thesis, PGFPlotsX, DifferentialEquations, DifferentialEquations.EnsembleAnalysis, Statistics, ProgressMeter
setup_pgfplotsx()

## Helper functions
"""the expected firing rate of an lIF neuron with threshold `δ` and time-constant `α` in response to a constant input `x`"""
r_est_lif(x,δ,α)::Float64 = x ≤ α*δ ? 0.0 : -α/log(1-δ*α/x)

"""the expected firing rate of an LNP neuron with gain `λ` and time-constant `α` in response to a constant input `x`"""
r_est_lnp(x,λ,α)::Float64 = x ≤ 0.0 ? 0.0 : λ*x

"""the expected membrane potential of the decoding lIF neuron for an encoding lIF neuron with threshold `δ` and time-constant `α` in response to a constant input `x`"""
z_est_lif(x,δ,α) = δ*r_est_lif(x,δ,α)

"""the expected threshold of an lIF neuron with time-constant `α` that leads to a firing rate `r` in response to a constant input `x`"""
δ_est(x,r, α) = x/α*(1-exp(-α/r))

"""the expected gain of an LNP neuron with time-constant `α` that leads to a firing rate `r` in response to a constant input `x`"""
λ_est(x,r, α) = r/x

"""the RMSE expected for a ternary value in a finite `range` """
RMSE_ternary(range) = √(1/12)*range/3

"""the RMSE expected for a binary-encoded value in a finite `range` with fixed number of `bits` """
RMSE_digital(bits, range) = √(1/12)*range*2.0^(-bits)

"""the RMSE expected for an lIF neuron with threshold δ and time-constant α """
RMSE_lif(δ, α) = √(1/12)*α*δ

"""the threshold δ expected for an lIF neuron with time-constant α required for a fixed `RMSE`"""
RMSE_lif_inverse(RMSE, α) = √(12)*RMSE/α

"""the RMSE expected for an LNP neuron with gain `λ` and time-constant α """
RMSE_lnp(λ, α) = √(α/(2*λ))

"""the gain `λ` expected for an LNP neuron with time-constant α required for a fixed `RMSE`"""
RMSE_lnp_inverse(RMSE, α) = α/(2*RMSE^2)

"""the estimated rate at which the expected RMSE will equal that of an `nbits` encoding """
r_nbit_lif(s₀, nbits, α) = r_est_lif(s₀, RMSE_lif_inverse(RMSE_digital(nbits, 1.0), α), α)

"""the estimated rate at which the expected RMSE will equal that of a ternary encoding """
r_ternary_lif(α) = r_est_lif(s₀, RMSE_lif_inverse(RMSE_ternary(1.0), α), α)

## Configuration types

abstract type SNNConfig end;

struct lIFConfig{A,D,II,W,IW} <: SNNConfig
    num_neurons::Int
    α::A
    δ::D
    input_neurons::II
    input_weights::IW
    weights::W
    can_spike::Vector{Bool}
end

struct LNPConfig{A,L,II,W,IW} <: SNNConfig
    num_neurons::Int
    α::A
    λ₀::L
    input_neurons::II
    input_weights::IW
    weights::W
    can_spike::Vector{Bool}
end

## ODE and problem definition for lIF and LNP neurons

"""dynamics of leaky neurons (lIF and LNP) with exponential filter"""
function f′(du, u, p, t)
    du .= -p.cfg.α*u

    # work-around for jump problems (broadcasting doesnt work)
    tmp = p.cfg.input_weights .* p.input(t)
    for idx ∈ p.cfg.input_neurons
        du[idx] += tmp[idx]
    end
end

function lIF_spike_condition(out, u, t, integrator)
    out .= u .- integrator.p.cfg.δ
    # ignore neurons that cant spike
    out[(~).(can_spike)] .= Inf
end

function lIF_spike!(integrator, idx)
    integrator.u .+= integrator.p.cfg.weights[:,idx]

    has_spiked = zeros(Bool, integrator.p.cfg.num_neurons)
    has_spiked[idx] = true
    push!(integrator.p.spike_train, (integrator.t, idx))
    while true
        new_spikes = findall(integrator.p.cfg.can_spike .* (~).(has_spiked) .* (integrator.u .>= integrator.p.cfg.δ))
        if isempty(new_spikes)
            break
        end
        for i ∈ new_spikes
            integrator.u .+= integrator.p.cfg.weights[:,i]
            has_spiked[i] = true
            push!(integrator.p.spike_train, (integrator.t, i))
        end
    end
    integrator.u[has_spiked] .= 0.0
end

function make_LNP_jump(i)
    function LNP_rate(u,p,t)
        return u[i]*p.cfg.λ₀[i]
    end

    function LNP_spike!(integrator)
        integrator.u.u .+= integrator.p.cfg.weights[:,i]
        push!(integrator.p.spike_train, (integrator.t, i))
    end

    VariableRateJump(LNP_rate, LNP_spike!)
end

"""run an ensemble of lIF encoder-decoder-pair simulations over timespan `tspan` for a given `input` (function of time) and an lIFConfig `p`
Returns a solution object that tracks, for each trajectory, the decoding neuron's membrane potential (first trace) and spike rate (second trace)
"""
function runEnsemble(tspan, input, p::lIFConfig; trajectories=100, kwargs...)
    p.weights[2,1]=p.δ[1] * p.α

    spike_trains = [Tuple{Float64,Int}[] for i ∈ 1:trajectories]

    prob_func(prob,i,repeat) = ODEProblem(f′, rand(p.num_neurons).*p.δ, tspan, (spike_train=spike_trains[i], cfg=p, input=input))
    prob = prob_func(nothing,1,0)
    cb = VectorContinuousCallback(lIF_spike_condition, lIF_spike!, nothing, p.num_neurons)
    ensemble_prob = EnsembleProblem(
        prob,
        prob_func=prob_func;
        output_func=(sol,i)->(
            [
                sol(tspan[2])[2],
                length(filter(((t,n),)->(t>(tspan[2]+tspan[1])/2), spike_trains[i]))
            ],
            false
        )
    )
    sol = solve(ensemble_prob, Tsit5(), scheduler(); callback=cb, trajectories=trajectories, kwargs...)

    return hcat(sol...)'
end

"""run an ensemble of LNP encoder-decoder-pair simulations over timespan `tspan` for a given `input` (function of time) and an LNPConfig `p`
Returns a solution object that tracks, for each trajectory, the decoding neuron's membrane potential (first trace) and spike rate (second trace)
"""
function runEnsemble(tspan, input, p::LNPConfig; trajectories=100, kwargs...)
    p.weights[2,1]=p.α / p.λ₀[1]

    spike_trains = [Tuple{Float64,Int}[] for i ∈ 1:trajectories]

    prob_func(prob,i,repeat) = JumpProblem(
        ODEProblem(f′, rand(p.num_neurons), tspan, (spike_train=spike_trains[i], cfg=p, input=input)),
        Direct(),
        JumpSet((make_LNP_jump(j) for j ∈ 1:p.num_neurons if p.can_spike[j])...)
    )
    prob = prob_func(nothing,1,0)
    ensemble_prob = EnsembleProblem(
        prob,
        prob_func=prob_func;
        output_func=(sol,i)->(
            [
                sol(tspan[2])[2],
                length(filter(((t,n),)->(t>(tspan[2]+tspan[1])/2), spike_trains[i]))
            ],
            false
        )
    )
    sol = solve(ensemble_prob, Tsit5(), scheduler(); trajectories=trajectories, kwargs...)

    return hcat(sol...)'
end

"""simulates an ensemble of encoder-decoder pairs of neurons (lIF and LNP) for constant inputs with different `signal_levels`.
returns the `signal_levels`, means `μ` of lIF and LNP neurons and RMSE `err` of lIF and LNP neurons
"""
function calculate_transfer(tspan, signal_levels, lIF_cfg, LNP_cfg; kwargs...)
    μ = (lif=Float64[],lnp=Float64[])
    err = (lif=(neg=Float64[], pos=Float64[]),lnp=(neg=Float64[], pos=Float64[]))
    @showprogress 1.0 "Calculating transfer functions..." for level ∈ signal_levels
        for (i,cfg) ∈ enumerate((lIF_cfg, LNP_cfg))
            sol=runEnsemble(tspan, t->[level], cfg; kwargs...)

            quant = quantile(sol[:,1],[0.25,0.75])
            m = mean(sol[:,1])

            push!(err[i].neg,m-quant[1])
            push!(μ[i],m)
            push!(err[i].pos,quant[2]-m)
        end
    end
    return (signal_levels=signal_levels, μ=μ, err=err)
end

"""simulates an ensemble of encoder-decoder pairs of neurons (lIF and LNP) as well as a digital encoding for an `input` function.
returns digital encoding as well as decodings of lIF and LNP neurons
"""
function compare_digital_lif_lnp(times, input, (bits, rate), lIF_cfg, LNP_cfg; kwargs...)
    tspan = (minimum(times),maximum(times))
    spike_trains = (digital=Tuple{Float64,Int}[], lif=Tuple{Float64,Int}[], lnp=Tuple{Float64,Int}[])
    sample_points = minimum(times):1.0/rate:maximum(times)
    samples = input.(sample_points)

    digital = map(t->round(samples[searchsortedlast(sample_points, t)]*2^bits)/2^bits, times)

    prob = ODEProblem(f′, rand(lIF_cfg.num_neurons).*lIF_cfg.δ, tspan, (spike_train=spike_trains.lif, cfg=lIF_cfg, input=input))
    cb = VectorContinuousCallback(lIF_spike_condition, lIF_spike!, nothing, lIF_cfg.num_neurons)
    sol_lIF = solve(prob, Tsit5(); callback=cb, kwargs...)

    prob = JumpProblem(
        ODEProblem(f′, rand(LNP_cfg.num_neurons), tspan, (spike_train=spike_trains.lnp, cfg=LNP_cfg, input=input)),
        Direct(),
        JumpSet((make_LNP_jump(i) for i ∈ 1:LNP_cfg.num_neurons if LNP_cfg.can_spike[i])...)
    )
    sol_LNP = solve(prob, Tsit5(); kwargs...)

    return (solution=(digital=digital, lif=sol_lIF(times, idxs=2)[:].+0.5*lIF_cfg.δ[1]*lIF_cfg.α, lnp=sol_LNP(times, idxs=2)[:]), spikes=spike_trains)
end

"""simulates an ensemble of encoder-decoder pairs of neurons (lIF and LNP) for a fixed `input` function with different `scales` of firing rates (adjust δ or λ accordingly).
returns the `scales`, mean firing rates `rate` of lIF and LNP neurons and RMSE `rmse` of lIF and LNP neurons
"""
function calculate_error_vs_spikes!(tspan, scales, input, lIF_cfg, LNP_cfg; kwargs...)
    rate = (lif=Float64[],lnp=Float64[])
    rmse = (lif=Float64[],lnp=Float64[])
    @showprogress 1.0 "Calculating error/rate relation..." for s ∈ scales
        for (i,cfg) ∈ enumerate((lIF_cfg,LNP_cfg))
            if i==1
                cfg.δ[1] = δ_est(input(0)[1],s,cfg.α)
            else
                cfg.λ₀[1] = λ_est(input(0)[1],s,cfg.α)
            end

            sol=runEnsemble(tspan, input, cfg; kwargs...)

            push!(rate[i], mean(sol[:,2]./(tspan[2]-tspan[1])*2))

            v = sol[:,1]

            if i==1
                v.+= cfg.α.*cfg.δ[1]./2
            end
            push!(rmse[i], √(mean((v .- input(tspan[2])).^2)))
        end
    end

    return (scales=scales, rate=rate, rmse=rmse)
end

## Set up
use_multithreading = false                           # whether to use multithreading or not
num_neurons = 2                                     # we always simulate 2 neurons here
α=1000.0/25.0                                       # α = 40 corresponds to a time-scale of 25ms
δ=[δ_est(0.5,r_nbit_lif(0.5,4,α),α), 2.0]               # δ = 0.0015625 gives similar RMSE at a signal value 0.5 as a 4-bit encoding would
λ₀=1.0 ./ δ                                         # λ = 1/δ results in the same mean firing rate as for the lIF neuron
can_spike = [true,false]                            # only the encoding neuron can spike
input_neurons=[1]                                   # only the encoding neuron receives input
input_weights_lIF=[1.0]                             # the input-weights for the lIF neuron are 1
input_weights_LNP=input_weights_lIF.*α              # the input-weights for the LNP neuron are scaled up by α
weights=[0 0; 1 0].*δ'.*α                           # the weight from encoder to decode is scaled by αδ
trajectories = 50                                  # we simulate an ensemble of multiple trajectories to determine the distribution
signal_levels = LinRange(-0.2,1.0, 100)             # we evaluate the input-output function of the neurons for constant input signals with various levels
scales = 10.0 .^ LinRange(0, 4, 10)                 # we evaluate the RMSE for different parameterizations by sweeping over multiple scales (changes either δ or λ, respectively)
s₀ = 0.5                                            # we use a constant signal s(t)=0.5 to evaluate different parameterizations
equivalent_sampling_rate = α/2                      # a rectangular filter with equivalent delay as the exponential filter with time-constant α has width 2/α -> rate α/2
reference_bits = [2,4,6,8]                          # plot refrence lines for digital encodings using different numbers of bits
tspan = (0.0,10.0)                                  # we simulate for 10 seconds and take the last half of that for calculating averages
sine(t) = 0.4*sin(0.8*π*t)+0.5                      # we visualize the encoding for a sine-wave
plot_times = LinRange(-1,2.5,1000)                  # we plot the sine wave in the interval 0 to 10

inset_xmin,inset_xmax = 0.3,0.5
inset_ymin,inset_ymax = 0.7,1.0
plot_times_zoom = LinRange(inset_xmin,inset_xmax,200)

# the decoded signal will approximate the filtered version of the input signal
#∫ᵗsin(ωτ)α exp(-α(t-τ))dτ
#   = α exp(-αt)/(2i)∫ᵗ(exp((α+iω)τ)-exp((α-iω)τ))dτ
#   = α exp(-αt)/(2i)[1/(α+iω)exp((α+iω)t)-1/(α-iω)exp((α-iω)t)]
#   = α 1/(2iα-2ω)exp(iωt)-1/(2iα+2ω)exp(-iωt)
#   = α (α sin(ωt) - ω cos(ωt))/(α² + ω²)
sine_filtered(t) = 0.4*(α*sin(0.8*π*t)-0.8*π*cos(0.8*π*t)) * α/(α^2+(0.8π)^2)+0.5

##

scheduler = use_multithreading ? EnsembleThreads : EnsembleSerial
lIF_cfg = lIFConfig(num_neurons, α, copy(δ), input_neurons, copy(input_weights_lIF), copy(weights), can_spike)
LNP_cfg = LNPConfig(num_neurons, α, copy(λ₀), input_neurons, copy(input_weights_LNP), copy(weights), can_spike)

# Calculate the input-output transfer function
(s, μ, err) = calculate_transfer(tspan, signal_levels, lIF_cfg, LNP_cfg; trajectories=trajectories)

# Calculate input and decoding of a sine-wave
(decoded,spikes) = compare_digital_lif_lnp(plot_times, sine, (4, equivalent_sampling_rate), lIF_cfg, LNP_cfg; trajectories=trajectories)

# Calculate the RMSE/firing rate relationship as a function of varying δ or λ
(scales, rate, rmse) = calculate_error_vs_spikes!(tspan, scales, t->[s₀], lIF_cfg, LNP_cfg; trajectories=trajectories)

ss = 10 .^ LinRange(log10(minimum(scales)), log10(maximum(scales)), 200)

reference_rmse = [RMSE_digital(bits, 1.0) for bits in reference_bits]
reference_rate = [r_nbit_lif(s₀, bits,α) for bits in reference_bits]
reference_text = ["\\approx $(bits)-bit @ \$$(equivalent_sampling_rate)\$Hz " for bits in reference_bits]

t0=searchsortedfirst(plot_times,0)

## Plot the monstrosity
@pgf gp = TikzPicture(
    {remember_picture},
    GroupPlot(
    {
        group_style={
            group_size="2 by 3",
            horizontal_sep="2cm",
            vertical_sep="2cm",
        },
        height="6cm",
        width="8cm",
    },
    {xlabel=raw"input value $c$", ylabel=raw"rate$(c)$", title=raw"decoded rate", legend_pos="north west", xmin=signal_levels[1], xmax=signal_levels[end]},
    Plot(
        {"name path=lnplow", draw="none", forget_plot},
        Table(signal_levels, μ.lnp .- err.lnp.neg)
    ),
    Plot(
        {"name path=lnphigh", draw="none", forget_plot},
        Table(signal_levels, μ.lnp .+ err.lnp.pos)
    ),
    PlotInc(
        {color="Set1-B", opacity=0.5, forget_plot},
        raw"fill between [of=lnplow and lnphigh]"
    ),
    PlotInc(
        {color="Set1-B", very_thick},
        Table(signal_levels, μ.lnp)
    ),
    LegendEntry("LNP response"),
    Plot(
        {"name path=liflow", draw="none", forget_plot},
        Table(signal_levels, μ.lif .- err.lif.neg)
    ),
    Plot(
        {"name path=lifhigh", draw="none", forget_plot},
        Table(signal_levels, μ.lif .+ err.lif.pos)
    ),
    PlotInc(
        {color="Set1-A", forget_plot},
        raw"fill between [of=liflow and lifhigh]"
    ),
    PlotInc(
        {color="Set1-A", very_thick},
        Table(signal_levels, μ.lif)
    ),
    LegendEntry("lIF response"),
    PlotInc(
        {dashed, very_thick},
        Table(signal_levels, max.(0,signal_levels.- 0.5*α*δ[1]))
    ),
    LegendEntry(raw"shifted ReLU"),
    PlotInc(
        {dotted, very_thick},
        Table(signal_levels, z_est_lif.(signal_levels, δ[1], α))
    ),
    LegendEntry(raw"expected lIF response"),
    {xlabel=raw"firing rate $r$ [Hz]", ylabel=raw"RMSE$(r)$", title="RMSE for constant input \$x(t) = $(s₀)\$", ymode="log", xmode="log",xmin=ss[1],xmax=ss[end]},
    (Plot(
        {forget_plot, thick, color= i==2 ? "Set1-C" : "black!30"},
        Coordinates([(ss[1],e),(ss[end],e)])
    ) for (i,e) ∈ enumerate(reference_rmse))...,
    PlotInc(
        {forget_plot, very_thick, color="Set1-B"},
        Table(ss, RMSE_lnp.(λ_est.(s₀,ss,α),α))
    ),
    PlotInc(
        {only_marks, very_thick, color="Set1-B"},
        Table(rate.lnp, rmse.lnp)
    ),
    LegendEntry(raw"LNP neuron"),
    PlotInc(
        {forget_plot, very_thick, color="Set1-A"},
        Table(ss, RMSE_lif.(δ_est.(s₀,ss,α),α))
    ),
    PlotInc(
        {only_marks, very_thick, color="Set1-A"},
        Table(rate.lif, rmse.lif)
    ),
    LegendEntry(raw"lIF neuron"),
    Plot(
        {only_marks, mark="+", color="black"},
        Table(reference_rate, reference_rmse)
    ),
    ["\\node[anchor=north east, color=$(i==2 ? "Set1-C" : "black!50")] at (axis cs:$(r),$(e)) {$(t)};" for (i,(r,e,t)) ∈ enumerate(zip(reference_rate, reference_rmse, reference_text))],
    {width="16.5cm", xshift="4.25cm", title="digital code vs. lIF neuron vs. LNP neuron", ylabel="decoded signal", xticklabels="\\empty", xmin=0,xmax=2.5,ymin=0, ymax=2.0},
    "\\coordinate (boxsw) at (axis cs:$(inset_xmin),$(inset_ymin));",
    "\\coordinate (boxse) at (axis cs:$(inset_xmax),$(inset_ymin));",
    "\\coordinate (boxne) at (axis cs:$(inset_xmax),$(inset_ymax));",
    "\\coordinate (boxnw) at (axis cs:$(inset_xmin),$(inset_ymax));",
    "\\filldraw[color=black!90, very thick, fill=black!10] (boxsw) rectangle (boxne);",
    PlotInc(
        {very_thick, color="Set1-B"},
        Table(plot_times,decoded.lnp)
    ),
    LegendEntry("LNP decoding ($(length(filter(>(0), first.(spikes.lnp)))) spikes)"),
    PlotInc(
        {very_thick, color="Set1-A"},
        Table(plot_times,decoded.lif)
    ),
    LegendEntry("lIF decoding ($(length(filter(>(0), first.(spikes.lif)))) spikes)"),
    PlotInc(
        {very_thick},
        Table(plot_times,decoded.digital)
    ),
    LegendEntry("digital decoding ($(sum(sum.(thesis.discretize.(decoded.digital, Ref(LinRange(0,1,2^4+1)), Ref(digits.(0:2^4,base=2,pad=4)))[t0:end])))) bits)"),
    Plot(
        {ultra_thick, dashed},
        Table(plot_times,sine_filtered.(plot_times))
    ),
    LegendEntry(raw"target signal $(s*\kappa_\alpha)(t)$"),
    raw"\coordinate (insetpos) at (axis cs:0.1,1.1);",
    {"group/empty plot", "axis background/.style={fill=none}", xmin=0,xmax=1,ymin=0,ymax=1},
    {yshift="1.75cm", width="16.5cm", height="3cm", xlabel="time [s]", xmin=0,xmax=2.5,ymin=0, ymax=2.0, ylabel="spikes", ytick="\\empty"},
    (Plot(
        {color="Set1-A", forget_plot, opacity=0.5},
        Coordinates([(t,0.05),(t,0.95)])
    ) for t ∈ first.(spikes.lif))...,
    (Plot(
        {color="Set1-B", forget_plot, opacity=0.5},
        Coordinates([(t,1.05),(t,1.95)])
    ) for t ∈ first.(spikes.lnp))...,
    ),
    Axis(
        {very_thick, height="1.75cm", width="5cm", scale_only_axis, "axis background/.style"={fill="black!10", opacity=0.9},  xmin=inset_xmin,xmax=inset_xmax, ymin=inset_ymin, ymax=inset_ymax, at={"(insetpos)"}, anchor={"south west"}, ticks="none"},
        PlotInc(
            {ultra_thick, color="Set1-B"},
            Table(plot_times,decoded.lnp)
        ),
        PlotInc(
            {ultra_thick, color="Set1-A"},
            Table(plot_times,decoded.lif)
        ),
        PlotInc(
            {ultra_thick},
            Table(plot_times,decoded.digital)
        ),
        Plot(
            {ultra_thick, dashed},
            Table(plot_times,sine_filtered.(plot_times))
        ),
        "\\coordinate (insw) at (axis cs:$(inset_xmin),$(inset_ymin));",
        "\\coordinate (inse) at (axis cs:$(inset_xmax),$(inset_ymin));",
        "\\coordinate (inne) at (axis cs:$(inset_xmax),$(inset_ymax));",
        "\\coordinate (innw) at (axis cs:$(inset_xmin),$(inset_ymax));",
    ),
    raw"\draw[color=black!75, very thick] (boxsw) -- (insw);",
    raw"\draw[color=black!75, very thick] (boxse) -- (inse);",
)
#pgfsave("fig/encoding_comparison.pdf",gp)

max_stars_repo_name: jleugeri/phd

max_issues_repo_name: jleugeri/phd

max_forks_repo_name: jleugeri/phd


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7cd44a22-91f0-11eb-056d-699740ca7da9
#Set up the working environment
begin
	using Pkg
	Pkg.activate("..")
	using Plots, PlutoUI,LaTeXStrings, DifferentialEquations
end

# ╔═╡ 3199bb4e-91ed-11eb-37fa-abed99588e6b
md"""
# Epidemics - The SIR Model
- developed by *Kermack* and *McKendric* in the 1930s to describe the spread of malaria
- describes the spread of a disease from a population of infected individuals over time
- The total population is constant and is divided into three groups: *susceptible (S), infected (I) and recovered (R)*
- immune and deceased individuals are "lumped" together into group R
- only susceptible individuals can get the disease
- the incubation period is negligible - infected individuals are infectious right away

$$S \xrightarrow{r} I \xrightarrow{a} R$$

This gives rise to the following set of differential equations:

$$\begin{aligned}
	\frac{dS}{dt} = -rSI\\
	\frac{dI}{dt} = rSI - aI\\
	\frac{dR}{dt} = aI
\end{aligned}$$

S, I & R represent parts of a constant total population, hence the following two relations hold

$$S(t) + I(t) + R(t) = N$$

$$\frac{dS}{dt} + \frac{dI}{dt} + \frac{dR}{dt} = 0$$

> If the latter relationship does not immediately make sense to you, pause for a second to figure it out.

In order to make things a bit more comparable, we'll consider the total population to be normalised. Hence $N=1$, and values for $S$, $I$, and $R$ correspond to fractions of the population.

## Finding the right parameters
Once we have formulated a model - establishing a way how things depend on each other - we need to find the right parameters, and thereby establish how strong these links are.

The best way to do this is to look at real data. Let us try to parametrise the model to roughly mimic the dynamics of COVID-19. We know, that the base reproduction number, i.e. the number of individuals infected by a single infected person, is on the order of $3$, if no interventions (masks, lockdowns, etc.) are in place. This base reproduction number falls right out of our model. (In fact it was first defined based on a similar model). It is given by the following relationship:

$$R_0 = \frac{rS_0}{a}$$

If we have a rough idea of the fractions of the population that is either infected ($I_0$) or resistant ($R_0$), we can easily determine the fraction of susceptible individuals ($S_0$).

The parameters $r$ and $a$ are harder to determine. However, if we have an idea of the amount of time infected individuals remain infectious on average, we can estimate the recovery rate $a$ using the following relationship:

$$t_{infectious} = \frac{1}{a}$$

>  By thinking about what the equations in the model mean, convince yourself that the relationship between our parameter $a$ and the average infectious time holds
"""

# ╔═╡ 1fc5f952-91ff-11eb-3366-9f477c640915
md"""
The code below implements the SIR-model and plots the output. Parameters can be adjusted using the sliders above the graphs. The simulation, as well as values for $R_0$ and $t_{infectious}$ are updated accordingly.

> **Exercise 1 - Understanding the model:** Try to paramterise the model to capture the dynamics of COVID-19 without any interventions. Start by changing $a$ and $r$ and observe the effect on the overall dynamics (as shown in the timecourse plot), as well as on the base reproduction rate and the infectious period. - Use an average infectious period of 10 days and a base reproduction rate of approximately 3.5 to guide your parametrisation.
>
> - How strong is the effect of $r$ and $a$ on the base reproduction number?
> - What is the effect of a higher proportion of infected individuals at the start on the dynamics of the outbreak and $R_0$?
> - Assuming a starting population of infectious individuals of 1%, how long does it take until the outbreak peaks? Roughly how large is the fraction of the population that has recovered at the end?
> - Assuming a case mortality (percentage of people with the disease that die) of 2%, what would be your estimate of the total number of fatalities at the end of the outbreak?
> - Research online to try to come up with reasonable estimates for the fraction of people currently infected and the number of recovered people (including vaccinations) to be used as $I(0)$ and $R(0)$ respectively
"""

# ╔═╡ bfe81088-9203-11eb-162b-f373131bc105
md"""
> **Exercise 2 - Lowering $R_0$:** The current base reproduction rate in Austria is estimated to be slightly above 1. Change the transmission rate accordingly.
>
> - How has the change affected the dynamics of the outbreak (time to peak, fraction of population that has recovered, fraction of the population infected at peak, estimated number of fatalities)?
"""

# ╔═╡ e574a398-9203-11eb-0006-c38430911ee1
md"""
>**Exercise 3 - Intensive Care:** Roughly 1.5% of cases need to be treated in intensive care. Roughly 2500 beds in intensive care units are available, of which roughly 1000 beds can be made available to treat COVID-19 patients. This translates to approximately 0.0001 ICU beds per capita.
>
> - Assuming that the recovery rate of ICU patients is unchanged, what will be the demand for ICU beds at the peak of the outbreak with a base reproduction number of 1.17? How does this compare to the ICU bed capacity estimated above?
> - How low would $R_0$ need to be in order not to overwhelm ICU capacity?
> - Give five reasons why this scenario might be an oversimplification? Make an educated guess, as to whether you're likely to over- or underestimate the number of patients needing ICU-care.
"""

# ╔═╡ ac0c0230-91f0-11eb-057d-2103908afe51
md""" ## Simulation
*Change the sliders to adjust the corresponding values*:

| Parameter         | Value                  |
|-------------------|------------------------|
 transmission rate $r$ | $(@bind r Slider(0:0.005:1, default=0.125, show_value=true)) |
| recovery rate $a$ | $(@bind a Slider(0:0.001:0.25, default=0.1, show_value=true)) |
| initial number of infected $I(0)$ | $(@bind I0 Slider(0.0001:0.001:0.1, default=0.01, show_value=true)) |
| initial number of resistant $R(0)$ | $(@bind R0 Slider(0:0.001:0.1, default=0,show_value=true)) |
| timespan $t_{max}$ | $(@bind tmax Slider(10.0:10.0:2000.0, default=100.0,show_value=true)) days |
| zoom y-axis | $(@bind ymax Slider(0.05:0.05:1.0, default=1.0,show_value=true)) |
| Population size             | $(@bind popsize NumberField(500000:10000000000, default=9044650)) |
"""

# ╔═╡ a5f6c71c-920a-11eb-09e9-098f52c01f11
md"""
!!! danger "Beware!"

    The blocks below are used to specify the model. Feel free to play around with them, but proceed with care!
"""

# ╔═╡ 618d69dc-91f1-11eb-211f-d7548bf4a306
function sir!(du,u,p,t)
	du[1] = -p[1] * u[1] * u[2]
	du[2] = p[1] * u[1] * u[2] - p[2] * u[2]
	du[3] = p[2] * u[2]
end

# ╔═╡ 9a4c2f80-91f1-11eb-2b14-1784007ca92c
p = [r, a]

# ╔═╡ a6c393c8-91f1-11eb-0be0-99d47c2b207f
u0 = [1.0 - I0 - R0, I0, R0]

# ╔═╡ eccbab9e-91f1-11eb-19b2-d733609880d0
tspan = (0.0, 2000)

# ╔═╡ f2abe15a-91f1-11eb-17c2-e99b8b67502f
prob = ODEProblem(sir!,u0,tspan,p)

# ╔═╡ 4c7f9b82-91f1-11eb-0064-e70bb182dae0
begin
	sol = solve(prob, dt=1.0, dense = true)

	simplot = plot(sol.t, [sol[3,:] + sol[2,:] + sol[1,:] sol[3,:] + sol[2,:] sol[2,:]],
		fill = true,
		color = ["gray89" "paleturquoise3" "lightcoral"],
		label = ["S(t)" "R(t)" "I(t)"],
		xlabel = "t / days",
		title = "Timecourse",
		xlims = (0, tmax),
		ylims = (0, ymax))

	phaseplot1 = plot(sol, vars = (1,2),
		xlabel = "S(t)",
		ylabel = "I(t)",
		label = "",
		title = "Phaseplanes",
		xlims = (0,ymax),
		ylims = (0,ymax),
		color = "skyblue4")

	phaseplot2 = plot(sol, vars = (3,2),
		xlabel = "R(t)",
		ylabel = "I(t)",
		label = "",
		xlims = (0,ymax),
		ylims = (0,ymax),
		color = "skyblue4")

	phaseplot = plot(phaseplot1, phaseplot2, layout = (2,1))

	#plot(simplot, phaseplot, layout = (1, 2),fontfamily="Computer Modern")
	plot(simplot,fontfamily="Computer Modern")
end

# ╔═╡ 0a8e0895-2de8-435e-827c-fb216ab07f10
function getres(solution, i, t, p)
	v = solution(t)[i]*p
	return round(maximum([0.0 v]))
end

# ╔═╡ 0b66bc58-91fa-11eb-24f1-7f9f839a6cbf
md"""
Base reprodution rate = $(round((r*(1-I0-R0)/a), sigdigits = 3)),
Infectious period = $(round(1/a, sigdigits = 3))

| t(days)   | I(t)   | R(t)    |
|-----------|--------|---------|
| 0         | $(getres(sol, 2, 0, popsize)) | $(getres(sol, 3, 0, popsize)) |
| $tmax     | $(getres(sol, 2, tmax, popsize)) | $(getres(sol, 3, tmax, popsize)) |
"""

# ╔═╡ Cell order:
# ╟─7cd44a22-91f0-11eb-056d-699740ca7da9
# ╟─3199bb4e-91ed-11eb-37fa-abed99588e6b
# ╟─1fc5f952-91ff-11eb-3366-9f477c640915
# ╟─bfe81088-9203-11eb-162b-f373131bc105
# ╟─e574a398-9203-11eb-0006-c38430911ee1
# ╟─ac0c0230-91f0-11eb-057d-2103908afe51
# ╟─0b66bc58-91fa-11eb-24f1-7f9f839a6cbf
# ╟─4c7f9b82-91f1-11eb-0064-e70bb182dae0
# ╟─a5f6c71c-920a-11eb-09e9-098f52c01f11
# ╟─618d69dc-91f1-11eb-211f-d7548bf4a306
# ╟─9a4c2f80-91f1-11eb-2b14-1784007ca92c
# ╟─a6c393c8-91f1-11eb-0be0-99d47c2b207f
# ╟─eccbab9e-91f1-11eb-19b2-d733609880d0
# ╟─f2abe15a-91f1-11eb-17c2-e99b8b67502f
# ╟─0a8e0895-2de8-435e-827c-fb216ab07f10

max_stars_repo_name: el-uhu/intro-to-modelling

max_issues_repo_name: el-uhu/intro-to-modelling

max_forks_repo_name: el-uhu/intro-to-modelling


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import JSON3
import Tables

using RAI
using Test

@testset "api.jl" begin
    include("api.jl")
end

# def output =
#     1, "foo", 3.4, :foo;
#     2, "bar", 5.6, :foo
const body = """{
    "output": [{
        "rel_key": {
            "name": "output",
            "values": [],
            "keys": ["Int64", "String", "Float64", ":foo"],
            "type": "RelKey"
        },
        "type": "Relation",
        "columns": [[1, 2], ["foo", "bar"], [3.4, 5.6]]
    }],
    "problems": [],
    "actions": [{
        "name": "action1",
        "type": "LabeledActionResult",
        "result": {
            "output": [],
            "type": "QueryActionResult"}
    }],
    "debug_level": 0,
    "aborted": false,
    "type": "TransactionResult"
}
"""

@testset "JSON Response" begin
    rsp = JSON3.read(body)
    @test rsp.type == "TransactionResult"
    @test rsp.aborted == false
    @test length(rsp.problems) == 0
    @test length(rsp.output) == 1
    @test rsp.output[1].type == "Relation"
    @test rsp.output[1].rel_key.type == "RelKey"
    @test length(rsp.output[1].columns) == 3
    @test length(rsp.output[1].columns[1]) == 2
    @test rsp.output[1].columns[1] == [1, 2]
    @test rsp.output[1].columns[2] == ["foo", "bar"]
    @test rsp.output[1].columns[3] == [3.4, 5.6]
end

@testset "TransactionResult" begin
    rsp = JSON3.read(body)
    result = TransactionResult(rsp)
    @test result.type == "TransactionResult"
    @test result.aborted == false
    @test length(result.problems) == 0
    @test length(result.output) == 1
    @test rsp.output[1].type == "Relation"
    @test rsp.output[1].rel_key.type == "RelKey"
    @test length(result.output[1].columns) == 3
    @test length(result.output[1].columns[1]) == 2
    @test result.output[1].columns[1] == [1, 2]
    @test result.output[1].columns[2] == ["foo", "bar"]
    @test result.output[1].columns[3] == [3.4, 5.6]
end

@testset "TransactionResult.relations" begin
    rsp = JSON3.read(body)
    result = TransactionResult(rsp)
    @test length(result.relations) == 1
    @test length(result.relations[1]) == 2
    @test schema(result.relations[1]) == [Int64, String, Float64, Symbol]
    @test getrow(result.relations[1], 1) == (1, "foo", 3.4, :foo)
    @test getrow(result.relations[1], 2) == (2, "bar", 5.6, :foo)
    @test result.relations[1][1] == (1, "foo", 3.4, :foo)
    @test result.relations[1][2] == (2, "bar", 5.6, :foo)
end

# Instantiate `Relation` on the selected output.
@testset "Relation" begin
    rsp = JSON3.read(body)
    rel = Relation(rsp.output[1])
    @test length(rel) == 2  # two rows
    @test schema(rel) == [Int64, String, Float64, Symbol]
    @test getrow(rel, 1) == (1, "foo", 3.4, :foo)
    @test getrow(rel, 2) == (2, "bar", 5.6, :foo)
    @test rel[1] == (1, "foo", 3.4, :foo)
    @test rel[2] == (2, "bar", 5.6, :foo)
end

@testset "Tables.AbstractRow" begin
    rsp = JSON3.read(body)
    rel = Relation(rsp.output[1])
    @test Tables.istable(typeof(rel))
    @test Tables.rowaccess(typeof(rel))
    @test Tables.rows(rel) === rel
    row = first(rel)
    @test eltype(rel) == typeof(row)
    @test Tables.getcolumn(row, :Column1) == 1
    @test Tables.getcolumn(row, 1) == 1
    @test row.values == (1, "foo", 3.4, :foo)
    @test row.Column1 == 1
    @test row.Column2 == "foo"
    @test row.Column3 == 3.4
    @test row.Column4 == :foo
    @test row[1] == 1
    @test row[2] == "foo"
    @test row[3] == 3.4
    @test row[4] == :foo
    @test Tables.columnnames(rel) == [:Column1, :Column2, :Column3, :Column4]
    @test Tables.columnnames(row) == [:Column1, :Column2, :Column3, :Column4]
end

@testset "Tables.dictrowtable" begin
    rsp = JSON3.read(body)
    rel = Relation(rsp.output[1])
    table = Tables.dictrowtable(rel)
    row = first(table)
    @test row.Column1 == 1
    @test row.Column2 == "foo"
    @test row.Column3 == 3.4
    @test row.Column4 == :foo
end

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
module RetroCap

import Pkg
import RegistryTools
import UUIDs

export CapStrategy, UpperBound, MonotonicUpperBound
export LatestVersionOption, CapLatestVersion, ExcludeLatestVersion
export add_caps

include("types.jl")

include("add_caps.jl")
include("assert.jl")
include("compat_entries.jl")
include("parse_registry.jl")
include("utils.jl")

end # end module RetroCap


text:
using Test
using Example2018

# Polynomial form a ring
# 1 addition is associative****
# 2 additive identity
# 3 additive inverse
# 4 additive commutes****
# 5 multiplication is...

@testset "Polynomials form a ring" begin

	z = Polynomial([0])

	for n in 1:100

		deg = rand(0:100)

		p = Polynomial(rand(-1000000:1000000, deg+1))

		@test p + z == p
		@test z + p == p
		@test p + (-p) == z
	end
end


text:
# Autogenerated wrapper script for libnl_jll for armv7l-linux-musleabihf
export libnl, libnl_route

JLLWrappers.@generate_wrapper_header("libnl")
JLLWrappers.@declare_library_product(libnl, "libnl-3.so.200")
JLLWrappers.@declare_library_product(libnl_route, "libnl-route-3.so.200")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnl,
        "lib/libnl-3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnl_route,
        "lib/libnl-route-3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("funcs.jl")
using Printf

let

    cd("../../erickson/output_files/")
    done = false
    dir_names = []
    titles = []
    dir_count = 1
    while !done
        @printf "Directory %d (-1 is done): " dir_count
        dir = chomp(readline())

        if dir == "-1"
            done = true
        else
            @printf "Title: "
            push!(titles, chomp(readline()))
            push!(dir_names, dir)
            dir_count += 1
        end
    end

    nb_arr= []
    bi_arr = []
    for dir_name in dir_names
        slip_file = open(string(dir_name, "/slip.dat"))
        slip_data = eachline(slip_file)
        num_breaks, break_indices = get_break_indices(slip_data)
        push!(nb_arr, num_breaks)
        push!(bi_arr, break_indices)
        close(slip_file)
    end

    @printf "The minimum number of cycles over the simulations is %d\n" floor(minimum(nb_arr)/2)

    @printf "Offset contours by # of cycles (or -1 is no offset): "
    cycle_offset = parse(Int64, chomp(readline()))
    @printf "Last cycle: "
    final_cycle = parse(Int64, chomp(readline()))

    plts = []
    count = 1
    for (i, dir_name) in enumerate(dir_names)

        break_indices = bi_arr[i]
        total_cycles = floor(nb_arr[i])/2
        final_index, index_offset = get_plot_indices(cycle_offset,
                                                     final_cycle,
                                                     break_indices)

        slip_file = open(string(dir_name, "/slip.dat"))
        slip_data = collect(eachline(slip_file))

        y, nn, slip_data = get_y_nn(slip_data)

        slip_data = @view slip_data[index_offset:final_index,:]

        push!(plts, plot_slip(slip_data, y, nn, titles[i]))

        close(slip_file)

        @printf "\rGot plot %d" count
        count += 1
    end

    plot(plts..., layout=(length(plts),1))
    gui()

end


text:
using Combinatorics, IterTools

function create_partition(lengths::Vector{Int})
    num = 0
    partition = UnitRange{Int}[]
    for i = 1:length(lengths)
        push!(partition, (1:lengths[i]) .+ num)
        num += lengths[i]
    end
    return partition
end

function create_partition(lengths::Vector{Int}, names::NTuple{N,Symbol} where N)
    N = length(names)
    M = names
    if length(lengths) == N
        partition = create_partition(lengths)
        named_part = NamedTuple{names}(partition)
        return named_part::NamedTuple{M,NTuple{N,UnitRange{Int}}}
    else
        throw(ArgumentError("lengths must be the same"))
    end
end

function create_partition(sizes::Vector{NTuple{N,Int}} where N)
    num = 0
    partition = NTuple{N,UnitRange{Int}}[]
    for i = 1:length(sizes)
        push!(partition, (1:lengths[1]) .+ num)
        num += lengths[1]
    end
    return partition
end

function create_partition2(lengths::Vector{Int})
    part1 = create_partition(lengths)
    partition = NTuple{2,UnitRange{Int}}[]
    N = sum(lengths)
    for (rng1,rng2) in Iterators.product(part1,part1)
        push!(partition, (rng1,rng2))
    end
    # linds = LinearIndices((N,N))
    # return [linds[id[1],id[2]] for id in partition]
    # return CartesianIndices.(partition)
    return partition
end

function create_partition2(lengths::Vector{Int}, names::NTuple{N,Symbol} where N)
    if length(lengths) == length(names)
        partition = create_partition2(lengths)
        names_all = vec(collect(Iterators.product(names,names)))
        names_all = Tuple([Symbol(string(a) * string(b)) for (a,b) in names_all])
        named_part = NamedTuple{names_all}(partition)
        return named_part
    else
        throw(ArgumentError("lengths must be the same"))
    end
end


text:
typealias Syllable Tuple{Array{Int, 1}, Array{Int, 1}}
typealias SyllableBox Tuple{Tuple{Int, Int}, Tuple{Int, Int}}

@doc """
 # Type aliases
 * `Syllable` is an alias for `Tuple{Array{Int, 1}, Array{Int, 1}}`
 * `SyllableBox` is an alias for `Tuple{Tuple{Int, Int}, Tuple{Int, Int}}`
""" aliased_types

const aliased_types = 0 #Never used -- this is a hack to make Documenter.jl work


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("ahmi/ahmi.jl")
include("ahmi_integration.jl")
include("bridge_sampling_integration.jl")


text:
using InstantiateFromURL
# optionally add arguments to force installation: instantiate = true, precompile = true
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.8.0")

# Exercise 1

# Part 1: compute inner product using zip()

A = [1 2 3]
B = [4 5 6]
@show sum(a * b for (a,b) ∈ zip(A,B))

# part 2: count even numbers between 0,99 using comprehension
@show even = sum(iseven, 1:99)

# part 3: count number of pairs where both elements are even using comprehension
A = [2 4 9 12]
B = [5 2 8 10]
@show sum( iseven(a)&&iseven(b) for (a,b) ∈ zip(A,B))

# or
pairs = ((2, 5), (4, 2), (9, 8), (12, 10))
@show sum(xy -> all(iseven, xy), pairs);

# Exercise 2

function p(x,coef)
    values = zeros(0)
    for (i, α) in enumerate(coef)
        val = α * x^(i-1)
        push!(values,val)
    end
    return sum(values)
end

coef = [12 0.5 3]
x = 2

@show sol = p(x,coef)

# or

p(x,coef) = sum( α * x^(i-1) for (i,α) in enumerate(coef))

@show sol2 = p(x,coef);

# Exercise 3

uppercase("foo")

function countuppercase(string)
    count = 0
    for letter in string
        if letter == uppercase(letter) && isletter(letter)
            count += 1
        end
    end
    return count
end

@show sol = countuppercase("Romulus and Remus founded Rome")

# or

countuppercase(string) = sum( letter == uppercase(letter) && isletter(letter) for letter in string)

@show sol2 = countuppercase("Romulus and Remus founded Rome");

# Exercise 4

sameset(seq_a,seq_b) = all( a ∈ seq_b for a in seq_a )

A = ["apple" , 4, 2]
B = [2, "apple", 4]

sameset(A,B)

# Exercise 5 piecewise linear interpolation

function linapprox(f,a,b,n,x)

    interval_length = b-a # total length of interval

    no_subintervals = n-1 # number of subintervals = gridpoints - 1

    subinterval_length = interval_length/no_subintervals

    step = subinterval_length

    # find grid point larger than x
    point = a
    while point < x
        point+= step
    end

    # x must lie between the gridpoints (point - step) and point
    u, v = point - step, point

    # computer interpolation polynomial
    return f(u) + (x-u) * (f(v)-f(u)) / (v-u)
end

f(x) = exp(x) + x^2
f_approx(x) = linapprox(f, -1, 1, 3, x)

using Plots

x_grid = range(-5.0, 5.0, length = 400)
y_vals = f.(x_grid)
y = f_approx.(x_grid)
plot(x_grid, y_vals, label = "true")
plot!(x_grid, y, label = "approximation")

# Exercise 6

open("us_cities.txt", "w") do f
  write(f,
"new york: 8244910
los angeles: 3819702
chicago: 2707120
houston: 2145146
philadelphia: 1536471
phoenix: 1469471
san antonio: 1359758
san diego: 1326179
dallas: 1223229")
end

pop = open("us_cities.txt", "r")

total_pop =0
for line in eachline(pop)
    city, population = split(line, ":")
    total_pop += parse(Int, population)
end

close(pop)
println("total population = $total_pop")

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
module BenchDictMigration

using BenchmarkTools
using ConcurrentCollections
using ConcurrentCollections.Implementations:
    LINEAR_PROBING_DICT_EXPAND_BASESIZE, migrate_serial!

pad16(x) = string(x; pad = 16)

function generate(f = pad16; datasize = LINEAR_PROBING_DICT_EXPAND_BASESIZE[])
    vs = UInt64.(1:datasize)
    ks = f.(vs)
    dict = ConcurrentDict{eltype(ks),eltype(vs)}(zip(ks, vs))
    return dict
end

const CACHE = Ref{Any}()

function setup(; generate_options...)
    CACHE[] = Dict(
        "String-UInt64" => generate(; generate_options...),
        "Uint32-UInt64" => generate(UInt32; generate_options...),
        "UInt64-UInt64" => generate(UInt64; generate_options...),
        # "Uint32-Nothing" => ???,
    )

    suite = BenchmarkGroup()
    for key in keys(CACHE[])
        CacheType = typeof(CACHE[][key])
        suite[key] = @benchmarkable(
            migrate_serial!(newslots, slots),
            setup = begin
                dict = CACHE[][$key]::$CacheType
                slots = copy(dict.slots)
                newslots = similar(slots, length(slots) * 2)
            end,
            evals = 1,
        )
    end
    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
facts("The ex1 example") do
    m = Model()
    @variable(m, x)
    @variable(m, y)
    @constraints(m, begin
        12 + 2x - y ≥ 0
        -6 - x + 2y ≥ 0
        -3 + x + y ≥ 0
        1 + x ≥ 0
    end)

    v, r = get_extrema(m)

    @fact length(r) --> 2
    @fact length(v) --> 3

    for ex in Vector[[ 0, 3],
        [-1, 4],
        [-1,10]]
        @fact is_approx_included(v, ex) --> true
    end

    for ex in Vector[[1,2],
        [2,1]]
        @fact is_approx_included(r, ex) --> true
    end
end