text:
module Uci

import Util

DEBUG_MODE = false
function boot()
    Moves.initSlidersAttack(true)
    Moves.initSlidersAttack(false)
    board = Util.startPositions()
    while(true)
        s = readline(stdin)
        if("uci" in s)
            id()
            option()
            uciok()
        elseif ("debug" in s)
            DEBUG_MODE = (DEBUG_MODE == true) ? false : true
        elseif ("isready" in s)
            println("readyok")
        elseif ("setoption name" in s)
            continue
        elseif ("register" in s)
            continue
        elseif ("ucinewgame" in s)
            continue
        elseif ("position" in s)
            input = split(s," ")
            if(input[2] == "startpos")
                board = Util.startPositions()
            end
            let colour = 'W'
            for index in 4:length(input)
                move = input[index]
                board = makeMoveUCI(move,board,colour)
                colour = (colour == 'W') ? 'B' : 'W'
            end
        elseif ("go" in s)
            move = pickMoveUCI(board,[],'B')
            println("bestmove ", move)
            board = makeMoveUCI(move,board,'B')
        elseif ("stop" in s)
        elseif ("ponderhit" in s)
        elseif ("quit" in s)

        end

end

function id()
    println("id name Europa")
    println("id name Kubrick")
end

function uciok()
    println("uciok")
end

function readyok()
    println("readyok")
end

function bestmove()
end

function copyprotection()
end

function registration()
end

function info()
end

function option()
end

boot()

end

max_stars_repo_name: grahamclyne/chessEngine

max_issues_repo_name: grahamclyne/chessEngine

max_forks_repo_name: grahamclyne/chessEngine


text:
# Based on
# U. Fjordholm, S. Mishra, E. Tadmor, Arbitrarly high-order accurate entropy
# stable essentially nonoscillatory schemes for systems of conservation laws.
# 2012. SIAM. vol. 50. No 2. pp. 544-573

struct FVTecnoScheme <: AbstractFVAlgorithm
  order :: Int
  Nflux :: Function #Entropy stable 2 point flux
  ve    :: Function #Entropy variable
end

function FVTecnoScheme(Nflux;order=2, ve = u -> u)
  FVTecnoScheme(order, Nflux, ve)
end

"""
compute_fluxes!(hh, Flux, u, mesh, dt, alg::FVTecnoScheme, nonscalar::Bool, ::Type{Val{false}})
Numerical flux of Tecno Scheme in 1D
"""
function compute_fluxes!(hh, Flux, u, mesh, dt, alg::FVTecnoScheme, nonscalar::Bool, ::Type{Val{false}})
    order = alg.order; Nflux = alg.Nflux; ve = alg.ve
    M,N = size(u)
    dx = cell_volume(mesh, 1)
    #Eno Reconstrucion
    @inbounds ul=cellval_at_left(1,u,mesh)
    @inbounds ur=cellval_at_right(1,u,mesh)
    RΛ1 = eigen(Flux.Dflux(0.5*(ul+ur)))
    MatR = Vector{typeof(RΛ1.vectors)}(undef,0)
    MatΛ = Vector{typeof(RΛ1.values)}(undef,0)
    push!(MatR,RΛ1.vectors)
    push!(MatΛ,RΛ1.values)
    for j = 2:(N+1)
      @inbounds ul=cellval_at_left(j,u,mesh)
      @inbounds ur=cellval_at_right(j,u,mesh)
      RΛj = eigen(Flux.Dflux(0.5*(ul+ur)))
      push!(MatR,RΛj.vectors); push!(MatΛ,RΛj.values)
    end
    dd = fill(0.0,M,N+1) #Extra numerical diffusion
    k = order - 1
    v = fill!(similar(u),zero(eltype(u)))
    for j in axes(v, 2)
      v[:,j] = ve(u[:,j]) #entropy variables
    end
    vminus = fill(0.0,M,N); vplus = fill(0.0,M,N)
    wminus = fill(0.0,M,N); wplus = fill(0.0,M,N)
    for j in cell_indices(mesh)
      for i in 1:M
        v_eno = get_cellvals(v,mesh,(i,j-k:j+k)...)
        vminus[i,j],vplus[i,j] = reconstruct(v_eno,dx,ENO_Reconstruction(2*order-1))
      end
      wminus[:,j] = MatR[j]*vminus[:,j]
      wplus[:,j] = MatR[j+1]*vplus[:,j]
    end
    wdiff = fill(0.0,M,N+1)
    for j = 2:N
      wdiff[:,j] = wminus[:,j] - wplus[:,j-1]
    end

    if mesh.left_boundary isa Periodic; wdiff[:,1] = wminus[:,1] - wplus[:,N];end
    if mesh.right_boundary isa Periodic; wdiff[:,N+1] = wdiff[:,1];end

    for j in node_indices(mesh)
      dd[:,j] = MatR[j]*[abs(MatΛ[j][i])*wdiff[i,j] for i in 1:M]
    end

    ff = fill(0.0,M,N+1)
    if order == 2
      for j in node_indices(mesh)
          @inbounds ul=cellval_at_left(j,u,mesh)
          @inbounds ur=cellval_at_right(j,u,mesh)
        ff[:,j] = Nflux(ul,ur)
      end
    elseif order == 3 || order == 4
      for j in node_indices(mesh)
        @inbounds ull=cellval_at_left(j-1,u,mesh)
        @inbounds ul=cellval_at_left(j,u,mesh)
        @inbounds ur=cellval_at_right(j,u,mesh)
        @inbounds urr=cellval_at_right(j+1,u,mesh)
        ff[:,j] = 4.0/3.0*Nflux(ul,ur)-1.0/6.0*(Nflux(ull,ur)+Nflux(ul,urr))
      end
    elseif order == 5
      for j in node_indices(mesh)
        @inbounds ulll=cellval_at_left(j-2,u,mesh)
        @inbounds ull=cellval_at_left(j-1,u,mesh)
        @inbounds ul=cellval_at_left(j,u,mesh)
        @inbounds ur=cellval_at_right(j,u,mesh)
        @inbounds urr=cellval_at_right(j+1,u,mesh)
        @inbounds urrr=cellval_at_right(j+2,u,mesh)
        ff[:,j] = 3.0/2.0*Nflux(ul,ur)-3.0/10.0*(Nflux(ull,ur)+Nflux(ul,urr))+
            1.0/30.0*(Nflux(ulll,ur)+Nflux(ull,urr)+Nflux(ul,urrr))
      end
    end
    hh[:,:] = ff - dd
end

ext: jl


text:
export  PlaneStressIncompressibleRivlinSaunders, getStress

@doc """
Pascon, João Paulo.
"Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements."
International Journal of Advanced Structural Engineering (2019): 1-20.
"""->
mutable struct PlaneStressIncompressibleRivlinSaunders
    ρ::Float64 # density
    C1::Float64 # Young's modulus
    C2::Float64 # Poisson's ratio

    # The classical incompressible PlaneStressIncompressibleRivlinSaunders strain energy function
    # I3 = 0
    # W = C1(I1 - 3) + C2(I2 - 2)

    σ0::Array{Float64} # stress at last time step
    σ0_::Array{Float64} # σ0 to be updated in `commitHistory`
    ε0::Array{Float64} # stress at last time step
    ε0_::Array{Float64} # σ0 to be updated in `commitHistory`
end

function PlaneStressIncompressibleRivlinSaundersJacobian(E11::T, E22::T, G12::T, C1::S, C2::S) where {S<:Number, T<:Number}
    f = x->PlaneStressIncompressibleRivlinSaundersStress(x[1],x[2],x[3], C1, C2)
    J = ForwardDiff.jacobian(f, [E11;E22;G12])
    σ = f([E11;E22;G12])
    σ, J
end

function PlaneStressIncompressibleRivlinSaundersStress(E11::T, E22::T, G12::T, C1::S, C2::S) where {S<:Number, T<:Number}
    C11, C22, C12 = 2*E11+1.0, 2*E22+1.0, G12
    det22C = C11*C22 - C12*C12
    C33 = 1/det22C

    I1 = C11 + C22 + C33

    C = [ C11    C12      0
          C12    C22      0
          0       0     C33]

    Cinv = [ C22/det22C    -C12/det22C        0
            -C12/det22C     C11/det22C        0
                  0             0         1.0/C33]
    svol = Cinv
    Siso = 2.0*C1*UniformScaling(1.0) + 2.0*C2*(C + UniformScaling(I1))

    p = -C33*Siso[3,3]

    σ = zeros(T, 3)
    σ[1] = Siso[1,1] + p*svol[1,1]
    σ[2] = Siso[2,2] + p*svol[2,2]
    σ[3] = Siso[1,2] + p*svol[1,2]
    return σ
end

function PlaneStressIncompressibleRivlinSaunders(prop::Dict{String, Any})
    C1 = prop["C1"]; C2 = prop["C2"]; ρ = prop["rho"];
    σ0 = zeros(3); σ0_ = zeros(3); ε0 = zeros(3); ε0_ = zeros(3)
    PlaneStressIncompressibleRivlinSaunders(ρ, C1, C2, σ0, σ0_, ε0, ε0_)
end

function getStress(self::PlaneStressIncompressibleRivlinSaunders,  strain::Array{Float64},  Dstrain::Array{Float64}, Δt::Float64 = 0.0)
    local dΔσdΔε
    C1, C2 = self.C1, self.C2
    E11, E22, G12 = strain[1], strain[2], strain[3]
    σ, dΔσdΔε = PlaneStressIncompressibleRivlinSaundersJacobian(E11, E22, G12, C1, C2)

    self.σ0_ = σ[:]
    self.ε0_ = strain

    #@show σ, dΔσdΔε
    return σ, dΔσdΔε
end

function getTangent(self::PlaneStressIncompressibleRivlinSaunders)
    error("Not implemented")
end

function commitHistory(self::PlaneStressIncompressibleRivlinSaunders)
    self.σ0 = self.σ0_
    self.ε0 = self.ε0_
end

# # W = C1(I1 - 3) + C2(I2 - 2)
# #
# # E11, E22, G12 = strain[1], strain[2], strain[3]
# # E12 = G12/2.0
# #
# #      E11 E12 0
# # E =  E21 E22 0
# #      0   0   E33
# # incompressibility
# # det(2E + 1) = det(C) = det(F^T F) = 1
# # E33 = (1.0/(1 + 2*E11 + G12 + 4*E11*E22 - G12*G12) - 1.0) /2.0
# # C = 2E + 1
# # I1(C) = lambda1 + lambda2 + lambda3 = 3 + 2(E11 + E22 + E33)
# # I2(C) = lambda1*lambda2 + lambda2*lambda3 + lambda3*lambda1
# #       = (1 + 2E11)*(1 + 2E22) + (1 + 2E22)*(1 + 2E33) + (1 + 2E33)*(1 + 2E11) - G12 * G21

# using SymPy
# E11, E22, G12, C1, C2 = @vars E11 E22 G12 C1 C2
# E33 = (1.0/(1 + 2*E11 + G12 + 4*E11*E22 - G12*G12) - 1.0) /2.0
# I1 = 3 + 2*(E11 + E22 + E33)
# I2 = 3 + 4*(E11 + E22 + E33) + 4*E11*E22 + 4*E11*E22 + 4*E11*E33 - G12*G12
# W = C1*(I1 - 3) + C2*(I2 - 3)
# σ = Array{Sym}(undef,3)
# σ[1] = simplify(diff(W, E11))
# σ[2] = simplify(diff(W, E22))
# σ[3] = simplify(diff(W, G12))
# J = Matrix{Sym}(undef, 3, 3)
# for i = 1:3
#     J[i,1] = simplify(diff(σ[i], E11))
#     J[i,2] = simplify(diff(σ[i], E22))
#     J[i,3] = simplify(diff(σ[i], G12))
# end
# s = "J = zeros(3,3);\n"
# for i = 1:3
#     for j = 1:3
#         global s *= "J[$i,$j]="*sympy.julia_code(J[i,j])*"\n"
#     end
# end
# t = "σ = zeros(3)\n"
# for i = 1:3
#     global t *= "σ[$i]="*sympy.julia_code(σ[i])*"\n"
# end


text:
using Test

function binary_search(arr,item)
    low = 0
    high = length(arr)

    while low <= high
        mid = (low + high) ÷ 2
        guess = 0
        try
            guess = arr[mid]
        catch
            return nothing
        end
        if guess == item
            return mid
        elseif guess > item
            high = mid - 1
        else
            low = mid + 1
        end
    end
end

arr = [1,3,5,7,9]
@test binary_search(arr,3) == 2
@test binary_search(arr,-1) == nothing


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

using Markdown
using InteractiveUtils

# ╔═╡ 95eaf04a-08e5-11eb-043f-27646f1bf207
using Dictionaries, NPZ, StatsPlots, StatsPlots.PlotMeasures, Statistics

# ╔═╡ 849adbca-08e5-11eb-1956-cd02f2a6be7a
md"""
# Figure: impact of imbalance
"""

# ╔═╡ eb7bec9e-5453-11eb-0876-c93037e57988
hypername((n,m,q)) = "HYPER (m=$m, q=$q)"

# ╔═╡ 210d331a-09c3-11eb-06ed-e1f4ddc023fc
md"""
## Load data
"""

# ╔═╡ fd643c2c-08e5-11eb-0d25-014798a505d1
md"""
### Utility functions
"""

# ╔═╡ fe7cf842-08e5-11eb-0d07-41ed2902921a
function loadresults(dir, pattern)
	regmatches = filter(!isnothing, match.(pattern, readdir(dir)))
	pairs = map(regmatches) do regmatch
		design = (n=parse(Int,regmatch[:n]), m=parse(Int,regmatch[:m]), q=parse(Int,regmatch[:q]))
		config = (t=parse(Int,regmatch[:t]), numpos=parse(Int,regmatch[:numpos]), numiters=parse(Int,regmatch[:numiters]), randseed=parse(Int,regmatch[:randseed]))
		result = npzread(joinpath(dir,regmatch.match))
		dictionary([design => dictionary([config=>result])])
	end
	return reduce((d1,d2)->mergewith(merge,d1,d2),pairs)
end

# ╔═╡ c4415a4a-ade7-11eb-16a4-1d071360686f
loadnumstage2(dir; pattern=r"^num_stage_two.n-(?<n>\d+)_m-(?<m>\d+)_q-(?<q>\d+).t-(?<t>\d+).numpos-(?<numpos>\d+).numiters-(?<numiters>\d+).randseed-(?<randseed>\d+).npy$") = loadresults(dir, pattern)

# ╔═╡ 981df61c-ade7-11eb-2524-c77bcf95ebd0
loadsens_perm(dir; pattern=r"^sens_avg_perm.n-(?<n>\d+)_m-(?<m>\d+)_q-(?<q>\d+).t-(?<t>\d+).numpos-(?<numpos>\d+).numiters-(?<numiters>\d+).randseed-(?<randseed>\d+).npy$") = loadresults(dir, pattern)

# ╔═╡ 122151e0-08e6-11eb-18f5-3b277145bc9d
md"""
### HYPER
"""

# ╔═╡ 4d02357e-ade9-11eb-3065-05222f3dfe21
hyp_numstage2 = loadnumstage2("designs/hypergraph")

# ╔═╡ 497b1ba0-ade9-11eb-39ac-3774907d79dd
hyp_sens_perm = loadsens_perm("designs/hypergraph")

# ╔═╡ 2e0cc554-08e6-11eb-0186-bba8c2aee7ec
md"""
### Lexicographic pooling
"""

# ╔═╡ dbc88970-ade9-11eb-2f5d-37beac0ae0de
lex_numstage2 = loadnumstage2("designs/lexicographic")

# ╔═╡ def24668-ade9-11eb-0e67-47bd5c7f64d7
lex_sens_perm = loadsens_perm("designs/lexicographic")

# ╔═╡ 0e4bfb8a-a21d-11eb-07d9-51acf99015b3
md"""
### Consecutive pooling
"""

# ╔═╡ fbb979f6-ade9-11eb-078e-cd247e41d9d0
con_numstage2 = loadnumstage2("designs/consecutive")

# ╔═╡ 03fae046-adea-11eb-1418-fbfb4abf42e0
con_sens_perm = loadsens_perm("designs/consecutive")

# ╔═╡ 1b3b908a-a21d-11eb-1ed1-318a1a5fd572
md"""
### Random assignment design
"""

# ╔═╡ b655856e-aeaa-11eb-2067-9140a3dd8094
rad_seeds = 0:2

# ╔═╡ 203f9d60-a21d-11eb-3d17-47daf91ce8db
rad_numstage2 = map(dictionary(rad_seeds)) do seed
	loadnumstage2("designs/rand,assign.seed-$seed")
end

# ╔═╡ 1ee354b0-adea-11eb-2226-8f78a0301e1c
rad_sens_perm = map(dictionary(rad_seeds)) do seed
	loadsens_perm("designs/rand,assign.seed-$seed")
end

# ╔═╡ 13866a4a-a21d-11eb-18f0-19efe13c34d1
md"""
### Double-pooling
"""

# ╔═╡ 0f4a798e-aead-11eb-05c9-db633b438e4b
dbl_seeds = 0:2

# ╔═╡ 2bb38bce-adea-11eb-3f9e-a1a101f25278
dbl_numstage2 = map(dictionary(dbl_seeds)) do seed
	loadnumstage2("designs/double,pooling.seed-$seed")
end

# ╔═╡ 2b5771e0-adea-11eb-32cd-99191a065b00
dbl_sens_perm = map(dictionary(dbl_seeds)) do seed
	loadsens_perm("designs/double,pooling.seed-$seed")
end

# ╔═╡ b25dd0a0-aeb0-11eb-07b6-2bcfce3c132d
md"""
## Collect results from designs
"""

# ╔═╡ 480a805a-aeaf-11eb-24fc-b7759b08b9a8
viz_numstage2 = dictionary([
	"Consecutive"    => con_numstage2,
	"Lexicographic"  => lex_numstage2,
	["Random assign. (draw $(i+1))" => rad_numstage2[i] for i in rad_seeds]...,
	["Double-pooling (draw $(i+1))" => dbl_numstage2[i] for i in dbl_seeds]...,
	"HYPER"          => hyp_numstage2,
])

# ╔═╡ 380150da-aeaf-11eb-0665-23d1cd1fede3
viz_sens_perm = dictionary([
	"Consecutive"    => con_sens_perm,
	"Lexicographic"  => lex_sens_perm,
	["Random assign. (draw $(i+1))" => rad_sens_perm[i] for i in rad_seeds]...,
	["Double-pooling (draw $(i+1))" => dbl_sens_perm[i] for i in dbl_seeds]...,
	"HYPER"          => hyp_sens_perm,
])

# ╔═╡ 0fd98c50-2203-11eb-0f58-d946dcf484b4
md"""
## Figures
"""

# ╔═╡ ca673ac2-adea-11eb-0289-07b30ddadda9
DES = (n=96,m=16,q=2)

# ╔═╡ 465ff718-aeb3-11eb-0a6b-dd01373ff488
CONFIG = (t=80,numpos=1,numiters=100000,randseed=1)

# ╔═╡ a2fd7cae-b7e9-11eb-32d9-21a05872a518
md"""
### Figs a-b
"""

# ╔═╡ f999303e-d7c0-11eb-034e-231a56f3e828
eff(numstage2) = DES.n/(DES.m+numstage2)

# ╔═╡ 99742f0c-b7e9-11eb-32a4-bd2e2de185e7
let
	plts = map(collect(pairs(viz_numstage2))) do (label,numstage2)
		plot(numstage2[DES][CONFIG] .|> eff;
			label="",title=replace(label," ("=>"\n("),titlefontsize=8,
			ylims=(3.7,5.8),yticks=4:0.5:5.5,
			xlims=(1,DES.n),xticks=[1; 32:32:96],
			xguidefontsize=7,
			xlabel="location of\npositive individual",ylabel="individuals/test",
			linewidth=1.5,
		)
	end
	plot!.(plts[2:end],ylabel="",yticks=(4:0.5:5.5,""),left_margin=-4mm)
	annotate!(plts[1],-50,6.21,text("a",14),left_margin=1.5mm)
	plot(plts...,layout=(1,9),size=(900,170),
		bottom_margin=5mm,top_margin=2mm,right_margin=1mm)
	savefig("fig-a.png"); plot!()
end

# ╔═╡ 9d432a66-b7e9-11eb-131f-1507c7391ef5
let
	numstage2list = map(viz_numstage2) do numstage2
		numstage2[DES][CONFIG]
	end
	boxplot(
		reduce(vcat,fill.(1:length(numstage2list),collect(length.(numstage2list)))),
		reduce(vcat,numstage2list) .|> eff;
		label=nothing,
		xlims=(0.5,9.5),
		xticks=(1:length(numstage2list),replace.(keys(numstage2list)," ("=>"\n(")),
		size=(900,170),right_margin=-4mm,bottom_margin=5mm,
		fillopacity=0.5,ylims=(3.7,5.8),
		ylabel="individuals/test",
		bar_width=0.5,
	)
	annotate!(0.060,6,text("b",14),top_margin=3mm,left_margin=1.5mm)
	savefig("fig-b.png"); plot!()
end

# ╔═╡ aaea1eba-b7e9-11eb-39d0-a5b9fd85c26b
md"""
### Figs c-d
"""

# ╔═╡ b0b8bb60-b7e9-11eb-3bb4-69c931bea667
let
	plts = map(collect(pairs(viz_sens_perm))) do (label,sens_perm)
		plot(sens_perm[DES][CONFIG];
			label="",title=replace(label," ("=>"\n("),titlefontsize=8,
			ylims=(0.72,0.76),yticks=0.72:0.01:0.76,
			xlims=(1,DES.n),xticks=[1; 32:32:96],
			xguidefontsize=7,
			xlabel="location of\npositive individual",ylabel="sensitivity",
			linewidth=1.5,
		)
	end
	plot!.(plts[2:end],ylabel="",yticks=(0.72:0.01:0.76,""),left_margin=-4mm)
	annotate!(plts[1],-58,0.77,text("c",14))
	plot(plts...,layout=(1,9),size=(900,170),
		bottom_margin=5mm,top_margin=3mm,right_margin=1mm)
	savefig("fig-c.png"); plot!()
end

# ╔═╡ eb6bae70-b7e9-11eb-22d7-af685cb504fb
let
	sens_permlist = map(viz_sens_perm) do sens_perm
		sens_perm[DES][CONFIG]
	end
	boxplot(
		reduce(vcat,fill.(1:length(sens_permlist),collect(length.(sens_permlist)))),
		reduce(vcat,sens_permlist);
		label=nothing,
		xlims=(0.5,9.5),
		xticks=(1:length(sens_permlist),replace.(keys(viz_sens_perm)," ("=>"\n(")),
		size=(900,170),right_margin=-4mm,bottom_margin=5mm,
		fillopacity=0.5,ylims=(0.72,0.76),
		ylabel="sensitivity",
		bar_width=0.5,
	)
	annotate!(-0.01,0.763,text("d",14),top_margin=3mm)
	savefig("fig-d.png"); plot!()
end

# ╔═╡ Cell order:
# ╟─849adbca-08e5-11eb-1956-cd02f2a6be7a
# ╠═95eaf04a-08e5-11eb-043f-27646f1bf207
# ╟─eb7bec9e-5453-11eb-0876-c93037e57988
# ╟─210d331a-09c3-11eb-06ed-e1f4ddc023fc
# ╟─fd643c2c-08e5-11eb-0d25-014798a505d1
# ╟─fe7cf842-08e5-11eb-0d07-41ed2902921a
# ╟─c4415a4a-ade7-11eb-16a4-1d071360686f
# ╟─981df61c-ade7-11eb-2524-c77bcf95ebd0
# ╟─122151e0-08e6-11eb-18f5-3b277145bc9d
# ╟─4d02357e-ade9-11eb-3065-05222f3dfe21
# ╟─497b1ba0-ade9-11eb-39ac-3774907d79dd
# ╟─2e0cc554-08e6-11eb-0186-bba8c2aee7ec
# ╟─dbc88970-ade9-11eb-2f5d-37beac0ae0de
# ╟─def24668-ade9-11eb-0e67-47bd5c7f64d7
# ╟─0e4bfb8a-a21d-11eb-07d9-51acf99015b3
# ╟─fbb979f6-ade9-11eb-078e-cd247e41d9d0
# ╟─03fae046-adea-11eb-1418-fbfb4abf42e0
# ╟─1b3b908a-a21d-11eb-1ed1-318a1a5fd572
# ╠═b655856e-aeaa-11eb-2067-9140a3dd8094
# ╟─203f9d60-a21d-11eb-3d17-47daf91ce8db
# ╟─1ee354b0-adea-11eb-2226-8f78a0301e1c
# ╟─13866a4a-a21d-11eb-18f0-19efe13c34d1
# ╠═0f4a798e-aead-11eb-05c9-db633b438e4b
# ╟─2bb38bce-adea-11eb-3f9e-a1a101f25278
# ╟─2b5771e0-adea-11eb-32cd-99191a065b00
# ╟─b25dd0a0-aeb0-11eb-07b6-2bcfce3c132d
# ╟─480a805a-aeaf-11eb-24fc-b7759b08b9a8
# ╟─380150da-aeaf-11eb-0665-23d1cd1fede3
# ╟─0fd98c50-2203-11eb-0f58-d946dcf484b4
# ╠═ca673ac2-adea-11eb-0289-07b30ddadda9
# ╠═465ff718-aeb3-11eb-0a6b-dd01373ff488
# ╟─a2fd7cae-b7e9-11eb-32d9-21a05872a518
# ╟─f999303e-d7c0-11eb-034e-231a56f3e828
# ╟─99742f0c-b7e9-11eb-32a4-bd2e2de185e7
# ╟─9d432a66-b7e9-11eb-131f-1507c7391ef5
# ╟─aaea1eba-b7e9-11eb-39d0-a5b9fd85c26b
# ╟─b0b8bb60-b7e9-11eb-3bb4-69c931bea667
# ╟─eb6bae70-b7e9-11eb-22d7-af685cb504fb


text:
using  jInv.Mesh
using DelimitedFiles
using  jInv.Utils
using jInvSeismic.FWI
using jInvSeismic.Utils

function readModelAndGenerateMeshMrefOverthrust(readModelFolder::String,modelFilename::String,dim::Int64,pad::Int64,domain::Vector{Float64},newSize::Vector=[],velBottom::Float64=1.75,velHigh::Float64=2.9)
########################## m,mref are in Velocity here. ###################################

if dim==2
	# SEGmodel2Deasy.dat
	m = readdlm(string(readModelFolder,"/",modelFilename));
	m = m*1e-3;
	m = Matrix(m');
	mref = getSimilarLinearModel(m,velBottom,velHigh);
else
	# 3D SEG slowness model
	file = matopen(string(readModelFolder,"/",modelFilename)); DICT = read(file); close(file);
	m = DICT["VELs"];
	m = m*1e-3;
	mref = getSimilarLinearModel(m,velBottom,velHigh);
end

sea_level = 1.5;
sea = abs.(m[:] .- 1.5) .< 1e-3;
mref[sea] = m[sea];
if newSize!=[]
	m    = expandModelNearest(m,   collect(size(m)),newSize);
	mref = expandModelNearest(mref,collect(size(mref)),newSize);
end

Minv = getRegularMesh(domain,collect(size(m)));

(mPadded,MinvPadded) = addAbsorbingLayer(m,Minv,pad);
(mrefPadded,MinvPadded) = addAbsorbingLayer(mref,Minv,pad);

mPaddedNoSalt = copy(mPadded);
mPaddedNoSalt[mPaddedNoSalt .> 5.5] .= 5.5;

sea_level = 1.5;
sea = abs.(mPadded[:] .- 1.5) .< 1e-2;

file = matopen(string("examples/mrefOverthrust3D.mat")); DICT = read(file); close(file);
mrefPadded = DICT["VELs"];

N = prod(MinvPadded.n);
boundsLow  = minimum(mPadded);
boundsHigh = maximum(mPadded);

boundsLow  = ones(N)*boundsLow;
boundsLow = convert(Array{Float32},boundsLow);
boundsHigh = ones(N)*boundsHigh;
boundsHigh = convert(Array{Float32},boundsHigh);

return (mPadded,MinvPadded,mrefPadded,boundsHigh,boundsLow);
end

dim     = 3;
pad     = 16;
jumpSrc    = 10;
jumpRcv = 4;
offset  = 1000;
domain = [0.0,7.5,0.0,7.5,0.0,4.65]; # without the pad for Marmousi 1
newSize = [172,172,108];
modelDir = pwd();

(m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMrefOverthrust(modelDir,"examples/overthrust3D_small.mat",dim,pad,domain,newSize,1.25,4.0);


text:
function available(interpreter::Interpreter, domain::Domain, state::State)
    # Ground all action definitions with arguments
    actions = Term[]
    for act in values(get_actions(domain))
        act_name = get_name(act)
        act_vars, act_types = get_argvars(act), get_argtypes(act)
        # Directly check precondition if action has no parameters
        if isempty(act_vars) && satisfy(domain, state, get_precond(act))
            push!(actions, Const(act_name))
            continue
        end
        # Include type conditions when necessary for correctness
        typecond = (@julog($ty(:v)) for (v, ty) in zip(act_vars, act_types))
        p = get_precond(act)
        if has_func(p, domain) || has_derived(p, domain) || has_quantifier(p)
            conds = prepend!(flatten_conjs(p), typecond)
        elseif get_requirements(domain)[:typing]
            conds = append!(flatten_conjs(p), typecond)
        else
            conds = flatten_conjs(p)
        end
        # Find all substitutions that satisfy preconditions
        subst = satisfiers(interpreter, domain, state, conds)
        if isempty(subst) continue end
        for s in subst
            args = [s[v] for v in act_vars if v in keys(s)]
            if any(!is_ground(a) for a in args) continue end
            term = isempty(args) ? Const(act_name) : Compound(act_name, args)
            push!(actions, term)
        end
    end
    return actions
end

function available(interpreter::Interpreter,
                   domain::Domain, state::State, action::Action, args)
    if any(!is_ground(a) for a in args)
       error("Not all arguments are ground.")
    end
    act_vars, act_types = get_argvars(action), get_argtypes(action)
    subst = Subst(var => val for (var, val) in zip(act_vars, args))
    # Construct type conditions of the form "type(val)"
    typecond = (all(ty == :object for ty in action.types) ? Term[] :
               [@julog($ty(:v)) for (v, ty) in zip(args, act_types)])
    # Check whether preconditions hold
    precond = substitute(get_precond(action), subst)
    conds = has_func(precond, domain) || has_quantifier(precond) ?
        [typecond; precond] : [precond; typecond]
    return satisfy(interpreter, domain, state, conds)
end


text:
module ARCHER
using ..AGFFileReader: Glass, GlassID, AGF
export S_LAL12M, S_BSL7M, S_LAL13M, S_NSL3M, S_BAL41M, S_LAM60M, S_BAL3M, L_BAL35M, S_LAH60M, S_TIM2M, S_FSL5M, L_LAM69M, K_VC79M, S_BAL42M, S_LAH53M

const S_LAL12M = Glass(GlassID(AGF, 914), 1, 2.753166, -0.01207074, 0.01625163, 0.001041721, -0.0001156154, 6.773623e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, -7.683e-7, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0047, 4.2, -1.0, 7.2, 2.0, 0, 53.0, [(0.27, 0.0, 10.0), (0.28, 0.06, 10.0), (0.29, 0.15, 10.0), (0.3, 0.29, 10.0), (0.31, 0.45, 10.0), (0.32, 0.61, 10.0), (0.334, 0.774, 10.0), (0.35, 0.89, 10.0), (0.365, 0.952, 10.0), (0.37, 0.962, 10.0), (0.38, 0.976, 10.0), (0.39, 0.984, 10.0), (0.4, 0.988, 10.0), (0.404, 0.989, 10.0), (0.42, 0.992, 10.0), (0.436, 0.994, 10.0), (0.46, 0.995, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.996, 10.0), (1.529, 0.991, 10.0), (1.97, 0.966, 10.0), (2.325, 0.809, 10.0)], 1.674061, -1.0, -1.0, 0, 55.017718, 0, 4.01, 0)
const S_BSL7M = Glass(GlassID(AGF, 915), 1, 2.262505, -0.01037921, 0.009341724, 0.0006453071, -7.463858e-5, 4.096215e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 3.095e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0024, 2.0, -1.0, 7.2, 1.5, 0, 1.0, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.08, 10.0), (0.3, 0.31, 10.0), (0.31, 0.58, 10.0), (0.32, 0.77, 10.0), (0.334, 0.909, 10.0), (0.35, 0.968, 10.0), (0.365, 0.988, 10.0), (0.37, 0.991, 10.0), (0.38, 0.991, 10.0), (0.39, 0.996, 10.0), (0.4, 0.997, 10.0), (0.404, 0.997, 10.0), (0.42, 0.996, 10.0), (0.436, 0.995, 10.0), (0.46, 0.995, 10.0), (0.5, 0.996, 10.0), (0.546, 0.998, 10.0), (0.58, 0.997, 10.0), (0.62, 0.997, 10.0), (0.66, 0.997, 10.0), (0.7, 0.998, 10.0), (1.06, 0.996, 10.0), (1.529, 0.989, 10.0), (1.97, 0.965, 10.0), (2.325, 0.862, 10.0)], 1.513232, -1.0, -1.0, 0, 63.709676, 0, 2.52, 0)
const S_LAL13M = Glass(GlassID(AGF, 916), 1, 2.800915, -0.0140473, 0.01703856, 0.001170426, -0.0001225957, 6.868454e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 7.171e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0081, 3.0, -1.0, 5.7, 1.5, 0, 52.0, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.01, 10.0), (0.3, 0.03, 10.0), (0.31, 0.08, 10.0), (0.32, 0.19, 10.0), (0.334, 0.411, 10.0), (0.35, 0.68, 10.0), (0.365, 0.844, 10.0), (0.37, 0.88, 10.0), (0.38, 0.932, 10.0), (0.39, 0.958, 10.0), (0.4, 0.972, 10.0), (0.404, 0.977, 10.0), (0.42, 0.986, 10.0), (0.436, 0.989, 10.0), (0.46, 0.993, 10.0), (0.5, 0.996, 10.0), (0.546, 0.997, 10.0), (0.58, 0.996, 10.0), (0.62, 0.995, 10.0), (0.66, 0.995, 10.0), (0.7, 0.996, 10.0), (1.06, 0.995, 10.0), (1.529, 0.99, 10.0), (1.97, 0.962, 10.0), (2.325, 0.762, 10.0)], 1.689007, -1.0, -1.0, 0, 52.782035, 0, 3.6, 0)
const S_NSL3M = Glass(GlassID(AGF, 917), 1, 2.267846, -0.008623029, 0.01087772, 0.0005394286, -5.529385e-5, 3.416585e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, -1.095e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0005, 1.0, -1.0, 9.0, 1.0, 0, 1.0, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.15, 10.0), (0.334, 0.655, 10.0), (0.35, 0.924, 10.0), (0.365, 0.978, 10.0), (0.37, 0.984, 10.0), (0.38, 0.99, 10.0), (0.39, 0.995, 10.0), (0.4, 0.997, 10.0), (0.404, 0.997, 10.0), (0.42, 0.997, 10.0), (0.436, 0.997, 10.0), (0.46, 0.997, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.999, 10.0), (0.62, 0.999, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.992, 10.0), (1.97, 0.937, 10.0), (2.325, 0.827, 10.0)], 1.516509, -1.0, -1.0, 0, 58.745047, 0, 2.48, 0)
const S_BAL41M = Glass(GlassID(AGF, 918), 1, 2.402809, -0.01041283, 0.01136728, 0.0007464327, -8.517315e-5, 4.764139e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 2.775e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0031, 3.0, -1.0, 7.5, 2.5, 0, 51.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.05, 10.0), (0.31, 0.25, 10.0), (0.32, 0.51, 10.0), (0.334, 0.782, 10.0), (0.35, 0.925, 10.0), (0.365, 0.971, 10.0), (0.37, 0.978, 10.0), (0.38, 0.985, 10.0), (0.39, 0.99, 10.0), (0.4, 0.993, 10.0), (0.404, 0.993, 10.0), (0.42, 0.994, 10.0), (0.436, 0.994, 10.0), (0.46, 0.995, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.997, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.992, 10.0), (1.97, 0.974, 10.0), (2.325, 0.861, 10.0)], 1.560983, -1.0, -1.0, 0, 60.308093, 0, 2.78, 0)
const S_LAM60M = Glass(GlassID(AGF, 919), 1, 2.958448, -0.01542036, 0.02030973, 0.001302646, -0.0001238708, 7.401598e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 8.217e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0085, 2.0, -1.0, 5.4, 2.0, 0, 52.0, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.06, 10.0), (0.31, 0.27, 10.0), (0.32, 0.53, 10.0), (0.334, 0.786, 10.0), (0.35, 0.922, 10.0), (0.365, 0.967, 10.0), (0.37, 0.975, 10.0), (0.38, 0.984, 10.0), (0.39, 0.989, 10.0), (0.4, 0.992, 10.0), (0.404, 0.992, 10.0), (0.42, 0.993, 10.0), (0.436, 0.993, 10.0), (0.46, 0.995, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.993, 10.0), (1.97, 0.98, 10.0), (2.325, 0.867, 10.0)], 1.737928, -1.0, -1.0, 0, 49.027571, 0, 4.06, 0)
const S_BAL3M = Glass(GlassID(AGF, 920), 1, 2.417473, -0.008685888, 0.01396835, 0.0006180845, -5.274288e-5, 3.679696e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, -3.881e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0005, 1.0, -1.0, 9.5, 2.0, 0, 1.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.334, 0.096, 10.0), (0.35, 0.63, 10.0), (0.365, 0.894, 10.0), (0.37, 0.928, 10.0), (0.38, 0.963, 10.0), (0.39, 0.979, 10.0), (0.4, 0.988, 10.0), (0.404, 0.99, 10.0), (0.42, 0.994, 10.0), (0.436, 0.994, 10.0), (0.46, 0.995, 10.0), (0.5, 0.997, 10.0), (0.546, 0.998, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.991, 10.0), (1.97, 0.951, 10.0), (2.325, 0.847, 10.0)], 1.568151, -1.0, -1.0, 0, 52.737315, 0, 2.98, 0)
const L_BAL35M = Glass(GlassID(AGF, 921), 1, 2.481845, -0.01213254, 0.01179439, 0.000804573, -9.005814e-5, 4.898641e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 6.78e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0043, 1.0, -1.0, 6.6, 2.5, 0, 4.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.06, 10.0), (0.31, 0.27, 10.0), (0.32, 0.53, 10.0), (0.334, 0.786, 10.0), (0.35, 0.922, 10.0), (0.365, 0.967, 10.0), (0.37, 0.975, 10.0), (0.38, 0.984, 10.0), (0.39, 0.989, 10.0), (0.4, 0.992, 10.0), (0.404, 0.992, 10.0), (0.42, 0.993, 10.0), (0.436, 0.993, 10.0), (0.46, 0.995, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.993, 10.0), (1.97, 0.98, 10.0), (2.325, 0.867, 10.0)], 1.58642, -1.0, -1.0, 0, 60.713872, 0, 2.82, 0)
const S_LAH60M = Glass(GlassID(AGF, 922), 1, 3.244481, -0.01528454, 0.03201574, 0.001865793, -0.000141547, 1.246734e-5, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 9.635e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0037, 1.0, -1.0, 5.6, 1.0, 0, 4.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.334, 0.0, 10.0), (0.35, 0.27, 10.0), (0.365, 0.643, 10.0), (0.37, 0.72, 10.0), (0.38, 0.83, 10.0), (0.39, 0.88, 10.0), (0.4, 0.924, 10.0), (0.404, 0.938, 10.0), (0.42, 0.957, 10.0), (0.436, 0.97, 10.0), (0.46, 0.98, 10.0), (0.5, 0.99, 10.0), (0.546, 0.996, 10.0), (0.58, 0.997, 10.0), (0.62, 0.997, 10.0), (0.66, 0.997, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.992, 10.0), (1.97, 0.968, 10.0), (2.325, 0.808, 10.0)], 1.828944, -1.0, -1.0, 0, 37.007308, 0, 4.43, 0)
const S_TIM2M = Glass(GlassID(AGF, 923), 1, 2.535373, -0.01040568, 0.02209869, 0.0009396727, -6.499788e-5, 9.894408e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 3.044e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.0051, 1.0, -1.0, 8.1, 2.0, 0, 1.0, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.334, 0.0, 10.0), (0.35, 0.0, 10.0), (0.365, 0.267, 10.0), (0.37, 0.44, 10.0), (0.38, 0.73, 10.0), (0.39, 0.87, 10.0), (0.4, 0.942, 10.0), (0.404, 0.955, 10.0), (0.42, 0.978, 10.0), (0.436, 0.986, 10.0), (0.46, 0.99, 10.0), (0.5, 0.994, 10.0), (0.546, 0.997, 10.0), (0.58, 0.997, 10.0), (0.62, 0.996, 10.0), (0.66, 0.996, 10.0), (0.7, 0.997, 10.0), (1.06, 0.999, 10.0), (1.529, 0.996, 10.0), (1.97, 0.973, 10.0), (2.325, 0.917, 10.0)], 1.613318, -1.0, -1.0, 0, 36.219537, 0, 2.69, 0)
const S_FSL5M = Glass(GlassID(AGF, 924), 1, 2.181284, -0.009793636, 0.007745517, 0.0005928784, -7.340057e-5, 4.012184e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, -5.246e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.0022, 2.0, -1.0, 9.0, 1.5, 0, 3.0, [(0.27, 0.0, 10.0), (0.28, 0.13, 10.0), (0.29, 0.43, 10.0), (0.3, 0.73, 10.0), (0.31, 0.89, 10.0), (0.32, 0.961, 10.0), (0.334, 0.988, 10.0), (0.35, 0.995, 10.0), (0.365, 0.997, 10.0), (0.37, 0.998, 10.0), (0.38, 0.998, 10.0), (0.39, 0.999, 10.0), (0.4, 0.999, 10.0), (0.404, 0.999, 10.0), (0.42, 0.999, 10.0), (0.436, 0.999, 10.0), (0.46, 0.999, 10.0), (0.5, 0.999, 10.0), (0.546, 0.999, 10.0), (0.58, 0.999, 10.0), (0.62, 0.999, 10.0), (0.66, 0.999, 10.0), (0.7, 0.999, 10.0), (1.06, 0.998, 10.0), (1.529, 0.989, 10.0), (1.97, 0.974, 10.0), (2.325, 0.864, 10.0)], 1.484524, -1.0, -1.0, 0, 69.75631, 0, 2.46, 0)
const L_LAM69M = Glass(GlassID(AGF, 925), 1, 2.909506, -0.01373872, 0.02464585, 0.001285182, -9.119209e-5, 8.193186e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 3.17e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0031, 3.1, -1.0, 8.6, 2.0, 0, 52.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.334, 0.0, 10.0), (0.35, 0.22, 10.0), (0.365, 0.612, 10.0), (0.37, 0.7, 10.0), (0.38, 0.81, 10.0), (0.39, 0.87, 10.0), (0.4, 0.912, 10.0), (0.404, 0.925, 10.0), (0.42, 0.95, 10.0), (0.436, 0.962, 10.0), (0.46, 0.974, 10.0), (0.5, 0.989, 10.0), (0.546, 0.996, 10.0), (0.58, 0.995, 10.0), (0.62, 0.994, 10.0), (0.66, 0.996, 10.0), (0.7, 0.998, 10.0), (1.06, 0.0, 10.0), (1.529, 0.0, 10.0), (1.97, 0.0, 10.0), (2.325, 0.0, 10.0)], 1.727801, -1.0, -1.0, 0, 40.364781, 0, 3.24, 0)
const K_VC79M = Glass(GlassID(AGF, 926), 1, 2.537275, -0.01000404, 0.01614788, -8.966925e-5, 4.332545e-5, -1.820481e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0076, -1.0, -1.0, 9.6, 2.0, 0, -1.0, [(0.27, 0.022, 10.0), (0.28, 0.103, 10.0), (0.29, 0.26, 10.0), (0.3, 0.378, 10.0), (0.31, 0.669, 10.0), (0.32, 0.799, 10.0), (0.334, 0.891, 10.0), (0.35, 0.964, 10.0), (0.365, 0.978, 10.0), (0.37, 0.982, 10.0), (0.38, 0.985, 10.0), (0.39, 0.985, 10.0), (0.4, 0.985, 10.0), (0.404, 0.985, 10.0), (0.42, 0.985, 10.0), (0.436, 0.988, 10.0), (0.46, 0.99, 10.0), (0.5, 0.992, 10.0), (0.546, 0.994, 10.0), (0.58, 0.994, 10.0), (0.62, 0.994, 10.0), (0.66, 0.997, 10.0), (0.7, 0.997, 10.0), (1.06, 0.998, 10.0), (1.529, 0.997, 10.0), (1.97, 0.978, 10.0), (2.325, 0.0, 10.0)], 1.606477, -1.0, -1.0, 0, 57.60776, 0, 3.09, 0)
const S_BAL42M = Glass(GlassID(AGF, 927), 1, 2.460903, -0.01005747, 0.01237562, 0.0007373934, -7.773189e-5, 4.239549e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 4.545e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.002, 1.0, -1.0, 6.6, 1.5, 0, 1.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.03, 10.0), (0.3, 0.15, 10.0), (0.31, 0.36, 10.0), (0.32, 0.58, 10.0), (0.334, 0.8, 10.0), (0.35, 0.932, 10.0), (0.365, 0.973, 10.0), (0.37, 0.979, 10.0), (0.38, 0.986, 10.0), (0.39, 0.991, 10.0), (0.4, 0.993, 10.0), (0.404, 0.994, 10.0), (0.42, 0.995, 10.0), (0.436, 0.995, 10.0), (0.46, 0.996, 10.0), (0.5, 0.998, 10.0), (0.546, 0.999, 10.0), (0.58, 0.998, 10.0), (0.62, 0.998, 10.0), (0.66, 0.998, 10.0), (0.7, 0.998, 10.0), (1.06, 0.997, 10.0), (1.529, 0.993, 10.0), (1.97, 0.978, 10.0), (2.325, 0.881, 10.0)], 1.580467, -1.0, -1.0, 0, 59.002542, 0, 3.19, 0)
const S_LAH53M = Glass(GlassID(AGF, 928), 1, 3.156971, -0.01480577, 0.02790549, 0.00158408, -0.00012293, 9.591655e-6, 0.0, 0.0, 0.0, 0.0, 0.365, 2.325, 8.688e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, -0.0052, 2.0, -1.0, 5.9, 1.0, 0, 4.2, [(0.27, 0.0, 10.0), (0.28, 0.0, 10.0), (0.29, 0.0, 10.0), (0.3, 0.0, 10.0), (0.31, 0.0, 10.0), (0.32, 0.0, 10.0), (0.334, 0.0, 10.0), (0.35, 0.38, 10.0), (0.365, 0.741, 10.0), (0.37, 0.8, 10.0), (0.38, 0.88, 10.0), (0.39, 0.925, 10.0), (0.4, 0.947, 10.0), (0.404, 0.954, 10.0), (0.42, 0.969, 10.0), (0.436, 0.977, 10.0), (0.46, 0.985, 10.0), (0.5, 0.993, 10.0), (0.546, 0.996, 10.0), (0.58, 0.996, 10.0), (0.62, 0.996, 10.0), (0.66, 0.997, 10.0), (0.7, 0.998, 10.0), (1.06, 0.999, 10.0), (1.529, 0.996, 10.0), (1.97, 0.973, 10.0), (2.325, 0.813, 10.0)], 1.801019, -1.0, -1.0, 0, 40.721498, 0, 4.43, 0)
end # module

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
mutable struct Tim2DDriver <: DriverModel{LatLonAccel}
    rec::SceneRecord
    mlon::LaneFollowingDriver
    mlat::LateralDriverModel
    mlane::LaneChangeModel

    function Tim2DDriver(
        timestep::Float64;
        mlon::LaneFollowingDriver=IntelligentDriverModel(),
        mlat::LateralDriverModel=ProportionalLaneTracker(),
        mlane::LaneChangeModel=TimLaneChanger(timestep),
        rec::SceneRecord = SceneRecord(1, timestep)
        )

        retval = new()

        retval.rec = rec
        retval.mlon = mlon
        retval.mlat = mlat
        retval.mlane = mlane

        retval
    end
end
get_name(::Tim2DDriver) = "Tim2DDriver"
function set_desired_speed!(model::Tim2DDriver, v_des::Float64)
    set_desired_speed!(model.mlon, v_des)
    set_desired_speed!(model.mlane, v_des)
    model
end
function track_longitudinal!(driver::LaneFollowingDriver, scene::Scene, roadway::Roadway, vehicle_index::Int, fore::NeighborLongitudinalResult)
    v_ego = scene[vehicle_index].state.v
    if fore.ind != nothing
        headway, v_oth = fore.Δs, scene[fore.ind].state.v
    else
        headway, v_oth = NaN, NaN
    end
    return track_longitudinal!(driver, v_ego, v_oth, headway)
end
function observe!(driver::Tim2DDriver, scene::Scene, roadway::Roadway, egoid::Int)

    update!(driver.rec, scene)
    observe!(driver.mlane, scene, roadway, egoid)

    vehicle_index = findfirst(egoid, scene)
    lane_change_action = rand(driver.mlane)
    laneoffset = get_lane_offset(lane_change_action, driver.rec, roadway, vehicle_index)
    lateral_speed = convert(Float64, get(VELFT, driver.rec, roadway, vehicle_index))

    if lane_change_action.dir == DIR_MIDDLE
        fore = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    elseif lane_change_action.dir == DIR_LEFT
        fore = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    else
        @assert(lane_change_action.dir == DIR_RIGHT)
        fore = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    end

    track_lateral!(driver.mlat, laneoffset, lateral_speed)
    track_longitudinal!(driver.mlon, scene, roadway, vehicle_index, fore)

    driver
end
Base.rand(driver::Tim2DDriver) = LatLonAccel(rand(driver.mlat), rand(driver.mlon).a)
Distributions.pdf(driver::Tim2DDriver, a::LatLonAccel) = pdf(driver.mlat, a.a_lat) * pdf(driver.mlon, a.a_lon)
Distributions.logpdf(driver::Tim2DDriver, a::LatLonAccel) = logpdf(driver.mlat, a.a_lat) * logpdf(driver.mlon, a.a_lon)

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
module AudioSchedules

import Base:
    close,
    getindex,
    iterate,
    IteratorEltype,
    IteratorSize,
    show,
    size,
    wait,
    write
using Base: EltypeUnknown, SizeUnknown
using Base.FastMath: sin_fast
using Base.Meta: ParseError
using Base.Threads: nthreads
using DataStructures: SortedSet
using InfiniteArrays: Fill, ∞
using Interpolations: CubicSplineInterpolation
using NLsolve: nlsolve
using PortAudio: Messanger, PortAudioStream, Scribe, write_buffer
import PortAudio: get_input_type, get_output_type
using RegularExpressions: capture, CONSTANTS, of, pattern, raw, short
import SampledSignals: nchannels, SampleBuf, samplerate
using Unitful: dB, Hz, μPa, s

const FREQUENCY = typeof(1.0Hz)
const RATE = typeof(1.0 / s)
const TIME = typeof(0.0s)

include("miscellaneous.jl")
include("series.jl")

export AudioSchedule

const ORCHESTRA = Vector{Tuple{Any, TIME}}
const ACTIVATED = Vector{Bool}
const TRIGGERS = SortedSet{Tuple{TIME, Int}}

struct AudioSchedule
    sample_rate::FREQUENCY
    orchestra::ORCHESTRA
    activated::ACTIVATED
    triggers::TRIGGERS
end

export AudioSchedule

"""
    AudioSchedule(; sample_rate = 44100Hz)

Create an empty `AudioSchedule`.
You can [`add!`](@ref) new synthesizers to the a_schedule.
Specify a `sample_rate` with units per time, like `1/s` or `Hz`.

```jldoctest audio_schedule
julia> using AudioSchedules

julia> using Unitful: s, Hz

julia> a_schedule = AudioSchedule()
0.0 s 44100.0 Hz AudioSchedule

julia> add!(
           a_schedule,
           Map(sin, Cycles(440Hz)),
           0s,
           0,
           Line => 1s,
           0.05,
           Line => 1s,
           0.05,
           Line => 1s,
           0,
       )

julia> add!(
            a_schedule,
            Map(sin, Cycles(660Hz)),
            2s,
            0,
            Line => 1s,
            0.05,
            Line => 1s,
            0.05,
            Line => 1s,
            0,
        )

julia> a_schedule
5.0 s 44100.0 Hz AudioSchedule
```

You can use write the a_schedule directly to a `PortAudioStream`.
However, to do so, the `PortAudioStream` must have a [`Weaver`](@ref) writer.
The `PortAudioStream` must have exactly 0 input channels, 1 output channel, and a matching
sample rate.

```jldoctest audio_schedule
julia> using PortAudio: PortAudioStream

julia> PortAudioStream(0, 1, writer = Weaver(), warn_xruns = false) do stream
           write(stream, a_schedule)
       end

```

You can save an `AudioSchedule` as a `SampledSignals.SampleBuf`.

```jldoctest audio_schedule
julia> using SampledSignals: SampleBuf

julia> saved = SampleBuf(a_schedule)
220500-frame, 1-channel SampleBuf{Float64, 2}
5.0s sampled at 44100.0Hz
▂▂▃▃▃▃▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▃▃▃▃▂▂

julia> PortAudioStream(0, 1, warn_xruns = false) do stream
           write(stream, saved)
       end
220500
```
"""
function AudioSchedule(; sample_rate = 44100Hz)
    AudioSchedule(sample_rate, ORCHESTRA(), ACTIVATED(), TRIGGERS())
end

function IteratorSize(::Type{AudioSchedule})
    SizeUnknown()
end
function IteratorEltype(::Type{AudioSchedule})
    EltypeUnknown()
end

function samplerate(a_schedule::AudioSchedule)
    a_schedule.sample_rate / Hz
end
function nchannels(a_schedule::AudioSchedule)
    1
end

function duration(a_schedule::AudioSchedule)
    triggers = a_schedule.triggers
    if length(triggers) > 0
        last(triggers)[1] - first(triggers)[1]
    else
        0.0s
    end
end

function show(io::IO, a_schedule::AudioSchedule)
    print(io, "$(duration(a_schedule)) $(a_schedule.sample_rate) AudioSchedule")
end

function make_envelope(start_level, (shape, duration), end_level, more_segments...)
    segments(shape, start_level, duration, end_level)...,
    make_envelope(end_level, more_segments...)...
end
function make_envelope(_)
    ()
end

function add_series!(a_schedule::AudioSchedule, series, start_time, duration)
    triggers = a_schedule.triggers
    orchestra = a_schedule.orchestra
    push!(orchestra, (series, start_time))
    push!(a_schedule.activated, false)
    instrument_index = length(orchestra)
    push!(triggers, (start_time, instrument_index))
    push!(triggers, (start_time + duration, instrument_index))
    nothing
end

function second(something)
    something[2]
end

function all_but_last(thing1, thing2, rest...)
    thing1, all_but_last(thing2, rest...)...
end

function all_but_last(_)
    ()
end

function skip(sample_rate, wave, skip_time)
    wave[(round(Int, skip_time * sample_rate)+1):end]
end

function make_wave(sample_rate, wave, shape, segment_start_time, segment_skip_time, duration)
    broadcast(
        *,
        # skip ahead to time in wave, then the one after
        skip(sample_rate, wave, segment_skip_time),
        make_series(shape, sample_rate),
    ),
    segment_start_time,
    duration
end

function triples(sample_rate, synthesizer, start_time, envelope...)
    wave = make_series(synthesizer, sample_rate)
    shapes_durations = make_envelope(envelope...)
    shapes = map(first, shapes_durations)
    durations = map(second, shapes_durations)
    segment_skip_times = all_but_last(cumsum((0.0s, durations...))...)
    segment_start_times = start_time .+ segment_skip_times
    map(
        let sample_rate = sample_rate, wave = wave
            function (shape, segment_start_time, segment_skip_time, duration)
                make_wave(sample_rate, wave, shape, segment_start_time, segment_skip_time, duration)
            end
        end,
        shapes,
        segment_start_times,
        segment_skip_times,
        durations,
    )
end

"""
    add!(a_schedule::AudioSchedule, synthesizer, start_time,
        start_level, shape => duration, end_level, more_segments...
    )

Add a synthesizer to a [`AudioSchedule`](@ref), where `synthesizer` is anything that supports [`make_series`](@ref),
`start_time` has units of time (like `s`), and the rest of the arguments specify the shape of the envelope.

For all envelope [`segments`](@ref), call

```
segments(shape, start_level, duration, end_level)
```

`duration` should have units of time (like `s`). For example,

```
add!(a_schedule, synthesizer, start_time, 0, Line => 1s, 1, Line => 1s, 0)
```

will call segments twice:

```
segments(Line, 0, 1s, 1)
segments(Line, 1, 1s, 0)
```
"""
function add!(a_schedule::AudioSchedule, synthesizer, start_time, piece_1, rest...)
    sample_rate = a_schedule.sample_rate
    map(
        let a_schedule = a_schedule
            function ((shaped, time, duration),)
                add_series!(a_schedule, shaped, time, duration)
            end
        end,
        triples(sample_rate, synthesizer, start_time, piece_1, rest...),
    )
    nothing
end

export add!

function iterate(a_schedule::AudioSchedule)
    iteration = iterate(a_schedule.triggers)
    # Done if no more triggers
    if iteration === nothing
        nothing
    else
        # Start at 0s
        trigger_iterate(a_schedule, 0.0s, iteration...)
    end
end

function iterate(a_schedule::AudioSchedule, (schedule_time, trigger_state))
    iteration = iterate(a_schedule.triggers, trigger_state)
    if iteration === nothing
        nothing
    else
        trigger_iterate(a_schedule, schedule_time, iteration...)
    end
end

function trigger_iterate(
    a_schedule,
    schedule_time,
    (trigger_time, instrument_index),
    trigger_state,
)
    activated = a_schedule.activated
    triggers = a_schedule.triggers
    # pull the rest of the triggers at the current time
    while schedule_time == trigger_time
        activated[instrument_index] = !activated[instrument_index]
        iteration = iterate(triggers, trigger_state)
        if iteration === nothing
            return nothing
        end
        (trigger_time, instrument_index), trigger_state = iteration
    end
    # adjust these based on how long ago they went on
    active = view(a_schedule.orchestra, activated)
    result = (
        if isempty(active)
            Fill(0.0, ∞)
        else
            broadcast(+, Iterators.map(
                let sample_rate = a_schedule.sample_rate, schedule_time = schedule_time
                    function ((instrument, instrument_start_time),)
                        skip(sample_rate, instrument, schedule_time - instrument_start_time)
                    end
                end,
                active
            )...)
        end,
        round(Int, (trigger_time - schedule_time) * a_schedule.sample_rate),
    ),
    (trigger_time, trigger_state)
    # pull the next trigger
    activated[instrument_index] = !activated[instrument_index]
    result
end

function trigger_iterate(_, __, ::Nothing)
    nothing
end

# TODO: make a SampledSignal without sacrificing performance
# TODO: figure out a way to avoid peaking !!!

"""
    struct Weaver <: Scribe end

A special `PortAudio` `Scribe`.
You can use a `Weaver` to write an [`AudioSchedule`](@ref) directly to your speakers.
"""
struct Weaver <: Scribe end

get_input_type(::Weaver, _) = Nothing
get_output_type(::Weaver, _) = Nothing

export Weaver

struct DividedRange <: AbstractVector{UnitRange{Int}}
    previous::Int
    share::Float64
    divisions::Int
end

function size(divided_range::DividedRange)
    (divided_range.divisions,)
end

function getindex(divided_range::DividedRange, index)
    share = divided_range.share
    divided_range.previous .+
    ((round(Int, (index - 1) * share)+1):round(Int, index * share))
end

function add_to(divided_range::DividedRange, increment)
    DividedRange(
        divided_range.previous + increment,
        divided_range.share,
        divided_range.divisions,
    )
end

function DividedRange(a_range::UnitRange, divisions)
    DividedRange(first(a_range) - 1, length(a_range) / divisions, divisions)
end

struct DualSubRanges <: AbstractVector{Tuple{UnitRange{Int},UnitRange{Int}}}
    frames_per_buffer::Int
    the_length::Int
    left_in_buffer::Int
    first_item::Tuple{DividedRange,DividedRange,Bool}
    last_item::Tuple{DividedRange,DividedRange,Bool}
    full_buffer_ranges::DividedRange
end

function DualSubRanges(
    first_buffer_at,
    frames_per_buffer,
    series_total,
    number_of_subsections,
)
    left_in_buffer = frames_per_buffer - first_buffer_at
    middle_buffers, last_buffer_at =
        fldmod(series_total - left_in_buffer, frames_per_buffer)
    is_last_complete = last_buffer_at == 0
    if is_last_complete
        middle_buffers = middle_buffers - 1
        last_buffer_at = frames_per_buffer
    end
    if series_total <= left_in_buffer
        first_series_sub_ranges = DividedRange(1:series_total, number_of_subsections)
        is_first_complete = is_last_complete
    else
        first_series_sub_ranges = DividedRange(1:left_in_buffer, number_of_subsections)
        is_first_complete = true
    end
    last_buffer_sub_ranges = DividedRange(1:last_buffer_at, number_of_subsections)
    DualSubRanges(
        frames_per_buffer,
        middle_buffers + 2,
        left_in_buffer,
        (
            first_series_sub_ranges,
            add_to(first_series_sub_ranges, first_buffer_at),
            is_first_complete,
        ),
        (
            add_to(
                last_buffer_sub_ranges,
                left_in_buffer + middle_buffers * frames_per_buffer,
            ),
            last_buffer_sub_ranges,
            is_last_complete,
        ),
        DividedRange(1:frames_per_buffer, number_of_subsections),
    )
end

function size(dual_sub_ranges::DualSubRanges)
    (dual_sub_ranges.the_length,)
end

function getindex(dual_sub_ranges::DualSubRanges, index::Int)
    the_length = dual_sub_ranges.the_length
    @boundscheck if index < 0 || index > the_length
        throw(BoundsError(dual_sub_ranges, index))
    end
    if index == 1
        dual_sub_ranges.first_item
    elseif index == the_length
        dual_sub_ranges.last_item
    else
        full_buffer_ranges = dual_sub_ranges.full_buffer_ranges
        (
            add_to(
                full_buffer_ranges,
                dual_sub_ranges.left_in_buffer +
                (index - 2) * dual_sub_ranges.frames_per_buffer,
            ),
            full_buffer_ranges,
            true,
        )
    end
end

struct TaskIO
    task::Task
    series_buffers_channel::Channel{Tuple{Any,Int}}
    sub_ranges_channel::Channel{Tuple{UnitRange{Int},UnitRange{Int},Bool}}
    output_channel::Channel{Nothing}
end

function close(task_io::TaskIO)
    close(task_io.series_buffers_channel)
    close(task_io.sub_ranges_channel)
    wait(task_io.task)
    close(task_io.output_channel)
end

function fill_buffer_signal!(
    buffer,
    series,
    sub_ranges_channel,
    output_channel,
)
    (series_sub_range, buffer_sub_range, is_complete) = take!(sub_ranges_channel)
    buffer.data[buffer_sub_range] .= view(series, series_sub_range)
    if is_complete
        put!(output_channel, nothing)
    end
end

function fill_series_signal!(buffer, (series, buffers), sub_ranges_channel, output_channel)
    foreach(
        let buffer = buffer, series = series, sub_ranges_channel = sub_ranges_channel, output_channel = output_channel
            function (_)
                fill_buffer_signal!(
                    buffer,
                    series,
                    sub_ranges_channel,
                    output_channel,
                )
            end
        end,
        1:buffers,
    )
end

function fill_all_signal!(
    buffer,
    series_buffers_channel,
    sub_ranges_channel,
    output_channel,
)
    foreach(
        let buffer = buffer,
            sub_ranges_channel = sub_ranges_channel,
            output_channel = output_channel

            function (series_buffers)
                fill_series_signal!(
                    buffer,
                    series_buffers,
                    sub_ranges_channel,
                    output_channel,
                )
            end
        end,
        series_buffers_channel,
    )
end

function fill_all_task_io(buffer)
    series_buffers_channel = Channel{Tuple{Any,Int}}(Inf)
    sub_ranges_channel = Channel{Tuple{UnitRange{Int},UnitRange{Int},Bool}}(0)
    output_channel = Channel{Nothing}(0)
    task = Task(
        let buffer = buffer,
            series_buffers_channel = series_buffers_channel,
            sub_ranges_channel = sub_ranges_channel,
            output_channel = output_channel

            function ()
                fill_all_signal!(
                    buffer,
                    series_buffers_channel,
                    sub_ranges_channel,
                    output_channel,
                )
            end
        end,
    )
    task.sticky = false
    schedule(task)
    TaskIO(task, series_buffers_channel, sub_ranges_channel, output_channel)
end

function fill_all_task_ios(buffer; number_of_tasks = nthreads() - 1)
    map(let buffer = buffer
        function (_)
            fill_all_task_io(buffer)
        end
    end, 1:number_of_tasks)
end

function feed_write!(task_ios, buffer, (series_sub_ranges, buffer_sub_ranges, is_complete))
    foreach(
        let is_complete = is_complete
            function ((task_io, series_sub_range, output_sub_range),)
                put!(
                    task_io.sub_ranges_channel,
                    (series_sub_range, output_sub_range, is_complete),
                )
            end
        end,
        zip(task_ios, series_sub_ranges, buffer_sub_ranges),
    )
    if is_complete
        foreach(function (task_io)
            take!(task_io.output_channel)
        end, task_ios)
        write_buffer(buffer)
    end
end

function write_series!(task_ios, series, series_total, buffer, buffer_at)
    dual_sub_ranges =
        DualSubRanges(buffer_at, buffer.frames_per_buffer, series_total, length(task_ios))
    foreach(let series_buffers = (series, length(dual_sub_ranges))
        function (task_io)
            put!(task_io.series_buffers_channel, series_buffers)
        end
    end, task_ios)
    foreach(let task_ios = task_ios, buffer = buffer
        function (ranges)
            feed_write!(task_ios, buffer, ranges)
        end
    end, dual_sub_ranges)
    last(last(last(dual_sub_ranges)[2]))
end

function write(
    stream::PortAudioStream{<:Messanger{<:Any,Weaver},<:Any},
    a_schedule::AudioSchedule;
    number_of_tasks = nthreads() - 1,
)
    if number_of_tasks == 0
        throw(ArgumentError("Must write with at least 1 task"))
    end
    if nchannels(stream.source_messanger) != 0
        throw(ArgumentError("$stream does not have 0 input channels"))
    end
    sink_messanger = stream.sink_messanger
    if nchannels(sink_messanger) != 1
        throw(ArgumentError("$stream does not have 1 output channel"))
    end
    if samplerate(stream) != samplerate(a_schedule)
        throw(ArgumentError("Sample rates of $stream and $a_schedule do not match"))
    end
    buffer = sink_messanger.buffer
    buffer_at = 0
    # TODO: Move task_ios inside of the Weaver object and reuse them
    task_ios = fill_all_task_ios(buffer; number_of_tasks = number_of_tasks)
    try
        for (series, series_total) in a_schedule
            buffer_at = write_series!(task_ios, series, series_total, buffer, buffer_at)
        end
    finally
        foreach(close, task_ios)
    end
    # This might make the previous TODO difficult. Do we need to pass-back a shorter-than-expected result?
    if buffer_at > 0
        write_buffer(buffer, buffer_at)
    end
    nothing
end

struct TaskInput{Input}
    task::Task
    input_channel::Channel{Input}
end

function close(task_input::TaskInput)
    close(task_input.input_channel)
    wait(task_input.task)
end

function feed!(task_inputs, series, series_total, output_at)
    number_of_tasks = length(task_inputs)
    series_ranges = DividedRange(1:series_total, number_of_tasks)
    foreach(
        let series = series
            function feed_task_input!((task_input, series_range, output_range))
                put!(task_input.input_channel, (series, series_range, output_range))
            end
        end,
        zip(task_inputs, series_ranges, add_to(series_ranges, output_at)),
    )
    output_at + series_total
end

function fill_output!(output_data, (series, series_sub_range, output_range))
    output_data[output_range] .= view(series, series_sub_range)
end

function fill_all!(input_channel, output_data)
    foreach(let output_data = output_data
        function (series_ranges)
            fill_output!(output_data, series_ranges)
        end
    end, input_channel)
end

function fill_all_task_input(output_data)
    input_channel = Channel{Tuple{Any,UnitRange{Int},UnitRange{Int}}}(Inf)
    task = Task(let input_channel = input_channel, output_data = output_data
        function ()
            fill_all!(input_channel, output_data)
        end
    end)
    task.sticky = false
    schedule(task)
    TaskInput(task, input_channel)
end

function SampleBuf(a_schedule::AudioSchedule; number_of_tasks = nthreads() - 1)
    if number_of_tasks == 0
        throw(ArgumentError("Must write with at least 1 task"))
    end
    number = sum(Iterators.map(second, a_schedule))
    output = SampleBuf(Float64, samplerate(a_schedule), number, 1)
    task_inputs = map(let output_data = output.data
        function (_)
            fill_all_task_input(output_data)
        end
    end, 1:number_of_tasks)
    try
        output_at = 0
        for (series, series_total) in a_schedule
            output_at = feed!(task_inputs, series, series_total, output_at)
        end
    finally
        foreach(close, task_inputs)
    end
    output
end

include("equal_loudness.jl")

end


text:
function random_vec(::Type{VT}, range; rng=Random.GLOBAL_RNG) where {VT}
    dim = length(VT)
    T = eltype(VT)
    (lb, up) = range
    p = VT(lb .+ rand(rng, T, dim) .* (up - lb))

    return p
end

function adjust!(opts::EnsembleOptions)
    if opts.nattempt % opts.naccept == 0
        ratio = opts.naccept / opts.nattempt
        if ratio > opts.ensemble.accept
            opts.ensemble.δr *= 1.05
        else
            opts.ensemble.δr *= 0.95
        end
    end

    return nothing
end

function _box_information(box::CellListMap.Box)
    box_size = box.unit_cell_max[begin]
    cutoff = box.cutoff^2

    return box_size, cutoff
end


text:
function BaremetalPi.gpio_set(::Nothing)
    nothing
end

function BaremetalPi.gpio_clear(::Nothing)
    nothing
end

function BaremetalPi.gpio_set_mode(::Nothing)
    nothing
end


text:
export kl_divergence

function kl_divergence(p::Var, q::Var)
    y = kl_divergence(p.data, q.data)
    function df(gy)
        hasgrad(p) && ∇kl_divergence_p!(p.data, p.grad, q.data)
        hasgrad(q) && ∇kl_divergence_q!(p.data, q.data, q.grad)
    end
    Var(y, [p,q], kl_divergence, df)
end

function kl_divergence{T}(p::Array{T}, q::Array{T})
    d = T(0)
    @inbounds @simd for i = 1:length(p)
        d += p[i] * (log(p[i]) - log(q[i]))
    end
    d
end

function ∇kl_divergence_p!{T}(p::Array{T}, gp::Array{T}, q::Array{T})
    @inbounds @simd for i = 1:length(p)
        gp[i] += log(p[i]) - log(q[i]) + 1
    end
end

function ∇kl_divergence_q!{T}(p::Array{T}, q::Array{T}, gq::Array{T})
    @inbounds @simd for i = 1:length(p)
        gq[i] += -p[i] / q[i]
    end
end


text:
using MetaGraphs
using Test
using LightGraphs

property_set=[
    (src="thirty",dst="inf",color="missing",penwidth="missing",style="missing",),
    (src="twenty",dst="mort",color="missing",penwidth="missing",style="missing",),
    (src="rate",dst="personper",color="missing",penwidth="missing",style="missing",),
    (src="person",dst="unit",color="missing",penwidth="missing",style="missing",),
    (src="personper",dst="unit",color="missing",penwidth="missing",style="missing",),
    (src="second",dst="unit",color="missing",penwidth="missing",style="missing",),
    (src="dollars",dst="unit",color="missing",penwidth="missing",style="missing",),
    (src="dGDP",dst="prices",color="missing",penwidth="missing",style="missing",),
    (src="prices",dst="regres",color="missing",penwidth="missing",style="missing",),
    (src="fed",dst="dGDP",color="missing",penwidth="missing",style="missing",),
    (src="ode",dst="epi",color="missing",penwidth="missing",style="missing",),
    (src="dollars",dst="cost",color="missing",penwidth="missing",style="missing",),
    (src="mort",dst="rate",color="black",penwidth="2.0",style="dashed",),
    (src="inf",dst="rate",color="black",penwidth="2.0",style="dashed",),
    (src="birth",dst="rate",color="black",penwidth="2.0",style="dashed",),
    (src="mort",dst="ind",color="orange",penwidth="4.0",style="solid",),
    (src="ind",dst="epi",color="orange",penwidth="4.0",style="solid",),
    (src="ind",dst="epi",color="orange",penwidth="4.0",style="solid",),
    (src="inf",dst="ind",color="orange",penwidth="4.0",style="solid",),
    (src="birth",dst="ind",color="orange",penwidth="4.0",style="solid",),
    (src="temp",dst="inf",color="orange",penwidth="4.0",style="solid",),
    (src="age",dst="mort",color="orange",penwidth="4.0",style="solid",),
    (src="demo",dst="birth",color="orange",penwidth="4.0",style="solid",),
    (src="epi",dst="cases",color="orange",penwidth="4.0",style="solid",),
    (src="cases",dst="regres",color="orange",penwidth="4.0",style="solid",),
    (src="weather",dst="temp",color="orange",penwidth="4.0",style="solid",),
    (src="demo",dst="age",color="orange",penwidth="4.0",style="solid",),
    (src="regres",dst="cost",color="orange",penwidth="4.0",style="solid",)
]

# name, label, color
vprops = [
    ("weather","NOAA\\nWeather","orange"),
    ("cost","Vaccination\\nCost","orange"),
    ("demo","Census\\nDemographics","orange"),
    ("fed", "Fed Forcast", "#5DADE2"),
    ("epi","SIR","#5DADE2"),
    ("ode", "ODE Solver", "#5DADE2"),
    ("rate","{Transition\\nRate}","#66AA55"),
    ("unit","Unit","#66AA55"),
    ("personper","Person/s","#66AA55"),
    ("person","Person","#66AA55"),
    ("second","second (s)","#66AA55"),
    ("dollars","\$","#66AA55"),
    ("inf","Infection\\nRate","#66AA55"),
    ("mort","Mortality\\nRate","#66AA55"),
    ("birth","Birth\\nRate","#66AA55"),
    ("twenty","0.2 Persons/s","#DD1133"),
    ("thirty","0.3 Persons/s","#DD1133"),
    ("ind","Individual\\nContact\\nModel","#5DADE2"),
    ("temp","Temperature","#5DADE2"),
    ("age","Age","#5DADE2"),
    ("dGDP","Economic Growth","#5DADE2"),
    ("cases","Flu\\nCases","#5DADE2"),
    ("prices","Vacc\\nPrice","#5DADE2"),
    ("regres", "Regression", "#5DADE2"),
]

g = MetaDiGraph(24)

set_prop!(g, :pack, true)

for (v, prop) in enumerate(vprops)
    set_prop!(g, v, :name,  prop[1])
    set_prop!(g, v, :label, prop[2])
    set_prop!(g, v, :color, prop[3])
end

# enable indexing with string names
set_indexing_prop!(g, :name)

# add edges
for prop in property_set
    src, dst = g[prop.src, :name], g[prop.dst, :name]
    add_edge!(g, src,dst)
    set_prop!(g, src, dst, :color, prop.color)
    set_prop!(g, src, dst, :penwidth, prop.penwidth)
    set_prop!(g, src, dst, :style, prop.style)
end

# set global edge properties
for e in edges(g)
    set_prop!(g, e, :dir, :none)
end
# set global vertex properties
for v in vertices(g)
    set_prop!(g, v, :shape, :record)
    set_prop!(g, v, :style, :filled)
    set_prop!(g, v, :fillcolor, "#dddddd")
    set_prop!(g, v, :penwidth, 2.0)
end

savegraph("diagram.dot", g, DOTFormat())
open("diagram_ref.dot") do fp
    s_ref = read(fp)
    open("diagram.dot") do fp
        s = read(fp)
        @test s == s_ref
    end
end


text:
export softmax, logsoftmax

doc"""
    softmax(x, dim::Int)

Softmax function over the given dimension.

```math
f(x) = \exp(x) \over \sum \exp(x)
```
"""
softmax(x::Var) = Var(softmax(x.data), ∇softmax!, (x,))
softmax(x::CuArray) = CUDNN.softmax(x)
softmax(x::Node) = Node(softmax, (x,))

function softmax(x::Matrix{T}) where T
    y = similar(x)
    @inbounds for j = 1:size(x,2)
        maxv = x[1,j]
        for i = 1:size(x,1)
            maxv = max(maxv, x[i,j])
        end
        z = T(0)
        for i = 1:size(x,1)
            y[i,j] = exp(x[i,j] - maxv)
            z += y[i,j]
        end
        z == T(0) && throw("z == 0")
        invz = 1 / z
        for i = 1:size(x,1)
            y[i,j] *= invz
        end
    end
    y
end

function ∇softmax!(y::Var, x::Var)
    isnothing(x.grad) && return
    ∇softmax!(y.data, y.grad, x.grad)
end

∇softmax!(y::CuArray, gy::CuArray, gx::CuArray) = CUDNN.∇softmax!(y, gy, gx)

function ∇softmax!(y::Matrix{T}, gy::Matrix{T}, gx::Matrix{T}) where T
    @inbounds for j = 1:size(y,2)
        sum = T(0)
        for i = 1:size(y,1)
            sum += gy[i,j] * y[i,j]
        end
        for i = 1:size(y,1)
            gx[i,j] += y[i,j] * (gy[i,j]-sum)
        end
    end
end

"""
    logsoftmax(x)

Logarithm of softmax function.
"""
logsoftmax(x::Var) = Var(logsoftmax(x.data), ∇logsoftmax!, (x,))
logsoftmax(x::CuArray) = CUDNN.softmax(x, CUDNN.CUDNN_SOFTMAX_LOG)
logsoftmax(x::Node) = Node(logsoftmax, (x,))

function logsoftmax(x::Matrix{T}) where T
    y = similar(x)
    max = maximum(x, dims=1)
    @inbounds for j = 1:size(x,2)
        sum = T(1e-10)
        for i = 1:size(x,1)
            sum += exp(x[i,j] - max[j])
        end
        logz = log(sum)
        for i = 1:size(x,1)
            y[i,j] = x[i,j] - max[j] - logz
        end
    end
    y
end

function ∇logsoftmax!(y::Var, x::Var)
    isnothing(x.grad) && return
    ∇logsoftmax!(y.data, y.grad, x.grad)
end

function ∇logsoftmax!(y::CuArray, gy::CuArray, gx::CuArray)
    CUDNN.∇softmax!(y, gy, gx, CUDNN.CUDNN_SOFTMAX_LOG)
end

function ∇logsoftmax!(y::Matrix{T}, gy::Matrix{T}, gx::Matrix{T}) where T
    @inbounds for j = 1:size(y,2)
        sum = T(0)
        for i = 1:size(y,1)
            sum += gy[i,j]
        end
        for i = 1:size(y,1)
            gx[i,j] += gy[i,j] - exp(y[i,j]) * sum
        end
    end
end


text:
#Module for either generating data or exporting from an existing dataset

module DataExamples

export get_Waveform, get_BreastCancer

function get_BreastCancer()
    data = readdlm("../data/Processed_BreastCancer.data",',')
    X = convert(Array{Float64,2},data[:,1:end-1])
    y = convert(Array{Float64,1},data[:,end])
    nSamples = size(X,1)
    shuffling = true
    DatasetName = "BreastCancer"
    Z = hcat(X,y)
    if shuffling
      Z = Z[shuffle(collect(1:nSamples)),:] #Shuffle the data
    end
    (Z[:,1:end-1],Z[:,end],DatasetName)
end

function get_Waveform()
    data = readdlm("../data/Processed_Waveform.data",',')
    X = convert(Array{Float64,2},data[:,1:end-1])
    y = convert(Array{Float64,1},data[:,end])
    nSamples = size(X,1)
    shuffling = false
    DatasetName = "Waveform"
    Z = hcat(X,y)
    if shuffling
      Z = Z[shuffle(collect(1:nSamples)),:] #Shuffle the data
    end
    (Z[:,1:end-1],Z[:,end],DatasetName)
end

end #end of module


text:
# Automatically-generated headers for cusparse
#
# DO NOT EDIT THIS FILE DIRECTLY.

# Automatically generated using Clang.jl

const CUSPARSE_VER_MAJOR = 10
const CUSPARSE_VER_MINOR = 3
const CUSPARSE_VER_PATCH = 1
const CUSPARSE_VER_BUILD = 89
const CUSPARSE_VERSION = CUSPARSE_VER_MAJOR * 1000 + CUSPARSE_VER_MINOR * 100 + CUSPARSE_VER_PATCH
const cusparseContext = Cvoid
const cusparseHandle_t = Ptr{cusparseContext}
const cusparseSolveAnalysisInfo = Cvoid
const cusparseSolveAnalysisInfo_t = Ptr{cusparseSolveAnalysisInfo}
const csrsv2Info = Cvoid
const csrsv2Info_t = Ptr{csrsv2Info}
const csrsm2Info = Cvoid
const csrsm2Info_t = Ptr{csrsm2Info}
const bsrsv2Info = Cvoid
const bsrsv2Info_t = Ptr{bsrsv2Info}
const bsrsm2Info = Cvoid
const bsrsm2Info_t = Ptr{bsrsm2Info}
const csric02Info = Cvoid
const csric02Info_t = Ptr{csric02Info}
const bsric02Info = Cvoid
const bsric02Info_t = Ptr{bsric02Info}
const csrilu02Info = Cvoid
const csrilu02Info_t = Ptr{csrilu02Info}
const bsrilu02Info = Cvoid
const bsrilu02Info_t = Ptr{bsrilu02Info}
const cusparseHybMat = Cvoid
const cusparseHybMat_t = Ptr{cusparseHybMat}
const csrgemm2Info = Cvoid
const csrgemm2Info_t = Ptr{csrgemm2Info}
const csru2csrInfo = Cvoid
const csru2csrInfo_t = Ptr{csru2csrInfo}
const cusparseColorInfo = Cvoid
const cusparseColorInfo_t = Ptr{cusparseColorInfo}
const pruneInfo = Cvoid
const pruneInfo_t = Ptr{pruneInfo}

@cenum cusparseStatus_t::UInt32 begin
    CUSPARSE_STATUS_SUCCESS = 0
    CUSPARSE_STATUS_NOT_INITIALIZED = 1
    CUSPARSE_STATUS_ALLOC_FAILED = 2
    CUSPARSE_STATUS_INVALID_VALUE = 3
    CUSPARSE_STATUS_ARCH_MISMATCH = 4
    CUSPARSE_STATUS_MAPPING_ERROR = 5
    CUSPARSE_STATUS_EXECUTION_FAILED = 6
    CUSPARSE_STATUS_INTERNAL_ERROR = 7
    CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED = 8
    CUSPARSE_STATUS_ZERO_PIVOT = 9
    CUSPARSE_STATUS_NOT_SUPPORTED = 10
end

@cenum cusparsePointerMode_t::UInt32 begin
    CUSPARSE_POINTER_MODE_HOST = 0
    CUSPARSE_POINTER_MODE_DEVICE = 1
end

@cenum cusparseAction_t::UInt32 begin
    CUSPARSE_ACTION_SYMBOLIC = 0
    CUSPARSE_ACTION_NUMERIC = 1
end

@cenum cusparseMatrixType_t::UInt32 begin
    CUSPARSE_MATRIX_TYPE_GENERAL = 0
    CUSPARSE_MATRIX_TYPE_SYMMETRIC = 1
    CUSPARSE_MATRIX_TYPE_HERMITIAN = 2
    CUSPARSE_MATRIX_TYPE_TRIANGULAR = 3
end

@cenum cusparseFillMode_t::UInt32 begin
    CUSPARSE_FILL_MODE_LOWER = 0
    CUSPARSE_FILL_MODE_UPPER = 1
end

@cenum cusparseDiagType_t::UInt32 begin
    CUSPARSE_DIAG_TYPE_NON_UNIT = 0
    CUSPARSE_DIAG_TYPE_UNIT = 1
end

@cenum cusparseIndexBase_t::UInt32 begin
    CUSPARSE_INDEX_BASE_ZERO = 0
    CUSPARSE_INDEX_BASE_ONE = 1
end

@cenum cusparseOperation_t::UInt32 begin
    CUSPARSE_OPERATION_NON_TRANSPOSE = 0
    CUSPARSE_OPERATION_TRANSPOSE = 1
    CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE = 2
end

@cenum cusparseDirection_t::UInt32 begin
    CUSPARSE_DIRECTION_ROW = 0
    CUSPARSE_DIRECTION_COLUMN = 1
end

@cenum cusparseHybPartition_t::UInt32 begin
    CUSPARSE_HYB_PARTITION_AUTO = 0
    CUSPARSE_HYB_PARTITION_USER = 1
    CUSPARSE_HYB_PARTITION_MAX = 2
end

@cenum cusparseSolvePolicy_t::UInt32 begin
    CUSPARSE_SOLVE_POLICY_NO_LEVEL = 0
    CUSPARSE_SOLVE_POLICY_USE_LEVEL = 1
end

@cenum cusparseSideMode_t::UInt32 begin
    CUSPARSE_SIDE_LEFT = 0
    CUSPARSE_SIDE_RIGHT = 1
end

@cenum cusparseColorAlg_t::UInt32 begin
    CUSPARSE_COLOR_ALG0 = 0
    CUSPARSE_COLOR_ALG1 = 1
end

@cenum cusparseAlgMode_t::UInt32 begin
    CUSPARSE_ALG0 = 0
    CUSPARSE_ALG1 = 1
    CUSPARSE_ALG_NAIVE = 0
    CUSPARSE_ALG_MERGE_PATH = 1
end

@cenum cusparseCsr2CscAlg_t::UInt32 begin
    CUSPARSE_CSR2CSC_ALG1 = 1
    CUSPARSE_CSR2CSC_ALG2 = 2
end

@cenum cusparseFormat_t::UInt32 begin
    CUSPARSE_FORMAT_CSR = 1
    CUSPARSE_FORMAT_CSC = 2
    CUSPARSE_FORMAT_COO = 3
    CUSPARSE_FORMAT_COO_AOS = 4
end

@cenum cusparseOrder_t::UInt32 begin
    CUSPARSE_ORDER_COL = 1
    CUSPARSE_ORDER_ROW = 2
end

@cenum cusparseSpMVAlg_t::UInt32 begin
    CUSPARSE_MV_ALG_DEFAULT = 0
    CUSPARSE_COOMV_ALG = 1
    CUSPARSE_CSRMV_ALG1 = 2
    CUSPARSE_CSRMV_ALG2 = 3
end

@cenum cusparseSpMMAlg_t::UInt32 begin
    CUSPARSE_MM_ALG_DEFAULT = 0
    CUSPARSE_COOMM_ALG1 = 1
    CUSPARSE_COOMM_ALG2 = 2
    CUSPARSE_COOMM_ALG3 = 3
    CUSPARSE_CSRMM_ALG1 = 4
end

@cenum cusparseIndexType_t::UInt32 begin
    CUSPARSE_INDEX_16U = 1
    CUSPARSE_INDEX_32I = 2
    CUSPARSE_INDEX_64I = 3
end

const cusparseSpVecDescr = Cvoid
const cusparseDnVecDescr = Cvoid
const cusparseSpMatDescr = Cvoid
const cusparseDnMatDescr = Cvoid
const cusparseSpVecDescr_t = Ptr{cusparseSpVecDescr}
const cusparseDnVecDescr_t = Ptr{cusparseDnVecDescr}
const cusparseSpMatDescr_t = Ptr{cusparseSpMatDescr}
const cusparseDnMatDescr_t = Ptr{cusparseDnMatDescr}

"""
Describes shape and properties of a CUSPARSE matrix. A convenience wrapper.

Contains:
* `MatrixType` - a [`cusparseMatrixType_t`](@ref)
* `FillMode` - a [`cusparseFillMode_t`](@ref)
* `DiagType` - a [`cusparseDiagType_t`](@ref)
* `IndexBase` - a [`cusparseIndexBase_t`](@ref)
"""
struct cusparseMatDescr
    MatrixType::cusparseMatrixType_t
    FillMode::cusparseFillMode_t
    DiagType::cusparseDiagType_t
    IndexBase::cusparseIndexBase_t
    function cusparseMatDescr(MatrixType,FillMode,DiagType,IndexBase)
        new(MatrixType,FillMode,DiagType,IndexBase)
    end
end

const cusparseMatDescr_t = Ptr{cusparseMatDescr}


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

abstract AbstractTime

"""
    Period
    Year
    Month
    Week
    Day
    Hour
    Minute
    Second
    Millisecond

`Period` types represent discrete, human representations of time.
"""
abstract Period     <: AbstractTime
abstract DatePeriod <: Period
abstract TimePeriod <: Period

for T in (:Year,:Month,:Week,:Day)
    @eval immutable $T <: DatePeriod
        value::Int64
        $T(v::Number) = new(v)
    end
end
for T in (:Hour,:Minute,:Second,:Millisecond)
    @eval immutable $T <: TimePeriod
        value::Int64
        $T(v::Number) = new(v)
    end
end

"""
    Year(v)
    Month(v)
    Week(v)
    Day(v)
    Hour(v)
    Minute(v)
    Second(v)
    Millisecond(v)

Construct a `Period` type with the given `v` value. Input must be losslessly convertible
to an `Int64`.
"""
Period(v)

"""
    Instant

`Instant` types represent integer-based, machine representations of time as continuous
timelines starting from an epoch.
"""
abstract Instant <: AbstractTime

"""
    UTInstant{T}

The `UTInstant` represents a machine timeline based on UT time (1 day = one revolution of
the earth). The `T` is a `Period` parameter that indicates the resolution or precision of
the instant.
"""
immutable UTInstant{P<:Period} <: Instant
    periods::P
end

# Convenience default constructors
UTM(x) = UTInstant(Millisecond(x))
UTD(x) = UTInstant(Day(x))

# Calendar types provide rules for interpretating instant
# timelines in human-readable form.
abstract Calendar <: AbstractTime

# ISOCalendar implements the ISO 8601 standard (en.wikipedia.org/wiki/ISO_8601)
# Notably based on the proleptic Gregorian calendar
# ISOCalendar provides interpretation rules for UTInstants to civil date and time parts
immutable ISOCalendar <: Calendar end

abstract TimeZone
immutable UTC <: TimeZone end

"""
    TimeType

`TimeType` types wrap `Instant` machine instances to provide human representations of the
machine instant. Both `DateTime` and `Date` are subtypes of `TimeType`.
"""
abstract TimeType <: AbstractTime

"""
    DateTime

`DateTime` wraps a `UTInstant{Millisecond}` and interprets it according to the proleptic
Gregorian calendar.
"""
immutable DateTime <: TimeType
    instant::UTInstant{Millisecond}
    DateTime(instant::UTInstant{Millisecond}) = new(instant)
end

"""
    Date

`Date` wraps a `UTInstant{Day}` and interprets it according to the proleptic Gregorian calendar.
"""
immutable Date <: TimeType
    instant::UTInstant{Day}
    Date(instant::UTInstant{Day}) = new(instant)
end

# Fallback constructors
_c(x) = convert(Int64,x)
DateTime(y,m=1,d=1,h=0,mi=0,s=0,ms=0) = DateTime(_c(y),_c(m),_c(d),_c(h),_c(mi),_c(s),_c(ms))
Date(y,m=1,d=1) = Date(_c(y),_c(m),_c(d))

# Convert y,m,d to # of Rata Die days
# Works by shifting the beginning of the year to March 1,
# so a leap day is the very last day of the year
const SHIFTEDMONTHDAYS = [306,337,0,31,61,92,122,153,184,214,245,275]
function totaldays(y,m,d)
    # If we're in Jan/Feb, shift the given year back one
    z = m < 3 ? y - 1 : y
    mdays = SHIFTEDMONTHDAYS[m]
    # days + month_days + year_days
    return d + mdays + 365z + fld(z,4) - fld(z,100) + fld(z,400) - 306
end

# If the year is divisible by 4, except for every 100 years, except for every 400 years
isleapyear(y) = ((y % 4 == 0) && (y % 100 != 0)) || (y % 400 == 0)

# Number of days in month
const DAYSINMONTH = [31,28,31,30,31,30,31,31,30,31,30,31]
daysinmonth(y,m) = DAYSINMONTH[m] + (m == 2 && isleapyear(y))

### CONSTRUCTORS ###
# Core constructors
"""
    DateTime(y, [m, d, h, mi, s, ms]) -> DateTime

Construct a `DateTime` type by parts. Arguments must be convertible to `Int64`.
"""
function DateTime(y::Int64,m::Int64=1,d::Int64=1,
                  h::Int64=0,mi::Int64=0,s::Int64=0,ms::Int64=0)
    0 < m < 13 || throw(ArgumentError("Month: $m out of range (1:12)"))
    0 < d < daysinmonth(y,m)+1 || throw(ArgumentError("Day: $d out of range (1:$(daysinmonth(y,m)))"))
    -1 < h < 24 || throw(ArgumentError("Hour: $h out of range (0:23)"))
    -1 < mi < 60 || throw(ArgumentError("Minute: $mi out of range (0:59)"))
    -1 < s < 60 || throw(ArgumentError("Second: $s out of range (0:59)"))
    -1 < ms < 1000 || throw(ArgumentError("Millisecond: $ms out of range (0:999)"))
    rata = ms + 1000*(s + 60mi + 3600h + 86400*totaldays(y,m,d))
    return DateTime(UTM(rata))
end

"""
    Date(y, [m, d]) -> Date

Construct a `Date` type by parts. Arguments must be convertible to `Int64`.
"""
function Date(y::Int64,m::Int64=1,d::Int64=1)
    0 < m < 13 || throw(ArgumentError("Month: $m out of range (1:12)"))
    0 < d < daysinmonth(y,m)+1 || throw(ArgumentError("Day: $d out of range (1:$(daysinmonth(y,m)))"))
    return Date(UTD(totaldays(y,m,d)))
end

# Convenience constructors from Periods
function DateTime(y::Year,m::Month=Month(1),d::Day=Day(1),
                  h::Hour=Hour(0),mi::Minute=Minute(0),
                  s::Second=Second(0),ms::Millisecond=Millisecond(0))
    return DateTime(value(y),value(m),value(d),
                        value(h),value(mi),value(s),value(ms))
end

Date(y::Year,m::Month=Month(1),d::Day=Day(1)) = Date(value(y),value(m),value(d))

# To allow any order/combination of Periods

"""
    DateTime(periods::Period...) -> DateTime

Construct a `DateTime` type by `Period` type parts. Arguments may be in any order. DateTime
parts not provided will default to the value of `Dates.default(period)`.
"""
function DateTime(periods::Period...)
    y = Year(1); m = Month(1); d = Day(1)
    h = Hour(0); mi = Minute(0); s = Second(0); ms = Millisecond(0)
    for p in periods
        isa(p, Year) && (y = p::Year)
        isa(p, Month) && (m = p::Month)
        isa(p, Day) && (d = p::Day)
        isa(p, Hour) && (h = p::Hour)
        isa(p, Minute) && (mi = p::Minute)
        isa(p, Second) && (s = p::Second)
        isa(p, Millisecond) && (ms = p::Millisecond)
    end
    return DateTime(y,m,d,h,mi,s,ms)
end

"""
    Date(period::Period...) -> Date

Construct a `Date` type by `Period` type parts. Arguments may be in any order. `Date` parts
not provided will default to the value of `Dates.default(period)`.
"""
function Date(periods::Period...)
    y = Year(1); m = Month(1); d = Day(1)
    for p in periods
        isa(p, Year) && (y = p::Year)
        isa(p, Month) && (m = p::Month)
        isa(p, Day) && (d = p::Day)
    end
    return Date(y,m,d)
end

# Traits, Equality
Base.isfinite{T<:TimeType}(::Union{Type{T},T}) = true
calendar(dt::DateTime) = ISOCalendar
calendar(dt::Date) = ISOCalendar

"""
    eps(::DateTime) -> Millisecond
    eps(::Date) -> Day

Returns `Millisecond(1)` for `DateTime` values and `Day(1)` for `Date` values.
"""
Base.eps

Base.eps(dt::DateTime) = Millisecond(1)
Base.eps(dt::Date) = Day(1)

Base.typemax(::Union{DateTime,Type{DateTime}}) = DateTime(146138512,12,31,23,59,59)
Base.typemin(::Union{DateTime,Type{DateTime}}) = DateTime(-146138511,1,1,0,0,0)
Base.typemax(::Union{Date,Type{Date}}) = Date(252522163911149,12,31)
Base.typemin(::Union{Date,Type{Date}}) = Date(-252522163911150,1,1)
# Date-DateTime promotion, isless, ==
Base.promote_rule(::Type{Date},x::Type{DateTime}) = DateTime
Base.isless(x::Date,y::Date) = isless(value(x),value(y))
Base.isless(x::DateTime,y::DateTime) = isless(value(x),value(y))
Base.isless(x::TimeType,y::TimeType) = isless(promote(x,y)...)
==(x::TimeType,y::TimeType) = ===(promote(x,y)...)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
export @ein_str, @ein, ein
"""
    ein"ij,jk -> ik"(A,B)

String macro interface which understands `numpy.einsum`'s notation.
Translates strings into `StaticEinCode`-structs that can be called to evaluate
an `einsum`.
To control evaluation order, use parentheses - instead of an `EinCode`,
a `NestedEinsum` is returned which evaluates the expression
according to parens.
The valid character ranges for index-labels are `a-z` and `α-ω`.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b, c = rand(10,10), rand(10,10), rand(10,1);

julia> ein"ij,jk,kl -> il"(a,b,c) ≈ ein"(ij,jk),kl -> il"(a,b,c) ≈ a * b * c
true
```
"""
macro ein_str(s::AbstractString)
    ein(s)
end

function ein(s::AbstractString)
    s = replace(replace(s, "\n" => ""), " "=>"")
    m = match(r"([\(\)a-z,α-ω]*)->([a-zα-ω]*)", s)
    m === nothing && throw(ArgumentError("invalid einsum specification $s"))
    sixs, siy = m.captures
    if '(' in sixs
        return parse_nested(sixs, collect(siy))
    else
        iy  = Tuple(siy)
        ixs = Tuple(Tuple(ix) for ix in split(sixs,','))
        return StaticEinCode{ixs, iy}()
    end
end

function (code::DynamicEinCode{LT})(@nospecialize(xs...); size_info=nothing) where LT
    size_dict = get_size_dict!(getixs(code), xs, size_info===nothing ? Dict{LT,Int}() : copy(size_info))
    einsum(code, xs, size_dict)
end

function (code::StaticEinCode)(xs...; size_info=nothing)
    size_dict = get_size_dict!(getixs(code), xs, size_info===nothing ? Dict{labeltype(code),Int}() : copy(size_info))
    einsum(code, xs, size_dict)
end

# 2us overheads if @nospecialize
@doc raw"
    get_size_dict!(ixs, xs, size_info)

return a dictionary that is used to get the size of an index-label
in the einsum-specification with input-indices `ixs` and tensors `xs` after
consistency within `ixs` and between `ixs` and `xs` has been verified.
"
@inline function get_size_dict!(ixs, xs, size_info::Dict{LT}) where LT
    if length(ixs) == 1
        get_size_dict_unary!(ixs[1], size(xs[1]), size_info)
    else
        get_size_dict_!(ixs, [collect(Int, size(x)) for x in xs], size_info)
    end
end

function get_size_dict_!(ixs, sizes::AbstractVector, size_info::Dict{LT}) where LT
    # check size of input tuples
    length(sizes)<1 && error("empty input tensors")
    length(ixs) != length(sizes) && throw(ArgumentError("$(length(sizes)) tensors labelled by $(length(ixs)) indices"))
    # check tensor orders
    @inbounds for i=1:length(sizes)
        ix, s = ixs[i], sizes[i]
        length(ix) == length(s) || throw(
            ArgumentError("indices $ix invalid for tensor with ndims = $(length(s))"))
        for j = 1:length(ix)
            k = ix[j]
            if haskey(size_info, k)
                s[j] == size_info[k] || throw(DimensionMismatch("$k = $(size_info[k]) or $(s[j]))?"))
            else
                size_info[k] = s[j]
            end
        end
    end
    return size_info
end
# to speed up unary operations
function get_size_dict_unary!(ix, s, size_info::Dict{LT}) where LT
    @inbounds for j = 1:length(ix)
        k = ix[j]
        if haskey(size_info, k)
            s[j] == size_info[k] || throw(DimensionMismatch("$k = $(size_info[k]) or $(s[j]))?"))
        else
            size_info[k] = s[j]
        end
    end
    return size_info
end

@inline function get_size_dict(ixs, xs, size_info=nothing)
    LT = promote_type(eltype.(ixs)...)
    return get_size_dict!(ixs, xs, size_info===nothing ? Dict{LT,Int}() : size_info)
end

using MacroTools
"""
    @ein A[i,k] := B[i,j] * C[j,k]     # A = B * C

Macro interface similar to that of other packages.

You may use numbers in place of letters for dummy indices, as in `@tensor`,
and need not name the output array. Thus `A = @ein [1,2] := B[1,ξ] * C[ξ,2]`
is equivalent to the above. This can also be written `A = ein"ij,jk -> ik"(B,C)`
using the numpy-style string macro.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> @ein c[i,k] := a[i,j] * b[j,k];

julia> c ≈ a * b
true

julia> c ≈ ein"ij,jk -> ik"(a,b)
true
```
"""
macro ein(exs...)
    _ein_macro(exs...)
end

primefix!(ind) = map!(i -> @capture(i, (j_)') ? Symbol(j, '′') : i, ind, ind)

function _ein_macro(ex; einsum=:einsum)
    @capture(ex, (left_ := right_)) || throw(ArgumentError("expected A[] := B[]... "))
    @capture(left, Z_[leftind__] | [leftind__] ) || throw(
        ArgumentError("can't understand LHS, expected A[i,j] etc."))
    if Z===nothing
        @gensym Z
    end
    primefix!(leftind)

    rightind, rightpairs = [], []
    @capture(right, *(factors__)) || (factors = Any[right])
    for fact in factors
        @capture(fact, A_[Aind__]) || return _nested_ein_macro(ex)
        primefix!(Aind)
        append!(rightind, Aind)
        push!(rightpairs, (A, Aind) )
    end
    unique!(rightind)
    isempty(setdiff(leftind, rightind)) || throw(
        ArgumentError("some indices appear only on the left"))

    lefttuple = Tuple(indexin(leftind, rightind))
    righttuples = [ Tuple(indexin(ind, rightind)) for (A, ind) in rightpairs ]
    rightnames = [ esc(A) for (A, ind) in rightpairs ]

    return :( $(esc(Z)) = $einsum( EinCode(($(righttuples...),), $lefttuple), ($(rightnames...),)) )
end

@doc raw"
    einsum(code::EinCode, xs, size_dict)
    einsum(rule, ixs, iy, xs, size_dict)

return the tensor that results from contracting the tensors `xs` according
to their indices `ixs` (`getixs(code)`), where all indices that do not appear in the output `iy` (`getiy(code)`) are
summed over.
The result is permuted according to `out`.

- `ixs` - tuple of tuples of index-labels of the input-tensors `xs`

- `iy` - tuple of index-labels of the output-tensor

- `xs` - tuple of tensors

- `size_dict` - a dictionary that maps index-labels to their sizes

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> einsum(EinCode((('i','j'),('j','k')),('i','k')), (a, b)) ≈ a * b
true

julia> einsum(EinCode((('i','j'),('j','k')),('k','i')), (a, b)) ≈ permutedims(a * b, (2,1))
true
```
"
@generated function einsum(code::StaticEinCode{ixs, iy}, xs::Tuple, size_dict::Dict{LT}) where {LT, ixs, iy}
    rule = match_rule(ixs, iy)
    :(einsum($rule, $ixs, $iy, xs, size_dict))
end

function einsum(code::DynamicEinCode, @nospecialize(xs::Tuple), size_dict::Dict)
    rule = match_rule(getixs(code), getiy(code))
    einsum(rule, getixs(code), getiy(code), xs, size_dict)
end

function einsum(code::EinCode, @nospecialize(xs::Tuple))
    einsum(code, xs, get_size_dict!(getixs(code), xs, Dict{labeltype(code),Int}()))
end

# the fallback
function einsum(::DefaultRule, ixs, iy, xs::Tuple, size_dict)
    @debug "DefaultRule loop_einsum" ixs => iy size.(xs)
    loop_einsum(EinCode(ixs, iy), (xs...,), size_dict)
end


text:
"""
    logsumexp(x::Vector{Complex{T}})::Complex{T}

Log-sum-exp with choosing the shift to be the maximum absolute value for
the real and imaginary parts of x's components.

# Example:
```
a = randn(3) + randn(3) .* im
u = exp.(a)
RHS = logsumexp(u)
LHS = log(sum(exp.(u)))
println("RHS = ", RHS)
println("LHS = ", LHS)
```
"""
function logsumexp(x::Vector{Complex{T}})::Complex{T} where T <: Real

    # use the maximum magnitude for real and imaginary parts to get the shift, a.
    val_unused, 𝑖_real = findmax( real.(x) )
    val_unused, 𝑖_imag = findmax( imag.(x) )

    a::Complex{T} = real(x[𝑖_real]) + im*imag(𝑖_imag)

    # log-sum-exp.
    running_sum::Complex{T} = zero(Complex{T})
    for i = 1:length(x)
        running_sum += exp( x[i] - a )
    end

    return a + log(running_sum)
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
# This Julia script is the code for the example of a piecewise-linear
# function.
using HomotopyContinuation
using LinearAlgebra
r = 10; # number of points
d = 3*r; # highest used moment
s = div(d+1, 2); # bound on Fourier coefficients
theseed = 123;

# We generate random input data.
import Distributions: Random, Uniform, Exponential, Normal
Random.seed!(theseed);

t = sort(rand(Uniform(-π,π), r)); # locations
t_ext = [t; 0];
f₀ = [0; rand(Normal(0,1), r-1); 0]; # jumps
f₁ = [0; rand(Normal(0,1), r-1); 0]; # slopes
# distribute slopes equally on both end points of segments
f₀ -= f₁ .* [0; t_ext[2:end] - t_ext[1:end-1]] / 2;
noise = rand(Normal(0,1e-12), 2*s+1);
f₀[8] -= 0.5;

fourierCoeffScaled = k -> ( # scaled by 2π(ik)^2
    @assert(-s ≤ k ≤ s);
    weights = (im * k * (f₀[2:end] - f₀[1:end-1]
                         + f₁[1:end-1].*(t_ext[1:end-1] - t_ext[2:end]))
               + f₁[2:end] - f₁[1:end-1]);
    transpose(weights) * exp.(-im * k * t) + noise[s+k+1]);

moment = k -> (@assert(0 ≤ k ≤ d ≤ 2*s);
               k == s ? 0.0*im : fourierCoeffScaled(k-s));

# We solve the quadratic system given by the moments via homotopy continuation.
@polyvar p[0:r-1];
P = [p;1];
X = [sum(P[i]*P[j] + 0.0 for i=1:(r+1) for j= 1:(r+1) if i+j == k) for k = 2:2*r+2];
momentmat = (a,b) -> hcat([[moment(i+j) for i = 0:a] for j = 0:b]...);
H₀ = momentmat(r-1, 2*r);
H₁ = momentmat(r, 2*r);
sols = solve(H₀ * X, seed=theseed);

# We pick the solution polynomial q that best solves H₁*q^2=0, using one additional moment.
computeErr(q) = (Q² = [map(xk -> xk(p => q), X[1:2*r]); 1];
                 norm(H₁ * Q²));
errs = map(q -> computeErr(q.solution), finite(sols));
@assert(!isempty(errs));
q = finite(sols)[findmin(errs)[2]].solution;

# We compute the roots of the solution polynomial.
using PolynomialRoots
ξ = map(z -> z/norm(z), sort(roots([q; 1]), by=(-) ∘ angle))

# We compute the weights using a confluent Vandermonde system
V = transpose(hcat([vcat([[ξ[j]^k, k * ξ[j]^(k-1)] for j = 1:r]...) for k = 0:d]...));
L = V \ map(moment, 0:d);
λ₀ = L[1:2:end-1];
λ₁ = L[2:2:end];

# We compare the estimated recovered parameters to the original parameters
t_est = real(im * log.(ξ));
f₁_est = real([0; cumsum((λ₀ + s * λ₁ ./ ξ) .* (ξ.^s))]);
f₀_est = real([0; cumsum(-im * λ₁ .* (ξ.^(s-1))
                         - f₁_est[1:end-1] .* (t_est - [t_est; 0][2:end])
                        )]);
# total error
[norm(t - t_est), norm(f₀ - f₀_est), norm(f₁ - f₁_est)]

max_stars_repo_name: mwageringel/momentproblems

max_issues_repo_name: mwageringel/momentproblems

max_forks_repo_name: mwageringel/momentproblems


text:
########################################################################
function ppca_given_b(X, b; Q=1)
########################################################################

    D, N = size(X)

    @assert(length(b) == D)

    JITTER = 1e-8

    numparam = D*Q        +   1     +  1
    #    linear mapping    scale    noise

    @printf("Running PPCA_given_b for %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    #---------------------------------------------
    function marginalLogLikelihood(W, μ, σ)
    #---------------------------------------------

        sum(logpdf(MvNormal(μ, W*W' + σ*σ*I + JITTER*I), X))

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local W = reshape(param[1:end-2], D, Q)

        local s = exp(param[end-1])

        local σ = exp(param[end])

        return W, s, σ

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == numparam)

        local W, s, σ = unpack(param)

        return -1.0 * marginalLogLikelihood(W, s * b, σ)

    end

    #---------------------------------------------
    # Run optimiser
    #---------------------------------------------

    opt    = Optim.Options(show_trace = true, iterations = 10_000)

    result = optimize(objective, [randn(D*Q); randn()*3; randn()*3], LBFGS(), opt, autodiff=:forward)

    W, s, σ = unpack(result.minimizer)

    #---------------------------------------------
    # Return results
    #---------------------------------------------

    return W, s, σ

end


text:
module TestToeplitz
using LinearAlgebra
using CovarianceFunctions: durbin, durbin!, trench, trench!, levinson, levinson!
using ToeplitzMatrices
using Test

@testset "toeplitz" begin
    n = 16
    x = range(-1, 1, length = n+1)
    a_exp = @. exp(-abs(x[1]-x))
    a_rbf = @. exp(-abs(x[1]-x)^2) / 2 # this one can be ill-conditioned unless we add a diagonal, i.e. scale kernel down while keeping diagonal at 1
    a_rand = rand(n+1) / n # ensures diagonal dominance

    a_tuple = (a_rand, a_exp, a_rbf)
    for a in a_tuple
        a[1] = 1
        r = a[2:end]

        y = durbin(r)

        K = SymmetricToeplitz(vcat(1, r[1:end-1]))
        KM = Matrix(K)
        b = - (KM \ r)
        @test b ≈ y

        B = trench(r[1:end-1])
        invKM = inv(KM)
        @test B ≈ invKM

        a = a[1:end-1]
        r = a[2:end]
        b = randn(n)
        y = levinson(r, b)

        K = SymmetricToeplitz(vcat(1, r[1:end]))
        KM = Matrix(K)

        Kb = KM \ b
        @test Kb ≈ y
        @test Kb ≈ levinson(K, b)
    end
end

end # module


text:
struct JointFlow{F,R}
    f::F
    λref::R
end
hasrefresh(F::JointFlow) = true

waiting_time_ref(rng, F::JointFlow) = poisson_time(rng, F.λref)/length(F.f)
waiting_time_ref(F::JointFlow) = poisson_time(F.λref)/length(F.f)

function grad_correct!(y, x, F::JointFlow)
    for i in eachindex(y)
        y[i] += grad_correct(x[i], F[i])
    end
    y
end
Base.getindex(F::JointFlow, args...) = getindex(F.f, args...)

function smove_forward!(G, i, t, x, θ, t′, F::JointFlow)
    nhd = neighbours(G, i)
    for i in nhd
        x[i], θ[i] = move_forward(t[i], x[i], θ[i], t′, F[i])
        t[i] = t′
    end
    return t, x, θ
end
function smove_forward!(i::Int, t, x, θ, t′, F::JointFlow)
    x[i], θ[i] = move_forward(t[i], x[i], θ[i], t′, F[i])
    t[i] = t′
    return t, x, θ
end

function move_forward!(τ, t, x, θ, F::JointFlow)
    t′ = t + τ
    for i in eachindex(x)
        x[i], θ[i] = move_forward(t, x[i], θ[i], t′, F[i])
    end
    t += τ
    t, x, θ
end