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

using Markdown
using InteractiveUtils

# ╔═╡ b2a8e540-19af-11eb-0214-9d5f04ac2486
md"""
# Arrays: Slices and views
"""

# ╔═╡ 3df9d1c0-164c-11eb-335f-edb4366d2625
v = ones(Int, 10)

# ╔═╡ 78a80050-173c-11eb-2fe8-6b70bb6a2507
v[3]

# ╔═╡ 06176472-19b0-11eb-21ec-0744703214b6
v[3] = 2

# ╔═╡ 1c31ef66-19b0-11eb-3c7b-87e5b3b253ac
v[3:5]

# ╔═╡ 1b45c1cc-19b0-11eb-04d2-937a0b44df8f
v[3:5] = 4

# ╔═╡ 1b274daa-19b0-11eb-3d39-63379843e8f5
v[3:5] .= 4

# ╔═╡ 1b099b70-19b0-11eb-2cd2-b535b3df5c37
v

# ╔═╡ 05bc6914-19b0-11eb-0857-1320b5d269fc
v[5:7] = [100, 101, 102]

# ╔═╡ 0572949c-19b0-11eb-3861-8bb389451ba7
v

# ╔═╡ 051ef344-19b0-11eb-0215-1b67608a6e35
w = v[3:5]

# ╔═╡ e9c25c40-19b0-11eb-0e7d-5b8230a3990f
w[1] = 444

# ╔═╡ e9a11bd4-19b0-11eb-385e-4f08b8c3e34d
w

# ╔═╡ e97fd756-19b0-11eb-395d-4fd72a6b4998
v

# ╔═╡ e961047c-19b0-11eb-37f5-bdcb55b8c5bd
z = view(v, 3:5)

# ╔═╡ e94283e6-19b0-11eb-2fe0-c1799fe420af
z

# ╔═╡ e923aa8c-19b0-11eb-03b4-63258474e7d8
z[1] = 555

# ╔═╡ e903b1fa-19b0-11eb-2f85-99101781fcac
z

# ╔═╡ e8e52762-19b0-11eb-11f4-d373924283e3
v

# ╔═╡ e8c15eea-19b0-11eb-018e-eda8a769bebb
z .= 9

# ╔═╡ e87f1062-19b0-11eb-3a5e-d764456e1a0e
v

# ╔═╡ e85d9e00-19b0-11eb-3925-99918cddcbd4
typeof(w), typeof(z)

# ╔═╡ 62c234a2-19b2-11eb-357e-51f15ad594f7
md"""
## Nicer syntax for views
"""

# ╔═╡ 61d29104-19b2-11eb-02fe-4bac18bb7c78
z2 = @view v[3:5]

# ╔═╡ 3b7d1600-19b2-11eb-1fcc-af0792b5f3ea
typeof(z2)

# ╔═╡ 3b9a61c4-19b2-11eb-01b3-a3d54b983b27
z2 .= 1888

# ╔═╡ 3bdb17fc-19b2-11eb-0611-9db04a1299a0
v

# ╔═╡ 3bf85cc0-19b2-11eb-3a6b-15570698ac18
md"""
## Matrices: slices and views
"""

# ╔═╡ 3c11fa04-19b2-11eb-0065-91f886ef2445
M = [10i + j for i in 0:4, j in 1:4]

# ╔═╡ 3c2a697c-19b2-11eb-31fd-3d8d978ac8ab
M2 = [10i + j for i in 0:4 for j in 1:4]

# ╔═╡ 3c467bda-19b2-11eb-3900-1f7a1c9c723d
M[3:4, 1:2]

# ╔═╡ 3c63cbd6-19b2-11eb-3b2c-c72cf3fd8d2f
view(M, 3:5, 1:2)

# ╔═╡ 3c7753b8-19b2-11eb-200d-7b7b2c0358c7
@view M[3:4, 1:2]

# ╔═╡ 3c90fe26-19b2-11eb-2f44-256156778b3e
md"""
## Reshaping matrices
"""

# ╔═╡ 5385781c-19b4-11eb-1c42-c312f8e12946
M4 = reshape(M, 2, 10)

# ╔═╡ 53a3df6e-19b4-11eb-1de1-632c4c6d90a3
vv = vec(M)

# ╔═╡ 53bf61d0-19b4-11eb-2212-d91065dae28f

# ╔═╡ 53d99a46-19b4-11eb-21fa-859486e6289b

# ╔═╡ 53f647ea-19b4-11eb-2356-755781371a3a

# ╔═╡ 5413928c-19b4-11eb-1890-230a8077f4ae

# ╔═╡ 542f1244-19b4-11eb-207c-131b5699ad24

# ╔═╡ Cell order:
# ╟─b2a8e540-19af-11eb-0214-9d5f04ac2486
# ╠═3df9d1c0-164c-11eb-335f-edb4366d2625
# ╠═78a80050-173c-11eb-2fe8-6b70bb6a2507
# ╠═06176472-19b0-11eb-21ec-0744703214b6
# ╠═1c31ef66-19b0-11eb-3c7b-87e5b3b253ac
# ╠═1b45c1cc-19b0-11eb-04d2-937a0b44df8f
# ╠═1b274daa-19b0-11eb-3d39-63379843e8f5
# ╠═1b099b70-19b0-11eb-2cd2-b535b3df5c37
# ╠═05bc6914-19b0-11eb-0857-1320b5d269fc
# ╠═0572949c-19b0-11eb-3861-8bb389451ba7
# ╠═051ef344-19b0-11eb-0215-1b67608a6e35
# ╠═e9c25c40-19b0-11eb-0e7d-5b8230a3990f
# ╠═e9a11bd4-19b0-11eb-385e-4f08b8c3e34d
# ╠═e97fd756-19b0-11eb-395d-4fd72a6b4998
# ╠═e961047c-19b0-11eb-37f5-bdcb55b8c5bd
# ╠═e94283e6-19b0-11eb-2fe0-c1799fe420af
# ╠═e923aa8c-19b0-11eb-03b4-63258474e7d8
# ╠═e903b1fa-19b0-11eb-2f85-99101781fcac
# ╠═e8e52762-19b0-11eb-11f4-d373924283e3
# ╠═e8c15eea-19b0-11eb-018e-eda8a769bebb
# ╠═e87f1062-19b0-11eb-3a5e-d764456e1a0e
# ╠═e85d9e00-19b0-11eb-3925-99918cddcbd4
# ╟─62c234a2-19b2-11eb-357e-51f15ad594f7
# ╠═61d29104-19b2-11eb-02fe-4bac18bb7c78
# ╠═3b7d1600-19b2-11eb-1fcc-af0792b5f3ea
# ╠═3b9a61c4-19b2-11eb-01b3-a3d54b983b27
# ╠═3bdb17fc-19b2-11eb-0611-9db04a1299a0
# ╟─3bf85cc0-19b2-11eb-3a6b-15570698ac18
# ╠═3c11fa04-19b2-11eb-0065-91f886ef2445
# ╠═3c2a697c-19b2-11eb-31fd-3d8d978ac8ab
# ╠═3c467bda-19b2-11eb-3900-1f7a1c9c723d
# ╠═3c63cbd6-19b2-11eb-3b2c-c72cf3fd8d2f
# ╠═3c7753b8-19b2-11eb-200d-7b7b2c0358c7
# ╟─3c90fe26-19b2-11eb-2f44-256156778b3e
# ╠═5385781c-19b4-11eb-1c42-c312f8e12946
# ╠═53a3df6e-19b4-11eb-1de1-632c4c6d90a3
# ╠═53bf61d0-19b4-11eb-2212-d91065dae28f
# ╠═53d99a46-19b4-11eb-21fa-859486e6289b
# ╠═53f647ea-19b4-11eb-2356-755781371a3a
# ╠═5413928c-19b4-11eb-1890-230a8077f4ae
# ╠═542f1244-19b4-11eb-207c-131b5699ad24

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
# inject plane wave source
using PyPlot, ElasticWave
nz = 150; nx = 250; ext= 20  ; iflag = 1 ;
dx = 20.; dz = 20.; dt = 2e-3; tmax  = 4.; f0=10.;
vp = 2500.*ones(nz, nx)   ; vs = 1450.*ones(nz, nx); rho=2.5*ones(nz,nx);
vp[26:100, :] = 3000.; vs[26:100, :] = 1752.;
for i = 1 : 25
    vp[i+50, 75-i:75+i] = 3500.;
    vs[i+50, 175-i:175+i] = 2020.;
end
vp[101:125, :] = 4000.; vs[101:125, :] = 2309.;
vp[126:nz, :] = 4500.; vs[126:nz, :] = 2600.;
SeisPlot(vp, wbox=6.7, hbox=4, vmax=maximum(vp), vmin=minimum(vp), cmap="afmhot", dx=dx, dy=dz, name="/Users/wenyue/Desktop/vp.pdf")
SeisPlot(vs, wbox=6.7, hbox=4, vmax=maximum(vs), vmin=minimum(vs), cmap="afmhot", dx=dx, dy=dz, name="/Users/wenyue/Desktop/vs.pdf")
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

vp = 2500.*ones(nz, nx); vs = 1450.*ones(nz, nx); rho=2.5*ones(nz,nx);
fidMtx_dp = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

# inject plane wave, output records with removing direct wave
isx = collect(1:3:nx); isz = ones(Int64, length(isx));
flags = [false, false, true, true, false];
thetal = -pi/6; thetau = pi/6; np = 21; vps = vp[1];
irx = collect(1:nx); irz = ones(Int64, length(irx));
root = "/Users/wenyue/Desktop/plane_wave/"
p = collect(linspace(thetal, thetau, np));
tdind = zeros(Int64, length(isx))
ot    = zeros(Float64, length(isx))
for ip = 1 : np
    for is = 1 : length(isx)
        lx = 0.0
        if p[ip] > 0
           lx = (isx[is]-1)*dx
        else p[ip] < 0
             lx = (isx[end]-isx[is])*dx
        end
        tdind[is] = floor(Int64, abs(lx)*sin(abs(p[ip]))/vps/dt) + 1
        ot[is] = (tdind[is]-1)*dt
    end
    srcs = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
    shotv = MultiStepForward(irz, irx, srcs, fidMtx, tmax=tmax)
    shot_tmp = MultiStepForward(irz, irx, srcs, fidMtx_dp, tmax=tmax)
    shotv.Vx[:] = shotv.Vx - shot_tmp.Vx
    shotv.Vz[:] = shotv.Vz - shot_tmp.Vz
    path = join([root "shotv" "$ip" ".bin"])
    WriteShotV(path, shotv)
end

# produce source side wave field(vxx, vxz, vzx, vzz)
vps = modSmooth(vp, 20); vss = modSmooth(vs,20);
root1 = "/Users/wenyue/Desktop/plane_wave/SWFD/"
fidMtx_SM = CreateFidMtx(nz, nx, ext, iflag, vps, vss, rho, dz, dx, dt, f0);
p = collect(linspace(thetal, thetau, np));
tdind = zeros(Int64, length(isx))
ot    = zeros(Float64, length(isx))
for ip = 1 : np
    for is = 1 : length(isx)
        lx = 0.0
        if p[ip] > 0
           lx = (isx[is]-1)*dx
        else p[ip] < 0
             lx = (isx[end]-isx[is])*dx
        end
        tdind[is] = floor(Int64, abs(lx)*sin(abs(p[ip]))/vps[1]/dt) + 1
        ot[is] = (tdind[is]-1)*dt
    end
    srcs = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
    path1 = join([root1 "pv" "$ip" ".bin"])
    MultiStepForward(path1, srcs, fidMtx_SM, dz, dx, tmax=tmax)
end

root = "/Users/wenyue/Desktop/plane_wave/Shotv/"
root1= "/Users/wenyue/Desktop/plane_wave/SWFD/"
root2= "/Users/wenyue/Desktop/plane_wave/Image/"

for ip = 1 : np
    path = join([root "shotv" "$ip" ".bin"])
    path1= join([root1 "pv" "$ip" ".bin"])
    path2= join([root2 "image" "$ip" ".bin"])
    shotv = ReadShotV(path)
    (dm, du) = MultiStepAdjoint(shotv, fidMtx_SM, path1)
    # dm = laplaceFilter(dm)
    # du = laplaceFilter(du)
    WriteImage(path2, dm, du)
end

Ip = zeros(nz, nx);
Is = zeros(nz, nx);
for ip = 1 : np
    path2= join([root2 "image" "$ip" ".bin"])
    (dm, du) = ReadImage(path2)
    # dm = gec(dm); du = gec(du)
    # (dvp, dvs) = lambdamu2vpvs(dm, du, vps, vss, rho)
    Ip = Ip + dm
    Is = Is + du
end

tmax = 4.; dt = 2e-3; p = -pi/4; xs = collect(1:3:nx); dx=20.; v1 = 2500.;
path = "/Users/wenyue/Desktop/shotv_s_-45.bin";
root = "/Users/wenyue/Desktop/shotv/shotv";
tdind = planeSourceWfd(tmax, path, root, p, xs, dx, vp)
ot = (tdind-1)*dt; isx=collect(1:3:200); isz=ones(isx);
srcs = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
path_pv = "/Users/wenyue/Desktop//pv_-45.bin"
vps = modSmooth(vp, 20); vss = modSmooth(vs,20);
fidMtx_SM = CreateFidMtx(nz, nx, ext, iflag, vps, vss, rho, dz, dx, dt, f0);
MultiStepForward(path_pv, srcs, fidMtx_SM, dz, dx, tmax=tmax)
shotv_S = ReadShotV(path)
(dm, du) = MultiStepAdjoint(shotv_S, fidMtx_SM, path_pv);
figure(); a = quantile(abs(dm[:]),0.98); imshow(dm, vmax=a, vmin=-a, cmap="PuOr");
figure(); a = quantile(abs(du[:]),0.98); imshow(du, vmax=a, vmin=-a, cmap="PuOr");

function planeSrcShotV(tmax::Float64, path::ASCIIString, root::ASCIIString, p::Float64, xs::Array{Int64,1}, dx::Float64, vp::Float64)
    pathin = join([root "1" ".bin"])
    (nr, nz, nx, ext, iflag, irz, irx, dt, nst) = InfoShotV(pathin)
    nt = ceil(Int64, tmax/dt) + 1
    shotv_S = InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    tdind = zeros(Int64, length(xs))
    for ix = 1 : length(xs)
        lx = 0.0
        if p > 0
           lx = (xs[ix]-1)*dx
        else p < 0
           lx = (xs[end]-xs[ix])*dx
        end
        tdind[ix] = floor(Int64, abs(lx)*sin(abs(p))/vp/dt) + 1
    end
    for ix = 1 : length(xs)
        vx = zeros(nt,nr); vz = zeros(nt,nr);
        pathin = join([root "$ix" ".bin"])
        shotv  = ReadShotV(pathin)
        if tdind[ix]+nst-1 > nt
           tdend = nt
           nsend = nt - tdind[ix] + 1
        else
           tdend = tdind[ix] + nst - 1
           nsend = nst
        end
        vx[tdind[ix]:tdend,1:end] = shotv.Vx[1:nsend,1:end];
        vz[tdind[ix]:tdend,1:end] = shotv.Vz[1:nsend,1:end];
        shotv_S.Vx = shotv_S.Vx + vx
        shotv_S.Vz = shotv_S.Vz + vz
        WriteShotV(path, shotv_S)
    end
    return tdind
end

irx   = collect(1:1:nx); irz = 1*ones(Int64, length(irx));
flags = vec([false false true true false]);
xs = collect(1:3:200); root = "/Users/wenlei/Desktop/shotv"
for is = 1 : length(xs)
    isx = xs[is]; isz = 1; ot = 0.0;
    src = InitSource(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
    shotv = MultiStepForward(irz, irx, src, fidMtx, tmax=tmax);
    shotv_dp = MultiStepForward(irz, irx, src, fidMtx_dp, tmax=tmax);
    shotv.Vx[:] =shotv.Vx - shotv_dp.Vx
    shotv.Vz[:] =shotv.Vz - shotv_dp.Vz
    path = join([root "$is" ".bin"])
    WriteShotV(path, shotv)
    println("finished $is shot")
end
path = "/Users/wenyue/Desktop/shotv/shotv34.bin"
# shotv = ReadShotV(path);
# tmp=quantile(abs(vec(shotv.Vx)),0.95); imshow(shotv.Vx, aspect=0.1, vmax=tmp, vmin=-tmp, cmap="PuOr"); figure();
# tmp=quantile(abs(vec(shotv.Vz)),0.95); imshow(shotv.Vz, aspect=0.1, vmax=tmp, vmin=-tmp, cmap="PuOr")

function planeSrcShotV(tmax::Float64, path::ASCIIString, root::ASCIIString, p::Float64, xs::Array{Int64,1}, dx::Float64, vp::Float64)
    pathin = join([root "1" ".bin"])
    (nr, nz, nx, ext, iflag, irz, irx, dt, nst) = InfoShotV(pathin)
    nt = ceil(Int64, tmax/dt) + 1
    shotv_S = InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    tdind = zeros(Int64, length(xs))
    for ix = 1 : length(xs)
        lx = 0.0
        if p > 0
           lx = (xs[ix]-1)*dx
        else p < 0
           lx = (xs[end]-xs[ix])*dx
        end
        tdind[ix] = floor(Int64, abs(lx)*sin(abs(p))/vp/dt) + 1
    end
    for ix = 1 : length(xs)
        vx = zeros(nt,nr); vz = zeros(nt,nr);
        pathin = join([root "$ix" ".bin"])
        shotv  = ReadShotV(pathin)
        if tdind[ix]+nst-1 > nt
           tdend = nt
           nsend = nt - tdind[ix] + 1
        else
           tdend = tdind[ix] + nst - 1
           nsend = nst
        end
        vx[tdind[ix]:tdend,1:end] = shotv.Vx[1:nsend,1:end];
        vz[tdind[ix]:tdend,1:end] = shotv.Vz[1:nsend,1:end];
        shotv_S.Vx = shotv_S.Vx + vx
        shotv_S.Vz = shotv_S.Vz + vz
        WriteShotV(path, shotv_S)
    end
    return tdind
end

path_shotv = "/Users/wenyue/Desktop/shotv/shotv33.bin";
path_pv    = "/Users/wenyue/Desktop/pv33.bin";
ot = 0.0; isx=97; isz=1;
src = InitSource(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
MultiStepForward(path_pv, src, fidMtx_SM, dz, dx, tmax=tmax)
shotv = ReadShotV(path_shotv)
(dm, du) = MultiStepAdjoint(shotv, fidMtx_SM, path_pv);
figure(); a = quantile(abs(dm[:]),0.98); imshow(dm, vmax=a, vmin=-a, cmap="PuOr");
figure(); a = quantile(abs(du[:]),0.98); imshow(du, vmax=a, vmin=-a, cmap="PuOr");

ext: jl

max_stars_repo_name: WneleiGao/ElasticWave


text:
import AxisArrays: AxisArray, axisdim, axisnames, axisvalues, CategoricalVector
using JuMP.Containers: DenseAxisArray

tofloat(x::String) = parse(Float64, x)
tofloat(x) = float(x)

const typeparsers = Dict("float"=>tofloat, "string"=>Symbol,
                         "int"=>x->parse(Int, x),
                         "bool"=>let d=Dict('t'=>true, 'f'=>false); x->d[first(x)] end)
const typenames = Dict("float"=>Float64, "string"=>Symbol, "int"=>Int64, "bool"=>Bool)
const Scalar = Union{Float64,Symbol,Int64,Bool}

astype(typ::String, x::String) = typeparsers[typ](x)
astype(typ::String, x) = astype(typenames[typ], x)
astype(::String, ::Missing) = missing
astype(::DataType, ::Missing) = missing
astype(::Type{Symbol}, x) = Symbol(x)
astype(::Type{T}, x) where T = convert(Union{T,Missing}, x)

# !!Type-piracy!! We should make the case to move these definitions to JuMP
# Will be difficult as long as they use their one JuMPArray methods!
axisnames(::JuMP.Containers.DenseAxisArray{T,N,Ax}) where {T,N,Ax} = AxisArrays._axisnames(Ax)
axisnames(::Type{JuMP.Containers.DenseAxisArray{T,N,Ax}}) where {T,N,Ax} = AxisArrays._axisnames(Ax)
axisvalues(A::JuMP.Containers.DenseAxisArray) = axisvalues(A.axes...)

function axisdim(::Type{JuMP.Containers.DenseAxisArray{T,N,Ax,D}}, ::Type{<:Axis{name}}) where {T,N,Ax,name,D}
    isa(name, Int) && return name <= N ? name : error("axis $name greater than array dimensionality $N")
    names = axisnames(Ax)
    idx = findfirst(isequal(name), names)
    isnothing(idx) && error("axis $name not found in array axes $names")
    idx
end

axisdim(A::JuMP.Containers.DenseAxisArray, ax::Axis) = axisdim(A, typeof(ax))
@generated function axisdim(A::JuMP.Containers.DenseAxisArray, ax::Type{Ax}) where Ax<:Axis
    dim = axisdim(A, Ax)
    :($dim)
end

# AxisArray(A::JuMPArray) = AxisArray(A.innerArray, A.indexsets...)
AxisArray(A::JuMP.Containers.DenseAxisArray) = AxisArray(A.data, A.axes...)
JuMP.Containers.DenseAxisArray(A::AxisArray) = DenseAxisArray(A.data, AxisArrays.axes(A)...)

function _shiftamt(A, shifts::Pair{Symbol,T}...) where T<:Integer
    amt = zeros(T, ndims(A))
    for (ax, s) = shifts
        amt[axisdim(A, Axis{ax})] += s
    end
    amt
end

Base.circshift(A::AxisArray, shifts::Pair{Symbol,<:Integer}...) = circshift(A, _shiftamt(A, shifts...))
function Base.circshift(A::DenseAxisArray, shifts::Pair{Symbol,<:Integer}...)
    JuMP.Containers.DenseAxisArray(circshift(A.data, _shiftamt(A, shifts...)), A.axes...)
end

Base.indexin(a::AbstractArray, b::Axis) = indexin(a, b.val)
Base.findall(pred::Base.Fix2{typeof(in), <:AbstractArray}, b::Axis) = findall(pred, b.val)
Base.findfirst(pred::Base.Fix2{typeof(isequal)}, b::Axis) = findfirst(pred, b.val)

# vcat definition for categorical vector
Base.vcat(As::CategoricalVector...) = CategoricalVector(vcat(map(A -> A.data, As)...))

# function JuMP.constructvariable!(m::Model, _error::Function, lowerbound::Number, upperbound::AffExpr, args...; kwargs...)
#     v = constructvariable!(m, _error, lowerbound, Inf, args...; kwargs...)
#     addconstraint(m, constructconstraint!(addtoexpr(upperbound, -1.0, v), :(>=)))
#     v
# end

# function JuMP.constructvariable!(m::Model, _error::Function, lowerbound::AffExpr, upperbound::Number, args...; kwargs...)
#     v = constructvariable!(m, _error, -Inf, upperbound, args...; kwargs...)
#     addconstraint(m, constructconstraint!(addtoexpr(lowerbound, -1.0, v), :(<=)))
#     v
# end

# function JuMP.constructvariable!(m::Model, _error::Function, lowerbound::AffExpr, upperbound::AffExpr, args...; kwargs...)
#     v = constructvariable!(m, _error, -Inf, Inf, args...; kwargs...)
#     addconstraint(m, constructconstraint!(addtoexpr(lowerbound, -1.0, v), :(<=)))
#     addconstraint(m, constructconstraint!(addtoexpr(upperbound, -1.0, v), :(>=)))
#     v
# end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
"""
    clusterdurations, nobservations =
        computeclusterinfo(clusterdata::Vector{Matrix{Float32}})

Compute the durations of and number of member localizations in clusters.

# Description
This method computes the duration of each cluster and the number of
localizations in each cluster present in the input `clusterdata`.
"""
function computeclusterinfo(clusterdata::Vector{Matrix{Float32}})
    # Compute some quantities needed from each of the preclusters.
    nclusters = size(clusterdata)[1]
    clusterdurations = ones(nclusters)
    nobservations = ones(nclusters)
    for nn = 1:nclusters
        # Compute the duration of this cluster.
        currentframes = clusterdata[nn][:, 5]
        clusterdurations[nn] = maximum(currentframes) - minimum(currentframes) + 1

        # Compute the number of observations of this cluster excluding the
        # multiple observations of localizations in the same frame (that is,
        # ignore the generous pre-clustering artifact of multiple localizations
        # per frame, since this calculation assumes we have just one emitter).
        nobservations[nn] = length(unique(currentframes))
    end

    return clusterdurations, nobservations
end


text:
# Variable definitions
# mdl.z_Ol[i,j] = z_OpenLoop, open loop prediction of the state, i = state, j = step

# States in path following mode:
# i = 1 -> s
# i = 2 -> ey
# i = 3 -> epsi
# i = 4 -> v

# States in LMPC and system ID mode:
# i = 1 -> xDot
# i = 2 -> yDot
# i = 3 -> psiDot
# i = 4 -> ePsi
# i = 5 -> eY
# i = 6 -> s

function solveMpcProblem(mdl::MpcModel,mpcSol::MpcSol,mpcCoeff::MpcCoeff,mpcParams::MpcParams,trackCoeff::TrackCoeff,lapStatus::LapStatus,posInfo::PosInfo,modelParams::ModelParams,zCurr::Array{Float64},uPrev::Array{Float64})

    # Load Parameters
    sol_status::Symbol
    sol_u::Array{Float64,2}
    sol_z::Array{Float64,2}

    # Update current initial condition, curvature and System ID coefficients
    setvalue(mdl.z0,zCurr)
    setvalue(mdl.uPrev,uPrev)
    setvalue(mdl.c_Vx,mpcCoeff.c_Vx)            # System ID coefficients
    setvalue(mdl.c_Vy,mpcCoeff.c_Vy)
    setvalue(mdl.c_Psi,mpcCoeff.c_Psi)
    setvalue(mdl.coeff,trackCoeff.coeffCurvature)       # Track curvature
    setvalue(mdl.coeffTermCost,mpcCoeff.coeffCost)      # Terminal cost
    setvalue(mdl.coeffTermConst,mpcCoeff.coeffConst)    # Terminal constraints

    # Solve Problem and return solution
    sol_status  = solve(mdl.mdl)
    sol_u       = getvalue(mdl.u_Ol)
    sol_z       = getvalue(mdl.z_Ol)

    # export data
    mpcSol.a_x = sol_u[1,1]
    mpcSol.d_f = sol_u[1,2]
    mpcSol.u   = sol_u
    mpcSol.z   = sol_z
    mpcSol.solverStatus = sol_status
    mpcSol.cost = zeros(6)
    mpcSol.cost = [0,getvalue(mdl.costZTerm),getvalue(mdl.constZTerm),getvalue(mdl.derivCost),0,getvalue(mdl.laneCost)]
    #mpcSol.cost = [getvalue(mdl.costZ),0,0,getvalue(mdl.derivCost),0,0]

    # Print information
    # println("--------------- MPC START -----------------------------------------------")
    # println("z0             = $(zCurr')")
    # #println("z_Ol           = $(sol_z)")
    # #println("u_Ol           = $(sol_u)")
    # println("c_Vx           = $(mpcCoeff.c_Vx)")
    # println("c_Vy           = $(mpcCoeff.c_Vy)")
    # println("c_Psi          = $(getvalue(mdl.c_Psi))")
     println("ParInt         = $(getvalue(mdl.ParInt))")
    # #println("u_prev         = $(getvalue(mdl.uPrev))")
    println("Solved, status = $sol_status")
    # println("Predict. to s  = $(sol_z[end,6])")
    # println("costZ          = $(mpcSol.cost[1])")
    # println("termCost       = $(mpcSol.cost[2])")
    # println("termConst      = $(mpcSol.cost[3])")
    # println("derivCost      = $(mpcSol.cost[4])")
    # println("controlCost    = $(mpcSol.cost[5])")
    # println("laneCost       = $(mpcSol.cost[6])")
    # println("costZ: epsi    = $(norm(sol_z[:,4]))")
    # println("costZ: ey      = $(norm(sol_z[:,5]))")
    # println("costZ: v       = $(norm(sol_z[:,1]-0.8))")
    # println("--------------- MPC END ------------------------------------------------")
    nothing
end

function solveMpcProblem_pathFollow(mdl::MpcModel_pF,mpcSol::MpcSol,mpcParams_pF::MpcParams,trackCoeff::TrackCoeff,posInfo::PosInfo,
                                    modelParams::ModelParams,zCurr::Array{Float64},uPrev::Array{Float64},lapStatus::LapStatus)

    # Load Parameters
    coeffCurvature  = trackCoeff.coeffCurvature::Array{Float64,1}
    v_ref       = mpcParams_pF.vPathFollowing

    z_ref1 = cat(2,zeros(mpcParams_pF.N+1,3),v_ref*ones(mpcParams_pF.N+1,1))
    z_ref2 = cat(2,zeros(mpcParams_pF.N+1,1),0.2*ones(mpcParams_pF.N+1,1),zeros(mpcParams_pF.N+1,1),v_ref*ones(mpcParams_pF.N+1,1))
    z_ref3 = cat(2,zeros(mpcParams_pF.N+1,1),-0.1*ones(mpcParams_pF.N+1,1),zeros(mpcParams_pF.N+1,1),v_ref*ones(mpcParams_pF.N+1,1))

    sol_status::Symbol
    sol_u::Array{Float64,2}
    sol_z::Array{Float64,2}

    # Update current initial condition, curvature and previous input
    setvalue(mdl.z0,zCurr)
    setvalue(mdl.uPrev,uPrev)
    setvalue(mdl.coeff,coeffCurvature)
    if lapStatus.currentLap == 1
        setvalue(mdl.z_Ref,z_ref1)
    elseif lapStatus.currentLap == 2
        setvalue(mdl.z_Ref,z_ref1)
    elseif lapStatus.currentLap == 3
        setvalue(mdl.z_Ref,z_ref1)
    end

    # Solve Problem and return solution
    sol_status  = solve(mdl.mdl)
    sol_u       = getvalue(mdl.u_Ol)
    sol_z       = getvalue(mdl.z_Ol)

    mpcSol.a_x = sol_u[1,1]
    mpcSol.d_f = sol_u[1,2]
    mpcSol.u   = sol_u
    mpcSol.z   = sol_z
    mpcSol.solverStatus = sol_status
    mpcSol.cost = zeros(6)
    #mpcSol.cost = [getvalue(mdl.costZ),0,0,getvalue(mdl.derivCost),getvalue(mdl.controlCost),0]

    # Print information
    # println("--------------- MPC PF START -----------------------------------------------")
    # println("z0             = $(zCurr')")
    println("Solved, status = $sol_status")
    # println("Predict. to s  = $(sol_z[end,1])")
    # #println("uPrev          = $(uPrev)")
    # println("--------------- MPC PF END ------------------------------------------------")
    nothing
end

function solveMpcProblem_convhull(m::MpcModel_convhull,mpcSol::MpcSol,mpcCoeff::MpcCoeff,mpcParams::MpcParams,trackCoeff::TrackCoeff,lapStatus::LapStatus,
                                  posInfo::PosInfo,modelParams::ModelParams,zCurr::Array{Float64},uPrev::Array{Float64},selectedStates::SelectedStates)

 # Load Parameters
    sol_status::Symbol
    sol_u::Array{Float64,2}
    sol_z::Array{Float64,2}

    selStates       = selectedStates.selStates::Array{Float64,2}
    statesCost      = selectedStates.statesCost::Array{Float64,1}

    # Update current initial condition, curvature and System ID coefficients
    setvalue(m.z0,zCurr)
    setvalue(m.uPrev,uPrev)
    setvalue(m.c_Vx,mpcCoeff.c_Vx)            # System ID coefficients
    setvalue(m.c_Vy,mpcCoeff.c_Vy)
    setvalue(m.c_Psi,mpcCoeff.c_Psi)
    setvalue(m.coeff,trackCoeff.coeffCurvature)       # Track curvature
    setvalue(m.selStates,selStates)
    setvalue(m.statesCost,statesCost)

     # Solve Problem and return solution
    sol_status  = solve(m.mdl)
    sol_u       = getvalue(m.u_Ol)
    sol_z       = getvalue(m.z_Ol)

    # export data
    mpcSol.a_x = sol_u[1,1]
    mpcSol.d_f = sol_u[1,2]
    mpcSol.u   = sol_u
    mpcSol.z   = sol_z
    #mpcSol.eps_alpha = getvalue(m.eps_alpha)
    mpcSol.solverStatus = sol_status
    mpcSol.cost = zeros(6)
    mpcSol.cost = [0,getvalue(m.terminalCost),getvalue(m.controlCost),getvalue(m.derivCost),0,getvalue(m.laneCost)]

    mpcSol.costSlack = zeros(6)
    mpcSol.costSlack = [getvalue(m.slackVx),getvalue(m.slackVy),getvalue(m.slackPsidot),getvalue(m.slackEpsi),getvalue(m.slackEy),getvalue(m.slackS)]

    println("Solved, status = $sol_status")

    nothing
end

function solveMpcProblem_test(m::MpcModel_test,mpcSol::MpcSol,mpcCoeff::MpcCoeff,mpcParams::MpcParams,trackCoeff::TrackCoeff,lapStatus::LapStatus,
                                  posInfo::PosInfo,modelParams::ModelParams,zCurr::Array{Float64},uPrev::Array{Float64},selectedStates::SelectedStates)

 # Load Parameters
    sol_status::Symbol
    sol_u::Array{Float64,2}
    sol_z::Array{Float64,2}

    selStates       = selectedStates.selStates::Array{Float64,2}
    statesCost      = selectedStates.statesCost::Array{Float64,1}

    # Update current initial condition, curvature and System ID coefficients
    setvalue(m.z0,zCurr)
    setvalue(m.uPrev,uPrev)
    setvalue(m.c_Vx,mpcCoeff.c_Vx)            # System ID coefficients
    setvalue(m.c_Vy,mpcCoeff.c_Vy)
    setvalue(m.c_Psi,mpcCoeff.c_Psi)
    setvalue(m.coeff,trackCoeff.coeffCurvature)       # Track curvature
    setvalue(m.selStates,selStates)
    setvalue(m.statesCost,statesCost)

     # Solve Problem and return solution
    sol_status  = solve(m.mdl)
    sol_u       = getvalue(m.u_Ol)
    sol_z       = getvalue(m.z_Ol)

    # export data
    mpcSol.a_x = sol_u[1,1]
    mpcSol.d_f = sol_u[1,2]
    mpcSol.u   = sol_u
    mpcSol.z   = sol_z
    #mpcSol.eps_alpha = getvalue(m.eps_alpha)
    mpcSol.solverStatus = sol_status
    mpcSol.cost = zeros(6)
    #mpcSol.cost = [0,getvalue(m.terminalCost),getvalue(m.controlCost),getvalue(m.derivCost),0,getvalue(m.laneCost)]
    #println("o,terminal,control,deriv,0,lane= ",mpcSol.cost)

    println("Solved, status = $sol_status")

    nothing
end

function solveMpcProblem_obstacle(m::MpcModel_obstacle,mpcSol::MpcSol,mpcCoeff::MpcCoeff,mpcParams::MpcParams,trackCoeff::TrackCoeff,lapStatus::LapStatus,
                                  posInfo::PosInfo,modelParams::ModelParams,zCurr::Array{Float64},uPrev::Array{Float64},selectedStates::SelectedStates,
                                  obs_now::Array{Float64},obstacle::Obstacle)

 # Load Parameters
    sol_status::Symbol
    sol_u::Array{Float64,2}
    sol_z::Array{Float64,2}

    selStates       = selectedStates.selStates::Array{Float64,2}
    statesCost      = selectedStates.statesCost::Array{Float64,1}

    Q_obs           = mpcParams.Q_obs::Array{Float64,1}

#    println("Q_obs= ",Q_obs)

    obs      = zeros(mpcParams.N+1,3)
    obs[1,:] = obs_now

    # Compute the position of the obstacle  in the whole prediction horizon

    for i = 1:mpcParams.N
        obs[i+1,1] = obs[i,1] + modelParams.dt*i*obs[i,3]
        obs[i+1,2] = obs[i,2]
        obs[i+1,3] = obs[i,3]
    end

   # println("obstacle= ",obs)

    # Update current initial condition, curvature and System ID coefficients
    setvalue(m.z0,zCurr)
    setvalue(m.uPrev,uPrev)
    setvalue(m.c_Vx,mpcCoeff.c_Vx)            # System ID coefficients
    setvalue(m.c_Vy,mpcCoeff.c_Vy)
    setvalue(m.c_Psi,mpcCoeff.c_Psi)
    setvalue(m.coeff,trackCoeff.coeffCurvature)       # Track curvature
    setvalue(m.selStates,selStates)
    setvalue(m.statesCost,statesCost)
    setvalue(m.Q_obs,Q_obs)
    setvalue(m.obs,obs)

     # Solve Problem and return solution
    sol_status  = solve(m.mdl)
    sol_u       = getvalue(m.u_Ol)
    sol_z       = getvalue(m.z_Ol)

    # export data
    mpcSol.a_x = sol_u[1,1]
    mpcSol.d_f = sol_u[1,2]
    mpcSol.u   = sol_u
    mpcSol.z   = sol_z
    #mpcSol.eps_alpha = getvalue(m.eps_alpha)
    mpcSol.solverStatus = sol_status
    mpcSol.cost = zeros(6)
    mpcSol.cost = [0,getvalue(m.terminalCost),getvalue(m.controlCost),getvalue(m.derivCost),0,getvalue(m.laneCost)]

    mpcSol.costSlack = zeros(6)
    mpcSol.costSlack = [getvalue(m.slackVx),getvalue(m.slackVy),getvalue(m.slackPsidot),getvalue(m.slackEpsi),getvalue(m.slackEy),getvalue(m.slackS)]

    #println("cost for obstacle= ",getvalue(m.obstacleSlackCost))

    println("Solved, status = $sol_status")

    nothing
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
export disp, showdict, log, @log, setfilelogging, setlogfile, errorstring, savemv

disp(x...) = println(join(map(string,x),", "))

showdict(a::AbstractDict, args...; kargs...) = showdict(stdout, a, args..., kargs...)
function showdict(io::IO, a::AbstractDict, desc = ""; indent::Int = 0)
    nindent = 2
    if !isempty(desc)
        println(io, desc*":")
        indent += nindent
    end
    sk = sort(collect(keys(a)))
    for k in sk
        print(io, repeat(" ", indent),"$k: ")
        if isa(a[k], Dict)
            println(io, "")
            showdict(io, a[k], indent = indent + nindent)
        else
            println(io, isnil(a[k]) ? repr(a[k]) : "$(a[k])")
        end
    end
end

function setfilelogging(a::Bool)
    global LOGTOFILE
    LOGTOFILE = a
end

function setlogfile(a::AbstractString)
    global LOGFILE
    LOGFILE = a
end

import Base.log
log(a::AbstractString, args...; kargs...) = log(stdout, a, args...; kargs...)
function log(io::IO, a::AbstractString, args...; indent = 0, tofile = [], toSTDOUT = true)
    buf = IOBuffer()
    println(buf, Libc.strftime("%Y-%m-%d %T %z %Z", time()), "  |  ", repeat("  ", indent), join([a, args...], " "))
    str = String(take!(buf))

    toSTDOUT && println(io, str[end] == '\n' ? str[1:end-1] : str)
    if (tofile == [] && LOGTOFILE) || tofile == true || isa(tofile, AbstractString)
        filename = isa(tofile, AbstractString) ? tofile : LOGFILE
        @spawnat 1 open(filename, "a") do fid
            write(fid, str)
        end
    end
    nothing
end

macro log(a...)
    quote
        log($(a...))
    end
end

function errorstring(e)
    buf = IOBuffer()
    showerror(buf, e, catch_backtrace())
    String(take!(buf))
end

savemv(a, filename::String, f::FD.Callable, args...) = savemv(a, f, filename, args...)
function savemv(a, f::FD.Callable, filename::String, args...)
    tmpfilename = @p concat filename "." randstring(20) ".savemv.temp"
    f(a, tmpfilename, args...)
    mv(tmpfilename, filename, remove_destination = true)
    filename
end


text:
function broydenmix1!(vin, vout, beta, df, dv, cdf, iter, mixdim, brank)

    # [vnew,df,dv,cdf] =

    if (brank > mixdim)
        @printf("brank = %d, mixdim = %d\n", brank, mixdim)
        @printf("The rank of the Broyden update must be larger than the mixing memory size")
        exit()
    end
    Npoints = size(vin)[1]
    #
    # Note that the sign convention used in this function is different
    # from that used in other mixing functions!!!
    #
    # function evaluation overwrites vout
    #
    vout = vin - vout
    bragg = zeros(Float64,Npoints)
    vupd  = beta*vout

    if (iter > 1)
        # compute the changes in function evaluations and steps (changes
        # in potentials)
        if iter < mixdim + 1
            dv[:,iter-1] = -dv[:,iter-1] + vin
        else
            dv[:,mixdim] = -dv[:,mixdim] + vin
        end
        #
        # the difference between the current and the previous function value
        #
        df = -df + vout
        #
        # Construct C*df
        #
        yvec = beta*df
        if (iter > brank+1)
            @printf("Rank-%d update...\n",brank);
            for jbr = 1:iter-brank-1
                jsel = jbr:jbr+brank-1;
                dvmcdf  = dv[:,jsel] - cdf[:,jbr];
                G = dv[:,jsel]'*cdf[:,jbr]
                yvec = yvec + dvmcdf*(G\(dv[:,jsel]'*yvec))
            end

            if (brank > 1)
                cdf0 = cdf[:,jbr]
                cdf0 = cdf0[:,2:brank]
                cdf0 = cdf0 + dvmcdf*(G\(dv[:,jsel]'*cdf0))
                #cdf[:,jbr+1] = reshape([cdf0 yvec],n123*brank,1) ???
            else
                cdf[:,jbr+1] = yvec
            end
            #
            # Compute the potential update
            #
            for jbr = 1:iter-brank
                jsel = jbr:jbr+brank-1
                dvmcdf  = dv[:,jsel] - cdf[:,jbr]
                G = dv[:,jsel]'*cdf[:,jbr]
                vupd = vupd + dvmcdf*(G\(dv[:,jsel]'*vupd));
            end
        else
            #
            # perform rank-1 update until there are enough potential vectors to mix
            #
            @printf("Rank-1 update...\n");
            #
            # Save all initial df vectors in the first column of cdf
            #
            for jbr = 1:iter-2
                ibeg = (jbr-1)*Npoints + 1
                iend = jbr*Npoints
                dvmcdf = dv[:,jbr] - cdf[ibeg:iend,1]
                G = dv[:,jbr]'*cdf[ibeg:iend,1]
                yvec = yvec + dvmcdf*(G\(dv[:,jbr]'*yvec))
            end
            ibeg = Npoints*(iter-2)+1;
            iend = Npoints*(iter-1);
            cdf[ibeg:iend,1] = yvec;
            for jbr = 1:iter-1
                ibeg = (jbr-1)*Npoints + 1
                iend = jbr*Npoints
                dvmcdf = dv[:,jbr] - cdf[ibeg:iend,1]
                G = dv[:,jbr]'*cdf[ibeg:iend,1]
                vupd = vupd + dvmcdf*(G\(dv[:,jbr]'*vupd))
            end
        end
    end

    if iter < mixdim+1
        dv[:,iter] = vin
    else
        # delete the first column of dv
        #dv(:,1)=[];
        # append a new column at the end;
        dv[:,mixdim] = vin
    end
    df = vout
    vnew = vin - vupd

    return vnew

end


text:
module TestTimesteps

using Mimi
using Test

import Mimi:
    AbstractTimestep, FixedTimestep, VariableTimestep, TimestepVector,
    TimestepMatrix, TimestepArray, next_timestep, hasvalue, getproperty, Clock,
    time_index, get_timestep_array

#------------------------------------------------------------------------------
#  Test basic timestep functions and Base functions for Fixed Timestep
#------------------------------------------------------------------------------

t1 = FixedTimestep{1850, 10, 3000}(1)
@test is_first(t1)
@test t1 == TimestepIndex(1)
@test t1 == TimestepValue(1850)
@test TimestepIndex(1) == t1    # test both ways because to test both method definitions
@test t1 == TimestepValue(1850)
@test TimestepValue(1850) == t1 # test both ways because to test both method definitions
@test_throws ErrorException t1_prev = t1-1 #first timestep, so cannot get previous

t2 = next_timestep(t1)
@test t2.t == 2
@test t2 == TimestepIndex(2)
@test t2 == TimestepValue(1860)
@test t2 > TimestepIndex(1)
@test t2 > TimestepValue(1850)
@test t2 < TimestepIndex(3)
@test t2 < TimestepValue(1870)
@test_throws ErrorException t2_prev = t2 - 2 #can't get before first timestep

@test t2 == t1 + 1
@test t1 == t2 - 1

t3 = FixedTimestep{2000, 1, 2050}(51)
@test is_last(t3)
@test_throws ErrorException t3_next = t3 + 2 #can't go beyond last timestep

t4 = next_timestep(t3)
@test t4 == TimestepIndex(52)
@test t4 == TimestepValue(2051)
@test_throws ErrorException t_next = t4 + 1
@test_throws ErrorException next_timestep(t4)

#------------------------------------------------------------------------------
#  Test basic timestep functions and Base functions for Variable Timestep
#------------------------------------------------------------------------------
years = Tuple([2000:1:2024; 2025:5:2105])

t1 = VariableTimestep{years}()
@test is_first(t1)
@test t1 == TimestepIndex(1)
@test t1 == TimestepValue(2000)
@test_throws ErrorException t1_prev = t1-1 #first timestep, so cannot get previous

t2 = next_timestep(t1)
@test t2.t == 2
@test t2 == TimestepIndex(2)
@test t2 == TimestepValue(2001)
@test TimestepIndex(1) < t2
@test TimestepValue(2000) < t2
@test TimestepIndex(3) > t2
@test TimestepValue(2002) > t2
@test_throws ErrorException t2_prev = t2 - 2 #can't get before first timestep

@test t2 == t1 + 1
@test t1 == t2 - 1

t3 = VariableTimestep{years}(42)
@test is_last(t3)
@test ! is_first(t3)
@test_throws ErrorException t3_next = t3 + 2 #can't go beyond last timestep

t4 = next_timestep(t3)
@test t4 == TimestepIndex(43)

# note here that this comes back to an assumption made in variable
# timesteps that we may assume the next step is 1 year for the final year in TIMES
@test t4 == TimestepValue(2106)
@test_throws ErrorException t_next = t4 + 1
@test_throws ErrorException next_timestep(t4)

#------------------------------------------------------------------------------
#  Test some basic functions for TimestepIndex and TimestepValue
#------------------------------------------------------------------------------
start = 1
stop = 10
step = 2

# TimestepValue
@test TimestepValue(2000, offset = 1) + 1 == TimestepValue(2000, offset = 2)
@test TimestepValue(2000) + 1 == TimestepValue(2000, offset = 1)
@test TimestepValue(2000, offset = 1) - 1 == TimestepValue(2000)

# TimestepIndex
@test TimestepIndex(start):TimestepIndex(stop) == TimestepIndex.([start:stop...])
@test TimestepIndex(start):TimestepIndex(stop) == TimestepIndex(start):1:TimestepIndex(stop)
@test TimestepIndex(start):step:TimestepIndex(stop) == TimestepIndex.([start:step:stop...])

@test TimestepIndex(1) + 1 == TimestepIndex(2)
@test TimestepIndex(2) - 1 == TimestepIndex(1)

#------------------------------------------------------------------------------
#  Test a model with components with different offsets
#------------------------------------------------------------------------------

# we'll have Bar run from 2000 to 2010
# and Foo from 2005 to 2010

@defcomp Foo begin
    inputF = Parameter()
    output = Variable(index=[time])

    function run_timestep(p, v, d, ts)
        v.output[ts] = p.inputF + ts.t
    end
end

@defcomp Bar begin
    inputB = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, ts)
        if ts < TimestepValue(2005)
            v.output[ts] = p.inputB[ts]
        else
            v.output[ts] = p.inputB[ts] * ts.t
        end
    end
end

years = 2000:2010
first_foo = 2005

m = Model()
set_dimension!(m, :time, years)

# test that you can only add components with first/last within model's time index range
@test_throws ErrorException add_comp!(m, Foo; first=1900)
@test_throws ErrorException add_comp!(m, Foo; last=2100)

foo = add_comp!(m, Foo, first=first_foo)
bar = add_comp!(m, Bar)

set_param!(m, :Foo, :inputF, 5.)
set_param!(m, :Bar, :inputB, collect(1:length(years)))

run(m)

@test length(m[:Foo, :output]) == length(years)
@test length(m[:Bar, :output]) == length(years)

yr_dim = Mimi.Dimension(years)
idxs = yr_dim[first_foo]:yr_dim[years[end]]
foo_output = m[:Foo, :output]

offset = first_foo - years[1]
 for i in idxs
     @test foo_output[i] == 5+(i-offset) # incorporate offset into i now because we set ts.t to match component not model
 end

for i in 1:5
    @test m[:Bar, :output][i] == i
end

for i in 6:11
    @test m[:Bar, :output][i] == i*i
end

#------------------------------------------------------------------------------
#  test get_timestep_array
#------------------------------------------------------------------------------

m = Model()

#fixed timestep to start
set_dimension!(m, :time, 2000:2009)

vector = ones(5)
matrix = ones(3,2)

t_vector = get_timestep_array(m.md, Float64, 1, 1, vector)
t_matrix = get_timestep_array(m.md, Float64, 2, 1, matrix)

@test typeof(t_vector) <: TimestepVector
@test typeof(t_matrix) <: TimestepMatrix

# try with variable timestep
set_dimension!(m, :time, [2000:1:2004; 2005:2:2009])

t_vector = get_timestep_array(m.md, Float64, 1, 1, vector)
t_matrix = get_timestep_array(m.md, Float64, 2, 2, matrix)

@test typeof(t_vector) <: TimestepVector
@test typeof(t_matrix) <: TimestepMatrix

#------------------------------------------------------------------------------
#  Now build a model with connecting components
#------------------------------------------------------------------------------

@defcomp Foo2 begin
    inputF = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, ts)
        v.output[ts] = p.inputF[ts]
    end
end

m2 = Model()
set_dimension!(m2, :time, years)
bar = add_comp!(m2, Bar)
foo2 = add_comp!(m2, Foo2, first = first_foo)

set_param!(m2, :Bar, :inputB, collect(1:length(years)))
connect_param!(m2, :Foo2, :inputF, :Bar, :output)

run(m2)

foo_output2 = m2[:Foo2, :output][yr_dim[first_foo]:yr_dim[years[end]]]
for i in 1:6
    @test foo_output2[i] == (i+5)^2
end

#------------------------------------------------------------------------------
#  Connect them in the other direction
#------------------------------------------------------------------------------

@defcomp Bar2 begin
    inputB = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, ts)
        v.output[ts] = p.inputB[ts] * ts.t
    end
end

years = 2000:2010

m3 = Model()

set_dimension!(m3, :time, years)
add_comp!(m3, Foo, first=2005)
add_comp!(m3, Bar2)

set_param!(m3, :Foo, :inputF, 5.)
connect_param!(m3, :Bar2, :inputB, :Foo, :output, zeros(length(years)))

run(m3)

@test length(m3[:Foo, :output]) == 11
@test length(m3[:Bar2, :inputB]) == 11
@test length(m3[:Bar2, :output]) == 11

end


text:
using Test, AdvancedHMC, LinearAlgebra
using Statistics: mean, var, cov
include("common.jl")

θ_init = randn(D)
ϵ = 0.02
n_steps = 20
n_samples = 100_000
n_adapts = 2_000

@testset "HMC and NUTS" begin
    @testset "$(typeof(metric))" for metric in [
        UnitEuclideanMetric(D),
        DiagEuclideanMetric(D),
        DenseEuclideanMetric(D),
    ]
        h = Hamiltonian(metric, logπ, ∂logπ∂θ)
        @testset "$(typeof(prop))" for prop in [
            TakeLastProposal(Leapfrog(ϵ), n_steps),
            NUTS(Leapfrog(find_good_eps(h, θ_init))),
        ]
            samples = sample(h, prop, θ_init, n_samples; verbose=false)
            @test mean(samples[n_adapts+1:end]) ≈ zeros(D) atol=RNDATOL
            @testset "$(typeof(adaptor))" for adaptor in [
                PreConditioner(metric),
                NesterovDualAveraging(0.8, prop.integrator.ϵ),
                NaiveCompAdaptor(
                    PreConditioner(metric),
                    NesterovDualAveraging(0.8, prop.integrator.ϵ),
                ),
                StanNUTSAdaptor(
                    n_adapts,
                    PreConditioner(metric),
                    NesterovDualAveraging(0.8, prop.integrator.ϵ),
                ),
            ]
                samples = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; verbose=false)
                @test mean(samples[n_adapts+1:end]) ≈ zeros(D) atol=RNDATOL
            end
        end
    end
end


text:
using Random
using StatsBase
using Distributions
using LinearAlgebra
function lognormpdf(Y; mu = 0,sigma = 1)
    l = -0.5*((Y .- mu)/sigma).^2 .- 0.5*log(2*pi) .- log.(sigma)
    return l
end

function logsumexp(v)
    vm = maximum(v)
    lse = log(sum(exp.(v .- vm))) + vm
    return lse
end

function hugo_lognormpdf(Y,#::Array{Float64},
                    mu,#::Array{Float64},
                    sigma)#::Array{Float64})

    -0.5 * (  (Y-mu) / sigma )^2   - 0.5 * log(2.0*pi) - log(sigma)

end
#end

struct hugo_logsumexp
    #v::Array{Float64}

    function hugo_logsumexp(v)#::Array{Float64})

    vm = max(v)
    log(sum(exp(v-vm))) + vm
    end
end

function sim_data(k,n, T; doplot = false)
    # true values
    μ = randn(k)*10
    σ = abs.(randn(k))
    p = rand(k)
    p = p / sum(p)
    K_draw = sample(1:k, Weights(p), n, replace = true)
    Y = Matrix{Float64}(undef, n, T)
    for t in 1:T
        Y[:, t] = rand.(Normal.(μ[K_draw], σ[K_draw]))
    end

    return Dict(:y => Y, :μ => μ, :σ => σ, :p => p)
end

ed_data = sdata(5, 1000, 5)
ed_data[:y][:, 1]
using Plots

ed_data[:μ]
histogram(ed_data[:y][:, 1], bins = 50)

function expectation(Y, p_k, μ_k, σ_k)
    # Y = rand(10, 3)
    # p_k = [1/3, 1/3, 1/3]
    # μ_k = [1, 1, 1]
    # σ_k = [1, 2, 3]
    N = size(Y, 1)
    T = size(Y, 2)
    K = size(p_k, 2)

    LL = Matrix{Float64}(undef, (N, K))
    posterior = Matrix{Float64}(undef, (N, K))
    for k = 1:K
        LL[:, k] = sum(log.(p_k[:, k]) .+ lognormpdf(Y, mu = μ_k[k], sigma = σ_k[k]), dims = 2)
    end
    posterior[:,:] = exp.(LL .- logsumexp(LL))
    # posterior[:,:] = exp.(LL) ./ sum(exp.(LL), dims = 2)
    return posterior
end

## Hugo code

struct fake_hugo_data
    Y1::Vector
    Y2::Vector
    Y3::Vector
    μ0::Array
    π0::Array
    σ0::Array
    μ::Array
    π::Array
    σ::Array
    function fake_hugo_data()
    N = 10
    nk = 5
    tau = Array{Float64, 2}(undef, N, nk)
    lpm = Array{Float64, 2}(undef, N, nk)
    # Now I'm going to randomly generate my Y's from a mixture of Gaussians with a
    π0 = fill(1/(nk), ((nk),1))
    # This gives the initial values for the std dev of each Gaussian model
    σ0 = repeat(fill(1, ((nk),1)), 3)
        σ0 = reshape(σ0,3,(nk))
    # This gives the initial values for the means of each Gaussian model
    μ0 = repeat(collect(1:(nk)),inner = 3)
    μ0 = reshape(μ0, 3,(nk))
    # Simulating the Sample from the true DGP constructed from π0, μ0 and σ0
    Y = rand(MvNormal(μ0[1,:], I), 3*N)
        Y = π0' * Y
        Y = Y'

    Y1 = Y[1:N]
    Y2 = Y[N+1:2*N]
    Y3 = Y[2*N+1:3*N]

    # This denotes the initial probability assigned to each Gaussian model
    π = fill(1/nk, (nk,1))
    # This gives the initial values for the std dev of each Gaussian model
    σ = Array{Float64}(repeat(fill(1, (nk,1)),3))
        σ = reshape(σ,3,nk)
    # This gives the initial values for the means of each Gaussian model
    μ = Array{Float64}(repeat(collect(1:nk),inner = 3))
        μ = reshape(μ, 3,nk)
        new(Y1, Y2, Y3, μ0, π0, σ0, μ, π, σ )
    end
end

function Expectation(Y1,
                     Y2,
                     Y3,
                     μ,
                     σ,
                     τ,
                     N,
                     nk)

                     lnorm1 = Array{Float64, 2}(undef, N, nk)
                     lnorm2 = Array{Float64, 2}(undef, N, nk)
                     lnorm3 = Array{Float64, 2}(undef, N, nk)
                     lall = Array{Float64, 2}(undef, N, nk)
                     lik = Array{Float64, 2}(undef, N, nk)
                     lpm = Array{Float64, 2}(undef, N, nk)

    lτ = log.(τ)
    for i in 1:N
        lnorm1[i,:] = hugo_lognormpdf.(Y1[i], μ[1,:], σ[1,:])
        lnorm2[i,:] = hugo_lognormpdf.(Y2[i],μ[1,:],σ[1,:])
        lnorm3[i,:] = hugo_lognormpdf.(Y3[i],μ[1,:],σ[1,:])
    end
    lall = lτ + lnorm1 + lnorm2 + lnorm3
    lik = lik + hugo_logsumexp.(lall)

    for i in 1:N
        τ[i,:] = exp.(lall[i,:]) ./ sum(exp.(lall[i,:]))
    end

    return(τ,lik)

end

# TODO: This step is probably wrong?
function maximisation(Y, posterior)
    N = size(Y, 1)
    p_k = vec(sum(posterior, dims = 1) ./ N)
    μ = vec(sum(posterior .* mean(Y,dims = 2), dims = 1) ./ sum(posterior, dims = 1))
    σ = vec(sqrt.(sum(posterior .* (mean(Y, dims = 2) .- μ').^2, dims = 1 ) ./ sum(posterior, dims = 1)))
    return p_k, μ, σ
end

Y_test = randn(10, 3)
p_test = fill(1/3, 10, 3)
μ_test = [1, 2, 3]
σ_test = [1, 2, 3]

e_1 = expectation(Y_test, p_test,μ_test, σ_test)
p_1 = maximisation(Y_test, e_1)
p_1

function expectation_maximisation(Y, K)
    Y = fake_data[:y]
    K = 3
    μ = rand(K)
    p = rand(K)
    p = p ./ sum(p)
    σ = abs.(rand(K))
    diff = Inf
    iter = 0
    while diff > 1e-6
        iter += 1
        println("iter: $iter")
        old_p, old_μ, old_σ = p, μ, σ
        e_step = expectation(Y, p, μ, σ )
        p, μ, σ = maximisation(Y, e_step)
        p_diff = maximum(abs.(p - old_p))
        μ_diff = maximum(abs.(μ - old_μ))
        σ_diff = maximum(abs.(σ - old_σ))
        diff = maximum([p_diff,μ_diff, σ_diff])
        p
        old_p
        μ
        σ
    end
    return Dict(:p_k => p, :μ => μ, :σ => σ)
end

fake_data = sdata(3, 100, 3)
expectation_maximisation(fake_data[:y], 3)
fake_data

fake_hugo_df = fake_hugo_data()

fieldnames(fake_hugo_data)
hugo_τ = reshape(repeat(fake_hugo_df.π, 1000), 1000, 5)
hugo_post, hugo_lik  = Expectation(
    fake_hugo_df.Y1,
    fake_hugo_df.Y2,
    fake_hugo_df.Y3,
    fake_hugo_df.μ,
    fake_hugo_df.σ,
    hugo_τ,
    1000,
    5
    )

ed_Y = hcat(fake_hugo_df.Y1, fake_hugo_df.Y2, fake_hugo_df.Y3)

ed_post = expectation(ed_Y, hugo_τ, fake_hugo_df.μ[1, :], fake_hugo_df.σ[1, :])

hugo_post
ed_post


text:
tests = ["test1", "test2"]

print_with_color(:blue, "Running tests:\n")

for t in tests
	test_fn = "$t.jl"
	print_with_color(:green, "* $test_fn\n")
	include(test_fn)
end

max_stars_repo_name: xrounder/optunity

max_issues_repo_name: xrounder/optunity

max_forks_repo_name: xrounder/optunity


text:
# Automatically generated using Clang.jl

# Skipping MacroDefinition: paMakeVersionNumber ( major , minor , subminor ) ( ( ( major ) & 0xFF ) << 16 | ( ( minor ) & 0xFF ) << 8 | ( ( subminor ) & 0xFF ) )
# Skipping MacroDefinition: paNoDevice ( ( PaDeviceIndex ) - 1 )
# Skipping MacroDefinition: paUseHostApiSpecificDeviceSpecification ( ( PaDeviceIndex ) - 2 )
# Skipping MacroDefinition: paFloat32 ( ( PaSampleFormat ) 0x00000001 )
# Skipping MacroDefinition: paInt32 ( ( PaSampleFormat ) 0x00000002 )
# Skipping MacroDefinition: paInt24 ( ( PaSampleFormat ) 0x00000004 )
# Skipping MacroDefinition: paInt16 ( ( PaSampleFormat ) 0x00000008 )
# Skipping MacroDefinition: paInt8 ( ( PaSampleFormat ) 0x00000010 )
# Skipping MacroDefinition: paUInt8 ( ( PaSampleFormat ) 0x00000020 )
# Skipping MacroDefinition: paCustomFormat ( ( PaSampleFormat ) 0x00010000 )
# Skipping MacroDefinition: paNonInterleaved ( ( PaSampleFormat ) 0x80000000 )

const paFormatIsSupported = 0
const paFramesPerBufferUnspecified = 0

# Skipping MacroDefinition: paNoFlag ( ( PaStreamFlags ) 0 )
# Skipping MacroDefinition: paClipOff ( ( PaStreamFlags ) 0x00000001 )
# Skipping MacroDefinition: paDitherOff ( ( PaStreamFlags ) 0x00000002 )
# Skipping MacroDefinition: paNeverDropInput ( ( PaStreamFlags ) 0x00000004 )
# Skipping MacroDefinition: paPrimeOutputBuffersUsingStreamCallback ( ( PaStreamFlags ) 0x00000008 )
# Skipping MacroDefinition: paPlatformSpecificFlags ( ( PaStreamFlags ) 0xFFFF0000 )
# Skipping MacroDefinition: paInputUnderflow ( ( PaStreamCallbackFlags ) 0x00000001 )
# Skipping MacroDefinition: paInputOverflow ( ( PaStreamCallbackFlags ) 0x00000002 )
# Skipping MacroDefinition: paOutputUnderflow ( ( PaStreamCallbackFlags ) 0x00000004 )
# Skipping MacroDefinition: paOutputOverflow ( ( PaStreamCallbackFlags ) 0x00000008 )
# Skipping MacroDefinition: paPrimingOutput ( ( PaStreamCallbackFlags ) 0x00000010 )

struct PaVersionInfo
    versionMajor::Cint
    versionMinor::Cint
    versionSubMinor::Cint
    versionControlRevision::Cstring
    versionText::Cstring
end

const PaError = Cint

@cenum PaErrorCode::Int32 begin
    paNoError = 0
    paNotInitialized = -10000
    paUnanticipatedHostError = -9999
    paInvalidChannelCount = -9998
    paInvalidSampleRate = -9997
    paInvalidDevice = -9996
    paInvalidFlag = -9995
    paSampleFormatNotSupported = -9994
    paBadIODeviceCombination = -9993
    paInsufficientMemory = -9992
    paBufferTooBig = -9991
    paBufferTooSmall = -9990
    paNullCallback = -9989
    paBadStreamPtr = -9988
    paTimedOut = -9987
    paInternalError = -9986
    paDeviceUnavailable = -9985
    paIncompatibleHostApiSpecificStreamInfo = -9984
    paStreamIsStopped = -9983
    paStreamIsNotStopped = -9982
    paInputOverflowed = -9981
    paOutputUnderflowed = -9980
    paHostApiNotFound = -9979
    paInvalidHostApi = -9978
    paCanNotReadFromACallbackStream = -9977
    paCanNotWriteToACallbackStream = -9976
    paCanNotReadFromAnOutputOnlyStream = -9975
    paCanNotWriteToAnInputOnlyStream = -9974
    paIncompatibleStreamHostApi = -9973
    paBadBufferPtr = -9972
end

const PaDeviceIndex = Cint
const PaHostApiIndex = Cint

@cenum PaHostApiTypeId::UInt32 begin
    paInDevelopment = 0
    paDirectSound = 1
    paMME = 2
    paASIO = 3
    paSoundManager = 4
    paCoreAudio = 5
    paOSS = 7
    paALSA = 8
    paAL = 9
    paBeOS = 10
    paWDMKS = 11
    paJACK = 12
    paWASAPI = 13
    paAudioScienceHPI = 14
end

struct PaHostApiInfo
    structVersion::Cint
    type::PaHostApiTypeId
    name::Cstring
    deviceCount::Cint
    defaultInputDevice::PaDeviceIndex
    defaultOutputDevice::PaDeviceIndex
end

struct PaHostErrorInfo
    hostApiType::PaHostApiTypeId
    errorCode::Clong
    errorText::Cstring
end

const PaTime = Cdouble
const PaSampleFormat = Culong

struct PaDeviceInfo
    structVersion::Cint
    name::Cstring
    hostApi::PaHostApiIndex
    maxInputChannels::Cint
    maxOutputChannels::Cint
    defaultLowInputLatency::PaTime
    defaultLowOutputLatency::PaTime
    defaultHighInputLatency::PaTime
    defaultHighOutputLatency::PaTime
    defaultSampleRate::Cdouble
end

struct PaStreamParameters
    device::PaDeviceIndex
    channelCount::Cint
    sampleFormat::PaSampleFormat
    suggestedLatency::PaTime
    hostApiSpecificStreamInfo::Ptr{Cvoid}
end

const PaStream = Cvoid
const PaStreamFlags = Culong

struct PaStreamCallbackTimeInfo
    inputBufferAdcTime::PaTime
    currentTime::PaTime
    outputBufferDacTime::PaTime
end

const PaStreamCallbackFlags = Culong

@cenum PaStreamCallbackResult::UInt32 begin
    paContinue = 0
    paComplete = 1
    paAbort = 2
end

# Skipping Typedef: CXType_FunctionProto PaStreamCallback
# Skipping Typedef: CXType_FunctionProto PaStreamFinishedCallback

struct PaStreamInfo
    structVersion::Cint
    inputLatency::PaTime
    outputLatency::PaTime
    sampleRate::Cdouble
end


text:
function corrcov(C::AbstractMatrix)
    #should check if C is positive semidefinite

    sigma = sqrt.(diag(C))
    return C ./ (sigma*sigma')
end


text:
function sample_lambda!(Factors::LatentFactors, Ytil::Array{Float64, 2},
    resid::Residuals, genetic_effects::GeneticEffects, eig_ZAZ::SVD)
    #Sample factor loadings (Factors.Lambda) while marginalizing over residual
    #genetic effects: Y - Z_2W = FL' + E, vec(E)~N(0,kron(Psi_E,In) + kron(Psi_U, ZAZ^T))
    # note: conditioning on W, but marginalizing over U.
    #  sampling is done separately by trait because each column of Lambda is
    #  independent in the conditional posterior
    # note: eig_ZAZ has parameters that diagonalize aI + bZAZ for fast
    #  inversion: inv(aI + bZAZ) = 1/b*Ur*diag(1./(eta+a/b))*Ur'
    p=resid.p
    k=Factors.k

    Ur = eig_ZAZ.U;
    eta = eig_ZAZ.S;
    FtU = Factors.scores*Ur;
    UtY = Ur' * Ytil';

    Zlams = rand(Normal(0,1),k,p);
    for j = 1:p
       FUDi  = genetic_effects.ps[j] * broadcast(*, FtU, 1. ./ (eta' .+ genetic_effects.ps[j]/resid.ps[j]));
       means = FUDi * UtY[:,j];
       Qlam  = FUDi*FtU' + Diagonal(Factors.Plam[j,:]);
       Llam  = cholesky(Hermitian(Qlam)).L
       vlam  = Llam  \ means;
       mlam  = Llam' \ vlam;
       ylam  = Llam' \ Zlams[:,j];
       Factors.Lambda[j,:] = (ylam + mlam);
    end
end

function sample_means(Ytil::Array{Float64, 2}, Qt_Design::Array{Float64, 2}, N::Int64,
  resid::Residuals, random_precision::Array{Float64, 1}, svd_Design_Ainv::Dict)
    # when used to sample [B;D]:
    # Y - FL' - Z_2W = XB + ZD + E, vec(E)~N(0,kron(Psi_E,In)).
    # Note: conditioning on F, L and W.
    #  The vector [b_j;d_j] is sampled simultaneously. Each trait is sampled separately because their
    #  conditional posteriors factor into independent MVNs.
    # Note:svd_Design_Ainv has parameters to diagonalize mixed model equations for fast inversion:
    #  inv(a*blkdiag(fixed_effects.cov,Ainv) + b*[X; Z_1][X; Z_1]') = Q*diag(1./(a.*s1+b.*s2))*Q'
    # Qt_Design = Q'*Design, which doesn't change each iteration. Design = [X;Z_1]
    #
    # function also used to sample W:
    #  Y - FL' - XB - ZD = Z_2W + E, vec(E) ~ N(0,kron(Psi_E,In)).
    #  Here, conditioning is on B and D.

    p=resid.p;

    Q = svd_Design_Ainv["Q"];
    s1 = svd_Design_Ainv["s1"];
    s2 = svd_Design_Ainv["s2"];

    means = broadcast(*, Qt_Design * Ytil', resid.ps');
    location_sample = zeros(N,p);
    Zlams = randn(N,p);
    for j = 1:p,
        d = s1 * random_precision[j] + s2*resid.ps[j];
        mlam = broadcast(*, means[:,j], 1. ./ d);
        location_sample[:,j] = Q * (mlam + broadcast(*, Zlams[:,j], 1. ./ sqrt.(d)));
    end
    location_sample=location_sample';
    location_sample
end

function sample_h2s_discrete!(Factors::LatentFactors, eig_ZAZ::SVD)
    # sample factor heritibilties from a discrete set on [0,1)
    # prior places 50% of the weight at h2=0
    # samples conditional on F, marginalizes over U.

    Ur = eig_ZAZ.U;
    eta = eig_ZAZ.S;

    r = Factors.r;
    k = Factors.k;
    s = Factors.h2_divisions;

    log_ps = zeros(k,s);
    std_scores_b = Factors.scores*Ur;
    for i=1:s
        h2 = (i-1)/s;
        std_scores = Factors.scores;
        if h2 > 0
            std_scores = 1/sqrt(h2) * broadcast(*, std_scores_b, 1. ./ sqrt.(eta .+ (1-h2)/h2)');
            det = sum(log.((eta .+ (1-h2)/h2) * h2) / 2.);
        else
            det = 0;
        end
        log_ps[:,i] = sum(logpdf(Normal(0, 1), std_scores), dims = 2) .- det; #Prior on h2
        if i==1
            log_ps = log_ps .+ log(s-1);
        end
    end
    for j=1:k
        norm_factor = max.(log_ps[j,:]) .+ log(sum(exp.(log_ps[j,:] - max.(log_ps[j,:]))));
        ps_j = exp.(log_ps[j,:] - norm_factor);
        log_ps[j,:] = ps_j;
        Factors.h2[j] = sum(rand() .> cumsum(ps_j)) / s;
    end
end

function sample_Us!(Factors::LatentFactors, genetic_effects::GeneticEffects,
                    svd_ZZ_Ainv::Dict, Z_1::Array{Float64, 2})
    #samples genetic effects (U) conditional on the factor scores F:
    # F_i = U_i + E_i, E_i~N(0,s2*(h2*ZAZ + (1-h2)*I)) for each latent trait i
    # U_i = zeros(r,1) if h2_i = 0
    # it is assumed that s2 = 1 because this scaling factor is absorbed in
    # Lambda
    # svd_ZZ_Ainv has parameters to diagonalize a*Z_1*Z_1' + b*I for fast
    # inversion:

    Q = svd_ZZ_Ainv["Q"];
    s1 = svd_ZZ_Ainv["s1"];
    s2 = svd_ZZ_Ainv["s2"];

    k = Factors.k;
    n = genetic_effects.n;
    tau_e = 1. ./ (1 .- Factors.h2);
    tau_u = 1. ./ Factors.h2;
    b = Q' * Z_1 * broadcast(*, Factors.scores, tau_e)';
    z = randn(n,k);
    for j=1:k
        if tau_e[j]==1
            genetic_effects.U[j,:] = zeros(1,n);
        elseif isinf(tau_e[j])
            genetic_effects.U[j,:] = Factors.scores[j,:];
        else
            d = s2 * tau_u[j] + s1 * tau_e[j];
            mlam = broadcast(*, b[:,j], 1. ./ d);
            genetic_effects.U[j,:] = (Q*(mlam + broadcast(*, z[:,j], 1. ./ sqrt.(d))))';
        end
    end
end

function sample_factors_scores!(Ytil::Array{Float64, 2}, Factors::LatentFactors,
    resid::Residuals, genetic_effects::GeneticEffects, Z_1::Array{Float64, 2})
#Sample factor scores given factor loadings, U, factor heritabilities and
#phenotype residuals

    k = Factors.k;
    n = Factors.n;
    Lambda = Factors.Lambda;
    Lmsg = broadcast(*, Lambda, resid.ps);
    tau_e = reshape(1. ./ (1. .- Factors.h2), k);
    S = cholesky(Hermitian(Lambda' * Lmsg + Diagonal(tau_e))).L;
    Meta = S' \ (S \ (Lmsg' * Ytil + broadcast(*, genetic_effects.U * Z_1 , tau_e)));
    Factors.scores = Meta + S' \ randn(k,n);
end

function sample_delta!( Factors::LatentFactors, Lambda2_resid )
    #sample delta and tauh parameters that control the magnitudes of higher
    #index factor loadings.

    ad1 = Factors.ad1;
    ad2 = Factors.ad2;
    bd1 = Factors.bd1;
    bd2 = Factors.bd2;
    k = Factors.k;
    psijh = Factors.psijh;

    mat = broadcast(*, psijh, Lambda2_resid);
    n_genes = size(mat,1);
    ad = ad1 + 0.5*n_genes*k;
    bd = bd1 + 0.5*(1/Factors.delta[1]) * sum(Factors.tauh .* sum(mat, dims = 1)');
    Factors.delta[1] = rand(Gamma(ad,1. / bd));
    Factors.tauh = cumprod(Factors.delta);

    for h = 2:k
        ad = ad2 + 0.5*n_genes*(k-h+1);
        bd = bd2 + 0.5*(1. / Factors.delta[h])*sum(Factors.tauh[h:end].*sum(mat[:,h:end], dims = 1)');
        Factors.delta[h] = rand(Gamma(ad,1/bd));
        Factors.tauh = cumprod(Factors.delta);
    end
end

function  update_k!(Factors::LatentFactors, genetic_effects::GeneticEffects,
    b0::Float64, b1::Float64, i::Int64, epsilon::Float64, prop::Float64, Z_1::Array{Float64, 2} )
#adapt the number of factors by dropping factors with only small loadings
#if they exist, or adding new factors sampled from the prior if all factors
#appear important. The rate of adaptation decreases through the chain,
#controlled by b0 and b1

    df = Factors.df;
    ad2 = Factors.ad2;
    bd2 = Factors.bd2;
    p = Factors.p;
    k = Factors.k;
    r = Factors.r;
    gene_rows = 1:p;
    Lambda = Factors.Lambda;

    # probability of adapting
    prob = 1/exp(b0 + b1*i);
    uu = rand();

    # proportion of elements in each column less than eps in magnitude
    lind = mean(abs.(Lambda[1:p,:]) .< epsilon, dims = 1);
    vec = lind .>= prop; num = sum(vec);       # number of redundant columns

    Factors.num = num;
    Factors.no_f[i] = k-num;

    if uu < prob && i>200
        if  i > 20 && num == 0 && all(lind .< 0.995) && k < 2*p
            #add a column
            k=k+1;
            Factors.k = k;
            Factors.psijh = [Factors.psijh rand(Gamma(df/2, 2/df), p, 1)]
            Factors.delta = [Factors.delta; rand(Gamma(ad2, 1/bd2))];
            Factors.tauh = cumprod(Factors.delta);
            Factors.Plam = broadcast(*, Factors.psijh, Factors.tauh');
            Factors.Lambda = [Factors.Lambda randn(p,1) .* sqrt.(1 ./ Factors.Plam[:, k])];
            Factors.h2 = [Factors.h2; rand()];
            genetic_effects.U = [genetic_effects.U; randn(1,genetic_effects.n)];
            Factors.scores = [Factors.scores; genetic_effects.U[k,:]' * Z_1 + randn(1,r) .* sqrt(1-Factors.h2[k])];
        elseif num > 0      # drop redundant columns
            nonred = Array(1:k)[vec[:] .== 0]; # non-redundant loadings columns
            k = max(k - num, 1);
            Factors.k = k;
            Factors.Lambda = Lambda[:,nonred];
            Factors.psijh = Factors.psijh[:,nonred];
            Factors.scores = Factors.scores[nonred,:];
            for red = setdiff(1:k-1,nonred)
                #combine deltas so that the shrinkage of kept columns doesnt
                #decrease after dropping redundant columns
                Factors.delta[red+1] = Factors.delta[red+1] * Factors.delta[red];
            end
            Factors.delta = Factors.delta[nonred];
            Factors.tauh = cumprod(Factors.delta);
            Factors.Plam = broadcast(*, Factors.psijh, Factors.tauh');
            Factors.h2 = Factors.h2[nonred];
            genetic_effects.U = genetic_effects.U[nonred,:];
        end
    end
    Factors.nofout[i+1]=k;
end


text:
module test_testjulia07_pull

using Test
using TestJulia07

@test TestJulia07.pull isa Function

end


text:
@testset "746.min-cost-climbing-stairs.jl" begin
    @test min_cost_climbing_stairs([12, 10]) == 10
    @test min_cost_climbing_stairs([10, 15, 20]) == 15
    @test min_cost_climbing_stairs([1, 100, 1, 1, 1, 100, 1, 1, 100, 1]) == 6
end


text:
name: format-pr

on:
  schedule:
    - cron: '0 0 * * *'
jobs:
  build:
    runs-on: ${{ matrix.os }}
    strategy:
      matrix:
        julia-version: [1.3.0]
        julia-arch: [x86]
        os: [ubuntu-latest]
    steps:
      - uses: julia-actions/setup-julia@latest
        with:
          version: ${{ matrix.julia-version }}

      - uses: actions/checkout@v2
      - name: Install JuliaFormatter and format
        run: |
          julia  -e 'import Pkg; Pkg.add("JuliaFormatter")'
          julia  -e 'using JuliaFormatter; format(".", BlueStyle();verbose=true)'

      # https://github.com/marketplace/actions/create-pull-request
      - name: Create Pull Request
        uses: peter-evans/create-pull-request@v2
        with:
          token: ${{ secrets.GITHUB_TOKEN }}
          commit-message: Format .jl files
          title: 'Automatic JuliaFormatter.jl run'
          branch: auto-juliaformatter-pr
          labels: formatting, automated pr, no changelog
      - name: Check outputs
        run: |
          echo 'Pull Request Number - ${{ env.PULL_REQUEST_NUMBER }}'
          echo 'Pull Request Number - ${{ steps.cpr.outputs.pr_number }}'


text:
"""
    bioclim(layer::Integer; left=nothing, right=nothing, bottom=nothing, top=nothing)

Download and prepare bioclim layers from the CHELSA database, and returns
them as an array of `SimpleSDMPredictor`s. Layers are called by their number,
from 1 to 19. The list of available layers is given in a table below.

The keyword argument is `path`, which refers to the path where the function
will look for the geotiff files.

Note that these files are *large* due the fine resolution of the data, and for
this reason this function will return the *integer* version of the layers. Also
note that the bioclim data are only available for the V1 of CHELSA, and are not
from the V2.

It is recommended to *keep* the content of the `path` folder, as it will
eliminate the need to download the tiff files (which are quite large). For
example, calling `bioclim(1:19)` will download and everything, and future
calls will be much faster.

| Variable | Description                                                |
| ------   | ------                                                     |
| 1        | Annual Mean Temperature                                    |
| 2        | Mean Diurnal Range (Mean of monthly (max temp - min temp)) |
| 3        | Isothermality (BIO2/BIO7) (* 100)                          |
| 4        | Temperature Seasonality (standard deviation *100)          |
| 5        | Max Temperature of Warmest Month                           |
| 6        | Min Temperature of Coldest Month                           |
| 7        | Temperature Annual Range (BIO5-BIO6)                       |
| 8        | Mean Temperature of Wettest Quarter                        |
| 9        | Mean Temperature of Driest Quarter                         |
| 10       | Mean Temperature of Warmest Quarter                        |
| 11       | Mean Temperature of Coldest Quarter                        |
| 12       | Annual Precipitation                                       |
| 13       | Precipitation of Wettest Month                             |
| 14       | Precipitation of Driest Month                              |
| 15       | Precipitation Seasonality (Coefficient of Variation)       |
| 16       | Precipitation of Wettest Quarter                           |
| 17       | Precipitation of Driest Quarter                            |
| 18       | Precipitation of Warmest Quarter                           |
| 19       | Precipitation of Coldest Quarter                           |

"""
function bioclim(layer::Integer; left=nothing, right=nothing, bottom=nothing, top=nothing)
    return raster(SimpleSDMPredictor, BioClim(), layer=layer, left=left, right=right, bottom=bottom, top=top)
end

bioclim(layers::Vector{T}; args...) where {T <: Integer} = [bioclim(l; args...) for l in layers]
bioclim(layers::UnitRange{T}; args...) where {T <: Integer} = [bioclim(l; args...) for l in layers]


text:
module Solver

export solve, complete, isconditioned

import OmegaCore

import ..Var
using ..Tagging: Solve, tag
using ..Traits: trait
using ..Conditioning: Conditional

function isconditioned end

function solve end

# function Var.

# function Var.prehook(::trait(Solve), f, ω)
#   # Assume h has no specialization
#   solve(f, ω)
#     # then just continue as normal
#   # either the recurse itself
#   solve(f, ω)
# end

# Doesn't have specialization we'll get here
hasspecial(f, ω) = solve(f, ω)  ## don't want to do a prehook

# doesn't have specialization
solve(f, ω) = nothing

function Var.recurse(::trait(Solve), f, ω)
  solve(f, ω)
end

# Default case, no specialization:
# solve(f, ω) = Var.recurse(f, ω)

# # Defaualt case
# solve(f, id, ω) = f(ω)

function complete(x, ω = OmegaCore.defω())
  x(tag(ω, (solve = true,)))
end

# States:
# NoExpansion, always used specialised, don't expand
# Expand

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
#==============================================================================#
# MediaTailor.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/mediatailor-2018-04-23.normal.json
#==============================================================================#

__precompile__()

module MediaTailor

using AWSCore

"""
    using AWSSDK.MediaTailor.delete_playback_configuration
    delete_playback_configuration([::AWSConfig], arguments::Dict)
    delete_playback_configuration([::AWSConfig]; Name=)

    using AWSCore.Services.mediatailor
    mediatailor([::AWSConfig], "DELETE", "/playbackConfiguration/{Name}", arguments::Dict)
    mediatailor([::AWSConfig], "DELETE", "/playbackConfiguration/{Name}", Name=)

# DeletePlaybackConfiguration Operation

Deletes the configuration for the specified name.

# Arguments

## `Name = ::String` -- *Required*
The identifier for the configuration.

# Returns

`DeletePlaybackConfigurationResponse`

# Exceptions

.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mediatailor-2018-04-23/DeletePlaybackConfiguration)
"""
@inline delete_playback_configuration(aws::AWSConfig=default_aws_config(); args...) = delete_playback_configuration(aws, args)

@inline delete_playback_configuration(aws::AWSConfig, args) = AWSCore.Services.mediatailor(aws, "DELETE", "/playbackConfiguration/{Name}", args)

@inline delete_playback_configuration(args) = delete_playback_configuration(default_aws_config(), args)

"""
    using AWSSDK.MediaTailor.get_playback_configuration
    get_playback_configuration([::AWSConfig], arguments::Dict)
    get_playback_configuration([::AWSConfig]; Name=)

    using AWSCore.Services.mediatailor
    mediatailor([::AWSConfig], "GET", "/playbackConfiguration/{Name}", arguments::Dict)
    mediatailor([::AWSConfig], "GET", "/playbackConfiguration/{Name}", Name=)

# GetPlaybackConfiguration Operation

Returns the configuration for the specified name.

# Arguments

## `Name = ::String` -- *Required*
The identifier for the configuration.

# Returns

`GetPlaybackConfigurationResponse`

# Exceptions

.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mediatailor-2018-04-23/GetPlaybackConfiguration)
"""
@inline get_playback_configuration(aws::AWSConfig=default_aws_config(); args...) = get_playback_configuration(aws, args)

@inline get_playback_configuration(aws::AWSConfig, args) = AWSCore.Services.mediatailor(aws, "GET", "/playbackConfiguration/{Name}", args)

@inline get_playback_configuration(args) = get_playback_configuration(default_aws_config(), args)

"""
    using AWSSDK.MediaTailor.list_playback_configurations
    list_playback_configurations([::AWSConfig], arguments::Dict)
    list_playback_configurations([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.mediatailor
    mediatailor([::AWSConfig], "GET", "/playbackConfigurations", arguments::Dict)
    mediatailor([::AWSConfig], "GET", "/playbackConfigurations", <keyword arguments>)

# ListPlaybackConfigurations Operation

Returns a list of the configurations defined in AWS Elemental MediaTailor. You can specify a max number of configurations to return at a time. The default max is 50. Results are returned in pagefuls. If AWS Elemental MediaTailor has more configurations than the specified max, it provides parameters in the response that you can use to retrieve the next pageful.

# Arguments

## `MaxResults = ::Int`
Maximum number of records to return.

## `NextToken = ::String`
Pagination token returned by the GET list request when results overrun the meximum allowed. Use the token to fetch the next page of results.

# Returns

`ListPlaybackConfigurationsResponse`

# Exceptions

.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mediatailor-2018-04-23/ListPlaybackConfigurations)
"""
@inline list_playback_configurations(aws::AWSConfig=default_aws_config(); args...) = list_playback_configurations(aws, args)

@inline list_playback_configurations(aws::AWSConfig, args) = AWSCore.Services.mediatailor(aws, "GET", "/playbackConfigurations", args)

@inline list_playback_configurations(args) = list_playback_configurations(default_aws_config(), args)

"""
    using AWSSDK.MediaTailor.put_playback_configuration
    put_playback_configuration([::AWSConfig], arguments::Dict)
    put_playback_configuration([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.mediatailor
    mediatailor([::AWSConfig], "PUT", "/playbackConfiguration", arguments::Dict)
    mediatailor([::AWSConfig], "PUT", "/playbackConfiguration", <keyword arguments>)

# PutPlaybackConfiguration Operation

Adds a new configuration to AWS Elemental MediaTailor.

# Arguments

## `AdDecisionServerUrl = ::String`
The URL for the ad decision server (ADS). This includes the specification of static parameters and placeholders for dynamic parameters. AWS Elemental MediaTailor substitutes player-specific and session-specific parameters as needed when calling the ADS. Alternately, for testing you can provide a static VAST URL. The maximum length is 25000 characters.

## `CdnConfiguration = [ ... ]`
The configuration for using a content delivery network (CDN), like Amazon CloudFront, for content and ad segment management.
```
 CdnConfiguration = [
        "AdSegmentUrlPrefix" =>  ::String,
        "ContentSegmentUrlPrefix" =>  ::String
    ]
```

## `Name = ::String`
The identifier for the configuration.

## `SlateAdUrl = ::String`
The URL for a high-quality video asset to transcode and use to fill in time that's not used by ads. AWS Elemental MediaTailor shows the slate to fill in gaps in media content. Configuring the slate is optional for non-VPAID configurations. For VPAID, the slate is required because AWS Elemental MediaTailor provides it in the slots that are designated for dynamic ad content. The slate must be a high-quality asset that contains both audio and video.

## `VideoContentSourceUrl = ::String`
The URL prefix for the master playlist for the stream, minus the asset ID. The maximum length is 512 characters.

# Returns

`PutPlaybackConfigurationResponse`

# Exceptions

.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/mediatailor-2018-04-23/PutPlaybackConfiguration)
"""
@inline put_playback_configuration(aws::AWSConfig=default_aws_config(); args...) = put_playback_configuration(aws, args)

@inline put_playback_configuration(aws::AWSConfig, args) = AWSCore.Services.mediatailor(aws, "PUT", "/playbackConfiguration", args)

@inline put_playback_configuration(args) = put_playback_configuration(default_aws_config(), args)

end # module MediaTailor

#==============================================================================#
# End of file
#==============================================================================#


text:
@testset "OPF basics" begin
    g1 = MiniOPF.Generator("1", 0.1)
    g2 = MiniOPF.Generator("2", 0.2)
    l1 = MiniOPF.Line("a", ("1", "2"), 1.0)
    branches = [l1]
    bus_ids = ["1", "2"]
    branch_ids = ["a"]
    generators = [g1, g2]
    g = MiniOPF.Grid(bus_ids, branch_ids, generators, branches)
    load = MiniOPF.Load(["1", "2"], [1, 1])
    @test isempty(symdiff(load.bus_ids, bus_ids))
    @testset "DCOPF" begin
        q = MiniOPF.DCOPF(load, g)
        println(q)

        @testset "Test run_opf options" begin
            out = run_opf(q, options=Dict("foo" => "bar"))
            @test isempty(symdiff(keys(out), ["Type", "Options", "Solution"]))
            @test isa(out["Options"], Associative)
            @test haskey(out["Options"], "foo")
        end

        @testset "Default" begin
            out = run_opf(q)
            @test out["Options"] == "none"
        end
    end
end


text:
export something

something = "happy"


text:
#------------------------------------------------------------------------------#
#                                    Module                                    #
#------------------------------------------------------------------------------#

"""
# Description

    module aaHTC

AbstractAlgebra backend (aa) module to the Hermite Polynomial Tensor package.

# Usage

This module is `include`d by the `HermitePolyTensor` package and is not  usually
used in isolation, or directly; rather through the  interfaces  defined  on  the
`HPT` submodule, which is loaded automatically by the package:

```julia-repl
julia> using HermitePolyTensor

julia> typeof(HPT.aaHTC)
Module

```
"""
module aaHTC

#------------------------------------------------------------------------------#
#                                   Imports                                    #
#------------------------------------------------------------------------------#

import AbstractAlgebra

#------------------------------------------------------------------------------#
#                                   Exports                                    #
#------------------------------------------------------------------------------#

# Auxiliary Functions
export SVHP

# Tensor Component Functions
export HTC

#------------------------------------------------------------------------------#
#                             Auxiliary Functions                              #
#------------------------------------------------------------------------------#

"""
# Description

    SVHP(
        n::Int; ax::Int = 1, st::Bool = true
    )::AbstractAlgebra.Generic.MPoly{Rational{BigInt}}

Returns a symbolic `n`-th order single-variable Hermite polynomial in  the  axis
`ax`,  axis  1,  2,  and  3  corresponding  to  the  x,  y,  and  z  directions,
respectively. A "probabilist's" polynomial is returned if `st`  is  `true`  (the
default), otherwise, a "physicist's".

This function is not used directly; rather through the interfaces defined on the
`HPT` submodule, which is loaded automatically by the package, by specifying the
AbstractAlgebra ("aa") backend:

# Usage

```julia-repl
julia> using HermitePolyTensor

julia> HPT.SVHP(3, be = "aa")
x^3-3//1*x

julia> HPT.SVHP(4, ax = 2, be = "aa")
y^4-6//1*y^2+3//1

julia> HPT.aaHTC.SVHP(5, ax = 3, st = false) # Direct access
32//1*z^5-160//1*z^3+120//1*z

```
"""
function SVHP(
    n::Int; ax::Int = 1, st::Bool = true
)::AbstractAlgebra.Generic.MPoly{Rational{BigInt}}
    # Validations
    if n < 0
        throw(DomainError("n = $n outside the valid domain [0, ∞)"))
    elseif ax <= 0 || ax >= 4
        throw(DomainError("ax = $ax outside the valid domain [1, 3]"))
    end
    # Execution
    R, (x, y, z) = AbstractAlgebra.PolynomialRing(
        AbstractAlgebra.QQ, ["x", "y", "z"]
    )
    v = ax > 1 ? ax > 2 ? z : y : x
    po = 0 * v
    if st
        for m in 0:Int(floor(n/2))
            ex = n - 2m
            nm = factorial(big(n)) * (-1)^m
            dn = factorial(big(m)) * factorial(big(ex)) * 2^m
            po += nm//dn * v^ex
        end
    else
        for m in 0:Int(floor(n/2))
            ex = n - 2m
            nm = factorial(big(n)) * (-1)^m
            dn = factorial(big(m)) * factorial(big(ex))
            po += nm//dn * (2v)^ex
        end
    end
    po
end

#------------------------------------------------------------------------------#
#                          Tensor Component Functions                          #
#------------------------------------------------------------------------------#

"""
# Description

    HTC(
        id::String; D::Int = 2, st::Bool = true
    )::AbstractAlgebra.Generic.MPoly{Rational{BigInt}}

Returns a single symbolic *component* of a Hermite Polynomial Tensor.

The indices are specified by `id`, a `String` containing x's, y's,  and/or  z's;
invalid characters on `id` are ignored, while `id`  without  x's,  y's,  or  z's
represent a zero-rank  tensor  (whose  only  component  is  `1`  for  all  space
dimensions).

The Euclidean space dimensionality is specified  by  `D`.  The  tensor  rank  is
implied by the provided `id`, and is always the total  count  of  *valid*  index
characters, considering `D`.

A "probabilist's" polynomial is  returned  if  `st`  is  `true`  (the  default),
otherwise, a "physicist's".

# Usage

This function is not used directly; rather through the interfaces defined on the
`HPT` submodule, which is loaded automatically by the package, by specifying the
AbstractAlgebra ("aa") backend:

```julia-repl
julia> using HermitePolyTensor

julia> HPT.aaHTC.HTC("xxy")
x^2*y-y

```
"""
function HTC(
    id::String; D::Int = 2, st::Bool = true
)::AbstractAlgebra.Generic.MPoly{Rational{BigInt}}
    # Validations
    if D <= 0 || D >= 4
        throw(DomainError("dimension D = $D is outside the valid domain [1, 3]"))
    end
    # Execution
    xc = count(i -> (i == 'x'), id)
    yc = D > 1 ? count(i -> (i == 'y'), id) : 0
    zc = D > 2 ? count(i -> (i == 'z'), id) : 0
    ret = SVHP(xc, ax = 1, st = st)
    if yc != 0
        ret *= SVHP(yc, ax = 2, st = st)
    end
    if zc != 0
        ret *= SVHP(zc, ax = 3, st = st)
    end
    return ret
end

end # module


text:
# This file is part of Kpax3. License is MIT.

function test_fasta_data_processing_exceptions()
  missdna = UInt8['?', '*', '#', '-', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's', 'v', 'w', 'x', 'y', 'j', 'z']

  @test_throws Kpax3.KFASTAError Kpax3.readfasta("data/read_empty_file.fasta", false, missdna, 100000000, false, 0)
  @test_throws Kpax3.KFASTAError Kpax3.readfasta("data/read_no_id_char.fasta", false, missdna, 100000000, false, 0)
  @test_throws Kpax3.KFASTAError Kpax3.readfasta("data/read_no_1st_seq.fasta", false, missdna, 100000000, false, 0)
  @test_throws Kpax3.KFASTAError Kpax3.readfasta("data/read_no_nth_seq.fasta", false, missdna, 100000000, false, 0)

  @test_throws Kpax3.KFASTAError Kpax3.readfasta("data/read_utf8_id.fasta", false, missdna, 100000000, false, 0)
  @test_throws Kpax3.KFASTAError Kpax3.readfasta("data/read_utf8_seq.fasta", false, missdna, 100000000, false, 0)

  nothing
end

test_fasta_data_processing_exceptions()

function test_fasta_data_processing_read_blanks()
  missdna = UInt8['?', '*', '#', '-', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's', 'v', 'w', 'x', 'y', 'j', 'z']

  (data, id, refseq) = Kpax3.readfasta("data/read_blanks.fasta", false, missdna, 100000000, false, 0)
  @test data == UInt8['g' 'a';
                      'a' 'g';
                      'g' 'c';
                      'c' 'g';
                      'a' 'g';
                      'g' 'c']
  @test id == ["ID1", "ID5"]
  @test refseq == UInt8['a', 't', 'g', '.', '.', '.', 'g', '.', '.', 'a', '.', 'a']

  @test_throws Kpax3.KDomainError Kpax3.categorical2binary(data, UInt8(1), UInt8('?'))

  (bindata, val, key) = Kpax3.categorical2binary(data, UInt8(127), UInt8('?'))
  @test bindata == UInt8[0 1;
                         1 0;
                         1 0;
                         0 1;
                         0 1;
                         1 0;
                         1 0;
                         0 1;
                         1 0;
                         0 1;
                         0 1;
                         1 0]
  @test val == UInt8['a', 'g', 'a', 'g', 'c', 'g', 'c', 'g', 'a', 'g', 'c', 'g']
  @test key == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6]

  (data, id, refseq) = Kpax3.readfasta("data/read_blanks.fasta", false, missdna, 1, false, 0)
  @test data == UInt8['g' 'a';
                      'a' 'g';
                      'g' 'c';
                      'c' 'g';
                      'a' 'g';
                      'g' 'c']
  @test id == ["ID1", "ID5"]
  @test refseq == UInt8['a', 't', 'g', '.', '.', '.', 'g', '.', '.', 'a', '.', 'a']

  nothing
end

test_fasta_data_processing_read_blanks()

function test_fasta_data_processing_proper_dna_file()
  missdna = UInt8['?', '*', '#', '-', 'b', 'd', 'h', 'k', 'm', 'n', 'r', 's', 'v', 'w', 'x', 'y', 'j', 'z']

  (data, id, refseq) = Kpax3.readfasta("data/read_proper_nt.fasta", false, missdna, 100000000, false, 0)

  @test data == UInt8['t' 't' 't' 't' 't' 'c';
                      'g' 'g' 'a' 'g' 'a' 'a';
                      'a' 'a' '?' 'g' 'g' '?';
                      'g' 'g' 'g' 'c' 'c' 'c';
                      'c' 'g' 'a' 'c' 'g' 'c';
                      'a' 'a' 'a' 'g' 'g' 't';
                      'a' 'a' 't' 't' 'a' 'a';
                      'g' 'c' 'g' 'g' 'c' 'g']
  @test id == ["ID1", "ID2", "ID3", "ID4", "ID5", "ID6"]
  @test refseq == UInt8['a', '.', 'g', '.', '.', '.', 'g', '.', '.', '.', '.', 'a']

  @test_throws Kpax3.KDomainError Kpax3.categorical2binary(data, UInt8(1), UInt8('?'))

  (bindata, val, key) = Kpax3.categorical2binary(data, UInt8(127), UInt8('?'))
  @test bindata == UInt8[0 0 0 0 0 1;
                         1 1 1 1 1 0;
                         0 0 1 0 1 1;
                         1 1 0 1 0 0;
                         1 1 0 0 0 0;
                         0 0 0 1 1 0;
                         0 0 0 1 1 1;
                         1 1 1 0 0 0;
                         0 0 1 0 0 0;
                         1 0 0 1 0 1;
                         0 1 0 0 1 0;
                         1 1 1 0 0 0;
                         0 0 0 1 1 0;
                         0 0 0 0 0 1;
                         1 1 0 0 1 1;
                         0 0 1 1 0 0;
                         0 1 0 0 1 0;
                         1 0 1 1 0 1]
  @test val == UInt8['c', 't', 'a', 'g', 'a', 'g', 'c', 'g', 'a', 'c', 'g', 'a', 'g', 't', 'a', 't', 'c', 'g']
  @test key == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 8, 8]

  (data, id, refseq) = Kpax3.readfasta("data/read_proper_nt.fasta", false, missdna, 1, false, 0)
  @test data == UInt8['t' 't' 't' 't' 't' 'c';
                      'g' 'g' 'a' 'g' 'a' 'a';
                      'a' 'a' '?' 'g' 'g' '?';
                      'g' 'g' 'g' 'c' 'c' 'c';
                      'c' 'g' 'a' 'c' 'g' 'c';
                      'a' 'a' 'a' 'g' 'g' 't';
                      'a' 'a' 't' 't' 'a' 'a';
                      'g' 'c' 'g' 'g' 'c' 'g']
  @test id == ["ID1", "ID2", "ID3", "ID4", "ID5", "ID6"]
  @test refseq == UInt8['a', '.', 'g', '.', '.', '.', 'g', '.', '.', '.', '.', 'a']

  # consider all characters
  (data, id, refseq) = Kpax3.readfasta("data/read_proper_nt.fasta", false, zeros(UInt8, 1), 100000000, false, 0)
  @test data == UInt8['t' 't' 't' 't' 't' 'c';
                      'g' 'g' 'a' 'g' 'a' 'a';
                      'a' 'a' 'x' 'g' 'g' 'x';
                      'g' 'g' 'g' 'c' 'c' 'c';
                      'c' 'g' 'a' 'c' 'g' 'c';
                      'a' 'a' 'a' 'g' 'g' 't';
                      'a' 'a' 't' 't' 'a' 'a';
                      'g' 'c' 'g' 'g' 'c' 'g';
                      'a' 'x' 'a' 'a' 'a' 'a']
  @test id == ["ID1", "ID2", "ID3", "ID4", "ID5", "ID6"]
  @test refseq == UInt8['a', '.', 'g', '.', '.', '.', 'g', '.', '.', '.', '.', '.']

  @test_throws Kpax3.KDomainError Kpax3.categorical2binary(data, UInt8(1), UInt8('\0'))

  (bindata, val, key) = Kpax3.categorical2binary(data, UInt8(127), UInt8('\0'))
  @test bindata == UInt8[0 0 0 0 0 1;
                         1 1 1 1 1 0;
                         0 0 1 0 1 1;
                         1 1 0 1 0 0;
                         1 1 0 0 0 0;
                         0 0 0 1 1 0;
                         0 0 1 0 0 1;
                         0 0 0 1 1 1;
                         1 1 1 0 0 0;
                         0 0 1 0 0 0;
                         1 0 0 1 0 1;
                         0 1 0 0 1 0;
                         1 1 1 0 0 0;
                         0 0 0 1 1 0;
                         0 0 0 0 0 1;
                         1 1 0 0 1 1;
                         0 0 1 1 0 0;
                         0 1 0 0 1 0;
                         1 0 1 1 0 1;
                         1 0 1 1 1 1;
                         0 1 0 0 0 0]
  @test val == UInt8['c', 't', 'a', 'g', 'a', 'g', 'x', 'c', 'g', 'a', 'c', 'g', 'a', 'g', 't', 'a', 't', 'c', 'g', 'a', 'x']
  @test key == [1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9]

  nothing
end

test_fasta_data_processing_proper_dna_file()

function test_fasta_data_processing_proper_protein_file()
  misspro = UInt8['?', '*', '#', '-', 'b', 'j', 'x', 'z']

  (data, id, refseq) = Kpax3.readfasta("data/read_proper_aa.fasta", true, misspro, 100000000, false, 0)

  @test data == UInt8['k' 'k' 'k' 'k' 'k' 'e';
                      'k' 'k' 'i' 'k' 'i' 'i';
                      'l' 'l' '?' 'v' 'v' '?';
                      'l' 'l' 'l' 'v' 'v' 'v';
                      't' 'l' 'c' 't' 'l' 't';
                      'l' 'l' 'l' 'm' 'm' 'e';
                      'c' 'c' 'y' 'y' 'c' 'c';
                      'a' 'a' 't' 'a' 't' 't']
  @test id == ["ID1", "ID2", "ID3", "ID4", "ID5", "ID6"]
  @test refseq == UInt8['m', '.', 'a', '.', '.', '.', 'v', '.', '.', '.', '.', 'f']

  (bindata, val, key) = Kpax3.categorical2binary(data, UInt8(127), UInt8('?'))
  @test bindata == UInt8[0 0 0 0 0 1;
                         1 1 1 1 1 0;
                         0 0 1 0 1 1;
                         1 1 0 1 0 0;
                         1 1 0 0 0 0;
                         0 0 0 1 1 0;
                         1 1 1 0 0 0;
                         0 0 0 1 1 1;
                         0 0 1 0 0 0;
                         0 1 0 0 1 0;
                         1 0 0 1 0 1;
                         0 0 0 0 0 1;
                         1 1 1 0 0 0;
                         0 0 0 1 1 0;
                         1 1 0 0 1 1;
                         0 0 1 1 0 0;
                         1 1 0 1 0 0;
                         0 0 1 0 1 1]
  @test val == UInt8['e', 'k', 'i', 'k', 'l', 'v', 'l', 'v', 'c', 'l', 't', 'e', 'l', 'm', 'c', 'y', 'a', 't']
  @test key == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 8, 8]

  (data, id, refseq) = Kpax3.readfasta("data/read_proper_aa.fasta", true, misspro, 1, false, 0)
  @test data == UInt8['k' 'k' 'k' 'k' 'k' 'e';
                      'k' 'k' 'i' 'k' 'i' 'i';
                      'l' 'l' '?' 'v' 'v' '?';
                      'l' 'l' 'l' 'v' 'v' 'v';
                      't' 'l' 'c' 't' 'l' 't';
                      'l' 'l' 'l' 'm' 'm' 'e';
                      'c' 'c' 'y' 'y' 'c' 'c';
                      'a' 'a' 't' 'a' 't' 't']
  @test id == ["ID1", "ID2", "ID3", "ID4", "ID5", "ID6"]
  @test refseq == UInt8['m', '.', 'a', '.', '.', '.', 'v', '.', '.', '.', '.', 'f']

  # consider all characters
  (data, id, refseq) = Kpax3.readfasta("data/read_proper_aa.fasta", true, zeros(UInt8, 1), 100000000, false, 0)
  @test data == UInt8['k' 'k' 'k' 'k' 'k' 'e';
                      'k' 'k' 'i' 'k' 'i' 'i';
                      'l' 'l' 'x' 'v' 'v' 'x';
                      'l' 'l' 'l' 'v' 'v' 'v';
                      't' 'l' 'c' 't' 'l' 't';
                      'l' 'l' 'l' 'm' 'm' 'e';
                      'c' 'c' 'y' 'y' 'c' 'c';
                      'a' 'a' 't' 'a' 't' 't';
                      'f' 'x' 'f' 'f' 'f' 'f']
  @test id == ["ID1", "ID2", "ID3", "ID4", "ID5", "ID6"]
  @test refseq == UInt8['m', '.', 'a', '.', '.', '.', 'v', '.', '.', '.', '.', '.']

  (bindata, val, key) = Kpax3.categorical2binary(data, UInt8(127), UInt8('\0'))
  @test bindata == UInt8[0 0 0 0 0 1;
                         1 1 1 1 1 0;
                         0 0 1 0 1 1;
                         1 1 0 1 0 0;
                         1 1 0 0 0 0;
                         0 0 0 1 1 0;
                         0 0 1 0 0 1;
                         1 1 1 0 0 0;
                         0 0 0 1 1 1;
                         0 0 1 0 0 0;
                         0 1 0 0 1 0;
                         1 0 0 1 0 1;
                         0 0 0 0 0 1;
                         1 1 1 0 0 0;
                         0 0 0 1 1 0;
                         1 1 0 0 1 1;
                         0 0 1 1 0 0;
                         1 1 0 1 0 0;
                         0 0 1 0 1 1;
                         1 0 1 1 1 1;
                         0 1 0 0 0 0]
  @test val == UInt8['e', 'k', 'i', 'k', 'l', 'v', 'x', 'l', 'v', 'c', 'l', 't', 'e', 'l', 'm', 'c', 'y', 'a', 't', 'f', 'x']
  @test key == [1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9]

  nothing
end

test_fasta_data_processing_proper_protein_file()


text:
# # Louder

struct Louder{C}
    control::C
    by::Int64
end

"""
    IterationControl.louder(control, by=1)

Wrap `control` to make in more (or less) verbose. The same as
`control`, but as if the global `verbosity` were increased by the value
`by`.

"""
louder(c; by=1) = Louder(c, by)

# api:
done(d::Louder, state) = done(d.control, state)
update!(d::Louder, model, verbosity, args...) =
    update!(d.control, model, verbosity + d.by, args...)
takedown(d::Louder, verbosity, state) =
    takedown(d.control, verbosity + d.by, state)

# # WithStateDo

# helper:
shout(control, state) = @info "$(typeof(control)) state: $(flat(state))"

struct WithStateDo{C,F}
    control::C
    f::F
end

"""
    IterationControl.with_state_do(control,
                                  f=x->@info "\$(typeof(control)) state: \$x")

Wrap `control` to give access to it's internal state. Acts exactly
like `control` except that `f` is called on the internal state of
`control`. If `f` is not specified, the control type and state are
logged to `Info` at every update (useful for debugging new controls).

**Warning.** The internal state of a control is not yet considered
part of the public interface and could change between in any pre 1.0
release of IterationControl.jl.

"""
with_state_do(c; f=state->shout(c, state)) = WithStateDo(c, f)

# api:
for op in [:done, :takedown]
    eval(:($op(d::WithStateDo, args...) = $op(d.control, args...)))
end
@inline function update!(d::WithStateDo, args...)
    state = update!(d.control, args...)
    d.f(state)
    return state
end

# ## Skip

struct Skip{C,F<:Function}
    control::C
    predicate::F
end

_pred(predicate) = predicate
_pred(predicate::Int) = t -> mod(t, predicate) == 0

"""
    IterationControl.skip(control, predicate=1)

An iteration control wrapper.

If `predicate` is an **integer**, `k`: Apply `control` on every `k`
calls to apply the wrapped control, starting with the `k`th call.

If `predicate` is a **function**: Apply `control` as usual when
`predicate(n + 1)` is `true` but otherwise skip. Here `n` is the
number of control cycles applied so far.

"""
skip(control; predicate::Int=1) = Skip(control, _pred(predicate))

_state(s, model, verbosity, n, atomic_state...) = if s.predicate(n)
    new_atomic_state =
        update!(s.control, model, verbosity + 1, n, atomic_state...)
    return (atomic_state = new_atomic_state, n = n)
else
    return nothing
end

function update!(s::Skip, model, verbosity, n)
    state_candidate = _state(s, model, verbosity, n)
    state_candidate isa Nothing && return (n = n, )
    return state_candidate
end

# in case atomic state is not initialized in first `update` call:
function update!(s::Skip, model, verbosity, n, state::NamedTuple{(:n,)})
    state_candidate = _state(s, model, verbosity, n)
    state_candidate isa Nothing && return (n = n, )
    return state_candidate
end

# regular update:
function update!(s::Skip, model, verbosity, n, state)
    state_candidate = _state(s, model, verbosity, n, state.atomic_state)
    state_candidate isa Nothing &&
        return (atomic_state = state.atomic_state, n = n)
    return state_candidate
end

done(s::Skip, state) = done(s.control, state.atomic_state)

# can't be done if atomic state never intialized:
done(s::Skip, state::NamedTuple{(:n,)}) = false

takedown(s::Skip, verbosity, state) =
    takedown(s.control, verbosity, state.atomic_state)
takedown(::Skip, ::Any, ::NamedTuple{(:n,)}) = NamedTuple()


text:
p = 0.0:0.01:3.0

anim = @animate for ip ∈ p
    println(ip)
    S = simulate(waveDrivenHarmonicSim(X0 = [0.1, 0.0], parameter_profile_parameters = ((π,), (ip*π,), (π,), (π,))))
    plot(S, downsample=3, size=(1000, 1000), title="Ω = $ip", border=:none, grid=:none)
end every 1
gif(anim, fps = 20)


text:
struct Vertex{T, U<:Real}
  position::AbstractVector{T}
  value::U
end
Vertex(x::AbstractVector{T}, f::F) where {T, F<:Function} = Vertex(x, f(x))
value(v::Vertex) = v.value
position(v::Vertex) = v.position
newposition(a, ϵ::Number, b::Vertex) = a + ϵ * (a - b)

# must explicitly use <= and >= because == can't overridden and will
# be used in conjunction with < to create a <=
import Base: isless, +, -, <=, >=, isequal, isnan, hash
Base.isless(a::Vertex, b::Vertex) = value(a) < value(b)
Base.:<=(a::Vertex, b::Vertex) = value(a) <= value(b)
Base.:>=(a::Vertex, b::Vertex) = value(a) >= value(b)
Base.:+(a::Vertex, b) = position(a) .+ b
Base.:-(a::Vertex, b::Vertex) = position(a) .- position(b)
Base.:-(a, b::Vertex) = a .- position(b)
function Base.isequal(a::Vertex, b::Vertex)
  values_equal = value(a) == value(b) || (isnan(a) && isnan(b))
  positions_equal = all(position(a) .== position(b))
  return values_equal && positions_equal
end
Base.isnan(a::Vertex) = isnan(value(a))
Base.hash(v::Vertex) = hash(v, hash(:Vertex))
Base.hash(v::Vertex, h::UInt64) = hash(hash.(v.position), hash(v.value, h))