text:
using SteadyState
using Test
using QuantumOptics, SparseArrays, LinearAlgebra
using IterativeSolvers

@testset "SteadyState.jl" begin
    ωc = 1.2
    ωa = 0.9
    g = 1.0
    γ = 0.5
    κ = 1.1

    T = Float64[0.,1.]

    fockbasis = FockBasis(10)
    spinbasis = SpinBasis(1//2)
    basis = tensor(spinbasis, fockbasis)

    sx = sigmax(spinbasis)
    sy = sigmay(spinbasis)
    sz = sigmaz(spinbasis)
    sp = sigmap(spinbasis)
    sm = sigmam(spinbasis)

    Ha = embed(basis, 1, 0.5*ωa*sz)
    Hc = embed(basis, 2, ωc*number(fockbasis))
    Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
    H = Ha + Hc + Hint

    Ja_unscaled = embed(basis, 1, sm)
    Jc_unscaled = embed(basis, 2, destroy(fockbasis))
    Junscaled = [Ja_unscaled, Jc_unscaled]

    Ja = embed(basis, 1, sqrt(γ)*sm)
    Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
    J = [Ja, Jc]

    Hdense = dense(H)
    Junscaled_dense = map(dense, Junscaled)
    Jdense = map(dense, J)

    Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 5)
    ρ₀ = dm(Ψ₀)

    tout, ρt = timeevolution.master([0,100], ρ₀, Hdense, Jdense; reltol=1e-7)

    #
    # Test defaults
    #

    ρ1 = SteadyState.iterative(ρ₀, Hdense, Jdense; tol=1e-7)
    @test tracedistance(ρ1, ρt[end]) < 1e-6
    ρ1_bis = SteadyState.iterative(Hdense, Jdense; tol=1e-7)
    @test tracedistance(ρ1, ρ1_bis) < 1e-6

    ρ2 = SteadyState.iterative(ρ₀, H, J; tol=1e-7)
    @test tracedistance(DenseOperator(ρ2), ρt[end]) < 1e-6
    ρ2_bis = SteadyState.iterative(H, J; tol=1e-7)
    @test tracedistance(DenseOperator(ρ2), DenseOperator(ρ2_bis)) < 1e-6

    ρ3 = SteadyState.iterative(Ψ₀, Hdense, Jdense; tol=1e-7)
    @test tracedistance(ρ3, ρt[end]) < 1e-6
    ρ3_bis = SteadyState.iterative(Hdense, Jdense; tol=1e-7)
    @test tracedistance(ρ3, ρ3_bis) < 1e-6

    ρ4 = SteadyState.iterative(Ψ₀, H, J; tol=1e-7)
    @test tracedistance(DenseOperator(ρ4), ρt[end]) < 1e-6
    ρ4_bis = SteadyState.iterative(H, J; tol=1e-7)
    @test tracedistance(DenseOperator(ρ4), DenseOperator(ρ4_bis)) < 1e-6

    # Test explicit call to iterative solvers
    ρ1, ch = SteadyState.iterative(ρ₀, Hdense, Jdense, bicgstabl!, 4; log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(ρ1, ρt[end]) < 1e-6
    ρ1_bis, ch = SteadyState.iterative(Hdense, Jdense, bicgstabl!, 4; log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(ρ1, ρ1_bis) < 1e-6

    ρ2, ch = SteadyState.iterative(ρ₀, H, J, idrs!; s=15, log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(DenseOperator(ρ2), ρt[end]) < 1e-6
    ρ2_bis, ch = SteadyState.iterative(H, J, idrs!; s=15, log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(DenseOperator(ρ2), DenseOperator(ρ2_bis)) < 1e-6

    ρ3, ch = SteadyState.iterative(Ψ₀, Hdense, Jdense, idrs!; s=15, log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(ρ3, ρt[end]) < 1e-6
    ρ3_bis, ch = SteadyState.iterative(Hdense, Jdense, idrs!; s=15, log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(ρ3, ρ3_bis) < 1e-6

    ρ4, ch = SteadyState.iterative(Ψ₀, H, J, idrs!; s=15, log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(DenseOperator(ρ4), ρt[end]) < 1e-6
    ρ4_bis, ch = SteadyState.iterative(H, J, idrs!; s=15, log=true, tol=1e-7)
    @test ch.isconverged
    @test tracedistance(DenseOperator(ρ4), DenseOperator(ρ4_bis)) < 1e-6

    #
    # Mirror to QuantumOptics.timeevolution.master tests
    #

    ρ = SteadyState.iterative(ρ₀, Hdense, Jdense; tol=1e-7)

    ρ2 = SteadyState.iterative(ρ₀, H, J; tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5

    ρ2 = SteadyState.iterative(Ψ₀, Hdense, Jdense; tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀, H, Jdense; tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5

    ρ2 = SteadyState.iterative(Ψ₀, Hdense, J; tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5

    # Test explicit gamma vector
    rates_vector = [γ, κ]

    ρ2 = SteadyState.iterative(ρ₀, Hdense, Junscaled_dense; rates=rates_vector, tol=1e-7)
    @test tracedistance(ρ, ρ2) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀, H, Junscaled_dense; rates=rates_vector, tol=1e-7)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀, H, Junscaled; rates=rates_vector, tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5

    # Test explicit gamma matrix
    alpha = 0.3
    R = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)]
    Rt = transpose(R)
    Jrotated_dense = [R[1,1]*Junscaled_dense[1] + R[1,2]*Junscaled_dense[2], R[2,1]*Junscaled_dense[1] + R[2,2]*Junscaled_dense[2]]
    Jrotated = [SparseOperator(j) for j=Jrotated_dense]
    rates_matrix = diagm(0 => rates_vector)
    rates_matrix_rotated = R * rates_matrix * Rt

    ρ2 = SteadyState.iterative(ρ₀, Hdense, Jrotated_dense; rates=rates_matrix_rotated, tol=1e-7)
    @test tracedistance(ρ, ρ2) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀, H, Jrotated; rates=rates_matrix_rotated, tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀, Hdense, Jrotated_dense; rates=rates_matrix_rotated, tol=1e-7)
    @test tracedistance(ρ, ρ2) < 1e-5

    #
    # Test genericity
    #

    f(O::AbstractOperator) = Complex{BigFloat}.(O.data)

    ρ2 = SteadyState.iterative(f(Hdense), f.(Jdense); tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,Matrix(ρ2))) < 1e-5

    ρ2 = SteadyState.iterative(f(Hdense), f.(Junscaled_dense); Jdagger=f.(dagger.(Junscaled_dense)), rates=rates_vector, tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,ρ2)) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀.data, f(Hdense), f.(Jdense); tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,Matrix(ρ2))) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀.data, f(Hdense), f.(Junscaled_dense); Jdagger=f.(dagger.(Junscaled_dense)), rates=rates_vector, tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,ρ2)) < 1e-5

    ρ2 = SteadyState.iterative(f(H), f.(J); tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,Matrix(ρ2))) < 1e-5

    ρ2 = SteadyState.iterative(f(H), f.(Junscaled); Jdagger=f.(dagger.(Junscaled)), rates=rates_vector, tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,Matrix(ρ2))) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀.data, f(H), f.(J); tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,Matrix(ρ2))) < 1e-5

    ρ2 = SteadyState.iterative(ρ₀.data, f(H), f.(Junscaled); Jdagger=f.(dagger.(Junscaled)), rates=rates_vector, tol=1e-7)
    @test tracedistance(ρ, DenseOperator(basis,Matrix(ρ2))) < 1e-5

    #
    # Test residual computation for any possible signature
    #

    # no BLAS
    # rates=nothing, Jdagger=nothing
    ρ2 = SteadyState.iterative(ρ₀, H, J; tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5
    # rates=Vector, Jdagger=nothing
    ρ2 = SteadyState.iterative(ρ₀, H, Junscaled; rates=rates_vector, tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5
    # rates=Matrix, Jdagger=nothing
    ρ2 = SteadyState.iterative(ρ₀, H, Jrotated; rates=rates_matrix_rotated, tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5
    # rates=nothing, Jdagger=dagger.(J)
    ρ2 = SteadyState.iterative(ρ₀, H, J; Jdagger=dagger.(J), tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5
    # rates=Vector, Jdagger=dagger.(J)
    ρ2 = SteadyState.iterative(ρ₀, H, Junscaled; Jdagger=dagger.(Junscaled), rates=rates_vector, tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5
    # rates=Matrix, Jdagger=dagger.(J)
    ρ2 = SteadyState.iterative(ρ₀, H, Jrotated; Jdagger=dagger.(Jrotated), rates=rates_matrix_rotated, tol=1e-6)
    @test tracedistance(ρ, DenseOperator(ρ2)) < 1e-5

    # BLAS
    # rates=nothing, Jdagger=nothing
    ρ2 = SteadyState.iterative(DenseOperator(ρ₀), Hdense, Jdense; tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5
    # rates=Vector, Jdagger=nothing
    ρ2 = SteadyState.iterative(DenseOperator(ρ₀), Hdense, Junscaled_dense; rates=rates_vector, tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5
    # rates=Matrix, Jdagger=nothing
    ρ2 = SteadyState.iterative(DenseOperator(ρ₀), Hdense, Jrotated_dense; rates=rates_matrix_rotated, tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5
    # rates=nothing, Jdagger=dagger.(J)
    ρ2 = SteadyState.iterative(DenseOperator(ρ₀), Hdense, Jdense; Jdagger=dagger.(Jdense), tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5
    # rates=Vector, Jdagger=dagger.(J)
    ρ2 = SteadyState.iterative(DenseOperator(ρ₀), Hdense, Junscaled_dense; Jdagger=dagger.(Junscaled_dense), rates=rates_vector, tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5
    # rates=Matrix, Jdagger=dagger.(J)
    ρ2 = SteadyState.iterative(DenseOperator(ρ₀), Hdense, Jrotated_dense; Jdagger=dagger.(Jrotated_dense), rates=rates_matrix_rotated, tol=1e-6)
    @test tracedistance(ρ, ρ2) < 1e-5
end


text:
export fdtdAdj, fdtdAdj!

#xr and xs, no IC
function fdtdAdj{T,N}(xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		      f::FDTD{T,N},
		      p::Array{T})
	s = fdtd(xs,xr,Nt,f,flipdim(p,1))
	return flipdim(s,1)
end

function fdtdAdj!{T}(s_out::Array{T},
		     xr::Vector{Tuple{Int,Int,Int}}, xs::Vector{Tuple{Int,Int,Int}},Nt::Int,
		     f::FDTD{T,3}, p::Array{T})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = Nt:-1:1
		fdtdKernel!(n,s_out,p0,p1,p2,xs,xr,p,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return s_out
end

#full and xs, no IC
function fdtdAdj{T,N}(xs::Vector{Tuple{Int,Int,Int}}, Nt::Int,
		      f::FDTD{T,N},
		      p::Array{T,4})
	s = fdtd(xs,Nt,f,flipdim(p,4))
	return flipdim(s,1)
end

function fdtdAdj!{T}(s_out::Array{T},
		     xs::Vector{Tuple{Int,Int,Int}},Nt::Int,
		     f::FDTD{T,3}, p::Array{T,4})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = Nt:-1:1
		fdtdKernel!(n,s_out,p0,p1,p2,xs,p,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return s_out
end

#full and source everywhere, no IC
function fdtdAdj{T,N}(Nt::Int,
		      f::FDTD{T,N},
		      p::Array{T,4})
	s = fdtd(Nt,f,flipdim(p,4))
	return flipdim(s,4)
end

function fdtdAdj!{T}(s_out::Array{T,4},
		     Nt::Int,
		     f::FDTD{T,3}, p::Array{T,4})
	resetIC!(f)
	p2 = f.p2
	p1 = f.p1
	p0 = f.p0
	for n = Nt:-1:1
		fdtdKernel!(n,s_out,p0,p1,p2,p,f.d[1],f.d[2],f.d[3],f.d[4],f.q,f.G)
		p2, p1, p0 = p0, p2, p1
	end
	return s_out
end

ext: jl

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


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

"""
    SimulationSolution

A solution to a spatial simulation problem.
"""
struct SimulationSolution{D<:AbstractDomain}
  domain::D
  realizations::Dict{Symbol,Vector{<:AbstractVector}}
end

SimulationSolution(domain, realizations) =
  SimulationSolution{typeof(domain)}(domain, realizations)

"""
    getindex(solution, var)

Return simulation solution for specific variable `var`
as a vector of realizations.
"""
function Base.getindex(solution::SimulationSolution, var::Symbol)
  solution.realizations[var]
end

function Base.getindex(solution::SimulationSolution{<:RegularGrid}, var::Symbol)
  sz = size(solution.domain)
  [reshape(real, sz) for real in solution.realizations[var]]
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, solution::SimulationSolution)
  dim = ndims(solution.domain)
  print(io, "$(dim)D SimulationSolution")
end

function Base.show(io::IO, ::MIME"text/plain", solution::SimulationSolution)
  println(io, solution)
  println(io, "  domain: ", solution.domain)
  print(  io, "  variables: ", join(keys(solution.realizations), ", ", " and "))
end


text:
module ImgProcMic

using ColorTypes

export complex_show, gray_show

"""
    complex_show(arr)

Displays a complex array. Color encodes phase, brightness encodes magnitude.
Works within Jupyter and Pluto.
"""
function complex_show(cpx::AbstractArray{<:Complex, N}) where N
	ac = abs.(cpx)
	HSV.(angle.(cpx)./2pi*256,ones(Float32,size(cpx)),ac./maximum(ac))
end

"""
    gray_show(arr; set_one=false, set_zero=false)
Displays a real gray color array. Brightness encodes magnitude.
Works within Jupyter and Pluto.

## Keyword args
* `set_one=false` divides by the maximum to set maximum to 1
* `set_zero=false` subtracts the minimum to set minimum to 1
"""
function gray_show(arr; set_one=true, set_zero=false)
    arr = set_zero ? arr .- minimum(arr) : arr
    arr = set_one ? arr ./ maximum(arr) : arr
    Gray.(arr)
end

end


text:
using POMDPTools
using Test

@testset "POMDPTools.jl" begin
    # Write your tests here.
end


text:
using Plots
using BeliefPropagation
using Erdos
using Random, Statistics

# 2D lattice

L = 8
net = Grid([L, L], Network, periodic=true)
# Set J and H has graph properties since they are constants
gprop!(net, "J", 1)
gprop!(net, "H", 0)
writenetwork("ferromagnet_2d_L$L.graphml", net)

Ts =  1:0.1:4
Ms = Float64[]
for T in Ts
    fg = Ising.run_bp(net, T=T, maxiters=1000, μ=1, verbose=false);
    m = abs(mean(fg.mags))
    push!(Ms, m)
end
Tc = 2 / log(1+√2) #2.269185314213022 from Onsager's solution
plot(Ts, Ms, label="BP")
title!("2D ferromagnet")
xlabel!("T")
ylabel!("M")
vline!([Tc], label="Tc Onsager", color="black")


text:
RealVector{T} = AbstractArray{T,1} where T<:Real
RealMatrix{T} = AbstractArray{T,2} where T<:Real
RVVector{T} = AbstractArray{T,1} where T<:RealVector
RVVVector{T} = AbstractArray{T,1} where T<:RVVector

"Time Unit set to millisecond by default."
SecondPerUnit = 0.001
function settimeunit(secondperunit)
    global SecondPerUnit = secondperunit
end
function timetounit(t,tunit=:ms)
    if tunit==:ms
        ts = t*0.001
    else
        ts = t
    end
    ts/SecondPerUnit
end


text:
module PSCOPFio

using ..PSCOPF
using ..Networks

using Dates
using Printf
using DataStructures

OUTPUT_PREFIX = "pscopf_out_"

##########################
#   Readers
##########################

function read_buses!(network::Network, data::String)
    buses_ids = Set{String}();
    open(joinpath(data, "pscopf_ptdf.txt"), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);
                push!(buses_ids, buffer[2])
            end
        end
    end

    open(joinpath(data, "pscopf_units.txt"), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);
                push!(buses_ids, buffer[3])
            end
        end
    end

    Networks.add_new_buses!(network, collect(buses_ids));
end

function read_branches!(network::Network, data::String)
    branches = Dict{String,Float64}();
    default_limit = 0.
    open(joinpath(data, "pscopf_ptdf.txt"), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);
                branch_id = buffer[1]
                push!(branches, branch_id => default_limit)
            end
        end
    end

    open(joinpath(data, "pscopf_limits.txt"), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);
                branch_id = buffer[1]
                limit = parse(Float64, buffer[2]);
                push!(branches, branch_id=>limit)
            end
        end
    end

    for (id,limit) in branches
        Networks.add_new_branch!(network, id, limit);
    end
end

function read_ptdf!(network::Network, data::String, filename="pscopf_ptdf.txt")
    open(joinpath(data, filename), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);
                branch_id = buffer[1]
                bus_id = buffer[2]
                ptdf_value = parse(Float64, buffer[3])
                Networks.add_ptdf_elt!(network, branch_id, bus_id, ptdf_value)
            end
        end
    end
end

function read_generators!(network, data)
    open(joinpath(data, "pscopf_units.txt"), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);

                generator_id = buffer[1]
                gen_type = parse(Networks.GeneratorType, buffer[2])
                gen_bus_id = buffer[3]
                pmin = parse(Float64, buffer[4])
                pmax = parse(Float64, buffer[5])
                start_cost = parse(Float64, buffer[6])
                prop_cost = parse(Float64, buffer[7])
                dmo = Dates.Second(parse(Float64, buffer[8]))
                dp = Dates.Second(parse(Float64, buffer[9]))

                Networks.add_new_generator_to_bus!(network, gen_bus_id,
                                        generator_id, gen_type, pmin, pmax, start_cost, prop_cost, dmo, dp)
            end
        end
    end
end

function read_uncertainties_distributions(data)
    result = PSCOPF.UncertaintiesDistribution()
    open(joinpath(data, "uncertainties_distribution.txt"), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln);

                id = buffer[1]
                min_value = parse(Float64, buffer[2])
                max_value = parse(Float64, buffer[3])
                mu = parse(Float64, buffer[4])
                sigma = parse(Float64, buffer[5])
                time_factor = parse(Float64, buffer[6])

                PSCOPF.add_uncertainty_distribution!(result,
                                        id, min_value, max_value, mu, sigma,
                                        time_factor)
            end
        end
    end

    return result
end

function read_uncertainties(data, filename="pscopf_uncertainties.txt")
    result = PSCOPF.Uncertainties()
    open(joinpath(data, filename), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln)
                # "name", "ts", "ech", "scenario", "value"))
                name = buffer[1]
                ts = Dates.DateTime(buffer[2])
                ech = Dates.DateTime(buffer[3])
                scenario = buffer[4]
                value = parse(Float64, buffer[5])
                PSCOPF.add_uncertainty!(result, ech, name, ts, scenario, value)
            end
        end
    end
    return result
end

function read_initial_state(data, filename="pscopf_init.txt")
    result = SortedDict{String, PSCOPF.GeneratorState}()
    open(joinpath(data, filename), "r") do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = PSCOPF.split_with_space(ln)
                gen_id = buffer[1]
                state = parse(PSCOPF.GeneratorState, buffer[2])
                result[gen_id] = state
            end
        end
    end
    return result
end

##########################
#   Writers
##########################

function write(dir_path::String, network::Networks.Network)
    mkpath(dir_path)
    # if !isdir(dir_path)
    #     mkpath(dir_path)
    # else
    #     msg = @sprintf("data folder `%s` already exists!", dir_path)
    #     error(msg)
    # end
    #units
    write(dir_path, network.generators)
    #limits
    write(dir_path, network.branches)
    #ptdf
    write(dir_path, network.ptdf)
end

function write(dir_path::String, generators::SortedDict{String, Networks.Generator})
    output_file_l = joinpath(dir_path, "pscopf_units.txt")
    open(output_file_l, "w") do file_l
        Base.write(file_l, @sprintf("#%24s%16s%25s%16s%16s%16s%16s%16s%16s\n",
                    "name", "type", "bus_id", "minP","maxP", "start", "prop", "dmo(s)", "dp(s)"))
        for (id_l, generator_l) in generators
            Base.write(file_l, @sprintf("%25s%16s%25s%16.8E%16.8E%16.8E%16.8E%16.8E%16.8E\n",
                                    Networks.get_id(generator_l),
                                    Networks.get_type(generator_l),
                                    Networks.get_bus_id(generator_l),
                                    Networks.get_p_min(generator_l),
                                    Networks.get_p_max(generator_l),
                                    Networks.get_start_cost(generator_l),
                                    Networks.get_prop_cost(generator_l),
                                    Dates.value(Networks.get_dmo(generator_l)),
                                    Dates.value(Networks.get_dp(generator_l))
                                    )
                    )
        end
    end
end

function write(dir_path::String, branches::SortedDict{String, Networks.Branch})
    output_file_l = joinpath(dir_path, "pscopf_limits.txt")
    open(output_file_l, "w") do file_l
        Base.write(file_l, @sprintf("#%24s%16s\n", "branch", "limit"))
        for (id_l, branch_l) in branches
            Base.write(file_l, @sprintf("%25s%16.8E\n",
                                    Networks.get_id(branch_l),
                                    Networks.get_limit(branch_l)
                                    )
                    )
        end
    end
end

function write(dir_path::String, ptdf::SortedDict{String,SortedDict{String, Float64}})
    output_file_l = joinpath(dir_path, "pscopf_ptdf.txt")
    open(output_file_l, "w") do file_l
        Base.write(file_l, @sprintf("#%24s%16s\n", "REF_BUS", "unknown"))
        Base.write(file_l, @sprintf("#%24s%25s%16s\n", "branch", "bus", "value"))
        for (branch_id_l, _) in ptdf
            for (bus_id_l, val_l) in ptdf[branch_id_l]
                Base.write(file_l, @sprintf("%25s%25s%16.8E\n",
                                        branch_id_l,
                                        bus_id_l,
                                        val_l
                                        )
                            )
            end
        end
    end
end

function write(dir_path::String, gen_init::SortedDict{String, PSCOPF.GeneratorState})
    mkpath(dir_path)
    output_file_l = joinpath(dir_path, "pscopf_init.txt")
    open(output_file_l, "w") do file_l
        Base.write(file_l, @sprintf("#%24s%10s\n", "name", "state"))
        for (gen_id, state) in gen_init
            Base.write(file_l, @sprintf("%25s%10s\n", gen_id, state))
        end
    end
end

function write(dir_path::String, uncertainties::PSCOPF.Uncertainties)
    mkpath(dir_path)
    output_file_l = joinpath(dir_path, "pscopf_uncertainties.txt")
    open(output_file_l, "w") do file_l

        Base.write(file_l, @sprintf("#%24s%20s%20s%10s%16s\n", "name", "ts", "ech", "scenario", "value"))
        for (ech, _) in uncertainties
            for (nodal_injection_name, _) in uncertainties[ech]
                for (ts, _) in uncertainties[ech][nodal_injection_name]
                    for (scenario, value_l) in uncertainties[ech][nodal_injection_name][ts]
                        Base.write(file_l, @sprintf("%25s%20s%20s%10s%16.8E\n",
                                        nodal_injection_name,
                                        ts,
                                        ech,
                                        scenario,
                                        value_l
                                        )
                            )
                    end
                end
            end
        end

    end
end

function write_commitment_schedule(dir_path::String, schedule::PSCOPF.Schedule, prefix="")
    ech = schedule.decision_time
    commitment_filename_l = joinpath(dir_path, OUTPUT_PREFIX*prefix*"commitment_schedule.txt")
    open(commitment_filename_l, "a") do commitment_file_l
        if filesize(commitment_file_l) == 0
            Base.write(commitment_file_l, @sprintf("#%19s%10s%25s%20s%10s%6s%10s\n", "ech", "decider", "name", "ts", "scenario", "value", "firmness"))
        end
        for (gen_id, gen_schedule) in schedule.generator_schedules
            for (ts, uncertain_value) in gen_schedule.commitment
                firmness = PSCOPF.is_definitive(uncertain_value) ? "FIRM" : "FREE"
                for (scenario, value_l) in uncertain_value.anticipated_value
                    value_l = ismissing(value_l) ? -1. : value_l
                    Base.write(commitment_file_l, @sprintf("%20s%10s%25s%20s%10s%6s%10s\n",
                                        ech,
                                        schedule.decider_type,
                                        gen_id,
                                        ts,
                                        scenario,
                                        value_l,
                                        firmness
                                        )
                            )
                end
            end
        end
    end
end

function write_production_schedule(dir_path::String, schedule::PSCOPF.Schedule, prefix="")
    ech = schedule.decision_time
    schedule_filename_l = joinpath(dir_path, OUTPUT_PREFIX*prefix*"schedule.txt")
    open(schedule_filename_l, "a") do schedule_file_l
        if filesize(schedule_file_l) == 0
            Base.write(schedule_file_l, @sprintf("#%19s%10s%25s%20s%10s%16s%10s\n", "ech", "decider", "name", "ts", "scenario", "value", "firmness"))
        end
        for (gen_id, gen_schedule) in schedule.generator_schedules
            for (ts, uncertain_value) in gen_schedule.production
                firmness = PSCOPF.is_definitive(uncertain_value) ? "FIRM" : "FREE"
                for (scenario, value_l) in uncertain_value.anticipated_value
                    value_l = ismissing(value_l) ? -1. : value_l
                    Base.write(schedule_file_l, @sprintf("%20s%10s%25s%20s%10s%16.8E%10s\n",
                                        ech,
                                        schedule.decider_type,
                                        gen_id,
                                        ts,
                                        scenario,
                                        value_l,
                                        firmness
                                        )
                            )
                end
            end
        end
    end
end

function write_flows(dir_path::String, context::PSCOPF.AbstractContext, schedule::PSCOPF.Schedule, prefix="")
    ech = schedule.decision_time
    flows_filename_l = joinpath(dir_path, OUTPUT_PREFIX*prefix*"flows.txt")

    flows = PSCOPF.compute_flows(context, schedule)
    open(flows_filename_l, "a") do flows_file_l
        if filesize(flows_file_l) == 0
            Base.write(flows_file_l, @sprintf("#%19s%10s%25s%20s%10s%16s\n", "ech", "decider", "branch_name", "ts", "scenario", "value"))
        end
        for ((branch_id, ts, scenario), flow_value) in flows
            Base.write(flows_file_l, @sprintf("%20s%10s%25s%20s%10s%16.8E\n",
                                ech,
                                schedule.decider_type,
                                branch_id,
                                ts,
                                scenario,
                                flow_value
                                )
                    )
        end
    end
end

function write(context::PSCOPF.AbstractContext, schedule::PSCOPF.Schedule, prefix="")
    dir_path = context.out_dir
    if !isnothing(dir_path)
        mkpath(dir_path)
        write_commitment_schedule(dir_path, schedule, prefix)
        write_production_schedule(dir_path, schedule, prefix)
        write_flows(dir_path, context, schedule, prefix)
    end
end

function write_full(context::PSCOPF.AbstractContext, schedule::PSCOPF.Schedule, prefix="")
    dir_path = context.out_dir
    if !isnothing(dir_path)
        mkpath(dir_path)
        _write(dir_path, schedule, context, prefix)
    end
end

function _write(dir_path::String, schedule::PSCOPF.Schedule, context::PSCOPF.AbstractContext, prefix="")
    ech = schedule.decision_time
    decider = schedule.decider_type
    schedule_filename_l = joinpath(dir_path, OUTPUT_PREFIX*prefix*"full_schedule.txt")
    open(schedule_filename_l, "a") do schedule_file_l
        if filesize(schedule_file_l) == 0
            Base.write(schedule_file_l, @sprintf("#%19s%10s%25s%10s%20s%10s%16s%16s%16s%16s%15s%8s%15s\n",
                                            "ech", "decider", "unit/bus", "PROD/LOAD", "ts", "scenario",
                                            "inject_value", "power_value", "load_shed", "power_capped", "DP_firmness",
                                            "ON/OFF", "DMO_firmness",
                                            ))
        end

        TS =  PSCOPF.get_target_timepoints(context)
        scenarios = PSCOPF.get_scenarios(context)
        uncertainties_at_ech = PSCOPF.get_uncertainties(context, ech)
        network = PSCOPF.get_network(context)

        for generator in Networks.get_generators(network)
            gen_id = Networks.get_id(generator)
            gen_schedule = PSCOPF.get_sub_schedule(schedule, gen_id)

            load_shed = 0. #this is for buses only
            for ts in TS
                power_uncertain_value = PSCOPF.get_prod_uncertain_value(gen_schedule, ts)
                power_firmness = PSCOPF.is_definitive(power_uncertain_value) ? "FIRM" : "FREE"
                if Networks.needs_commitment(generator)
                    commitment_uncertain_value = PSCOPF.get_commitment_uncertain_value(gen_schedule, ts)
                    commitment_firmness = PSCOPF.is_definitive(commitment_uncertain_value) ? "FIRM" : "FREE"
                else
                    commitment_firmness = "X"
                end
                for scenario in scenarios
                    prod_value = PSCOPF.get_prod_value(gen_schedule, ts, scenario)
                    if Networks.is_limitable(generator)
                        power_value = PSCOPF.get_uncertainties(uncertainties_at_ech, gen_id, ts, scenario)
                        power_capped = PSCOPF.get_capping(schedule, gen_id, ts, scenario)
                        inject_value = prod_value - power_capped #TODO : add this to schedule and compute it by each Model
                                                                 #FIXME : for TSOModel the prod_value take account of capping already
                    else
                        inject_value = power_value = prod_value
                        power_capped = 0.
                    end
                    if Networks.needs_commitment(generator)
                        commitment_value = PSCOPF.get_commitment_value(gen_schedule, ts, scenario)
                    else
                        commitment_value = "X"
                    end

                    Base.write(schedule_file_l, @sprintf("%20s%8s%25s%10s%20s%10s%16.8E%16.8E%16.8E%16.8E%15s%8s%15s\n",
                                                    ech, decider, gen_id, "PROD", ts, scenario,
                                                    inject_value, power_value, load_shed, power_capped, power_firmness,
                                                    commitment_value, commitment_firmness)
                                )
                end
            end
        end

        for bus_id in Networks.get_id.(Networks.get_buses(network))
            for ts in TS
                for scenario in scenarios
                    power_value = PSCOPF.get_uncertainties(uncertainties_at_ech, bus_id, ts, scenario)
                    load_shed = PSCOPF.get_loss_of_load(schedule, bus_id, ts, scenario)
                    inject_value = power_value - load_shed

                    Base.write(schedule_file_l, @sprintf("%20s%8s%25s%8s%20s%10s%16.8E%16.8E%16.8E%16.8E%8s%6s%8s\n",
                                                        ech, decider, bus_id, "LOAD", ts, scenario,
                                                        inject_value, power_value, load_shed, 0., "X",
                                                        "X", "X")
                        )
                end
            end
        end
    end
end

end #module PSCOPFio

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# For CI, make sure the downloads do not hang
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

using Test

@testset "quasi-Rosenbrock derivative" begin
    include("rosenbrock.jl")
end

@testset "AIBECS test" begin
    include("AIBECS_test.jl")
end


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

name = "libpng"
version = v"1.6.37"

# Collection of sources required to build libpng
sources = [
    ArchiveSource("https://downloads.sourceforge.net/libpng/libpng-$(version).tar.gz",
               "daeb2620d829575513e35fecc83f0d3791a620b9b93d800b763542ece9390fb4"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libpng-*/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" ..
make -j${ncore}
make install
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    "Zlib_jll",
]

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


text:
"""
```
circle!(builder, center, radius; n=20)
```

Add points and facets approximating a circle.
"""
function circle!(builder::SimplexGridBuilder, center, radius; n=20)
    points=[point!(builder, center[1]+radius*sin(t),center[2]+radius*cos(t)) for t in range(0,2π,length=n)]
    for i=1:n-1
        facet!(builder,points[i],points[i+1])
    end
    facet!(builder,points[end],points[1])
    builder
end

function bregions!(builder::SimplexGridBuilder,grid,regionlist;facetregions=nothing)
    save_facetregion=builder.current_facetregion
    if isnothing(facetregions)
        facetregions=fill(builder.current_facetregion,length(regionlist))
    end
    if isa(facetregions,Number)
        facetregions=fill(facetregions,length(regionlist))
    end
    coord=grid[Coordinates]
    bfnodes=grid[BFaceNodes]
    bfregions=grid[BFaceRegions]

    for ibface=1:size(bfnodes,2)
        ireg=findfirst(i->i==bfregions[ibface],regionlist)
        if ireg!=nothing
            facetregion!(builder,facetregions[ireg])
            @views p1=point!(builder, coord[:,bfnodes[1,ibface]])
            @views p2=point!(builder, coord[:,bfnodes[2,ibface]])
            facet!(builder,p1,p2)
        end
    end
    facetregion!(builder,save_facetregion)
end

"""
```
rect2d!(builder, sw, ne; facetregions=nothing)
```

Add points and facets describing a rectangle via points describing its
south-west and north-east corners. On default, the corresponding facet
regions are deduced from the current facetregion. Alternatively,
a 4-vector of facetregions can be passed.
"""
function rect2d!(builder::SimplexGridBuilder,PA,PB;facetregions=nothing, nx=1, ny=1)
    save_facetregion=builder.current_facetregion
    if isnothing(facetregions)
        facetregions=fill(builder.current_facetregion,4)
    end
    if isa(facetregions,Number)
        facetregions=fill(facetregions,4)
    end
    p00=point!(builder,PA[1],PA[2])
    p10=point!(builder,PB[1],PA[2])
    p11=point!(builder,PB[1],PB[2])
    p01=point!(builder,PA[1],PB[2])

    x=range(PA[1],PB[1],length=nx+1)
    for i=1:nx
        facetregion!(builder,facetregions[1])
        p1=point!(builder,x[i],PA[2])
        p2=point!(builder,x[i+1],PA[2])
        facet!(builder,p1,p2)

        facetregion!(builder,facetregions[3])
        p1=point!(builder,x[i],PB[2])
        p2=point!(builder,x[i+1],PB[2])
        facet!(builder,p1,p2)
    end

    y=range(PA[2],PB[2],length=ny+1)
    for i=1:ny
        facetregion!(builder,facetregions[2])
        p1=point!(builder,PB[1],y[i])
        p2=point!(builder,PB[1],y[i+1])
        facet!(builder,p1,p2)

        facetregion!(builder,facetregions[4])
        p1=point!(builder,PA[1],y[i])
        p2=point!(builder,PA[1],y[i+1])
        facet!(builder,p1,p2)
    end

    facetregion!(builder,save_facetregion)
    builder
end

"""
```
rect3d!(builder, bsw, tne; facetregions=nothing)
```

Add points and facets describing a qudrilateral via points describing its
bottom south-west and top north-east corners. On default, the corresponding facet
regions are deduced from the current facetregion. Alternatively,
a 6-vector of facetregions can be passed (in the sequence s-e-n-w-b-t)
"""
function rect3d!(builder::SimplexGridBuilder,PA,PB;facetregions=nothing)
    save_facetregion=builder.current_facetregion
    if isnothing(facetregions)
        facetregions=fill(builder.current_facetregion,6)
    end
    if isa(facetregions,Number)
        facetregions=fill(facetregions,6)
    end
    p1=point!(builder,PA[1],PA[2],PA[3])
    p2=point!(builder,PB[1],PA[2],PA[3])
    p3=point!(builder,PB[1],PB[2],PA[3])
    p4=point!(builder,PA[1],PB[2],PA[3])
    p5=point!(builder,PA[1],PA[2],PB[3])
    p6=point!(builder,PB[1],PA[2],PB[3])
    p7=point!(builder,PB[1],PB[2],PB[3])
    p8=point!(builder,PA[1],PB[2],PB[3])

    facetregion!(builder,facetregions[1])
    facet!(builder,p1 ,p2 ,p3 ,p4)
    facetregion!(builder,facetregions[2])
    facet!(builder,p5 ,p6 ,p7 ,p8)
    facetregion!(builder,facetregions[3])
    facet!(builder,p1 ,p2 ,p6 ,p5)
    facetregion!(builder,facetregions[4])
    facet!(builder,p2 ,p3 ,p7 ,p6)
    facetregion!(builder,facetregions[5])
    facet!(builder,p3 ,p4 ,p8 ,p7)
    facetregion!(builder,facetregions[6])
    facet!(builder,p4 ,p1 ,p5 ,p8)
    facetregion!(builder,save_facetregion)
    builder
end

function refine(coord,tri)
    # Short and trivial method, without creating edges,
    # and with doubling points (relying on binned point list
    # upon insertion into builder

    function pinsert(p,istop)
        append!(coord,p)
        return size(coord,2)
    end

    rscale(p)=@. p/sqrt(p[1]^2+p[2]^2+p[3]^2)

    newtri=ElasticArray{Cint}(undef,3,0)
    istop=size(coord,2)
    ntri=size(tri,2)
    @views for itri=1:ntri
        i1=tri[1,itri]
        i2=tri[2,itri]
        i3=tri[3,itri]
        p1=coord[:,i1]
        p2=coord[:,i2]
        p3=coord[:,i3]

        # new points inserted on the spere
        i12=pinsert(rscale((p1+p2)/2),istop)
        i13=pinsert(rscale((p1+p3)/2),istop)
        i23=pinsert(rscale((p2+p3)/2),istop)

        append!(newtri,(i1,i12,i13))
        append!(newtri,(i2,i12,i23))
        append!(newtri,(i3,i23,i13))
        append!(newtri,(i12,i13,i23))
    end
    coord,newtri
end

"""
```
sphere!(builder, center, radius; nref=3)
```

Add points and facets approximating a sphere. `nref` is a refinement level.
"""
function sphere!(builder::SimplexGridBuilder, center, radius; nref=3)
    # Initial octahedron
    q=1.0/sqrt(2)

    coord=ElasticArray([-q -q  0;
                        -q  q  0;
                        q  q  0;
                        q -q  0;
                        0   0 -1;
                        0   0  1]')

    tri=[1 2 5;
         2 3 5;
         3 4 5;
         4 1 5;
         1 2 6;
         2 3 6;
         3 4 6;
         4 1 6]'

    for iref=1:nref
        coord,tri=refine(coord,tri)
    end

    @views pts=[point!(builder,(radius*coord[:,i].+center)) for i=1:size(coord,2) ]

    for i=1:size(tri,2)
        facet!(builder, pts[tri[1,i]],pts[tri[2,i]],pts[tri[3,i]])
    end

    builder
end


text:
"""
Representation of a state for the SingleOCcludedCrosswalk POMDP,
depends on ADM VehicleState type
"""
mutable struct SingleOCState
    crash ::Bool
    ego::VehicleState
    ped::VehicleState
end

const TERMINAL_STATE = SingleOCState(false, VehicleState(), VehicleState())

# copy b to a
function Base.copyto!(a::SingleOCState, b::SingleOCState)
    a.crash = b.crash
    a.ego = b.ego
    a.ped = b.ped
end

function Base.hash(s::SingleOCState, h::UInt64 = zero(UInt64))
    return hash(s.crash, hash(s.ego, hash(s.ped, h)))
end

function Base.:(==)(a::SingleOCState, b::SingleOCState)
    return a.crash == b.crash && a.ego == b.ego && a.ped == b.ped
end

function Base.isapprox(veh1::VehicleState, veh2::VehicleState; kwargs...)
    isapprox(veh1.posG, veh2.posG, kwargs...) && isapprox(veh1.posG, veh2.posG, kwargs...) && isapprox(veh1.v, veh2.v, kwargs...)
end

#### Observaton type

const SingleOCObs = SingleOCState

#### Action type

mutable struct SingleOCAction
    acc::Float64
end

function Base.copyto!(a::SingleOCAction, b::SingleOCAction)
    a.acc = b.acc
end

function Base.hash(a::SingleOCAction, h::UInt64 = zero(UInt64))
    return hash(a.acc, h)
end

function Base.:(==)(a::SingleOCAction, b::SingleOCAction)
    return a.acc == b.acc
end

#### POMDP type

mutable struct SingleOCPOMDP <: POMDP{SingleOCState, SingleOCAction, SingleOCObs}
    env::CrosswalkEnv
    ego_type::VehicleDef
    ped_type::VehicleDef
    max_acc::Float64
    pos_res::Float64
    vel_res::Float64
    x_start::Float64
    y_start::Float64
    x_goal::Float64
    y_goal::Float64
    ΔT::Float64 # decision frequency
    p_birth::Float64
    no_ped_prob::Float64
    no_ped::Bool
    v_noise::Float64
    pos_obs_noise::Float64
    vel_obs_noise::Float64
    collision_cost::Float64
    action_cost::Float64
    goal_reward::Float64
    γ::Float64 # discount factor
end

function SingleOCPOMDP(; env::CrosswalkEnv = CrosswalkEnv(),
                   ego_type::VehicleDef = VehicleDef(),
                   ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0),#,env.params.crosswalk_width),
                   max_acc::Float64 = 2.0,
                   pos_res::Float64 = 1.,
                   vel_res::Float64 = 1.,
                   x_start::Float64 = 5.,
                   y_start::Float64 = -env.params.crosswalk_length/4,
                   x_goal::Float64 = env.params.roadway_length/2 + 2*env.params.crosswalk_width,
                   y_goal::Float64 = env.params.crosswalk_length/4,
                   ΔT::Float64  = 0.5,
                   p_birth::Float64 = 0.3,
                   no_ped_prob::Float64 = 0.3,
                   no_ped::Bool = false,
                   v_noise::Float64 = 1.,
                   pos_obs_noise::Float64 = 1.0,
                   vel_obs_noise::Float64 = 1.0,
                   collision_cost::Float64 = -1.,
                   action_cost::Float64 = 0.0,
                   goal_reward::Float64 = 1.,
                   γ::Float64  = 0.9)
    return SingleOCPOMDP(env,
                   ego_type,
                   ped_type,
                   max_acc,
                   pos_res,
                   vel_res,
                   x_start,
                   y_start,
                   x_goal,
                   y_goal,
                   ΔT,
                   p_birth,
                   no_ped_prob,
                   no_ped,
                   v_noise,
                   pos_obs_noise,
                   vel_obs_noise,
                   collision_cost,
                   action_cost,
                   goal_reward,
                   γ)
end

### REWARD MODEL ##################################################################################

# function POMDPs.reward(pomdp::SingleOCPOMDP, s::SingleOCState, a::SingleOCAction, sp::SingleOCState)
#     r = 0.
#     if sp == TERMINAL_STATE
#         return r
#     end
#     if sp.crash
#         r += pomdp.collision_cost
#     end
#     if sp.ego.posG.x >= pomdp.x_goal
#         r += pomdp.goal_reward
#     elseif a.acc > 0.
#         r += pomdp.action_cost
#     else
#         r += pomdp.action_cost
#     end
#     return r
# end

# other method for SARSOP
function POMDPs.reward(pomdp::SingleOCPOMDP, s::SingleOCState, a::SingleOCAction)
    r = 0.
    if s.crash
        r += pomdp.collision_cost
    end
    if s.ego.posG.x >= pomdp.x_goal
        r += pomdp.goal_reward
    elseif a.acc > 0.
        r += pomdp.action_cost
    else
        r += pomdp.action_cost
    end
    if isterminal(pomdp, s)
        return 0.0
    end
    return r
end

function POMDPs.isterminal(pomdp::SingleOCPOMDP, s::SingleOCState)
    return s == TERMINAL_STATE
    # return s.ego.posG.x >= pomdp.x_goal || s.crash
end

### HELPERS

"""
    most_likely_state(d::SparseCat)
returns the most likely state given distribution d
"""
function most_likely_state(d::SparseCat)
    val, ind = argmax(d.probs)
    return d.vals[ind]
end

function POMDPs.discount(pomdp::SingleOCPOMDP)
    return pomdp.γ
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
nsoundings = 24
ncores = 48
nchainspersounding = 4
@assert mod(ncores,nchainspersounding) == 0
nparallelsoundings = Int(ncores/nchainspersounding)
nsequentialiters = ceil(Int, nsoundings/nparallelsoundings)
@info "starting"
for iter = 1:nsequentialiters
    if iter<nsequentialiters
        @show s = (iter-1)*nparallelsoundings+1:iter*nparallelsoundings
    else
        @show s = (iter-1)*nparallelsoundings+1:nsoundings
    end
    @show pids = [(i-1)*nchainspersounding+2:i*nchainspersounding+1 for i = 1:length(s)]
    # set up dBzdt operators for sounding numbers in s, at pids, with fnames
end

## trying parallel, workers calling other workers
using Distributed
addprocs(4)
@info "workers are $(workers())"
##
@everywhere function sleeper(t)
    @info "sleeping at $(myid())"
    sleep(t)
end
@everywhere function trycall(t, pids)
    @sync for pid in pids
        @async remotecall_fetch(sleeper, pid, t)
    end
end
## now use them
ntimes = 2
tsleep = 2
@time for i = 1:ntimes
    @sync for (jpid, pids) in enumerate([[2,3],[4,5]])
        @async remotecall_fetch(trycall, jpid*2, tsleep, pids)
        println()
    end
end


text:
module VariableTemplates

export varsize, Vars, Grad, @vars, varsindex, varsindices

using StaticArrays
using LinearAlgebra

"""
    varsindex(S, p::Symbol, [sp::Symbol...])

Return a range of indices corresponding to the property `p` and
(optionally) its subproperties `sp` based on the template type `S`.

# Examples
```julia-repl
julia> S = @vars(x::Float64, y::Float64)
julia> varsindex(S, :y)
2:2

julia> S = @vars(x::Float64, y::@vars(α::Float64, β::SVector{3, Float64}))
julia> varsindex(S, :y, :β)
3:5
```
"""
function varsindex(::Type{S}, insym::Symbol) where {S <: NamedTuple}
    offset = 0
    for varsym in fieldnames(S)
        T = fieldtype(S, varsym)
        if T <: Real
            offset += 1
            varrange = offset:offset
        elseif T <: SHermitianCompact
            LT = StaticArrays.lowertriangletype(T)
            N = length(LT)
            varrange = offset .+ (1:N)
            offset += N
        elseif T <: StaticArray
            N = length(T)
            varrange = offset .+ (1:N)
            offset += N
        else
            varrange = offset .+ (1:varsize(T))
            offset += varsize(T)
        end
        if insym == varsym
            return varrange
        end
    end
    error("symbol '$insym' not found")
end

# return `Symbol`s unchanged.
wrap_val(sym) = sym
# wrap integer in `Val`
wrap_val(i::Int) = Val(i)

# We enforce that calls to `varsindex` on
# an `NTuple` must be unrapped in `Val`.
# This is enforced to synchronize failures
# on the CPU and GPU, rather than allowing
# CPU-working and GPU-breaking versions.
# This means that users _must_ wrap `sym`
# in `Val`, which can be done with `wrap_val`
# above.
Base.@propagate_inbounds function varsindex(
    ::Type{S},
    sym::Symbol,
    rest...,
) where {S <: Union{NamedTuple, Tuple}}
    vi = varsindex(fieldtype(S, sym), rest...)
    return varsindex(S, sym)[vi]
end
Base.@propagate_inbounds function varsindex(
    ::Type{S},
    ::Val{i},
    rest...,
) where {i, S <: Union{NamedTuple, Tuple}}
    et = eltype(S)
    offset = (i - 1) * varsize(et)
    vi = varsindex(et, rest...)
    return (vi.start + offset):(vi.stop + offset)
end

Base.@propagate_inbounds function varsindex(
    ::Type{S},
    ::Val{i},
) where {i, S <: SArray}
    return i:i
end

"""
    varsindices(S, ps::Tuple)
    varsindices(S, ps...)

Return a tuple of indices corresponding to the properties
specified by `ps` based on the template type `S`. Properties
can be specified using either symbols or strings.

# Examples
```julia-repl
julia> S = @vars(x::Float64, y::Float64, z::Float64)
julia> varsindices(S, (:x, :z))
(1, 3)

julia> S = @vars(x::Float64, y::@vars(α::Float64, β::SVector{3, Float64}))
julia> varsindices(S, "x", "y.β")
(1, 3, 4, 5)
```
"""
function varsindices(::Type{S}, vars::Tuple) where {S <: NamedTuple}
    indices = Int[]
    for var in vars
        splitvar = split(string(var), '.')
        append!(indices, collect(varsindex(S, map(Symbol, splitvar)...)))
    end
    Tuple(indices)
end
varsindices(::Type{S}, vars...) where {S <: NamedTuple} = varsindices(S, vars)

"""
    varsize(S)

The number of elements specified by the template type `S`.
"""
varsize(::Type{T}) where {T <: Real} = 1
varsize(::Type{Tuple{}}) = 0
varsize(::Type{NamedTuple{(), Tuple{}}}) = 0
varsize(::Type{SArray{S, T, N} where L}) where {S, T, N} = prod(S.parameters)

include("var_names.jl")

# TODO: should be possible to get rid of @generated
@generated function varsize(::Type{S}) where {S}
    types = fieldtypes(S)
    isempty(types) ? 0 : sum(varsize, types)
end

function process_vars!(syms, typs, expr)
    if expr isa LineNumberNode
        return
    elseif expr isa Expr && expr.head == :block
        for arg in expr.args
            process_vars!(syms, typs, arg)
        end
        return
    elseif expr.head == :(::)
        push!(syms, expr.args[1])
        push!(typs, expr.args[2])
        return
    else
        error("Invalid expression")
    end
end

"""
    @vars(var1::Type1, var2::Type2)

A convenient syntax for describing a `NamedTuple` type.

# Example
```julia
julia> @vars(a::Float64, b::Float64)
NamedTuple{(:a, :b),Tuple{Float64,Float64}}
```
"""
macro vars(args...)
    syms = Any[]
    typs = Any[]
    for arg in args
        process_vars!(syms, typs, arg)
    end
    :(NamedTuple{$(tuple(syms...)), Tuple{$(esc.(typs)...)}})
end

struct GetVarError <: Exception
    sym::Symbol
end
struct SetVarError <: Exception
    sym::Symbol
end

abstract type AbstractVars{S, A, offset} end

"""
    Vars{S,A,offset}(array::A)

Defines property overloading for `array` using the type `S` as a template. `offset` is used to shift the starting element of the array.
"""
struct Vars{S, A, offset} <: AbstractVars{S, A, offset}
    array::A
end
Vars{S}(array) where {S} = Vars{S, typeof(array), 0}(array)

Base.parent(v::AbstractVars) = getfield(v, :array)
Base.eltype(v::AbstractVars) = eltype(parent(v))
Base.propertynames(::AbstractVars{S}) where {S} = fieldnames(S)
Base.similar(v::AbstractVars) = typeof(v)(similar(parent(v)))

@generated function Base.getproperty(
    v::Vars{S, A, offset},
    sym::Symbol,
) where {S, A, offset}
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :($T(array[$(offset + 1)]))
            offset += 1
        elseif T <: SHermitianCompact
            LT = StaticArrays.lowertriangletype(T)
            N = length(LT)
            retexpr = :($T($LT($([:(array[$(offset + i)]) for i in 1:N]...))))
            offset += N
        elseif T <: StaticArray
            N = length(T)
            retexpr = :($T($([:(array[$(offset + i)]) for i in 1:N]...)))
            offset += N
        else
            retexpr = :(Vars{$T, A, $offset}(array))
            offset += varsize(T)
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(GetVarError(sym))))
    expr
end

@generated function Base.setproperty!(
    v::Vars{S, A, offset},
    sym::Symbol,
    val,
) where {S, A, offset}
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :(array[$(offset + 1)] = val)
            offset += 1
        elseif T <: SHermitianCompact
            LT = StaticArrays.lowertriangletype(T)
            N = length(LT)
            retexpr = :(
                array[($(offset + 1)):($(offset + N))] .=
                    $T(val).lowertriangle
            )
            offset += N
        elseif T <: StaticArray
            N = length(T)
            retexpr = :(array[($(offset + 1)):($(offset + N))] .= val[:])
            offset += N
        else
            offset += varsize(T)
            continue
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(SetVarError(sym))))
    expr
end

"""
    Grad{S,A,offset}(array::A)

Defines property overloading along slices of the second dimension of `array` using the type `S` as a template. `offset` is used to shift the starting element of the array.
"""
struct Grad{S, A, offset} <: AbstractVars{S, A, offset}
    array::A
end
Grad{S}(array) where {S} = Grad{S, typeof(array), 0}(array)

@generated function Base.getproperty(
    v::Grad{S, A, offset},
    sym::Symbol,
) where {S, A, offset}
    if A <: SubArray
        M = size(fieldtype(A, 1), 1)
    else
        M = size(A, 1)
    end
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :(SVector{$M, $T}(
                $([:(array[$i, $(offset + 1)]) for i in 1:M]...),
            ))
            offset += 1
        elseif T <: StaticArray
            N = length(T)
            retexpr = :(SMatrix{$M, $N, $(eltype(T))}(
                $([:(array[$i, $(offset + j)]) for i in 1:M, j in 1:N]...),
            ))
            offset += N
        else
            retexpr = :(Grad{$T, A, $offset}(array))
            offset += varsize(T)
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(GetVarError(sym))))
    expr
end

@generated function Base.setproperty!(
    v::Grad{S, A, offset},
    sym::Symbol,
    val::AbstractArray,
) where {S, A, offset}
    if A <: SubArray
        M = size(fieldtype(A, 1), 1)
    else
        M = size(A, 1)
    end
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :(array[:, $(offset + 1)] = val)
            offset += 1
        elseif T <: StaticArray
            N = length(T)
            retexpr = :(
                array[
                    :,
                    # static range is used here to force dispatch to
                    # StaticArrays setindex! because generic setindex! is slow
                    StaticArrays.SUnitRange($(offset + 1), $(offset + N)),
                ] = val
            )
            offset += N
        else
            offset += varsize(T)
            continue
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(SetVarError(sym))))
    expr
end

export unroll_map, @unroll_map
"""
    @unroll_map(f::F, N::Int, args...) where {F}
    unroll_map(f::F, N::Int, args...) where {F}

Unroll N-expressions and wrap arguments in `Val`.
"""
@generated function unroll_map(f::F, ::Val{N}, args...) where {F, N}
    quote
        Base.@_inline_meta
        Base.Cartesian.@nexprs $N i -> f(Val(i), args...)
    end
end
macro unroll_map(func, N, args...)
    @assert func.head == :(->)
    body = func.args[2]
    pushfirst!(body.args, :(Base.@_inline_meta))
    quote
        $unroll_map($(esc(func)), Val($(esc(N))), $(esc(args))...)
    end
end

export vuntuple
"""
    vuntuple(f::F, N::Int)

Val-Unroll ntuple: wrap `ntuple`
arguments in `Val` for unrolling.
"""
vuntuple(f::F, N::Int) where {F} = ntuple(i -> f(Val(i)), Val(N))

# Inside unroll_map expressions, all indexes `i`
# are wrapped in `Val`, so we must redirect
# these methods:
Base.@propagate_inbounds Base.getindex(t::Tuple, ::Val{i}) where {i} =
    Base.getindex(t, i)
Base.@propagate_inbounds Base.getindex(a::SArray, ::Val{i}) where {i} =
    Base.getindex(a, i)

Base.@propagate_inbounds function Base.getindex(
    v::Vars{NTuple{N, T}, A, offset},
    ::Val{i},
) where {N, T, A, offset, i} # 1 <= i <= N
    return Vars{T, A, offset + (i - 1) * varsize(T)}(parent(v))
end

Base.@propagate_inbounds function Base.getindex(
    v::Grad{NTuple{N, T}, A, offset},
    ::Val{i},
) where {N, T, A, offset, i} # 1 <= i <= N
    return Grad{T, A, offset + (i - 1) * varsize(T)}(parent(v))
end

"""
    getpropertyorindex

An interchangeably and nested-friendly
`getproperty`/`getindex`.
"""
function getpropertyorindex end

# Redirect to Base getproperty/getindex:
Base.@propagate_inbounds getpropertyorindex(t::Tuple, ::Val{i}) where {i} =
    Base.getindex(t, i)
Base.@propagate_inbounds getpropertyorindex(
    a::AbstractArray,
    ::Val{i},
) where {i} = Base.getindex(a, i)
Base.@propagate_inbounds getpropertyorindex(v::AbstractVars, s::Symbol) =
    Base.getproperty(v, s)
Base.@propagate_inbounds getpropertyorindex(
    v::AbstractVars,
    ::Val{i},
) where {i} = Base.getindex(v, Val(i))

# Only one element left:
Base.@propagate_inbounds getpropertyorindex(
    v::AbstractVars,
    t::Tuple{A},
) where {A} = getpropertyorindex(v, t[1])
Base.@propagate_inbounds getpropertyorindex(
    a::AbstractArray,
    t::Tuple{A},
) where {A} = getpropertyorindex(a, t[1])

# Peel first element from tuple and recurse:
Base.@propagate_inbounds getpropertyorindex(v::AbstractVars, t::Tuple) =
    getpropertyorindex(getpropertyorindex(v, t[1]), Tuple(t[2:end]))

# Redirect to getpropertyorindex:
Base.@propagate_inbounds Base.getproperty(v::AbstractVars, tup_chain::Tuple) =
    getpropertyorindex(v, tup_chain)
Base.@propagate_inbounds Base.getindex(v::AbstractVars, tup_chain::Tuple) =
    getpropertyorindex(v, tup_chain)

include("flattened_tup_chain.jl")

end # module


text:
#= ATTENTION:
   your own folder is considered as the current working directory
   for running your solver
=#

include("../../../libSPP/librarySPP.jl")

function main()
    println("Etudiant(e)s : nom1 et nom2")

    # Collecting the names of instances to solve located in the folder Data ----
    target = "../../../Data"
    fnames = getfname(target)

    fres = splitdir(splitdir(pwd())[end-1])[end]
    io = open("../res/"*fres*".res", "w")
    for instance = 1:length(fnames)

        # Load one numerical instance ------------------------------------------
        C, A = loadSPP(string(target,"/",fnames[instance]))

        zInit = 0 ; zBest = 0 ; t1 =0.0 ; t2 = 0.0

        #=
         votre code :
         t1 = @elapsed x, zInit = GreedyConstruction(C, A)
         t2 = @elapsed xbest, zBest = GreedyImprovement(C, A, x, zInit)
        =#

        # Saving results -------------------------------------------------------
        println(io, fnames[instance], " ", zInit, " ", zBest, " ", t1, " ", t2, " ", t1+t2)
    end
    close(io)

end

main()

max_stars_repo_name: xgandibleux/solveSPP

max_issues_repo_name: xgandibleux/solveSPP

max_forks_repo_name: xgandibleux/solveSPP


text:
"Value and index tuple: enables sortandperm_radix"
Valindex{T, S<:Integer} = Tuple{T, S}

isbits(::Type{Valindex{T,S}}) where {T,S} = Base.isbits(T)
uint_mapping(o, vi::Valindex{T,S}) where {T,S} = uint_mapping(o, vi[1]) # enable sorting
ordtype(o, vs::Valindex{T,S}) where {T,S} = Base.ordtype(o, vs[1])
sizeof(::Type{Valindex{T,S}}) where {T,S} = Base.sizeof(T)

"""
	sortandperm(v, alg, o)

returns both the sort(v) as well as sortperm(v)
"""
function fsortandperm(v, alg::RadixSortAlg; order::Ordering = Base.ForwardOrdering())
	fsortandperm_radix(v, order = order)
end

function f_sortandperm_radix(v::AbstractVector{T}; order::Ordering= Base.ForwardOrdering()) where T
	vv = Valindex{T,Int}[(vv,i) for (i,vv) in enumerate(v)]
	sort!(vv, alg=RadixSort, order = order)
	vv
end

function fsortandperm_radix(v::AbstractVector{T}; order::Ordering= Base.ForwardOrdering()) where T
	vv = f_sortandperm_radix(v, order = order)
	res = Int[v[2] for v in vv]
	val = T[v[1] for v in vv]
	return (val, res)
end

"""
	sortperm_radix(v, o)

sortperm using the radixsort algorithm
"""
function fsortperm_radix(v; order::Ordering = Base.ForwardOrdering())
	vv = f_sortandperm_radix(v,order = order)
	Int[v[2] for v in vv]
end


text:
using .JSort
using DelimitedFiles: readdlm

function combinepeaks(calib::AbstractVector{V} where V<:ParticleCalibrator, fnames::AbstractVector{T} where T<:AbstractString)
    # We want to backtransform the peaks so that all peaks
    # are given in original channel numbers, which are the same for all
    # experimental data.
    calib = invert.(calib)
    # Load file
    peaks_exp = [Float64[] for f in 1:8, b in 1:8]
    peaks_kinz = [Float64[] for f in 1:8, b in 1:8]
    for (i, file) in readdlm.(fnames, ',') |> enumerate
        for line in eachrow(file)
            front, e, qkinz = line
            front = Int(front) + 1
            for back in 1:8
              push!(peaks_exp[front, back], calib[i][front, back](e))
              push!(peaks_kinz[front, back], qkinz)
            end
        end
    end

    newcalib = ParticleCalibrator()
    for front in 1:8, back in 1:8
        @assert length(peaks_exp[front, back]) > 1 "Provide at least 2 peaks"
        coeff = leastsquares(peaks_exp[front, back], peaks_kinz[front, back]; order=1)
        newcalib[front, back] = Poly(coeff, calib[1][1,1].var)
    end
    return newcalib
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


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

name = "tblis"
version = v"1.2.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/devinamatthews/tblis.git", "3e4c4b82943726c443b6f408c9c9791dcad7a847")
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd tblis/

for i in ./Makefile.* ./configure*; do

    # Building in container forbids -march options
    sed -i "s/-march[^ ]*//g" $i

done

case ${target} in
    # Unlike stated in Wiki,
    # TBLIS automatically detects threading model.
    *"x86_64"*"linux"*"gnu"*)
        export BLI_CONFIG=x86
        export BLI_THREAD=openmp
        ;;
    *"x86_64"*"linux"*"musl"*)
        export BLI_CONFIG=x86
        export BLI_THREAD=pthreads
        export CC=clang
        export CXX=clang++
        ;;
    *"x86_64"*"w64"*)
        # Windows lacks support for some instructions.
        # Building only for AMD processors.
        export BLI_CONFIG=amd
        export BLI_THREAD=openmp
        # Wrapper for posix_memalign calls.
        patch src/memory/aligned_allocator.hpp < ${WORKSPACE}/srcdir/patches/aligned_allocator.hpp.mingw.patch
        # Additional linking parameter needed for MinGW Autoconf.
        # Update Autoconf parameters and refresh.
        cd src/external/tci
        echo "AM_LDFLAGS=-no-undefined" >> Makefile.am
        echo "lib_libtci_la_LDFLAGS=-no-undefined" >> Makefile.am
        update_configure_scripts --reconf
        cd ../../..
        sed -i "s/.*AM_LDFLAGS.*/AM_LDFLAGS=-lpthread -no-undefined/" Makefile.am
        echo "lib_libtblis_la_LDFLAGS=-no-undefined" >> Makefile.am
        update_configure_scripts --reconf
        ;;
    *"x86_64"*"apple"*)
        export BLI_CONFIG=x86
        export BLI_THREAD=openmp
        export CC=gcc
        export CXX=g++
        ;;
    *"x86_64"*"freebsd"*)
        export BLI_CONFIG=x86
        export BLI_THREAD=openmp
        export CC=gcc
        export CXX=g++
        ;;
    *)
        ;;

esac

CFG_OPTION_POSIX="--prefix=${prefix} --build=${MACHTYPE} --host=${target}"
# ./configure will warn about --enable-thread-model but it is actually effective.
CFG_OPTION_TBLIS="--enable-config=${BLI_CONFIG} --enable-thread-model=${BLI_THREAD}"

./configure ${CFG_OPTION_TBLIS} ${CFG_OPTION_POSIX}
make -j${nproc}
make install

# Copy license file
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(cxxstring_abi=:cxx11)),
    Platform("x86_64", "macos"),
    Platform("x86_64", "freebsd"),
    Platform("x86_64", "windows")
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libtci", :tci),
    LibraryProduct("libtblis", :tblis)
]

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

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

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
using BenchmarkConfigSweeps
using DataFrames
using VegaLite

sweepresult = BenchmarkConfigSweeps.load(joinpath(@__DIR__, "build"))

df_raw = DataFrame(sweepresult)

begin
    df = select(df_raw, Not(:trial))
    df[:, :time_ns] = map(t -> minimum(t).time, df_raw.trial)
    rename!(df, :OPENBLAS_NUM_THREADS => :nthreads)
    df
end

df_speedup = combine(groupby(df, [:n])) do g
    t0 = only(g[g.nthreads.==1, :time_ns])
    (; g.nthreads, speedup = t0 ./ g.time_ns)
end

@vlplot(
    data = df_speedup,
    layer = [
        {
            mark = {:line, point = true},
            x = :nthreads,
            y = :speedup,
            color = {:n, type = :ordinal},
        },
        {mark = {type = :rule}, encoding = {y = {datum = 1}}},
    ]
)


text:
using BayesianNonparametrics

img_size = 3
(X, trueWeights) = bars(img_size = img_size)

D = img_size^2
H = DirichletMultinomial(D, 1.0)

model = HDP(H)
initialisation = RandomInitialisation(k = 20)

modelBuffer = init(X, model, initialisation)
model0 = BayesianNonparametrics.extractpointestimate(modelBuffer)

models = train(modelBuffer, HDPHyperparam(), Gibbs(maxiter = 50))

densities = Float64[m.energy for m in models]
activeComponents = Int[length(m.distributions) for m in models]

println(densities)


text:
module mwe

function func()
  println("hello again")
end

end # module

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Serialize "Hello World!" on julia 0.5.2 (not supported)
const img = JLSO._image()
const pkgs = JLSO._pkgs()
const project, manifest = Pkg.TOML.parse.(JLSO._env())
const hw_5 = UInt8[
    0x26, 0x15, 0x87, 0x48, 0x65,
    0x6c, 0x6c, 0x6f, 0x20, 0x57,
    0x6f, 0x72, 0x6c, 0x64, 0x21
]

const datas = Dict(
    :String => "Hello World!",
    :Vector => [0.867244, 0.711437, 0.512452, 0.863122, 0.907903],
    :Matrix => [
        0.400348 0.892196 0.848164;
        0.0183529 0.755449 0.397538;
        0.870458 0.0441878 0.170899
    ],
    :DateTime => DateTime(2018, 1, 28),
    :ZonedDateTime => ZonedDateTime(2018, 1, 28, tz"America/Chicago"),
    :DataFrame => DataFrame(
        :a => collect(1:5),
        :b => [0.867244, 0.711437, 0.512452, 0.863122, 0.907903],
        :c => ["a", "b", "c", "d", "e"],
        :d => [true, true, false, false, true],
    ),
    :Distribution => Normal(50.2, 4.3),
)

#==
for format in (:bson, :julia_serialize)
    for compression in (:none, :gzip, :gzip_fastest, :gzip_smallest)
        fn = "specimens/v2_$(format)_$(compression).jlso"
        time = @elapsed JLSO.save(fn, datas; format=format, compression=compression);
        time = @elapsed JLSO.save(fn, datas; format=format, compression=compression);
        @info "" format compression time size=filesize(fn)
    end
end
==#


text:
export UnitGroup, solvemod, gen_mod_pk,
       disc_log_bs_gs, disc_log_ph, disc_log_mod

function order(x::GenRes{fmpz}, fp::Dict{fmpz, Int64})
  error("missing")
end

@doc """
  is_primitive_root(x::GenRes{fmpz}, M::fmpz, fM::Dict{fmpz, Int64}) -> Bool

>  Given x in Z/MZ, the factorisation of M (in fM), decide if x is primitive.
>  Intrinsically, only makes sense if the units of Z/MZ are cyclic.
""" ->
function is_primitive_root(x::GenRes{fmpz}, M::fmpz, fM::Dict{fmpz, Int64})
  for p=keys(fM)
    if x^divexact(M, p) == 1
      return false
    end
  end
  return true
end

#=
  for p = 2 this is trivial, as <-1, 5> are generators independently of
  everything.
  Assume p>2 is odd.
  Then (Z/p^k)^* = <g> with some unknown g, independently of k
  We want
  <g>/<g^m> = <A> for A = g^a, independently of k, depending on m
  (to avoid having to factor p-1, assuming that m is small)

  Let G_k = <g+ p^kZ>/<g^m + p^k Z> and ord_k(A) the order of A in G_k

  Then ord_k(A) = |G_k|/gcd(|G_k|, a),
          |G_k| = gcd(phi(p^k), m) = gcd((p-1)p^(k-1), m)

  ord_k    (A) = |G_k|      iff gcd(|G_k|    , a) = 1
  ord_(k+1)(A) = |G_(k+1)|  iff gcd(|G_(k+1)|, a) = 1
  Since either |G_k| = |G_(k+1)| or
              p|G_k| = |G_(k+1)|
  we get (from the constant gcd) that wither a is coprime to p, hence
  the gcd will be stable for all subsequent k's as well, or
  |G_k| is stable (other gcd) hence will be stable. In either case,
  gcd(|G_(k+l)|, a) = 1 for all l, thus a is a generator
=#

@doc """
  gen_mod_pk(p::fmpz, mod::fmpz=0) -> fmpz

>  Find an integer x s.th. x is a primtive root for all powers of the (odd) prime p. If mod is non zero, it finds a generator for Z/p^kZ modulo mod powers only.
""" ->
function gen_mod_pk(p::fmpz, mod::fmpz=fmpz(0))
  @assert isodd(p)
  @assert isprime(p)
  gc = gcd(p-1, mod)
  mi = divexact(p-1, gc)
  fp = factor(gc)
  Rp = ResidueRing(FlintZZ, p)
  Rpp = ResidueRing(FlintZZ, p*p)

  g = fmpz(2)
  if is_primitive_root(Rp(g)^mi, gc, fp)
    if Rpp(g)^gc != 1
      return g
    else
      return g+p
    end
  end

  while true
#    g = rand(3:p-2)
    g += 1
    if is_primitive_root(Rp(g)^mi, gc, fp)
      if Rpp(g)^gc != 1
        return g
      else
        return g+p
      end
    end
  end
end

type MapUnitGroupModM{T} <: Map{T, GenResRing{fmpz}}
  header::Hecke.MapHeader

  function MapUnitGroupModM(G::T, R::GenResRing{fmpz}, dexp::Function, dlog::Function)
    r = new()
    r.header = Hecke.MapHeader(G, R, dexp, dlog)
    return r
  end
end

@doc """
  UnitGroup(R::GenResRing{fmpz}) -> FinGenGrpAb, Map

>  The unit group of R = Z/nZ together with the apropriate map.
""" ->
function UnitGroup(R::GenResRing{fmpz}, mod::fmpz=fmpz(0))
  m = R.modulus
  fm = factor(m)

  r = Array(fmpz, 0)
  g = Array(fmpz, 0)
  mi = Array(fmpz, 0)
  for p=keys(fm)
    k = fm[p]
    if gcd(mod, (p-1)*p^(max(0, k-1))) == 1
      continue
    end
    pk = p^k
    if p==2  && iseven(mod)
      if k==1
        continue
      elseif k==2
        push!(r, 2)
        push!(mi, pk)
        gg = fmpz(-1)
        if m == pk
          push(g, gg)
        else
          push!(g, crt(fmpz(-1), pk, fmpz(1), divexact(m, pk)))
        end
      else
        mpk = divexact(m, pk)
        push!(r, 2)
        push!(r, gcd(p^(k-2), mod))  # cannot be trivial since mod is even
        push!(mi, fmpz(4))
        push!(mi, pk)
        if mpk == 1
          push!(g, fmpz(-1))
          push!(g, fmpz(5))
        else
          push!(g, crt(fmpz(-1), pk, fmpz(1), mpk))
          push!(g, crt(fmpz(5), pk, fmpz(1), mpk))
        end
      end
    else
      mpk = divexact(m, pk)
      s = gcd((p-1)*p^(fm[p]-1), mod)
      if s==1
        continue
      end
      push!(r, s)
      push!(mi, pk)
      gg = gen_mod_pk(p, mod)
      gg = powmod(gg, divexact(p-1, gcd(p-1, mod)), m)
      if mpk == 1
        push!(g, gg)
      else
        push!(g, crt(gg, pk, fmpz(1), mpk))
      end
    end
  end

  G = DiagonalGroup(r)
  function dexp(x::FinGenGrpAbElem)
    return prod([R(g[i])^x[i] for i=1:ngens(G)])
  end
  function dlog(x::GenRes{fmpz})
    return G([disc_log_mod(g[i], lift(x), mi[i]) for i=1:ngens(G)])
  end
  return G, MapUnitGroupModM{typeof(G)}(G, R, dexp, dlog)
end

@doc """
  solvemod(a::fmpz, b::fmpz, M::fmpz)

>  Finds x s.th. ax == b mod M.
""" ->
function solvemod(a::fmpz, b::fmpz, M::fmpz)
  #solve ax = b (mod M)
  g = gcd(a, M)
  if g==1
    return invmod(a, M)*b % M
  end
  a = divexact(a, g)
  @assert b%g == 0
  Mg = divexact(M, g)
  return invmod(a, Mg)* divexact(b, g) %Mg
end

@doc """
  disc_log_mod(a::fmpz, b::fmpz, M::fmpz)

>  Computes g s.th. a^g == b mod M. M has to be a power of an odd prime
>  and a a generator for the multiplicative group mod M
""" ->
#solves a^x = b (mod M) for M a prime power
function disc_log_mod(a::fmpz, b::fmpz, M::fmpz)
  fM = factor(M)
  @assert length(keys(fM)) == 1
  p = first(keys(fM))
  if p==2
    if (a+1) % 4 == 0
      if b%4 == 3
        return fmpz(1)
      else
        return fmpz(0)
      end
    elseif a % M ==5
      if b%4 == 3
        b = -b
      end
      g = fmpz(0)
      if (b-5) % 8 == 0
        g += 1
        b = b* invmod(a, M) % M
      end
      @assert (b-1) % 8 == 0
      @assert (a^2-1) % 8 == 0
      F = FlintPadicField(p, fM[p])
      g += 2*lift(divexact(log(F(b)), log(F(a^2))))
      return g
    else
      error("illegal generator mod 2^$(fM[p])")
    end
  end

  ## we could do the odd priems case using the same p-adic
  #  log approach...

  @assert isodd(p)

  Fp = ResidueRing(FlintZZ, p)
  g = disc_log_bs_gs(Fp(a), Fp(b), p-1)
#  println("1st level ", g)

  #so b*a^-r = 1 mod p, a^(p-1) = 1 mod p
  # in fact, a^(p-1) should be a multiplicative generator for
  # (1+pZ) mod p^k for all k
  #= Plan:
    a^x = b, x = r + (p-1)*y, 0<= r < p-1 and y>=0
    a^(p-1)y = ba^-r and now both a^(p-1) and ba^-r are =1 (p)
    if A, B = 1 mod p and A^y = B, then
    A = 1*px_1, B = 1+pb_1, A^y = 1+pyx_1 + p^2...
    so yx_1 = b_1 mod p, solve this,

    By induction:
    A = 1+p^lx_l, B = 1+p^lb_l, then
    A^y = 1+p^lyx_l + p^(2l)..., thus
    yx_l = b_l mod p^l
  =#

  A = powmod(a, p-1, M)
  B = b*powmod(a, -g, M) %M
  @assert B%p == 1
  @assert A%p == 1
  lp = [fM[p]]
  while lp[end] > 1
    push!(lp, div(lp[end]+1, 2))
  end
  push!(lp, 0)
#  println("lp: $lp")
#TODO: too much powering and fiddling. It works and is correct
#      but needs sorting
  #we are in 1+p^lp[i+1]X mod p^lp[i]
  #so X is defined mod p^(lp[i]-lp[i+1])
  #since A=1+p^lp[i+1] = a^(p-1)^lp[i+1], g is adjusted by x*(p-1)*p^lp[i+1]
  for i=length(lp)-2:-1:1
    pim = p^lp[i]
    pim1 = p^lp[i+1]
    pd = p^(lp[i]-lp[i+1])
#    println(typeof(pim1), typeof(g), typeof(p), "pim1=$pim1")
    @assert A %pim1 == 1
    @assert B %pim1 == 1
    Ai = divexact(A-1, pim1) % pd
    Bi = divexact(B-1, pim1) % pd
    #need to solve x Ai = Bi (pd), however Ai might not be coprime
#    println("solve $Ai x X = $Bi mod $pd")
    yi = solvemod(Ai, Bi, pd)
#    println("Ai=$Ai, Bi=$Bi, yi=$yi, pim1 = $pim1")
    g += yi*(p-1)*p^(lp[i+1]-1)
#    println("for pim1=$pim1 yi=$yi g=$g")
#    println(valuation(b*powmod(a, -g, M)-1, p))
    B = B*powmod(A, -yi, M) % M
    A = powmod(a, (p-1)*divexact(pim, p), M)
  end
  return g
end

@doc """
  disc_log_bs_gs{T}(a::GenRes{T}, b::GenRes{T}, o::fmpz)

>  Tries to find g s.th. a^g == b under the assumption that g <= o.
>  Uses Baby-Step-Giant-Step
""" ->
function disc_log_bs_gs{T <: Union{PolyElem, fmpz, fq_nmod_poly, fq_poly, nmod_poly}}(a::GenRes{T}, b::GenRes{T}, o::fmpz)
  b==1 && return fmpz(0)
  b==a && return fmpz(1)
  if o < 100
    ai = inv(a)
    for g=1:Int(o)
      b *= ai
      b==1 && return fmpz(g)
    end
    throw("disc_log failed")
  else
    r = root(o, 2)
    r = Int(r)
    baby = Array{typeof(a), 1}(r)
    baby[1] = parent(a)(1)
    baby[2] = a
    for i=3:r
      baby[i] = baby[i-1]*a
      baby[i] == b && return fmpz(i-1)
    end
    giant = baby[end]*a
    @assert giant == a^r
    b == giant && return fmpz(r)
    giant = inv(giant)
    g = fmpz(0)
    for i=1:r
      b *= giant
      g += r
      f = findfirst(baby, b)
      f >0 && return fmpz(g+f-1)
    end
    throw("disc_log failed")
  end
end

@doc """
  disc_log_ph{T <:PolyElem}(a::Residue{T}, b::Residue{T}, o::fmpz, r::Int)
  disc_log_ph(a::Residue{fmpz}, b::Residue{fmpz}, o::fmpz, r::Int)
  disc_log_ph(a::Residue{fq_nmod_poly}, b::Residue{fq_nmod_poly}, o::fmpz, r::Int)
  disc_log_ph(a::Residue{fq_poly}, b::Residue{fq_poly}, o::fmpz, r::Int)
  disc_log_ph(a::Residue{nmod_poly}, b::Residue{nmod_poly}, o::fmpz, r::Int)

>  Tries to find g s.th. a^g == b under the assumption that ord(a) | o^r
>  Uses Pohlig-Hellmann and Baby-Step-Giant-Step for the size(o) steps.
  """ ->
function disc_log_ph{T <: Union{PolyElem, fmpz, fq_nmod_poly, fq_poly, nmod_poly}}(a::GenRes{T}, b::GenRes{T}, o::fmpz, r::Int)
  #searches for g sth. a^g = b
  # a is of order o^r
  # Pohlig-Hellmann a^g = b => (a^o)^g = b^g
  g = 0
  aa = a^(o^(r-1))
  for s=r:-1:1
    bb = b*inv(a)^g
    gg = disc_log_bs_gs(aa, bb^(o^(s-1)), o)
    g = g+o^(r-s)*gg
  end
  return g
end


text:
using Pkg

Pkg.add(PackageSpec(url="https://github.com/LANIMFE/ApproximationGrids.jl.git"))


text:
####
#### Resampling schemes for particle filters
####

# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

struct ResampleWithESSThreshold{R,T<:Real}
    resampler::R
    threshold::T
end

function ResampleWithESSThreshold(resampler=resample_systematic)
    return ResampleWithESSThreshold(resampler, 0.5)
end

# More stable, faster version of rand(Categorical)
function randcat(rng::Random.AbstractRNG, p::AbstractVector{<:Real})
    T = eltype(p)
    r = rand(rng, T)
    cp = p[1]
    s = 1
    n = length(p)
    while cp <= r && s < n
        @inbounds cp += p[s += 1]
    end
    return s
end

function resample_multinomial(
    rng::Random.AbstractRNG, w::AbstractVector{<:Real}, num_particles::Integer=length(w)
)
    return rand(rng, Distributions.sampler(Distributions.Categorical(w)), num_particles)
end

function resample_residual(
    rng::Random.AbstractRNG,
    w::AbstractVector{<:Real},
    num_particles::Integer=length(weights),
)
    # Pre-allocate array for resampled particles
    indices = Vector{Int}(undef, num_particles)

    # deterministic assignment
    residuals = similar(w)
    i = 1
    @inbounds for j in 1:length(w)
        x = num_particles * w[j]
        floor_x = floor(Int, x)
        for k in 1:floor_x
            indices[i] = j
            i += 1
        end
        residuals[j] = x - floor_x
    end

    # sampling from residuals
    if i <= num_particles
        residuals ./= sum(residuals)
        rand!(rng, Distributions.Categorical(residuals), view(indices, i:num_particles))
    end

    return indices
end

"""
    resample_stratified(rng, weights, n)

Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
generated by stratified resampling.

In stratified resampling `n` ordered random numbers `u₁`, ..., `uₙ` are generated, where
``uₖ \\sim U[(k - 1) / n, k / n)``. Based on these numbers the samples `x₁`, ..., `xₙ`
are selected according to the multinomial distribution defined by the normalized `weights`,
i.e., `xᵢ = j` if and only if
``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
"""
function resample_stratified(
    rng::Random.AbstractRNG, weights::AbstractVector{<:Real}, n::Integer=length(weights)
)
    # check input
    m = length(weights)
    m > 0 || error("weight vector is empty")

    # pre-calculations
    @inbounds v = n * weights[1]

    # generate all samples
    samples = Array{Int}(undef, n)
    sample = 1
    @inbounds for i in 1:n
        # sample next `u` (scaled by `n`)
        u = oftype(v, i - 1 + rand(rng))

        # as long as we have not found the next sample
        while v < u
            # increase and check the sample
            sample += 1
            sample > m &&
                error("sample could not be selected (are the weights normalized?)")

            # update the cumulative sum of weights (scaled by `n`)
            v += n * weights[sample]
        end

        # save the next sample
        samples[i] = sample
    end

    return samples
end

"""
    resample_systematic(rng, weights, n)

Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
generated by systematic resampling.

In systematic resampling a random number ``u \\sim U[0, 1)`` is used to generate `n` ordered
numbers `u₁`, ..., `uₙ` where ``uₖ = (u + k − 1) / n``. Based on these numbers the samples
`x₁`, ..., `xₙ` are selected according to the multinomial distribution defined by the
normalized `weights`, i.e., `xᵢ = j` if and only if
``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
"""
function resample_systematic(
    rng::Random.AbstractRNG, weights::AbstractVector{<:Real}, n::Integer=length(weights)
)
    # check input
    m = length(weights)
    m > 0 || error("weight vector is empty")

    # pre-calculations
    @inbounds v = n * weights[1]
    u = oftype(v, rand(rng))

    # find all samples
    samples = Array{Int}(undef, n)
    sample = 1
    @inbounds for i in 1:n
        # as long as we have not found the next sample
        while v < u
            # increase and check the sample
            sample += 1
            sample > m &&
                error("sample could not be selected (are the weights normalized?)")

            # update the cumulative sum of weights (scaled by `n`)
            v += n * weights[sample]
        end

        # save the next sample
        samples[i] = sample

        # update `u`
        u += one(u)
    end

    return samples
end


text:
mutable struct LeadTime{F, D<:Distribution}
    holding_cost::F
    onorder::Vector{Float64}
    leadtime::Int
    onorder_reset::D
    pull_orders::IdDict{Any,Float64}
    cost_log::Vector{Float64}
end

function LeadTime(holding_cost, leadtime::Number, onorder::NumDist)
    @assert hasmethod(holding_cost, Tuple{LeadTime}) "holding cost must have a method with (::LeadTime) signature"
    doo = parametrify(onorder)
    oo = max.(0.0, rand(doo, leadtime))
    LeadTime(holding_cost, oo, Int(leadtime), doo, IdDict{Any, Float64}(), zeros(0))
end

function LeadTime(holding_cost::Number, leadtime::Number, onorder::NumDist)
    LeadTime(LinearHoldingCost(holding_cost), Int(leadtime), onorder)
end

function LeadTime(leadtime::Number, onorder::NumDist)
    LeadTime(it-> 0.0, leadtime, onorder)
end

RLBase.state(lt::LeadTime) = vec(lt.onorder)
state_size(lt::LeadTime) = lt.leadtime
print_state(lt::LeadTime) = [" on order t+$i" => oo for (i, oo) in enumerate(lt.onorder)]

function Base.push!(lt::LeadTime, quantity, source)
    lt.pull_orders[source] = quantity
    push!(lt.onorder, quantity)
    nothing
end

function dispatch!(lt::LeadTime)
    destination, quantity = only(lt.pull_orders)
    push!(destination, popfirst!(lt.onorder), lt)
    nothing
end

function reward!(lt::LeadTime)
    cost = lt.holding_cost(lt)
    push!(lt.cost_log, cost)
    empty!(lt.pull_orders)
    return -cost
end

inventory_position(lt::LeadTime) = sum(lt.onorder)

function RLBase.reset!(lt::LeadTime)
    empty!(lt.cost_log)
    lt.onorder = max.(0., rand(lt.onorder_reset, lt.leadtime))
end

#Defined in item.jl
(f::LinearHoldingCost)(lt::LeadTime) = f.h*sum(lt.onorder)

Base.show(io::IO, lt::LeadTime{F,D}) where {F,D} = print(io, "LeadTime($(lt.leadtime) periods, $F)")

max_stars_repo_path: src/elementary blocks/leadtime.jl

max_issues_repo_path: src/elementary blocks/leadtime.jl

max_forks_repo_path: src/elementary blocks/leadtime.jl


text:
struct Grid{TF <: Union{Float32, Float64}}
    # Specified by user.
    itot::Int64
    jtot::Int64
    ktot::Int64
    ktoth::Int64

    xsize::TF
    ysize::TF
    zsize::TF

    igc::Int64
    jgc::Int64
    kgc::Int64

    # Calculated.
    icells::Int64
    jcells::Int64
    kcells::Int64

    is::Int64
    js::Int64
    ks::Int64

    ie::Int64
    je::Int64
    ke::Int64
    keh::Int64

    x::Vector{TF}
    xh::Vector{TF}
    y::Vector{TF}
    yh::Vector{TF}
    z::Vector{TF}
    zh::Vector{TF}

    dx::TF
    dy::TF
    dz::Vector{TF}
    dzh::Vector{TF}

    dxi::TF
    dyi::TF
    dzi::Vector{TF}
    dzhi::Vector{TF}
end

function Grid(d::Dict, TF)
    itot = d["itot"]
    jtot = d["jtot"]
    ktot = d["ktot"]
    ktoth = ktot+1

    xsize = d["xsize"]
    ysize = d["ysize"]
    zsize = d["zsize"]

    z_nogc::Vector{Real} = d["z"]

    igc = 1
    jgc = 1
    kgc = 1

    icells = itot + 2*igc
    jcells = jtot + 2*jgc
    kcells = ktot + 2*kgc

    is = igc + 1
    js = jgc + 1
    ks = kgc + 1

    ie = igc + itot
    je = jgc + jtot
    ke = kgc + ktot
    keh = ke + 1

    dx = xsize / itot
    dy = ysize / jtot

    dxi = 1. / dx
    dyi = 1. / dy

    x = collect(range(0.5-igc, length=icells)) .* dx
    y = collect(range(0.5-jgc, length=jcells)) .* dy
    xh = collect(range(-igc, length=icells)) .* dx
    yh = collect(range(-jgc, length=jcells)) .* dy

    z = zeros(kcells)
    z[ks:ke] = z_nogc[:]
    z[ks-1] = -z[ks]
    z[ke+1] = 2*zsize - z[ke]

    zh = zeros(kcells)
    zh[ks] = 0.
    zh[keh] = zsize
    zh[ks+1:ke] .= 0.5*(z[ks:ke-1] + z[ks+1:ke])
    zh[ks-1] = -zh[ks+1]

    dzh = zeros(kcells)
    dzh[2:kcells] .= z[2:kcells] .- z[1:kcells-1]
    dzhi = 1. ./ dzh[:]
    dzh[ks-1] = dzh[ks+1]
    dzhi[ks-1] = dzhi[ks+1]

    dz = zeros(kcells)
    dz[ks:ke] .= zh[ks+1:keh] .- zh[ks:ke]
    dzi = 1. ./ dz[:]
    dz[ks-1] = dz[ks]
    dz[keh] = dz[ke]
    dzi[ks-1] = dzi[ks]
    dzi[keh] = dzi[ke]

    dzi = 1. ./ dz[:]

    g = Grid{TF}(
        itot, jtot, ktot, ktoth,
        xsize, ysize, zsize,
        igc, jgc, kgc,
        icells, jcells, kcells,
        is, js, ks, ie, je, ke, keh,
        x, xh, y, yh, z, zh,
        dx, dy, dz, dzh,
        dxi, dyi, dzi, dzhi)
end


text:
using Random
using VegaLite
using StatsFuns
using DataFrames
using Distributions
using TripletEmbeddings
using GreenVideoExperiments

Random.seed!(4)

# n = 100
# Y = 2*rand(1,n) .- 1
Y = GreenVideoExperiments.load("TaskA", "ground_truth", "30Hz").data[begin:30:end]
Y = Matrix(transpose(filter(x -> !isnan(x), Y)))

triplets = Triplets(Y; f = x -> logistic(2x))
loss = STE(σ = 1/sqrt(2))
X = Embedding(size(Y))

@time violations = fit!(loss, triplets, X; verbose=true, max_iterations=1000)
procrustes!(X, Y)

dfX = DataFrame(embedding = "X", time = 1:length(X), value = X[1,:])
dfY = DataFrame(embedding = "Y", time = 1:length(Y), value = Y[1,:])

vcat(dfX, dfY) |> @vlplot(:line, x = :time, y = :value, color=:embedding, width=600, height=400)

# savefig("figures/1D.svg")


text:
## Disjuncive Domain A1 ∨ A2 ∨ A3 ∨ ... ∨ An
## =================================================================

# This is a simple non-relational domain of disjunctive domains (paramaterised by T)

struct SimpleDisjunctive{T} <: Domain{T}
  values::Set{T}
end

isrelational(::Union{SimpleDisjunctive, Type{SimpleDisjunctive}}) = false

## Set Operations
## ==============
issubset(x::SimpleDisjunctive, y::SimpleDisjunctive) = issubset(x.values, y.values)
isintersect(x::SimpleDisjunctive, y::SimpleDisjunctive) = intersect(x.values, y.values)
domaineq(x::SimpleDisjunctive, y::SimpleDisjunctive) = x.values == y.values
⊔{T}(x::SimpleDisjunctive{T},y::T) = push!(x.values, y)

# Apply f to every abstract element in disjunction
function setmap{T}(f::Function, s::Set{T})
  out = Set{T}()
  for x in s
    push!(out,f(x))
  end
  out
end

# Operations on Disjunctive Domains (to be sound) do a cartesian product
# Real × Real ->  Rea
for op = (:+, :-, :*, :/)
  @eval begin
    function ($op){T}(x::SimpleDisjunctive{T}, y::SimpleDisjunctive{T})
      let op = $op
        fcartproduct($op,T,x,y)
      end
    end
  end
end

# Boolean Valued Functions
# Real × Real ->  Rea
for op = (:>, :>=, :<, :<=, :&, :|)
  @eval begin
    function ($op){T}(x::SimpleDisjunctive{T}, y::SimpleDisjunctive{T})
      let op = $op
        fcartproduct($op,Bool,x,y)
      end
    end
  end
end

# Real ->  Rea
for op = (:sqr,)
  @eval begin
    function ($op){T}(x::SimpleDisjunctive{T})
      let op = $op
        SimpleDisjunctive{T}(setmap(op,x.values))
      end
    end
  end
end

function !(x::SimpleDisjunctive{Bool})
  if x.values == Set{Bool}(true) SimpleDisjunctive{Bool}(Set(false))
  elseif x.values == Set{Bool}(false) SimpleDisjunctive{Bool}(Set(true))
  else x end
end

# Apply f to the cartesian product of values in x and y
function fcartproduct(f::Function, T::DataType,
                      x::SimpleDisjunctive, y::SimpleDisjunctive)
  result = SimpleDisjunctive{T}(Set{T}())
  for args in Iterators.product(x.values,y.values)
    ⊔(result, f(args...))
  end
  result
end


text:
module aoc2019

using BenchmarkTools
import DataFrames: DataFrame
import Dates

export get_answer, print_day_answer, print_all_answers

answerfuncs = Dict{Int, Tuple{Function, Function}}()
inputreads = Dict{Int, Function}()
inputpaths = Dict{Int, String}()

for day = 1:25
    day_string = "d$(lpad(string(day), 2, "0"))"
    file = "$day_string/$day_string.jl"
    if isfile(joinpath(@__DIR__, file))
        include(file)
        funcs = ()
        answerfuncs[day] = (
            getfield(aoc2019, Symbol("$(day_string)_answer1")),
            getfield(aoc2019, Symbol("$(day_string)_answer2"))
        )
        inputreads[day] = getfield(aoc2019, Symbol("$(day_string)_inputread"))
        inputpaths[day] = getfield(aoc2019, Symbol("$(day_string)_inputpath"))
    end
end

function get_answers(; benchmark=false)
    df = DataFrame(
        day = Int[],
        answer1 = String[],
        answer2 = String[],
        benchmarki = String[],
        benchmark1 = String[],
        benchmark2 = String[],
    )
    for (day, (answer1func, answer2func)) in answerfuncs
        inputread = inputreads[day]
        inputpath = inputpaths[day]
        input = inputread(inputpath)
        if benchmark
            benchmarki = @belapsed($inputread($inputpath))
            benchmark1 = @belapsed($answer1func($input))
            benchmark2 = @belapsed($answer2func($input))
        else
            benchmarki = NaN
            benchmark1 = NaN
            benchmark2 = NaN
        end
        push!(df, Dict(
            :day=>day,
            :answer1=>string(answer1func(input)),
            :answer2=>string(answer2func(input)),
            :benchmarki=>BenchmarkTools.prettytime(1e9 * benchmarki),
            :benchmark1=>BenchmarkTools.prettytime(1e9 * benchmark1),
            :benchmark2=>BenchmarkTools.prettytime(1e9 * benchmark2),
        ))
    end
    df
end

end # module


text:
#=
using JuliaDB

t=table([1.0,2.0,3.0], [1.0,1.0,1.0], [1.0,1.0,1.0],[1.0,1.0,1.0],names=[:x,:y,:z,:w]);
d_t = descriptive_stats(t,"x",(:y,:z))
d_t_w = descriptive_stats(t,"x",(:y,:z),weight="w")

d = descriptive_stats(t,"x",:z)

d_w = descriptive_stats(t,"x",:z,weight="w")

@test select(d_t,:MEAN) == [2.0]

@test select(d_t_w,:MEAN) == [2.0]

@test select(d,:MEAN) == [2.0]

@test select(d_w,:MEAN) == [2.0]
=#


text:
using DashBootstrapComponents, DashHtmlComponents

carousel = html_div([
   dbc_carousel(
        id="carousel",
        items=[
            Dict("key" => "1", "src" => "/static/images/slide1.svg"),
            Dict("key" => "2", "src" => "/static/images/slide2.svg"),
            Dict("key" => "3", "src" => "/static/images/slide3.svg"),
        ],
        controls=false,
        indicators=false,
        interval=false,
    ),

    dbc_radioitems(
        id="slide-number",
        options=[
            Dict("label" => "Slide 1", "value" => 0),
            Dict("label" => "Slide 2", "value" => 1),
            Dict("label" => "Slide 3", "value" => 2),
        ],
        value=0,
        inline=true,
    ),
]);

callback!(
    app,
    Output("carousel", "active_index"),
    Input("slide-number", "value"),
) do idx
    return idx
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module VulkanShaders

using Vulkan
using SPIRV
using MLStyle
using Parameters
using DataStructures

import glslang_jll
import Vulkan: ShaderModule, DescriptorSetLayout, PipelineShaderStageCreateInfo
import SPIRV: IR

const glslangValidator = glslang_jll.glslangValidator(x -> x)

include("formats.jl")
include("bindings.jl")
include("compilation.jl")
include("reflection.jl")
include("shaders.jl")

export
        # formats
        ShaderFormat,
        BinaryFormat, TextFormat,
        FormatGLSL, FormatHLSL, FormatSPIRV,
        ShaderFile,
        bytes,

        # shaders
        Shader,
        compile,
        reflect,

        # descriptors
        DescriptorBinding,
        get_descriptor_sets,
        get_bindings,
        create_descriptor_set_layouts

end


text:
"""
    lagrangian_lower_bound(C, c, d, Γ, X, l, dg)

Computes Lagrangian lower bound.

Check section 5.3 _Lagrangian lower bound_ of [publication](https://arxiv.org/abs/1811.06719)
for more information about this algorithm.

# Arguments
- `C`: is a vector of nonnegative first stage costs.
- `c`: is a vector of a nonnegative nominal second stage costs.
- `d`: is a vector of maximal deviations of the costs from their nominal values.
- `Γ`: is a budget, or the amount of uncertainty, which can be allocated to the second stage costs.
- `X`: is a set of feasible solutions represented as a list functions, each of which accepts a list of JuMP variables as an argument and returns a JuMP linear constraint.
- `l`: value of parameter ``l=\\lceil m(1-\\alpha) \\rceil``
- `pd`: an instance of [`ProblemDescriptor`](@ref)
"""
function lagrangianLowerBound(C, c, d, Γ, X, l, dg)
    @assert hasEqualCardinalityProperty(dg)

    tₗ = getProperty("lagrangianLowerBound.overallTimeLimit", parameterType = Int)

    Δt = @elapsed begin
        ϵ = getProperty("lagrangianLowerBound.epsilon", parameterType = Float64)
        invϕ = (sqrt(5) - 1) / 2 # 1/ϕ
        invϕ² = (3 - sqrt(5)) / 2 # 1/ϕ²

        α = 0
        β = 100 # TODO: check this assumption

        h = β - α

        γ = α + invϕ² * h
        fγ, status = relaxedIncrementalProblem(C, c, d, Γ, X, γ, l)
        if status != :Optimal
            @warn "Lagrangian lower bound solved non-optimaly"
            return missing
        end

        δ = α + invϕ * h
        fδ, status = relaxedIncrementalProblem(C, c, d, Γ, X, δ, l)
        if status != :Optimal
            @warn "Lagrangian lower bound solved non-optimaly"
            return missing
        end

        n = ceil(log(ϵ/h)/log(invϕ))
    end
    for i in 0:n-2
        if Δt > tₗ
            @warn "Lagrangian lower bound time limit exceeded (current time is $Δt seconds with time limit $tₗ seconds)"
            break
        end

        Δt += @elapsed begin
            if fγ > fδ
                β, δ, fδ = δ, γ, fγ
                h *=  invϕ
                γ = α + invϕ² * h

                fγ, status = relaxedIncrementalProblem(C, c, d, Γ, X, γ, l)
                if status != :Optimal
                    @warn "Lagrangian lower bound solved non-optimaly"
                    return missing
                end
            else
                α, γ, fγ = γ, δ, fδ
                h *=  invϕ
                δ = α + invϕ * h  # TODO: Why `invphi` here

                fδ, status = relaxedIncrementalProblem(C, c, d, Γ, X, δ, l)
                if status != :Optimal
                    @warn "Lagrangian lower bound solved non-optimaly"
                    return missing
                end
            end
        end
    end

    if (fγ > fδ)
        μ = (α + δ) / 2
        t, status = relaxedIncrementalProblem(C, c, d, Γ, X, μ, l)
        if status != :Optimal
            @warn "Lagrangian lower bound solved non-optimaly"
            return missing
        else
            return t
        end
    else
        μ = (γ + β) / 2
        t, status = relaxedIncrementalProblem(C, c, d, Γ, X, μ, l)
        if status != :Optimal
            @warn "Lagrangian lower bound solved non-optimaly"
            return missing
        else
            return t
        end
    end
end

function relaxedIncrementalProblem(C, c, d, Γ, X, μ, l)
    @assert size(C) == size(c) == size(d)

    n = size(C, 1)

    model = Model(solver = CplexSolver(CPXPARAM_ScreenOutput = getProperty("lagrangianLowerBound.cplexLogging"),
        CPX_PARAM_TILIM = getProperty("lagrangianLowerBound.subproblemTimeLimit")))

    @variable(model, π >= 0)
    if ndims(C) == 1
        @variable(model, x[1:n], Bin)
        @variable(model, 0 <= y[1:n] <= 1)
        @variable(model, u[1:n] >= 0)
        @variable(model, z[1:n] >= 0)
    elseif ndims(C) == 2
        m = size(C, 2)
        @variable(model, x[1:n, 1:m], Bin)
        @variable(model, 0 <= y[1:n, 1:m] <= 1)
        @variable(model, u[1:n, 1:m] >= 0)
        @variable(model, z[1:n, 1:m] >= 0)
    else
        throw(ArgumentError("Unsupported number of dimensions ($(ndims(c)))"))
    end

    @objective(model, Min, vecdot(C, x) + π * Γ + vecdot(y, c) - μ * sum(z) + vecdot(u, d) + μ * l)

    # Constraints defining set of feasible solutions
    for constraint in X
        JuMP.addconstraint(model, constraint(x)) # Assuming X = {Ax=b, x∈{0, 1}ⁿ}
    end

    @constraint(model, -y + π + u .>= 0)

    # Constraints defining set of feasible solutions
    for constraint in X
        JuMP.addconstraint(model, constraint(y)) # Assuming X = {Ax=b, x∈{0, 1}ⁿ}
    end

    @constraint(model, z .<= x)
    @constraint(model, z .<= y)

    status = solve(model)
    if status == :Optimal
        return getobjectivevalue(model), status
    else
        return getobjectivebound(model), status
    end
end