text:
using Revise
using OptiMode,GeometryPrimitives, ChainRules, Zygote, FiniteDifferences, DataFrames, CSV, BenchmarkTools
using Zygote: @adjoint, ignore, dropgrad #,StaticArrays,FFTW
##
w           =   1.7
t_core      =   0.7
θ           =   π/14.
edge_gap    =   0.5               # μm
n_core      =   2.4
n_subs      =   1.4
λ           =   1.55                  # μm
Δx          =   6.                    # μm
Δy          =   4.                    # μm
Δz          =   1.
Nx          =   128
Ny          =   128
Nz          =   1
ω           =   1 / λ
##
using Plots, BenchmarkTools
g = MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
s1 = ridge_wg(w,t_core,θ,edge_gap,n_core,n_subs,Δx,Δy)
s2 = circ_wg(0.7,0.2,edge_gap,n_core,n_subs,Δx,Δy)
s3 = test_shapes(1.05)

tanθ = tan(θ)
tcore_tanθ = t_core*tanθ
w_bottom = w + 2*tcore_tanθ
verts = 0.5.*   [   w     -w     -w_bottom    w_bottom
                    t_core   t_core    -t_core      -t_core    ]'
size(verts)

plot_ε(εₛ(s1,g;npix_sm=2),g.x,g.y)
plot_ε(εₛ(s2,g),g.x,g.y)
plot_ε(εₛ(s3,g;npix_sm=2),g.x,g.y)
# plot_ε(εₛ⁻¹(s3,g),g.x,g.y)

esm = εₛ(s5,g)
minimum(esm[2,1,:,:])

H,k = solve_k(ω,s1::Vector{<:GeometryPrimitives.Shape},g::MaxwellGrid)
H_ω²,ω² = solve_ω²(k,s1::Vector{<:GeometryPrimitives.Shape},g::MaxwellGrid)
@show λ_ω² = 1.0 / sqrt(ω²) # check that solve_ω² finds same neff as solve_k
plot_d⃗(H,k,g)
plot_d⃗(H_ω²,k,g)

H,ω² = solve_ω²(1.6,s1,g)
H,ω² = solve_ω²(1.6,s1,Δx,Δy,Δz,Nx,Ny,Nz)
ds = MaxwellData(1.6,Δx,Δy,Δz,Nx,Ny,Nz)
epsi = make_εₛ⁻¹(s1,ds.grid)
solve_ω²(1.6,epsi,ds)
solve_ω²(1.6,make_εₛ⁻¹(s1,g),g)
gg = make_MG(Δx,Δy,Δz,Nx,Ny,Nz)
solve_ω²(1.6,make_εₛ⁻¹(s1,gg),gg)
solve_nω(1.6,s1,Δx,Δy,Δz,Nx,Ny,Nz)

ks = collect(1.2:0.03:1.8)
function ngom_rwg(kz,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_nω(kz,ridge_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))
end
function ngoms_rwg(kz,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    [solve_nω(kk,ridge_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))[2] for kk in kz]
end
ngom_rwg1(kk,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) = ngom_rwg(kk,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1]
ngom_rwg2(kk,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) = ngom_rwg(kk,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[2]
ngom_rwg(1.6,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
ngom_rwg2(1.6,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
Zygote.gradient(ngom_rwg2,1.6,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
Zygote.gradient(ngom_rwg2,k,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
ns_k, ngs_k = ngom_rwg(ks,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)

#### Manual benchmark of ridge waveguide dispersion modeling @ 128 x 128 x 1 (6μm x 4μm) grid #####
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=2)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 1.099 s (2477656 allocations: 260.41 MiB)
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=5)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 2.674 s (2981906 allocations: 493.96 MiB)
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=10)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 5.590 s (3823588 allocations: 895.45 MiB)
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=15)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 7.024 s (4660434 allocations: 1.22 GiB)
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=20)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 9.324 s (5500140 allocations: 1.59 GiB)
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=25)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 11.875 s (6340834 allocations: 1.98 GiB)
# @btime ns_k, ngs_k = ngom_rwg($collect(range(1.2,1.8,length=30)),$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz) # 15.313 s (7191928 allocations: 2.46 GiB)
# Nk_ngom_rwg = [2,5,10,15,20,25,30]
# t_ngom_rwg = [1.099,2.674,5.590,7.024,9.324,11.875,15.313]
# mem_ngom_rwg = [0.26041,0.49396,0.89545,1.22,1.59,1.98,2.46]
# allocs_ngom_rwg = [2477656,2981906,3823588,4660434,5500140,6340834,7191928]
# plt_t = plot(Nk_ngom_rwg,t_ngom_rwg,label="time (sec)",legend=:bottomright,xlabel="# k values");
# plt_mem = plot(Nk_ngom_rwg,mem_ngom_rwg,label="memory (GB)",legend=:bottomright,xlabel="# k values");
# plt_allocs = plot(Nk_ngom_rwg,allocs_ngom_rwg,label="# allocations",legend=:bottomright,xlabel="# k values");
# sct_t = scatter!(plt_t,Nk_ngom_rwg,t_ngom_rwg,label=nothing); sct_mem = scatter!(plt_mem,Nk_ngom_rwg,mem_ngom_rwg,label=nothing); sct_allocs = scatter!(plt_allocs,Nk_ngom_rwg,allocs_ngom_rwg,label=nothing);
# plt_ngom_rwg = plot(plt_t,plt_mem,plt_allocs,sct_t,sct_mem,sct_allocs)

"""
################################################################################
#																			   #
#					check gradients of ridge waveguide solve				   #
#																			   #
################################################################################
"""
## generate ridge waveguide sweep + grad data
∂n_k = [gradient(ngom_rwg1,kk,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for kk in ks]
∂ng_k = [gradient(ngom_rwg2,kk,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for kk in ks]
∂n∂k_AD = [dnk[1] for dnk in ∂n_k]
∂ng∂k_AD = [dngk[1] for dngk in ∂ng_k]
# @btime n_ng_rwg($ωs,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 998.123 ms (2063797 allocations: 482.05 MiB) for 11 frequencies, 64×64 spatial grid

n_ng_w = [ngom_rwg(k,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for ww in ws]
n_ng_t = [ngom_rwg(k,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for tt in ts]
n_w = [ nng[1] for nng in n_ng_w]
ng_w = [ nng[2] for nng in n_ng_w]
n_t = [ nng[1] for nng in n_ng_t]
ng_t = [ nng[2] for nng in n_ng_t]
∂ns_w = [gradient(ngom_rwg1,k,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for ww in ws]
∂ns_t = [gradient(ngom_rwg1,k,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for tt in ts]
∂ngs_w = [gradient(ngom_rwg2,k,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for ww in ws]
∂ngs_t = [gradient(ngom_rwg2,k,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for tt in ts]
∂n∂ω_AD = [ dnng[1] for dnng in ∂ns_ω ]
∂n∂w_AD = [ dnng[2] for dnng in ∂ns_w ]
∂n∂t_AD = [ dnng[3] for dnng in ∂ns_t ]
∂ng∂ω_AD = [ dnng[1] for dnng in ∂ngs_ω ]
∂ng∂w_AD = [ dnng[2] for dnng in ∂ngs_w ]
∂ng∂t_AD = [ dnng[3] for dnng in ∂ngs_t ]
∂n_ng∂k_FD(kk) = central_fdm(5, 1)(x->[ngom_rwg(x,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],kk)
∂n_ng∂w_FD(ww) = central_fdm(5, 1)(x->[ngom_rwg(k,x,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],ww)
∂n_ng∂t_FD(tt) = central_fdm(5, 1)(x->[ngom_rwg(k,w,x,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],tt)
∂nng∂k_FD = hcat([∂n_ng∂k_FD(kk) for kk in ks]...)
∂nng∂t_FD = hcat([∂n_ng∂t_FD(tt) for tt in ts]...)
∂nng∂w_FD = hcat([∂n_ng∂w_FD(ww) for ww in ws]...)
∂n∂k_FD = ∂nng∂k_FD[1,:]
∂ng∂k_FD = ∂nng∂k_FD[2,:]
∂n∂t_FD = ∂nng∂t_FD[1,:]
∂ng∂t_FD = ∂nng∂t_FD[2,:]
∂n∂w_FD = ∂nng∂w_FD[1,:]
∂ng∂w_FD = ∂nng∂w_FD[2,:]

## plot results
p_∂n∂k = plot(ks,∂n∂k_FD,marker=:dot,markersize=3,label="∂n∂k,FD")
plot!(p_∂n∂k,ks,real(∂n∂k_AD),marker=:dot,markersize=3,label="∂n∂k,AD")

p_∂n∂w = plot(ws,∂n∂w_FD,marker=:dot,markersize=3,label="∂n∂w,FD")
plot!(p_∂n∂w,ws,real(∂n∂w_AD),marker=:dot,markersize=3,label="∂n∂w,AD")

p_∂n∂t = plot(ts,∂n∂t_FD,marker=:dot,markersize=3,label="∂n∂t,FD")
plot!(p_∂n∂t,ts,real(∂n∂t_AD),marker=:dot,markersize=3,label="∂n∂t,AD")

p_∂ng∂k = plot(ks,∂ng∂k_FD,marker=:dot,markersize=3,label="∂ng∂k,FD")
plot!(p_∂ng∂k,ks,real(∂ng∂k_AD),marker=:dot,markersize=3,label="∂ng∂k,AD")

p_∂ng∂t = plot(ts,∂ng∂t_FD,marker=:dot,markersize=3,label="∂ng∂t,FD")
plot!(p_∂ng∂t,ts,∂ng∂t_AD/1e4,marker=:dot,markersize=3,label="∂ng∂t,AD")

p_∂ng∂w = plot(ws,∂ng∂w_FD,marker=:dot,markersize=3,label="∂ng∂w,FD")
plot!(p_∂ng∂w,ws,∂ng∂w_AD/10,marker=:dot,markersize=3,label="∂ng∂w,AD")

# plot(p_∂n∂k, p_∂ng∂k, p_∂n∂w, p_∂ng∂w, p_∂n∂t, p_∂ng∂t )
plot(p_∂n∂k, p_∂n∂w, p_∂n∂t, p_∂ng∂k, p_∂ng∂w, p_∂ng∂t )
##
n_rwg, ng_rwg = ngom_rwg(ks,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
ω_rwg = ks ./ n_rwg
λ_rwg = 1.0 ./ ω_rwg
plot(λ_rwg,ng_rwg)
scatter!(λ_rwg,ng_rwg)
plot!(λ_rwg,n_rwg)
scatter!(λ_rwg,n_rwg)

plot(ks,ngs_rwg)
# using ChainRulesCore
# ChainRulesCore.refresh_rules()
# Zygote.refresh()

################################################################################
#                                                                              #
#                            Ridge Waveguide Tests                             #
#                                                                              #
################################################################################
ts = collect(0.4:0.04:0.9)
ws = collect(1.2:0.05:2.0)
ωs = collect(0.645*(0.9:0.02:1.1))

function solve_rwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_k(om,ridge_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))
end

using BenchmarkTools
# plot_d⃗(solve_cwg(1/1.55,0.7,0.2,edge_gap,1.7,1.5,Δx,Δy,Δz,Nx,Ny,Nz)...,g)
plot_d⃗(solve_rwg(1/1.55,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...,g)

@btime solve_cwg(1/1.55,0.7,0.2,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 357.346 ms (569288 allocations: 96.75 MiB) for 1 value, 64×48 spatial grid
# 3.121 s (2935840 allocations: 686.20 MiB) for 1 value, 128×128 spatial grid

@benchmark solve_rwg(1/1.55,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 58.9 ms (169302 allocations: 21.02 MiB) for 1 ω values, 32 x 32 spatial grid
# 731.258 ms (1780443 allocations: 252.59 MiB) for 1 ω values, 128 x 96 spatial grid
# 3.232 s (3276403 allocations: 1.16 GiB)   # for 11 ω values, 128 x 96 spatial grid
# 175.092 ms (456693 allocations: 62.81 MiB) for 1 value, 64×48 spatial grid
# 1.001 s (2362625 allocations: 322.87 MiB) for 1 value, 128×128 spatial grid
Ns = Int.(round.((2).^(5:0.2:5.8)))
outs = Array{Tuple{Float64,Float64}}(undef,length(Ns))
#ngom_rwg(1.6,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
function plot_bm_single(df::DataFrame)
    plt = plot(Int.(df.Nsq),
                    df.time,
                    xlabel="grid dimension N",
                    ylabel="time (s)",
                    labelfontsize=14,
                    scale=:log10,
                    )
    scatter!(plt,Int.(df.Nsq),df.time,label=nothing)
    savefig(plt,*(path,"time_",name,".png"))
    savefig(plt,*(path,"time_",name,".pdf"))
    n_err_norm = (df.n[begin:(end-1)] .- df.n[end]) ./  df.n[end]
    ng_err_norm = (df.ng[begin:(end-1)] .- df.ng[end]) ./  df.ng[end]
    plt2 = plot(Int.(df.Nsq[begin:(end-1)]),
                    n_err_norm,
                    xlabel="grid dimension N",
                    label="n",
                    labelfontsize=14,
                    xscale=:log10,
                    )
    scatter!(plt2,Int.(df.Nsq[begin:(end-1)]),n_err_norm,label=nothing)
    plot!(plt2,
            Int.(df.Nsq[begin:(end-1)]),
            ng_err_norm,
            xlabel="grid dimension N",
            label="ng",
            labelfontsize=14,
            xscale=:log10,
            )
    scatter!(plt2,Int.(df.Nsq[begin:(end-1)]),ng_err_norm,label=nothing)
    savefig(plt2,*(path,"convrg_",name,".png"))
    savefig(plt2,*(path,"convrg_",name,".pdf"))
    plt3 = plot(Int.(df.Nsq[begin:(end-1)]),
                    abs.(n_err_norm),
                    xlabel="grid dimension N",
                    label="n",
                    labelfontsize=14,
                    xscale=:log10,
                    )
    scatter!(plt3,Int.(df.Nsq[begin:(end-1)]),abs.(n_err_norm),label=nothing)
    plot!(plt3,
            Int.(df.Nsq[begin:(end-1)]),
            abs.(ng_err_norm),
            xlabel="grid dimension N",
            label="ng",
            labelfontsize=14,
            scale=:log10,
            )
    scatter!(plt3,Int.(df.Nsq[begin:(end-1)]),abs.(ng_err_norm),label=nothing)
    savefig(plt3,*(path,"absconvrg_",name,".png"))
    savefig(plt3,*(path,"absconvrg_",name,".pdf"))
end

plot_bm_single(df)

function rwg_ng_single_bm_N(Ns;name="rwg_ng_single_bm_N",path="/home/dodd/github/OptiMode/test/")
    bm = [ @benchmark(outs[$i] = ngom_rwg(1.6,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Ns[$i],$Ns[$i],1))  for i=1:length(Ns) ]
    ns = [o[1] for o in outs]
    ngs = [o[2] for o in outs]
    # bm = [ @benchmark(ngom_rwg(1.6,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$nn,$nn,1))  for nn=Ns ]
    df = DataFrame( N = Ns,
                    Nsq = (Ns.^2),
                    n = ns,
                    ng = ngs,
                    time = (time.(bm).*1e-9),
                    gctime_frac = (gctime.(bm) ./ time.(bm)),
                    mem = memory.(bm),
                    allocs = allocs.(bm),
                   )
    CSV.write(*(path,name,".csv"), df)
    plot_bm_single(df)
    df
end
rwg_ng_single_bm_N(Ns)

function plot_bm(df::DataFrame)
    plt = plot(Int.(df.Nsq),
                    df.time,
                    xlabel="grid dimension N",
                    ylabel="time (s)",
                    labelfontsize=14,
                    scale=:log10,
                    )
    scatter!(plt,Int.(df.Nsq),df.time,label=nothing)
    savefig(plt,*(path,"time_",name,".png"))
    savefig(plt,*(path,"time_",name,".pdf"))
    n_err_norm = (df.n[begin:(end-1)] .- df.n[end]) ./  df.n[end]
    ng_err_norm = (df.ng[begin:(end-1)] .- df.ng[end]) ./  df.ng[end]
    plt2 = plot(Int.(df.Nsq[begin:(end-1)]),
                    n_err_norm,
                    xlabel="grid dimension N",
                    label="n",
                    labelfontsize=14,
                    xscale=:log10,
                    )
    scatter!(plt2,Int.(df.Nsq[begin:(end-1)]),n_err_norm,label=nothing)
    plot!(plt2,
            Int.(df.Nsq[begin:(end-1)]),
            ng_err_norm,
            xlabel="grid dimension N",
            label="ng",
            labelfontsize=14,
            xscale=:log10,
            )
    scatter!(plt2,Int.(df.Nsq[begin:(end-1)]),ng_err_norm,label=nothing)
    savefig(plt2,*(path,"convrg_",name,".png"))
    savefig(plt2,*(path,"convrg_",name,".pdf"))
    plt3 = plot(Int.(df.Nsq[begin:(end-1)]),
                    abs.(n_err_norm),
                    xlabel="grid dimension N",
                    label="n",
                    labelfontsize=14,
                    xscale=:log10,
                    )
    scatter!(plt3,Int.(df.Nsq[begin:(end-1)]),abs.(n_err_norm),label=nothing)
    plot!(plt3,
            Int.(df.Nsq[begin:(end-1)]),
            abs.(ng_err_norm),
            xlabel="grid dimension N",
            label="ng",
            labelfontsize=14,
            scale=:log10,
            )
    scatter!(plt3,Int.(df.Nsq[begin:(end-1)]),abs.(ng_err_norm),label=nothing)
    savefig(plt3,*(path,"absconvrg_",name,".png"))
    savefig(plt3,*(path,"absconvrg_",name,".pdf"))
end

Ns = Int.(round.((2).^(5:0.2:5.8)))
ks = collect(1.2:0.3:1.8)
outs = Array{Tuple{Vector{Float64},Vector{Float64}}}(undef,length(Ns))
function rwg_ng_bm_N(ks,Ns;name="rwg_ng_bm_N",path="/home/dodd/github/OptiMode/test/")
    bm = [ @benchmark(outs[$i] = ngom_rwg($ks,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Ns[$i],$Ns[$i],1))  for i=1:length(Ns) ]
    ns = [o[1] for o in outs]
    ngs = [o[2] for o in outs]
    # bm = [ @benchmark(ngom_rwg(1.6,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$nn,$nn,1))  for nn=Ns ]
    df = DataFrame( N = Ns,
                    Nsq = (Ns.^2),
                    n = ns,
                    ng = ngs,
                    time = (time.(bm).*1e-9),
                    gctime_frac = (gctime.(bm) ./ time.(bm)),
                    mem = memory.(bm),
                    allocs = allocs.(bm),
                   )
    CSV.write(*(path,name,".csv"), df)
    # plot_bm_single(df)
    df
end
df = rwg_ng_bm_N(ks,Ns)

Ns = Int.(round.((2).^(5:0.5:7)))
ks = collect(1.2:0.03:1.8)
ngs = Array{Float64}(undef,length(Ns))
[ @benchmark(ngs[$i] = ngom_rwg(1.6,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Ns[$i],$Ns[$i],1))  for i=1:length(Ns) ]

ngs

name="rwg_ng_single_bm_N"
path="/home/dodd/github/OptiMode/test/"
df = CSV.read(*(path,name,".csv"),DataFrame)
plt = plot(Int.(df.Nsq),
                df.time,
                xlabel="grid dimension N",
                ylabel="time (s)",
                labelfontsize=14,
                scale=:log10,
                )
scatter!(plt,Int.(df.Nsq),df.time,label=nothing)
savefig(plt,*(path,name,".png"))
savefig(plt,*(path,name,".pdf"))

plt2 = plot(Int.(df.Nsq),
                df.n .- df.n[end],
                xlabel="grid dimension N",
                label="n",
                labelfontsize=14,
                xscale=:log10,
                )
scatter!(plt2,Int.(df.Nsq),df.n .- df.n[end],label=nothing)
plot!(plt2,
        Int.(df.Nsq),
        df.ng .- df.ng[end],
        xlabel="grid dimension N",
        label="ng",
        labelfontsize=14,
        xscale=:log10,
        )
scatter!(plt2,Int.(df.Nsq),df.ng .- df.ng[end],label=nothing)

plot!()

function rwg_bm_N(Ns;name="rwg_bm_N",path="/home/dodd/github/OptiMode/test/")
    bm = [ @benchmark(solve_rwg(1/1.55,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$nn,$nn,1))  for nn=Ns ]
    df = DataFrame( N = Ns,
                    Nsq = (Ns.^2),
                    time = (time.(bm).*1e-9),
                    gctime_frac = (gctime.(bm) ./ time.(bm)),
                    mem = memory.(bm),
                    allocs = allocs.(bm),
                   )
    CSV.write(*(path,name,".csv"), df)
    plt = plot(Int.(df.Nsq),
                    df.time,
                    xlabel="grid dimension N",
                    ylabel="time (s)",
                    labelfontsize=14,
                    scale=:log10,
                    )
    scatter!(plt,df.time,label=nothing)
    savefig(plt,*(path,name,".png"))
    savefig(plt,*(path,name,".pdf"))
end

function ng_rwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_n(om,ridge_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))[2]
end

function n_rwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_n(om,ridge_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))[1]
end

function n_ng_rwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_n(om,ridge_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))
end

# ng_rwg(1/1.55,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
# ng_rwg(collect(0.645*(0.9:0.02:1.1)),w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)

# @btime ng_rwg(0.645,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 50.866 ms (173716 allocations: 22.07 MiB) for 1 frequency, 32×32 spatial grid
# 198.690 ms (630899 allocations: 75.48 MiB) for 1 frequency, 64×64 spatial grid
# 956.631 ms (2527945 allocations: 322.88 MiB) for 1 frequency, 128×128 spatial grid
# 9.661 s (9835867 allocations: 1.23 GiB) for 1 frequency, 256×256 spatial grid

gradient(ng_rwg,0.645,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
gradient(ng_rwg,ωs,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)

# @btime gradient(ng_rwg,0.645,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 476.851 ms (1772737 allocations: 118.51 MiB) for 1 frequency, 32×32 spatial grid
# 2.087 s (6093797 allocations: 499.66 MiB) for 1 frequency, 64×64 spatial grid
# 9.447 s (22503426 allocations: 2.55 GiB) for 1 frequency, 128×128 spatial grid
# 60.865 s (86090986 allocations: 15.07 GiB) for 1 frequency, 256×256 spatial grid

n_ng_rwg(ωs,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
ns_ω, ngs_ω = n_ng_rwg(ωs,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
∂ns_ω = [gradient(n_rwg,om,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for om in ωs]
∂ngs_ω = [gradient(ng_rwg,om,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for om in ωs]
# @btime n_ng_rwg($ωs,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 998.123 ms (2063797 allocations: 482.05 MiB) for 11 frequencies, 64×64 spatial grid

n_ng_w = [n_ng_rwg(ω,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for ww in ws]
n_ng_t = [n_ng_rwg(ω,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for tt in ts]
n_w = [ nng[1] for nng in n_ng_w]
ng_w = [ nng[2] for nng in n_ng_w]
n_t = [ nng[1] for nng in n_ng_t]
ng_t = [ nng[2] for nng in n_ng_t]
∂ns_w = [gradient(n_rwg,ω,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for ww in ws]
∂ns_t = [gradient(n_rwg,ω,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for tt in ts]
∂ngs_w = [gradient(ng_rwg,ω,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for ww in ws]
∂ngs_t = [gradient(ng_rwg,ω,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for tt in ts]

∂n∂ω_AD = [ dnng[1] for dnng in ∂ns_ω ]
∂n∂w_AD = [ dnng[2] for dnng in ∂ns_w ]
∂n∂t_AD = [ dnng[3] for dnng in ∂ns_t ]

∂ng∂ω_AD = [ dnng[1] for dnng in ∂ngs_ω ]
∂ng∂w_AD = [ dnng[2] for dnng in ∂ngs_w ]
∂ng∂t_AD = [ dnng[3] for dnng in ∂ngs_t ]

∂n_ng∂ω_FD(om) = central_fdm(5, 1)(x->[n_ng_rwg(x,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],om)
∂n_ng∂w_FD(ww) = central_fdm(5, 1)(x->[n_ng_rwg(ω,x,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],ww)
∂n_ng∂t_FD(tt) = central_fdm(5, 1)(x->[n_ng_rwg(ω,w,x,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],tt)

∂nng∂ω_FD = hcat([∂n_ng∂ω_FD(om) for om in ωs]...)
∂nng∂t_FD = hcat([∂n_ng∂t_FD(tt) for tt in ts]...)
∂nng∂w_FD = hcat([∂n_ng∂w_FD(ww) for ww in ws]...)
∂n∂ω_FD = ∂nng∂ω_FD[1,:]
∂ng∂ω_FD = ∂nng∂ω_FD[2,:]
∂n∂t_FD = ∂nng∂t_FD[1,:]
∂ng∂t_FD = ∂nng∂t_FD[2,:]
∂n∂w_FD = ∂nng∂w_FD[1,:]
∂ng∂w_FD = ∂nng∂w_FD[2,:]

p_∂n∂ω = plot(ωs,∂n∂ω_FD,label="∂n∂ω,FD")
plot!(p_∂n∂ω,ωs,real(∂n∂ω_AD),label="∂n∂ω,AD")

p_∂n∂w = plot(ws,∂n∂w_FD,label="∂n∂w,FD")
plot!(p_∂n∂w,ws,real(∂n∂w_AD),label="∂n∂w,AD")

p_∂ng∂ω = plot(ωs,∂ng∂ω_FD,label="∂ng∂ω,FD")
plot!(p_∂ng∂ω,ωs,real(∂ng∂ω_AD),label="∂ng∂ω,AD")

p_∂ng∂t = plot(ts,∂ng∂t_FD,label="∂ng∂t,FD")
plot!(p_∂ng∂t,ts,∂ng∂t_AD/1e4,label="∂ng∂t,AD")

p_∂ng∂w = plot(ws,∂ng∂w_FD,label="∂ng∂w,FD")
plot!(p_∂ng∂w,ws,∂ng∂w_AD/1e4,label="∂ng∂w,AD")

################################################################################
#                                                                              #
#                            Circle Waveguide Tests                            #
#                                                                              #
################################################################################
t = 0.2
w = 0.7
ts = collect(0.4:0.04:0.9)
ws = collect(0.6:0.02:0.8)
ωs = collect(0.645*(0.9:0.02:1.1))
# ε

ε_pt_w_cwg(ww) = εₛ(circ_wg(ww,t,edge_gap,n_core,n_subs,Δx,Δy),g)[1,1,40,40,1]
∂ε_pt_w_cwg_AD = [gradient(ε_pt_w_cwg,ww) for ww in ws]
∂ε_pt_w_cwg_FD
plot(ws,ε_w_cwg)

function solve_cwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_k(om,circ_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))
end

function ng_cwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_n(om,circ_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))[2]
end

function n_cwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_n(om,circ_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))[1]
end

function n_ng_cwg(om,ww,tt_core,eedge_gap,nn_core,nn_subs,Dx,Dy,Dz,nx,ny,nz)
    solve_n(om,circ_wg(ww,tt_core,eedge_gap,nn_core,nn_subs,dropgrad(Dx),dropgrad(Dy)),dropgrad(Dx),dropgrad(Dy),dropgrad(Dz),dropgrad(nx),dropgrad(ny),dropgrad(nz))
end

# ng_cwg(1/1.55,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
# ng_cwg(collect(0.645*(0.9:0.02:1.1)),w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)

# @btime ng_cwg(0.645,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 50.866 ms (173716 allocations: 22.07 MiB) for 1 frequency, 32×32 spatial grid
# 198.690 ms (630899 allocations: 75.48 MiB) for 1 frequency, 64×64 spatial grid
# 956.631 ms (2527945 allocations: 322.88 MiB) for 1 frequency, 128×128 spatial grid
# 9.661 s (9835867 allocations: 1.23 GiB) for 1 frequency, 256×256 spatial grid

# @btime gradient(ng_cwg,0.645,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 476.851 ms (1772737 allocations: 118.51 MiB) for 1 frequency, 32×32 spatial grid
# 2.087 s (6093797 allocations: 499.66 MiB) for 1 frequency, 64×64 spatial grid
# 9.447 s (22503426 allocations: 2.55 GiB) for 1 frequency, 128×128 spatial grid
# 60.865 s (86090986 allocations: 15.07 GiB) for 1 frequency, 256×256 spatial grid

n_ng_cwg(ωs,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
ns_ω, ngs_ω = n_ng_cwg(ωs,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
∂ns_ω = [gradient(n_cwg,om,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for om in ωs]
∂ngs_ω = [gradient(ng_cwg,om,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for om in ωs]
# @btime n_ng_cwg($ωs,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Nx,$Ny,$Nz)
# 998.123 ms (2063797 allocations: 482.05 MiB) for 11 frequencies, 64×64 spatial grid

n_ng_w = [n_ng_cwg(ω,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for ww in ws]
n_ng_t = [n_ng_cwg(ω,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for tt in ts]
n_w = [ nng[1] for nng in n_ng_w]
ng_w = [ nng[2] for nng in n_ng_w]
n_t = [ nng[1] for nng in n_ng_t]
ng_t = [ nng[2] for nng in n_ng_t]
∂ns_w = [gradient(n_cwg,ω,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for ww in ws]
∂ns_t = [gradient(n_cwg,ω,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for tt in ts]
∂ngs_w = [gradient(ng_cwg,ω,ww,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for ww in ws]
∂ngs_t = [gradient(ng_cwg,ω,w,tt,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)[1:6] for tt in ts]

∂n∂ω_AD = [ dnng[1] for dnng in ∂ns_ω ]
∂n∂w_AD = [ dnng[2] for dnng in ∂ns_w ]
∂n∂t_AD = [ dnng[3] for dnng in ∂ns_t ]

∂ng∂ω_AD = [ dnng[1] for dnng in ∂ngs_ω ]
∂ng∂w_AD = [ dnng[2] for dnng in ∂ngs_w ]
∂ng∂t_AD = [ dnng[3] for dnng in ∂ngs_t ]

∂n_ng∂ω_FD(om) = central_fdm(5, 1)(x->[n_ng_cwg(x,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],om)
∂n_ng∂w_FD(ww) = central_fdm(5, 1)(x->[n_ng_cwg(ω,x,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],ww)
∂n_ng∂t_FD(tt) = central_fdm(5, 1)(x->[n_ng_cwg(ω,w,x,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)...],tt)

∂nng∂ω_FD = hcat([∂n_ng∂ω_FD(om) for om in ωs]...)
∂nng∂t_FD = hcat([∂n_ng∂t_FD(tt) for tt in ts]...)
∂nng∂w_FD = hcat([∂n_ng∂w_FD(ww) for ww in ws]...)
∂n∂ω_FD = ∂nng∂ω_FD[1,:]
∂ng∂ω_FD = ∂nng∂ω_FD[2,:]
∂n∂t_FD = ∂nng∂t_FD[1,:]
∂ng∂t_FD = ∂nng∂t_FD[2,:]
∂n∂w_FD = ∂nng∂w_FD[1,:]
∂ng∂w_FD = ∂nng∂w_FD[2,:]

p_∂n∂ω = plot(ωs,∂n∂ω_FD,label="∂n∂ω,FD")
plot!(p_∂n∂ω,ωs,real(∂n∂ω_AD),label="∂n∂ω,AD")

N = 64 * 64 * 2

p_∂n∂w = plot(ws,∂n∂w_FD,label="∂n∂w,FD")
plot!(p_∂n∂w,ws,real(∂n∂w_AD/N),label="∂n∂w,AD")

p_∂n∂t = plot(ts,∂n∂t_FD,label="∂n∂t,FD")
plot!(p_∂n∂t,ts,real(∂n∂t_AD/N),label="∂n∂t,AD")

p_∂ng∂ω = plot(ωs,∂ng∂ω_FD,label="∂ng∂ω,FD")
plot!(p_∂ng∂ω,ωs,real(∂ng∂ω_AD),label="∂ng∂ω,AD")

p_∂ng∂t = plot(ts,∂ng∂t_FD,label="∂ng∂t,FD")
plot!(p_∂ng∂t,ts,∂ng∂t_AD/1e4,label="∂ng∂t,AD")

p_∂ng∂w = plot(ws,∂ng∂w_FD,label="∂ng∂w,FD")
plot!(p_∂ng∂w,ws,∂ng∂w_AD/1e4,label="∂ng∂w,AD")

################################################################################
################################################################################
ng_rwg(collect(ωs,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
grads_ng_rwg = [gradient(ng_rwg,om,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz) for om in ωs]

ng_rwg_soln, ng_rwg_pb = Zygote.pullback(ng_rwg,1/1.55,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
∂ng = ng_rwg_pb(1)
∂ng = gradient(ng_rwg,1/1.55,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)
∂ng[1:6]
Zygote.gradient(ng_rwg,collect(0.645*(0.9:0.02:1.1)),w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)

gradient(n_rwg,1/1.55,w,t_core,edge_gap,n_core,n_subs,Δx,Δy,Δz,Nx,Ny,Nz)

ωs = 1/1.55 * collect(0.9:0.02:1.1)
loo(ωs,w,t_core,edge_gap,n_core,n_subs)
@btime loo($ωs,$w,$t_core,$edge_gap,$n_core,$n_subs)

Zygote.gradient(loo,1/1.55,w,t_core,edge_gap,n_core,n_subs)

# plot_ε(hoo(1/1.55,w,t_core,edge_gap,n_core,n_subs,g),g.x,g.y)
# H,k = hoo(1/1.55,w,t_core,edge_gap,n_core,n_subs,g)
# plot_d⃗(H,k,g)

sum(eps_hoo)

function goo(w,Δx,Δy,Δz,Nx,Ny,Nz)
        g = MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
        sum(εₛ⁻¹(test_shapes2(w),g))
end
foo2(w,Δx,Δy,Δz,Nx,Ny,Nz)
Zygote.gradient(foo2,w,Δx,Δy,Δz,Nx,Ny,Nz)

ε⁻¹ = [ε⁻¹_SA[ix,iy,iz][a,b] for a=1:3,b=1:3,ix=1:Nx,iy=1:Ny,iz=1:Nz]
H,k = solve_k(ω,ε⁻¹,Δx,Δy,Δz)

compare_fields(d,d,g)

g = MaxwellGrid(Δx,Δy,Δz,Nx,Ny,Nz)
tree = KDTree(s)
# ε_sm_inv = copy(reshape( [inv(εₛ(shapes,tree,g.x[i],g.y[j],g.δx,g.δy)) for i=1:g.Nx,j=1:g.Ny], (g.Nx,g.Ny,1)) )
# ε_sm_inv = [inv(εₛ(shapes,tree,g.x[ix],g.y[iy],g.δx,g.δy))[a,b] for a=1:3,b=1:3,ix=1:g.Nx,iy=1:g.Ny,iz=1:g.Nz]
ε_sm_inv = [inv(εₛ(s,tree,g.x[ix],g.y[iy],g.δx,g.δy))[a,b] for a=1:3,b=1:3,ix=1:g.Nx,iy=1:g.Ny,iz=1:g.Nz]

@tullio einv[a,b]

##
using LinearAlgebra
ε = rand(3,3)
[   -1/ε[1,1]                  ε[1,2]/ε[1,1]                                 ε[1,3]/ε[1,1]
    ε[2,1]/ε[1,1]              ε[2, ]- ε[2,1]*ε[1,2]/ε[1,1]                  ε[2,3] - ε[2,1]*ε[1,3]/ε[1,1]
    ε[3,1]/ε[1,1]              ε[3,2] - ε[3,1]*ε[1,2]/ε[1,1]                 ε[3,3] - ε[3,1]*ε[1,3]/ε[1,1]       ]

a⃗ = SVector(1., 2., 3.)
b⃗ = SVector(4.,5.,6.)
normalize(a⃗×b⃗)

gradient(collect(1.:1.:6.)...) do a,b,c,d,e,f
    a⃗ = SVector(a, b, c)
    b⃗ = SVector(d, e, f)
    normalize(a⃗×b⃗)[1]
end

const Float = typeof(0.0)  # use Float = Float128 for quadruple precision in the future
const CFloat = Complex{Float}
const SVec3 = SVector{3}
const SVec3Complex = SVec3{CFloat}
const SVec3Number = SVec3{<:Number}
const SMat3Complex = SMatrix{3,3,CFloat,9}
const SVec2Float = SVec2{Float}
const SVec3Float = SVec3{Float}

ix,iy,iz = 6,9,1
tree = KDTree(s)
εₛ(s,tree,g.x[ix],g.y[iy],g.δx,g.δy)
[(println("$ix,$iy,$iz");inv(εₛ(s,tree,g.x[ix],g.y[iy],g.δx,g.δy))[a,b]) for a=1:3,b=1:3,ix=1:g.Nx,iy=1:g.Ny,iz=1:g.Nz]

# Equation (4) of the paper by Kottke et al.
# @inline function τ_trans(ε₁₁, ε₂₁, ε₃₁, ε₁₂, ε₂₂, ε₃₂, ε₁₃, ε₂₃, ε₃₃)
#     return SMat3Complex(
#         -1/ε₁₁, ε₂₁/ε₁₁, ε₃₁/ε₁₁,
#         ε₁₂/ε₁₁, ε₂₂ - ε₂₁*ε₁₂/ε₁₁, ε₃₂ - ε₃₁*ε₁₂/ε₁₁,
#         ε₁₃/ε₁₁, ε₂₃ - ε₂₁*ε₁₃/ε₁₁, ε₃₃ - ε₃₁*ε₁₃/ε₁₁
#     )
#     # t = zeros(3,3)
#     # t[2:3, 2:3] = s[2:3, 2:3]
#     # s11 = s[1,1]
#     # s[1,1] = -1
#     #
#     # t = t - (s(:,1) * s(1,:)) ./ s11
#     # return t
# end

A = [   11  21  31
        12  22  32
        13  23  33 ]

B = [   11  12  13
        21  22  23
        31  32  33 ]

function ftest(x1,x2,x3,x4,x5,x6,x7,x8,x9)
    println("x1: $x1")
    println("x2: $x2")
    println("x3: $x3")
    println("x4: $x4")
    println("x5: $x5")
    println("x6: $x6")
    println("x7: $x7")
    println("x8: $x8")
    println("x9: $x9")
end

ftest(A...)
# x1: 11
# x2: 12
# x3: 13
# x4: 21
# x5: 22
# x6: 23
# x7: 31
# x8: 32
# x9: 33

ftest(B...)
# x1: 11
# x2: 21
# x3: 31
# x4: 12
# x5: 22
# x6: 32
# x7: 13
# x8: 23
# x9: 33

ftest(vec(A)...)
ftest(vec(B)...)

ins = randn(9)
real(Array(τ_trans(ins...)))   ≈ τ_trans2(ins...)

τ_trans2(ins...)

Array(τ⁻¹_trans(collect(1:9)))   ≈ τ⁻¹_trans2(collect(1:9)...)

τ⁻¹_trans2(collect(1:9)...)

n12 = [1.; 0.; 0.]
n = n12 / sqrt(sum2(abs2,n12))
# Pick a vector that is not along n.
if any(n .== 0.)
    h_temp2 = Float64.(n .== 0.)
    println("n: $n")
    println("h_temp2: $h_temp2")
else
    h_temp2 = [1.; 1.; 1.]
end
h_temp2 = [1.; 1.; 1.]
# h_temp = [0., 0., 1.]
# Create two vectors that are normal to n and normal to each other.
# h_temp2 = n×[0.;0.;1.]
h = h_temp2 / sqrt(sum2(abs2,h_temp2))
v_temp = n×h
v = v_temp / sqrt(sum2(abs2,v_temp))
# Create a local Cartesian coordinate system.
S = [n h v]  # unitary

@inline function τ_trans3(εin)
    ε = εin
    [   -1/ε[1,1]                  ε[1,2]/ε[1,1]                                 ε[1,3]/ε[1,1]
        ε[2,1]/ε[1,1]              ε[2, ]- ε[2,1]*ε[1,2]/ε[1,1]                  ε[2,3] - ε[2,1]*ε[1,3]/ε[1,1]
        ε[3,1]/ε[1,1]              ε[3,2] - ε[3,1]*ε[1,2]/ε[1,1]                 ε[3,3] - ε[3,1]*ε[1,3]/ε[1,1]       ]
end

# Equation (23) of the paper by Kottke et al.
@inline function τ⁻¹_trans3(τin)
    τ = τin
    [   -1/τ[1,1]       -τ[1,2]/τ[1,1]                       -τ[1,3]/τ[1,1]
        -τ[2,1]/τ[1,1]    τ[2,2] - τ[2,1]*τ[1,2]/τ[1,1]       τ[2,3] - τ[2,1]*τ[1,3]/τ[1,1]
        -τ[3,1]/τ[1,1]    τ[3,2] - τ[3,1]*τ[1,2]/τ[1,1]       τ[3,3] - τ[3,1]*τ[1,3]/τ[1,1]      ]
end

param1 = ε_tensor(3.5)
param2 = ε_tensor(1.4)
rvol1 = 0.3
  # unitary
τ1 = τ_trans3(transpose(S) * param1 * S)  # express param1 in S coordinates, and apply τ transform
τ2 = τ_trans3(transpose(S) * param2 * S)  # express param2 in S coordinates, and apply τ transform
τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume
S * τ⁻¹_trans(τavg) * transpose(S)

@inline function τ_trans(ε₁₁, ε₂₁, ε₃₁, ε₁₂, ε₂₂, ε₃₂, ε₁₃, ε₂₃, ε₃₃)
    [   -1/ε₁₁               ε₁₂/ε₁₁                            ε₁₃/ε₁₁
        ε₂₁/ε₁₁              ε₂₂ - ε₂₁*ε₁₂/ε₁₁                  ε₂₃ - ε₂₁*ε₁₃/ε₁₁
        ε₃₁/ε₁₁              ε₃₂ - ε₃₁*ε₁₂/ε₁₁                  ε₃₃ - ε₃₁*ε₁₃/ε₁₁       ]
end

# Equation (23) of the paper by Kottke et al.
@inline function τ⁻¹_trans(τ₁₁, τ₂₁, τ₃₁, τ₁₂, τ₂₂, τ₃₂, τ₁₃, τ₂₃, τ₃₃)
    [   -1/τ₁₁      -τ₁₂/τ₁₁                -τ₁₃/τ₁₁
        -τ₂₁/τ₁₁    τ₂₂ - τ₂₁*τ₁₂/τ₁₁       τ₂₃ - τ₂₁*τ₁₃/τ₁₁
        -τ₃₁/τ₁₁    τ₃₂ - τ₃₁*τ₁₂/τ₁₁       τ₃₃ - τ₃₁*τ₁₃/τ₁₁      ]
end

function avg_param(param1, param2, n12, rvol1)
    n = n12 / sqrt(sum2(abs2,n12))
    # Pick a vector that is not along n.
    # if any(n .== 0)
    # 	h_temp = (n .== 0)
    # else
    # 	h_temp = [1., 0., 0.]
    # end
    # h_temp = [0., 0., 1.]
    # Create two vectors that are normal to n and normal to each other.
    h_temp2 = n×[0.;0.;1.]
    h = h_temp2 / sqrt(sum2(abs2,h_temp2))
    v_temp = n×h
    v = v_temp / sqrt(sum2(abs2,v_temp))
    # Create a local Cartesian coordinate system.
    S = [n h v]  # unitary
    τ1 = τ_trans(transpose(S) * param1 * S)  # express param1 in S coordinates, and apply τ transform
    τ2 = τ_trans(transpose(S) * param2 * S)  # express param2 in S coordinates, and apply τ transform
    τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume-weighted average
    return S * τ⁻¹_trans(τavg) * transpose(S)  # apply τ⁻¹ and transform back to global coordinates
end

ε_tensor(3.5)

gradient(3.5,1.4,1.,2.,3.,0.3) do n1,n2,a,b,c,r
    abs2(sum(avg_param(ε_tensor(n1),ε_tensor(n2),[a;b;c],r)))
end

function fgrad(p1,p2,p3,p4,p5,p6)
    gradient(p1,p2,p3,p4,p5,p6) do n1,n2,a,b,c,r
        abs2(sum(kottke_avg_param2(ε_tensor(n1),ε_tensor(n2),[a;b;c],r)))
    end
end

function f(n1,n2,a,b,c,r)
    # abs2(sum(kottke_avg_param2(ε_tensor(n1),ε_tensor(n2),[a,b,c],r)))
    kottke_avg_param2(ε_tensor(n1),ε_tensor(n2),[a,b,c],r)
end

@inline function τ_trans2(ε₁₁, ε₂₁, ε₃₁, ε₁₂, ε₂₂, ε₃₂, ε₁₃, ε₂₃, ε₃₃)
    [   -1/ε₁₁               ε₁₂/ε₁₁                            ε₁₃/ε₁₁
        ε₂₁/ε₁₁              ε₂₂ - ε₂₁*ε₁₂/ε₁₁                  ε₂₃ - ε₂₁*ε₁₃/ε₁₁
        ε₃₁/ε₁₁              ε₃₂ - ε₃₁*ε₁₂/ε₁₁                  ε₃₃ - ε₃₁*ε₁₃/ε₁₁       ]
end

# Equation (23) of the paper by Kottke et al.
@inline function τ⁻¹_trans2(τ₁₁, τ₂₁, τ₃₁, τ₁₂, τ₂₂, τ₃₂, τ₁₃, τ₂₃, τ₃₃)
    [   -1/τ₁₁      -τ₁₂/τ₁₁                -τ₁₃/τ₁₁
        -τ₂₁/τ₁₁    τ₂₂ - τ₂₁*τ₁₂/τ₁₁       τ₂₃ - τ₂₁*τ₁₃/τ₁₁
        -τ₃₁/τ₁₁    τ₃₂ - τ₃₁*τ₁₂/τ₁₁       τ₃₃ - τ₃₁*τ₁₃/τ₁₁      ]
end

function kottke_avg_param2(param1, param2, n12, rvol1)
    n = n12 / sqrt(sum2(abs2,n12))
    # Pick a vector that is not along n.
    # if any(n .== 0)
    # 	h_temp = (n .== 0)
    # else
    # 	h_temp = [1., 0., 0.]
    # end
    # h_temp = [0., 0., 1.]
    # Create two vectors that are normal to n and normal to each other.
    h_temp2 = n×[0.;0.;1.]
    h = h_temp2 / sqrt(sum2(abs2,h_temp2))
    v_temp = n×h
    v = v_temp / sqrt(sum2(abs2,v_temp))
    # Create a local Cartesian coordinate system.
    S = [n h v]  # unitary
    τ1 = τ_trans((transpose(S) * param1 * S)...)  # express param1 in S coordinates, and apply τ transform
    τ2 = τ_trans((transpose(S) * param2 * S)...)  # express param2 in S coordinates, and apply τ transform
    τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume-weighted average
    return S * τ⁻¹_trans(τavg...) * transpose(S)  # apply τ⁻¹ and transform back to global coordinates
end

f(3.5,1.4,1.,2.,3.,0.3)
fgrad(3.5,1.4,1.,2.,3.,0.3)
@btime f(3.5,1.4,1.,2.,3.,0.3)
@btime fgrad(3.5,1.4,1.,2.,3.,0.3)

##
tree = KDTree(s)
eps_bg = εₛ(s,g)[1,1,1]

x,y = g.x[50], g.y[50]
x1,y1 = x+g.δx/2.,y+g.δy/2
x2,y2 = x+g.δx/2.,y-g.δy/2
x3,y3 = x-g.δx/2.,y-g.δy/2
x4,y4 = x-g.δx/2.,y+g.δy/2

using MaxwellFDM: kottke_avg_param

vxl = (SVector{2,Float64}(x3,y3), SVector{2,Float64}(x1,y1))
rvol = volfrac(vxl,nout,r₀)
r₀,nout = surfpt_nearby([x, y], s[2])
e0 = kottke_avg_param(
        SHM3(s[2].data),
        SHM3(eps_bg),
        SVector{3,Float64}(nout[1],nout[2],0),
        rvol,)

@tullio eps[a,b,ix,iy,iz] :=

Array{Float64}(real.(e0))

ε_sm_inv = [inv(εₛ(s,tree,g.x[ix],g.y[iy],g.δx,g.δy))[a,b] for a=1:3,b=1:3,ix=1:g.Nx,iy=1:g.Ny,iz=1:g.Nz]

using Revise, StaticArrays, GeometryPrimitives, Zygote
##
Zygote.@adjoint (T::Type{<:SArray})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
Zygote.@adjoint (T::Type{<:SArray})(x::AbstractArray) = T(x), dv -> (nothing, dv)
Zygote.@adjoint (T::Type{<:SMatrix})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
Zygote.@adjoint (T::Type{<:SMatrix})(x::AbstractMatrix) = T(x), dv -> (nothing, dv)
Zygote.@adjoint (T::Type{<:SVector})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
Zygote.@adjoint (T::Type{<:SVector})(x::AbstractVector) = T(x), dv -> (nothing, dv)

Zygote.@adjoint enumerate(xs) = enumerate(xs), diys -> (map(last, diys),)
_ndims(::Base.HasShape{d}) where {d} = d
_ndims(x) = Base.IteratorSize(x) isa Base.HasShape ? _ndims(Base.IteratorSize(x)) : 1
Zygote.@adjoint function Iterators.product(xs...)
                    d = 1
                    Iterators.product(xs...), dy -> ntuple(length(xs)) do n
                        nd = _ndims(xs[n])
                        dims = ntuple(i -> i<d ? i : i+nd, ndims(dy)-nd)
                        d += nd
                        func = sum(y->y[n], dy; dims=dims)
                        ax = axes(xs[n])
                        reshape(func, ax)
                    end
                end

function sum2(op,arr)
    return sum(op,arr)
end

function sum2adj( Δ, op, arr )
    n = length(arr)
    g = x->Δ*Zygote.gradient(op,x)[1]
    return ( nothing, map(g,arr))
end

Zygote.@adjoint function sum2(op,arr)
    return sum2(op,arr),Δ->sum2adj(Δ,op,arr)
end
##
Zygote.refresh()

b = Box([0.,0.],[1.,2.]) # ([center],[size along axes])

sum(abs2,sum.(bounds(b)))
bounds(b)
gradient(0.5,0.3,1.0,2.0) do c1,c2,r1,r2
        sum(abs2,sum.(bounds(Box([c1,c2],[r1,r2]))))
end

gradient(0.5,0.6,1.4,2.8) do c1,c2,r1,r2
        sum(abs2,sum.(bounds2(Box([c1,c2],[r1,r2]))))
end

signmatrix(b::Box{1}) = SMatrix{1,1}(1)
signmatrix(b::Box{2}) = SMatrix{2,2}(1,1, -1,1)
signmatrix(b::Box{3}) = SMatrix{3,4}(1,1,1, -1,1,1, 1,-1,1, 1,1,-1)

function bounds2(b::Box)
    A = inv(b.p) .* b.r'
    # m = maximum(abs.(A * signmatrix(b)), dims=2)[:,1] # extrema of all 2^N corners of the box
    m = maximum(abs.(Array((inv(b.p) .* b.r') * signmatrix(b))), dims=2)[:,1] # extrema of all 2^N corners of the box
    return (b.c-m,b.c+m)
end

bounds(b)
bounds2(b)
bounds(b) == bounds2(b)
A = inv(b.p) .* b.r'

abs.(A * signmatrix(b))

x1 = maximum(abs.((inv(b.p) .* b.r') * signmatrix(b)), dims=2)[:,1]

SVector{2,Float64}(maximum(abs.(Array((inv(b.p) .* b.r') * signmatrix(b))), dims=2)[:,1])
x2 = maximum(abs.( Array((inv(b.p) .* b.r') * signmatrix(b)) ), dims=2)

"""
################################################################################
#																			   #
#					check gradients of polygon shape functions				   #
#																			   #
################################################################################
"""
##
ε_core = ε_tensor(n_core)
tanθ = tan(θ)
tcore_tanθ = t_core*tanθ
w_bottom = w + 2*tcore_tanθ
verts = 0.5.*   [   w     -w     -w_bottom    w_bottom
                    t_core   t_core    -t_core      -t_core    ]'
s = GeometryPrimitives.Polygon(					                        # Instantiate 2D polygon, here a trapazoid
                verts,			                                            # v: polygon vertices in counter-clockwise order
                ε_core,					                                    # data: any type, data associated with box shape
            )

##

x = 0.6:0.003:0.9
xy = [[xx,0.1] for xx in x]
spt_n = [surfpt_nearby(xxyy,s) for xxyy in xy]
sptx = [a[1][1] for a in spt_n]
spty = [a[1][2] for a in spt_n]

spt_n2 = [surfpt_nearby2(xxyy,s) for xxyy in xy]
sptx2 = [a[1][1] for a in spt_n2]
spty2 = [a[1][2] for a in spt_n2]

@assert sptx2 ≈ sptx
@assert spty2 ≈ spty

sptx_grad = [Zygote.gradient(v->surfpt_nearby2([v,0.1],s)[1][1],a)[1] for a in x]
spty_grad = [Zygote.gradient(v->surfpt_nearby2([v,0.1],s)[1][2],a)[1] for a in x]

##
plt = plot(x,sptx,label="sptx")
plot!(plt,x,spty,label="spty")
plt_grad = plot(x,sptx_grad,label="sptx_grad")
plot!(plt_grad,x,spty_grad,label="spty_grad")
l = @layout [   a
                b   ]
plot(plt,plt_grad,layout=l)
##
abs.((-)(bounds(s)...))

bounds(s)

((-)(bounds(s)...))

##

plot_ε(εₛ(s5,g;npix_sm=2),g.x,g.y)
esm = εₛ(s5,g)
minimum(esm[2,1,:,:])


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module MockGMM2

include("../../../src/julia/shared/load.jl")
using ADPerfTest
using GMMData

export get_gmm_test

mutable struct MockGMMContext
    waitcoef::Float64
end

get_gmm_test() = Test{GMMInput, GMMOutput}(
    MockGMMContext(2.0),
    (ctx, input) -> nothing,
    (ctx, times) -> nothing,
    (ctx, times) -> nothing,
    (out, ctx) -> begin out.objective = 2.0 end
)

end

ext: jl

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
using BinaryBuilder, Pkg

name = "SDPA"
version = v"7.3.8"

# Collection of sources required to build SDPABuilder
sources = [
    ArchiveSource("https://sourceforge.net/projects/sdpa/files/sdpa/sdpa_$(version).tar.gz",
                  "c7541333da2f0bb2d18e90dbf758ac7cc099f3f7da3f256b284b0725f96d4117")
    DirectorySource("./bundled")
]

MUMPS_seq_version = v"4.10.0"
MUMPS_seq_packagespec = PackageSpec(; name = "MUMPS_seq_jll",
                                    uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d",
                                    version = MUMPS_seq_version)

METIS_version = v"4.0.3"
METIS_packagespec = PackageSpec(; name = "METIS_jll",
                                uuid = "d00139f3-1899-568f-a2f0-47f597d42d70",
                                version = METIS_version)

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

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la
rm -f /opt/${target}/${target}/lib*/*.la
update_configure_scripts

# Apply patches
atomic_patch -p1 $WORKSPACE/srcdir/patches/shared.diff
mv configure.in configure.ac
atomic_patch -p1 $WORKSPACE/srcdir/patches/lt_init.diff
autoreconf -vi

export CPPFLAGS="${CPPFLAGS} -I${prefix}/include -I$prefix/include/coin"
export CXXFLAGS="${CXXFLAGS} -std=c++11"
if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L$prefix/bin"
elif [[ ${target} == *linux* ]]; then
    export LDFLAGS="-ldl -lrt"
fi

./configure --prefix=$prefix --with-pic --disable-pkg-config --build=${MACHTYPE} --host=${target} \
--enable-shared lt_cv_deplibs_check_method=pass_all \
--with-blas="-lopenblas" --with-lapack="-lopenblas" \
--with-mumps-libs="-L${prefix}/lib -ldmumps -lzmumps -lcmumps -lsmumps -lmumps_common -lmpiseq -lpord -lmetis -lopenblas -lgfortran -lpthread" \
--with-mumps-include="-I${prefix}/include/mumps_seq"

make -j${nproc}
make install

## Then build the libcxxwrap-julia wrapper
cd $WORKSPACE/srcdir/sdpawrap

mkdir build
cd build/

if [[ $target == i686-* ]] || [[ $target == arm-* ]]; then
    export processor=pentium4
else
    export processor=x86-64
fi

# Override compiler ID to silence the horrible "No features found" cmake error
if [[ $target == *"apple-darwin"* ]]; then
  macos_extra_flags="-DCMAKE_CXX_COMPILER_ID=AppleClang -DCMAKE_CXX_COMPILER_VERSION=10.0.0 -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
fi

cmake $macos_extra_flags -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DSDPA_DIR=$prefix \
      -DMUMPS_INCLUDE_DIR="../../destdir/include/mumps_seq" \
      -DCMAKE_FIND_ROOT_PATH=${prefix} \
      -DJulia_PREFIX=${prefix} \
      -DSDPA_LIBRARY="-lsdpa" \
      -D_GLIBCXX_USE_CXX11_ABI=1 \
      -DJlCxx_DIR=${prefix}/lib/cmake/JlCxx \
      ..
cmake --build . --config Release --target install

if [[ $target == *w64-mingw32* ]] ; then
    cp $WORKSPACE/destdir/lib/libsdpawrap.dll ${libdir}
fi

install_license $WORKSPACE/srcdir/sdpa-*/COPYING
"""

# The products that we will ensure are always built
products = [
    ExecutableProduct("sdpa", :sdpa),
    LibraryProduct("libsdpa", :libsdpa),
    LibraryProduct("libsdpawrap", :libsdpawrap)
]

# Pick platforms from L/libcxxwrap-julia/build_tarballs.jl
platforms = [
    FreeBSD(:x86_64; compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    Linux(:armv7l; libc=:glibc, compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    Linux(:aarch64; libc=:glibc, compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    Linux(:x86_64; libc=:glibc, compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    Linux(:i686; libc=:glibc, compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    MacOS(:x86_64; compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    Windows(:x86_64; compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
    Windows(:i686; compiler_abi=CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11)),
]
#platforms = expand_gfortran_versions(platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    Dependency("OpenBLAS32_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    BuildDependency(PackageSpec(name="Julia_jll", version="v1.4.1")),
    BuildDependency(MUMPS_seq_packagespec),
    BuildDependency(METIS_packagespec),
]

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

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
using DSGE

path = dirname(@__FILE__)
# Set parameters for testing
custom_settings = Dict{Symbol, Setting}(
    :date_forecast_start  => Setting(:date_forecast_start, quartertodate("2015-Q4")))
m = AnSchorfheide(custom_settings = custom_settings, testing = true)
m <= Setting(:use_parallel_workers, true)
parallel = get_setting(m, :use_parallel_workers)

# Set number of draws
draws = 10000

# Remove previous testing file
rm("$path/../reference/systematic_resampling_test.csv")

# Append weights vectors for testing
weights_vec=Any[]
# Increasing weight
append!(weights_vec, [collect(0.1:0.1:0.9)])
# Even distribution between indices
append!(weights_vec,[[0.4, 0.1, 0.1, 0.2, 0.3]])
# Index in the middle should appear far more frequently than others
append!(weights_vec,[[0.1, 0.8, 0.1]])
# Weights whose sum far exceeds 0
append!(weights_vec,[[9.0, 20.0, 2.0, 7.0, 14.0, 23.0, 4.0]])
# Index with 0 probability
append!(weights_vec,[[0.0, 0.5, 0.5]])

open("$path/../reference/systematic_resampling_test.csv","w") do file
   write(file,"---------------------")
end

for weights in weights_vec

    n_parts = length(weights)

    count = zeros(n_parts)
    out = zeros(n_parts)

    for i=1:draws
        out = resample(weights; method = :systematic, parallel = parallel)
        for idx in out
            count[idx] += 1
        end
    end

    count = count./(draws*n_parts)
    count = round(count,3)
    act_weights = round(weights./sum(weights),3)

    open("$path/../reference/systematic_resampling_test.csv","a") do file
        write(file,"\nActual Probabilities: ")
        writecsv(file,act_weights')
        write(file,"Tested Probabilities: ")
        writecsv(file,count')
        write(file,"---------------------")
    end

end

nothing


text:
using OrdinaryDiffEq, DiffEqSensitivity, Zygote, Test

function f!(du,u,p,t)
    du[1] = -p[1]'*u
    du[2] = (p[2].a + p[2].b)u[2]
    du[3] = p[3](u,t)
    return nothing
end

struct mystruct
    a
    b
end

function control(u,t)
    return -exp(-t)*u[3]
end

u0 = [10,15,20]
p = [[1;2;3], mystruct(-1,-2), control]
tspan = (0.0,10.0)

prob = ODEProblem(f!,u0, tspan, p)

sol = solve(prob, Tsit5()) # Solves without errors

function loss(p1)
    sol = solve(prob, Tsit5(), p=[p1, mystruct(-1,-2), control])
    return sum(abs2, sol)
end

grad(p) = Zygote.gradient(loss, p)

p2 = [4;5;6]
@test_throws DiffEqSensitivity.ForwardDiffSensitivityParameterCompatibilityError grad(p2)

function loss(p1)
    sol = solve(prob, Tsit5(), p=[p1, mystruct(-1,-2), control], sensealg = InterpolatingAdjoint())
    return sum(abs2, sol)
end

@test_throws DiffEqSensitivity.AdjointSensitivityParameterCompatibilityError grad(p2)

function loss(p1)
    sol = solve(prob, Tsit5(), p=[p1, mystruct(-1,-2), control], sensealg = ForwardSensitivity())
    return sum(abs2, sol)
end

@test_throws DiffEqSensitivity.ForwardSensitivityParameterCompatibilityError grad(p2)
@test_throws DiffEqSensitivity.ForwardSensitivityParameterCompatibilityError ODEForwardSensitivityProblem(f!,u0, tspan, p)


text:
const pi = 3.14;
println(pi);

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
"Infer fluents that are never modified by some action in a domain."
function infer_static_fluents(domain::Domain)
    affected = Set(infer_affected_fluents(domain))
    static = setdiff(keys(get_fluents(domain)), affected)
    return collect(static)
end

"Infer fluents that are modified by some action in a domain."
function infer_affected_fluents(domain::Domain)
    affected = Symbol[]
    for action in values(get_actions(domain))
        append!(affected, get_affected(action))
    end
    return unique!(affected)
end

const AFFECTED_FUNCS = Dict{Symbol, Function}()

"Return the names of all fluents affected by an action."
function get_affected(action::Action)
    return get_affected(get_effect(action))
end

function get_affected(effect::Term)
    affected_fn! = get(AFFECTED_FUNCS, effect.name, builtin_affected!)
    return affected_fn!(Symbol[], effect)
end

function builtin_affected!(affected::Vector{Symbol}, effect::Term)
    if effect.name == :and
        for eff in effect.args append!(affected, get_affected(eff)) end
    elseif effect.name == :when
        append!(affected, get_affected(effect.args[2]))
    elseif effect.name == :forall
        append!(affected, get_affected(effect.args[2]))
    elseif effect.name == :assign
        push!(affected, effect.args[1].name)
    elseif effect.name in keys(GLOBAL_MODIFIERS)
        push!(affected, effect.args[1].name)
    elseif effect.name == :not
        push!(affected, effect.args[1].name)
    else
        push!(affected, effect.name)
    end
    return affected
end

"Check if term is static or composed of static subterms."
function is_static(term::Term, domain::Domain,
                   statics=infer_static_fluents(domain))
    if is_external_func(term, domain)
        return all(is_static(a, domain, statics) for a in term.args)
    else
        return term.name in statics
    end
end

function is_static(term::Const, domain::Domain,
                   statics=infer_static_fluents(domain))
    !is_fluent(term, domain) || term.name in statics
end

"Simplify away static fluents within a `term`."
function simplify_statics(term::Term, domain::Domain, state::State,
                          statics=infer_static_fluents(domain))
    # Simplify predicates if they are static and ground
    if is_static(term, domain, statics) && is_ground(term)
        return Const(evaluate(domain, state, term))
    elseif !(term.name in (:and, :or, :imply, :not))
        return term
    end
    # Simplify logical compounds
    args = Term[simplify_statics(a, domain, state, statics) for a in term.args]
    if term.name == :and
        true_idxs = Int[]
        for (i, a) in enumerate(args)
            a.name == false && return Const(false)
            a.name == true && push!(true_idxs, i)
        end
        length(true_idxs) == length(args) && return Const(true)
        deleteat!(args, true_idxs)
        return length(args) == 1 ? args[1] : Compound(:and, args)
    elseif term.name == :or
        false_idxs = Int[]
        for (i, a) in enumerate(args)
            a.name == true && return Const(true)
            a.name == false && push!(false_idxs, i)
        end
        length(false_idxs) == length(args) && return Const(false)
        deleteat!(args, false_idxs)
        return length(args) == 1 ? args[1] : Compound(:or, args)
    elseif term.name == :imply
        cond, query = args
        cond.name == true && return query
        cond.name == false && return Const(true)
        return Compound(:imply, args)
    elseif term.name == :not
        val = args[1].name
        return val isa Bool ? Const(!val) : Compound(:not, args)
    else
        error("Unrecognized logical operator: $(term.name)")
    end
end


text:
using RecipesBase
using UnitfulRecipes

#include("Arrays.jl")
#include("RadiationPattern.jl")
#include("Arrays.jl")

function grid_pattern(ϕ, θ, r, min, max)
    shape = (length(ϕ), length(θ))
    # Preallocate matricies
    x = zeros(Float64, shape)
    y = zeros(Float64, shape)
    z = zeros(Float64, shape)
    # Reshape radius data
    r = [(data + abs(min)) / (max + abs(min)) for data in r]
    # Eliminate negative numbers
    replace!(x -> x < 0 ? 0 : x, r)
    for (i, phi) in enumerate(ϕ), (j, theta) in enumerate(θ)
        x[i,j] = r[i,j] * sind(theta) * cosd(phi)
        y[i,j] = r[i,j] * sind(theta) * sind(phi)
        z[i,j] = r[i,j] * cosd(theta)
    end
    return x, y, z
end

@userplot PolarPattern
@recipe function f(pp::PolarPattern)
    if !(((length(pp.args) == 1) && (pp.args[1] isa RadiationPattern)) ||
         ((length(pp.args) == 2) && (pp.args[1] isa AbstractVector &&
                                     pp.args[2] isa AbstractVector)))
       error("Polar pattern should be given two vectors of θ and gain, or a `RadiationPattern` .
              Got: $(typeof(pp.args))")
   end
    if pp.args[1] isa RadiationPattern
        if length(pp.args[1].dims) == 2
            θ = deg2rad.(pp.args[1].dims[2].val)
            r = gain(pp.args[1]) |> ustrip |> transpose
        else
            θ = deg2rad.(pp.args[1].dims[1].val)
            r = gain(pp.args[1]) |> ustrip
        end
    else
        θ, r = pp.args
    end
    # R is log-scalled, θ was in degrees
    if haskey(plotattributes, :lims)
        if plotattributes[:lims][1] == :auto
            min = minimum(r)
            max = plotattributes[:lims][2]
        elseif plotattributes[:lims][2] == :auto
            min = plotattributes[:lims][1]
            max = maximum(r)
        else
            min, max = plotattributes[:lims]
        end
    else
        min = minimum(r)
        max = maximum(r)
    end
    # Clamp all values less then the minimum to the minimum
    r = [x < min ? min : x for x in r]
    @series begin
        projection := :polar
        lims --> (min, max)
        θ, ustrip(r)
    end
end

@userplot SphericalPattern
@recipe function f(sp::SphericalPattern)
    if !(((length(sp.args) == 1) && (sp.args[1] isa RadiationPattern)) ||
         ((length(sp.args) == 3) && (sp.args[1] isa AbstractVector &&
                                     sp.args[2] isa AbstractVector &&
                                     sp.args[3] isa AbstractMatrix)))
        error("Spherical pattern should be given two vectors and a matrix, or a `RadiationPattern` .
               Got: $(typeof(sp.args))")
    end
    if sp.args[1] isa RadiationPattern
        ϕ = sp.args[1].dims[1].val
        θ = sp.args[1].dims[2].val
        r = gain(sp.args[1])
    else
        ϕ, θ, r = sp.args
    end
    r = ustrip(r)
    if haskey(plotattributes, :lims)
        if plotattributes[:lims][1] == :auto
            min = minimum(r)
            max = plotattributes[:lims][2]
        elseif plotattributes[:lims][2] == :auto
            min = plotattributes[:lims][1]
            max = maximum(r)
        else
            min, max = plotattributes[:lims]
        end
    else
        min = minimum(r)
        max = maximum(r)
    end
    x, y, z = grid_pattern(ϕ, θ, r, min, max)
    c = [x < min ? min : x for x in r]
    @series begin
        grid := false
        seriestype := :surface
        xaxis := false
        yaxis := false
        zaxis := false
        lims := :auto
        fill_z := c
        xlims := extrema(x)
        ylims := extrema(y)
        zlims := extrema(z)
        x, y, z
    end
end


text:
module DataGridReader

export Excel, Csv

include("Excel.jl")
using .Excel

include("Csv.jl")
using .Csv

end


text:
module NewPkgEval

    using BinaryBuilder
    using BinaryProvider
    using LightGraphs
    import Pkg.TOML
    using Pkg
    import Base: UUID
    using Dates
    using DataStructures
    import LibGit2

    downloads_dir(name) = joinpath(@__DIR__, "..", "deps", "downloads", name)
    julia_path(ver) = joinpath(@__DIR__, "..", "deps", "julia-$ver")
    versions_file() = joinpath(@__DIR__, "..", "deps", "Versions.toml")

    function read_versions()
        vers = TOML.parse(read(versions_file(), String))
    end

    function obtain_julia(the_ver)
        vers = read_versions()
        for (ver, data) in vers
            ver = VersionNumber(ver)
            ver == the_ver || continue
            if haskey(data, "url")
                file = get(data, "file", "julia-$ver.tar.gz")
                @assert !isabspath(file)
                download_verify_unpack(
                    data["url"],
                    data["sha"],
                    julia_path(ver);
                    tarball_path=downloads_dir(file),
                    force=true
                )
            else
                file = data["file"]
                !isabspath(file) && (file = downloads_dir(file))
                BinaryProvider.verify(file, data["sha"])
                isdir(julia_path(ver)) || BinaryProvider.unpack(file, julia_path(ver))
            end
            return
        end
        error("Requested Julia version not found")
    end

    function run_sandboxed_julia(args=``; ver=v"1.0", do_obtain=true, kwargs...)
        if do_obtain
            obtain_julia(ver)
        else
            @assert ispath(julia_path(ver))
        end
        #ispath(registry_path()) || error("Please run `NewPkgEval.get_registry()` first")
        jp = joinpath(julia_path(ver), first(readdir(julia_path(ver))))
        runner = BinaryBuilder.UserNSRunner(pwd(),
            workspaces=[
                jp => "/maps/julia",
                #registry_path() => "/maps/registries/Uncurated"
            ])
        BinaryBuilder.run_interactive(runner, `/maps/julia/bin/julia --color=yes $args`; kwargs...)
    end

    log_path(ver) = joinpath(@__DIR__, "..", "logs-$ver")
    function run_sandboxed_test(pkg; ver=v"1.0", do_depwarns=false, kwargs...)
        isdir(log_path(ver)) || mkdir(log_path(ver))
        log = joinpath(log_path(ver), "$pkg.log")
        c = quote
            mkpath("/root/.julia/registries")
            open("/etc/hosts", "w") do f
                println(f, "127.0.0.1\tlocalhost")
            end
            run(`mount -t devpts -o newinstance jrunpts /dev/pts`)
            run(`mount -o bind /dev/pts/ptmx /dev/ptmx`)
            run(`mount -t tmpfs tempfs /dev/shm`)
            #run(`ln -s /opt/General/ /root/.julia/registries/General`)
            Pkg.add($pkg)
            Pkg.test($pkg)
        end
        arg = "using Pkg; eval($(repr(c)))"
        try
            open(log, "w") do f
                cmd = ``
                if do_depwarns
                    cmd = `--depwarn=error`
                end
                cmd = `$cmd -e $arg`
                run_sandboxed_julia(cmd; ver=ver, kwargs..., stdout=f, stderr=f)
            end
            return true
        catch e
            return false
        end
    end

    const CSI = "\e["

    function show_result(io::IO, dg, queue, results)
        print(io, "Success: ")
        o, f, s = count(==(:ok), values(results)), count(==(:fail), values(results)),
            count(==(:skipped), values(results))
        q = length(queue)
        printstyled(io, o; color = :green)
        print(io, "\tFailed: ")
        printstyled(io, f; color = :red)
        print(io, "\tSkipped: ")
        printstyled(io, s; color = :yellow)
        print(io, "\tCurrent Frontier/Remaining: ")
        print(io, q, '/')
        print(io, length(vertices(dg.g)) - (o+f+s))
        println(io)
    end

    struct PkgEntry
        degree::Int
        idx::Int
    end
    Base.:<(a::PkgEntry, b::PkgEntry) = a.degree < b.degree

    function skip_pkg!(result, dg, to_skip)
        result[dg.names[to_skip]] = :skipped
        for revdep in inneighbors(dg.g, to_skip)
            if haskey(result, dg.names[revdep]) &&
                result[dg.names[revdep]] == :skipped
                continue
            end
            skip_pkg!(result, dg, revdep)
        end
    end

    const skip_list = [
        "AbstractAlgebra", # Hangs forever
        "DiscretePredictors", # Hangs forever
        "LinearLeastSquares", # Hangs forever
        "SLEEF", # Hangs forever
        "OrthogonalPolynomials", # Hangs forever
        "IndexableBitVectors",
        "LatinHypercubeSampling", # Hangs forever
        "DynamicalBilliards", # Hangs forever
        "ChangePrecision", # Hangs forever
        "Rectangle", # Hangs forever
        "Parts", # Hangs forever
        "ZippedArrays", # Hangs forever
        "Chunks", # Hangs forever
        "Electron",
        "DotOverloading",
        "ValuedTuples",
        "HCubature",
        "SequentialMonteCarlo",
        "RequirementVersions",
        "NumberedLines",
        "LazyContext",
        "RecurUnroll", # deleted, hangs
        "TypedBools", # deleted, hangs
        "LazyCall", # deleted, hangs
        "MeshCatMechanisms",
        "SessionHacker",
        "Embeddings",
        "GeoStatsDevTools",
    ]

    const ok_list = [
        "BinDeps", # Not really ok, but packages may list it just as a fallback
        "InteractiveUtils", # We rely on LD_LIBRARY_PATH working for the moment
        "Homebrew",
        "WinRPM",
        "NamedTuples", # As requested by quinnj
        "Compat",
        "LinearAlgebra"
    ]

    function run_all(dg, ninstances, ver, result = Dict{String, Symbol}(); do_depwarns=false)
        obtain_julia(ver)
        frontier = BitSet()
        pkgs = copy(dg.names)
        running = Vector{Union{Nothing, Symbol}}(nothing, ninstances)
        times = DateTime[now() for i = 1:ninstances]
        processed = BitSet()
        cond = Condition()
        queue = binary_maxheap(PkgEntry)
        completed = Channel(Inf)
        # Add packages without dependencies to the frontier
        for x in ok_list
            result[x] = :ok
            put!(completed, findfirst(==(x), dg.names))
        end
        for p in filter(v->length(outneighbors(dg.g, v)) == 0, vertices(dg.g))
            (dg.names[p] in skip_list) && continue
            (haskey(result, dg.names[p])) && continue
            push!(queue, PkgEntry(length(inneighbors(dg.g, p)), p))
        end
        for x in skip_list
            pkg = findfirst(==(x), dg.names)
            if pkg === nothing
                @warn "Package $(x) in skip list, but not found"
                continue
            end
            skip_pkg!(result, dg, pkg)
        end
        done = false
        processing = false
        all_workers = Task[]
        did_signal_workers = false
        function stop_work()
            if !done
                (done = true; notify(cond); put!(completed, -1))
                if !did_signal_workers
                    for task in all_workers
                        task == current_task() && continue
                        !Base.istaskdone(task) || continue
                        try; schedule(task, InterruptException(); error=true); catch; end
                    end
                    did_signal_workers = true
                end
            end
        end
        @sync begin
            @async begin
                try
                    buf = IOBuffer()
                    io = IOContext(buf, :color => true)
                    while !isempty(queue) || !all(==(nothing), running) || isready(completed) || processing
                        show_result(io, dg, queue, result)
                        for i = 1:ninstances
                            r = running[i]
                            if r === nothing
                                println(io, "Worker $i: -------")
                            else
                                println(io, "Worker $i: $(r) running for ", Dates.canonicalize(Dates.CompoundPeriod(now() - times[i])))
                            end
                        end
                        print(String(take!(buf)))
                        sleep(1)
                        print(io, "$(CSI)$(ninstances+1)A$(CSI)1G$(CSI)0J")
                    end
                    stop_work()
                    println("done")
                catch e
                    @Base.show e
                    stop_work()
                    !isa(e, InterruptException) && rethrow(e)
                end
            end
            # Scheduler
            @async begin
                try
                    while !done
                        pkgno = take!(completed)
                        pkgno == -1 && break
                        processing = true
                        push!(processed, pkgno)
                        for revdep in inneighbors(dg.g, pkgno)
                            if result[dg.names[pkgno]] != :ok
                                skip_pkg!(result, dg, revdep)
                            else
                                (revdep in processed) && continue
                                # Last dependency to finish adds it to
                                # the frontier
                                all_processed = true
                                for dep in outneighbors(dg.g, revdep)
                                    if !(dep in processed) || result[dg.names[dep]] != :ok
                                        all_processed = false
                                        break
                                    end
                                end
                                all_processed || continue
                                haskey(result, dg.names[revdep]) && continue
                                push!(queue, PkgEntry(length(inneighbors(dg.g, revdep)), revdep))
                            end
                        end
                        notify(cond)
                        processing = false
                    end
                catch e
                    @Base.show e
                    stop_work()
                    !isa(e, InterruptException) && rethrow(e)
                end
            end
            # Workers
            for i = 1:ninstances
                push!(all_workers, @async begin
                    try
                        while !done
                            if isempty(queue)
                                wait(cond)
                                continue
                            end
                            pkgno = pop!(queue).idx
                            pkg = dg.names[pkgno]
                            running[i] = Symbol(pkg)
                            times[i] = now()
                            result[pkg] = run_sandboxed_test(pkg, do_depwarns=do_depwarns, ver=ver, do_obtain=false) ? :ok : :fail
                            running[i] = nothing
                            put!(completed, pkgno)
                        end
                    catch e
                        @Base.show e
                        stop_work()
                        !isa(e, InterruptException) && rethrow(e)
                    end
                end)
            end
        end
        result
    end

    function run_pkg_deps(pkg, dg, ninstances)
        result = Dict{String, Symbol}()
        nbs = unique(LightGraphs.neighborhood(dg.g, findfirst(==(pkg), dg.names), typemax(Int)))
        for pkg in dg.names[filter(x->!(x in nbs), 1:length(dg.names))]
            result[pkg] = :skipped
        end
        run_all(dg, ninstances, result)
    end

    function read_all_pkgs()
        pkgs = Tuple{String, UUID, String}[]
        for registry in (registry_path(),)
            open(joinpath(registry, "Registry.toml")) do io
                for (_uuid, pkgdata) in Pkg.Types.read_registry(joinpath(registry, "Registry.toml"))["packages"]
                    uuid = UUID(_uuid)
                    name = pkgdata["name"]
                    path = abspath(registry, pkgdata["path"])
                    push!(pkgs, (name, uuid, path))
                end
            end
        end
        pkgs
    end

    registry_path() = joinpath(@__DIR__, "..", "work", "registry")

    function get_registry()
        if !isdir(registry_path())
            creds = LibGit2.CachedCredentials()
            url = Pkg.Types.DEFAULT_REGISTRIES["General"]
            repo = Pkg.GitTools.clone(url, registry_path(); header = "registry Uncurated from $(repr(url))", credentials = creds)
            close(repo)
        else
            errors = Pkg.API.do_update_registry!(registry_path())
            Pkg.API.print_errors(errors)
        end
    end

    struct PkgDepGraph
        vertex_map::Dict{UUID, Int}
        uuid::Vector{UUID}
        names::Vector{String}
        g::LightGraphs.SimpleGraphs.SimpleDiGraph
    end

    function read_stdlib(ver)
        obtain_julia(ver)
        stdlib_path = joinpath(julia_path(ver), first(readdir(julia_path(ver))),
                "share/julia/stdlib/v$(ver.major).$(ver.minor)")
        ret = Tuple{String, UUID, Vector{UUID}}[]
        stdlibs = readdir(stdlib_path)
        for stdlib in stdlibs
            proj = Pkg.Types.read_project(joinpath(stdlib_path, stdlib, "Project.toml"))
            deps = UUID.(collect(values(proj["deps"])))
            push!(ret, (proj["name"], UUID(proj["uuid"]), deps))
        end
        ret
    end

    function PkgDepGraph(pkgs, ver)
        # Add packages
        vertex_map = Dict(pkg[2] => i for (i, pkg) in enumerate(pkgs))
        uuids = map(x->x[2], pkgs)
        names = map(x->x[1], pkgs)
        # Add stdlibs to vertex map
        stdlibs = read_stdlib(ver)
        stdlib_uuids = Set(x[2] for x in stdlibs)
        for (i, (name, uuid, _)) in enumerate(stdlibs)
            x = findfirst(==(uuid), uuids)
            if x === nothing
                push!(names, name)
                push!(uuids, uuid)
                vertex_map[uuid] = length(names)
            end
        end
        g = LightGraphs.SimpleGraphs.SimpleDiGraph(length(names))
        for (_, uuid, deps) in stdlibs
            for dep in deps
                add_edge!(g, vertex_map[uuid], vertex_map[dep])
            end
        end
        for (name, uuid, path) in pkgs
            (uuid in stdlib_uuids) && continue
            vers = Pkg.Operations.load_versions(path)
            max_ver = maximum(keys(vers))
            data = Pkg.Operations.load_package_data(UUID, joinpath(path, "Deps.toml"), max_ver)
            data === nothing && continue
            for (k,v) in data
                if !haskey(vertex_map, v)
                    @warn "Dependency $k ($v) not in registry"
                    continue
                end
                add_edge!(g, vertex_map[uuid], vertex_map[v])
            end
        end
        # Arbitrarily break dependency cycles
        cycs = LightGraphs.simplecycles(g)
        for cyc in cycs
            rem_edge!(g, cyc[end], cyc[1])
        end
        PkgDepGraph(vertex_map, uuids, names, g)
    end

    function no_dep_pkgs(g::PkgDepGraph)
        filter(v->length(outneighbors(g.g, v)) == 0, vertices(g.g))
    end

    function analyze_results(dg, result)
        # Print the top packages by recursive benefit if fixed
        failed_pkgs = filter(((k,v),)->v == :fail, result)
        r = Dict(begin
            id = findfirst(==(k), dg.names)
            # Avoid cycles
            visited = BitSet()
            stack = Int[id]
            push!(visited, id)
            n = 0
            while !isempty(stack)
                x = pop!(stack)
                for y in inneighbors(dg.g, x)
                    (y in visited) && continue
                    push!(visited, y)
                    push!(stack, y)
                    n += 1
                end
            end
            k => n
        end for (k, v) in failed_pkgs)
        sort(collect(r), by = x->x[2], rev = true)
    end

    function analyze_results_pkg(pkg, dg, result)
        nbs = unique(LightGraphs.neighborhood(dg.g, findfirst(==(pkg), dg.names), typemax(Int)))
        nbs_names = dg.names[nbs]
        nresult = filter(((k,v),)->k in nbs_names, result)
        analyze_results(dg, nresult)
    end

end # module


text:
using DashBootstrapComponents

row = html_div([
    dbc_row(
        dbc_col(
            html_div("A single, half-width column"),
            width = Dict("size" => 6, "offset" => 3),
        ),
    ),
    dbc_row([
        dbc_col(
            html_div("The last of three columns"),
            width = Dict("size" => 3, "order" => "last", "offset" => 1),
        ),
        dbc_col(
            html_div("The first of three columns"),
            width = Dict("size" => 3, "order" => 1, "offset" => 2),
        ),
        dbc_col(
            html_div("The second of three columns"),
            width = Dict("size" => 3, "order" => 5),
        ),
    ]),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
function get_block2d(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    radius=0.0,
    side=0.5,
    contact=true,
    contact_type=:nonlinear,
    mode=:block2d,
    T=Float64)

    # Parameters
    axis = [1.0, 0.0, 0.0]

    origin = Origin{T}()
    pbody = Box(side, side, side, 1.;
        color=RGBA(1., 1., 0.))
    joint1 = JointConstraint(PlanarAxis(origin, pbody, axis))
    bodies = [pbody]
    joints = [joint1]

    if contact
        # Corner vectors
        if mode == :particle
            corners = [[0.0, 0.0, 0.0]]
        elseif mode == :block2d
            corners = [
                [[0.0,  side / 2.0,  side / 2.0]]
                [[0.0,  side / 2.0, -side / 2.0]]
                [[0.0, -side / 2.0,  side / 2.0]]
                [[0.0, -side / 2.0, -side / 2.0]]
            ]
        else
            @error "incorrect mode specified, try :particle or :block2d"
        end
        n = length(corners)
        normal = [[0.0, 0.0, 1.0] for i = 1:n]
        contact_radius = [radius for i = 1:n]
        friction_coefficient = friction_coefficient * ones(n)

        contacts = contact_constraint(pbody, normal;
            friction_coefficient,
            contact_origins=corners,
            contact_radius,
            contact_type)

        mech = Mechanism(origin, bodies, joints, contacts;
            gravity,
            timestep)
    else
        mech = Mechanism(origin, bodies, joints;
            gravity,
            timestep)
    end
    return mech
end

function initialize_block2d!(mechanism::Mechanism{T};
    position=[0.0, 1.0],
    linear_velocity=[0.0, 0.0],
    orientation=0.0,
    angular_velocity=0.0) where T

    if length(mechanism.contacts) > 0
        model = mechanism.contacts[1].model
        side = model.collision.contact_origin[2]
        offset = model.collision.contact_radius
        z = side + offset
    else
        z = 0.0
    end

    body = mechanism.bodies[1]

    set_maximal_configurations!(body,
        x=[0.0; position] + [0.0, 0.0 , z],
        q=RotX(orientation))
    set_maximal_velocities!(body,
        v=[0.0; linear_velocity],
        ω=[angular_velocity, 0.0, 0.0])
end


text:
"""
    bisect(a, b, f, e, n)

Compute the root of the given function using Bisection method.
Uses midpoint of the interval `[a, b]`, replacing one of the endpoinds with
midpoint depending on the location of the root in relation to the midpoint and repeting the process
untill the approximation `p` of the root of the function `f` is found within the specified error margin.

    bisect(a::Integer, b:Integer, f::Function, e::Float64, n::Integer)

### OUTPUT:
  * `p` - approximation of the root of the function `f` (Float64) or,
  * `-1` - if `n` the maximum number of iterations was exceeded.
"""
function bisect(a::Integer, b:Integer, f::Function, e::Float64, n::Integer)

    Ap = []

    for i in 1:n

        p = a + 0.5*(b-a)
        ei = (b-a)/2
        push!(Ap, p)

        if sign(f(a))*sign(f(p)) < 0
            b = p
        else
            a = p
        end

        if ei < e
            return Ap
        end
    end

    return -1
end


text:
"""
Moore-Penrose inverse control allocator.
"""
struct PseudoInverseAllocator <: StaticAllocator
    B_pinv
    function PseudoInverseAllocator(B)
        B_pinv = pinv(B)
        new(B_pinv)
    end
end

"""
# Variables
ν: virtual input
# Notes
ν = B*u where u: control input
"""
function Command(allocator::PseudoInverseAllocator)
    @unpack B_pinv = allocator
    return function (ν, Λ=Diagonal(ones(size(B_pinv)[1])))
        (pinv(Λ) * B_pinv) * ν  # pinv(B*Λ) = pinv(Λ) * pinv(B)
    end
end


text:
using Statistics
using Dates
using FileIO: load
using ImageCore: channelview
using ColorTypes

LearningTarget{N} = Union{CodeBook, NTuple{N, CodeBook}, Multiscale}

function train!(target::LearningTarget, trainingSet::AbstractArray;
    epochs::Integer=1, shape=nothing,
    sparsecoder=SparseCoders.IHT,
    sparsecoder_options=(),
    optimizer=Optimizers.Steepest,
    optimizer_options=(),
    verbose::Union{Integer,Symbol}=1, logdir=Union{Nothing,AbstractString}=nothing,
    plot_function=t->nothing)

    vlevel = verboselevel(verbose)

    savesettings(logdir, target, trainingSet;
        vlevel=vlevel,
        epochs=epochs)

    log_configs = initializelogs(logdir, trainingSet)

    vlevel >= 1 && println("beginning dictionary training...")

    params_dic = getParamsDictionary(target)
    for itr = 1:epochs
        K = length(trainingSet)
        loss_sps = fill(Inf, K)
        loss_dus = fill(Inf, K)
        for k = 1:K
            x = gettrainingdata(trainingSet[k])
            shapek = getvalidshape(shape, target, sparsecoder, optimizer, x)

            aresult = trainPerData!(target, x, params_dic, itr;
                                    shape=shapek, sparsecoder=sparsecoder,
                                    sparsecoder_options=sparsecoder_options,
                                    optimizer=optimizer,
                                    optimizer_options=optimizer_options,
                                    vlevel=vlevel)

            loss_sps[k] = aresult.loss_sparse_coding
            loss_dus[k] = aresult.loss_dictionary_update

            params_dic = aresult.params_dictionary
        end
        if vlevel >= 1
            println("--- epoch #$itr, sum(loss) = $(sum(loss_sps)), var(loss) = $(var(loss_sps))")
        end
        # plot_function(target)

        savelogs(logdir, target, itr, log_configs...;
            vlevel=vlevel,
            loss_sparse_coding=loss_sps,
            loss_dictionary_update=loss_dus)
    end
    vlevel >= 1 && println("training finished.")
    return setParamsDictionary!(target, params_dic)
end

function trainPerData!(target, x, params_dic_init, epoch;
                        sparsecoder=SparseCoders.IHT,
                        sparsecoder_options=(),
                        optimizer=Optimizers.Steepest,
                        optimizer_options=(),
                        shape, vlevel)

    vlevel >= 3 && println("start Sparse Coding Stage.")
    sparse_coefs, loss_sp = stepSparseCoding(sparsecoder, sparsecoder_options, target, x; epoch=epoch, shape=shape, vlevel=vlevel)
    vlevel >= 3 && println("end Sparse Coding Stage.")

    vlevel >= 3 && println("start Dictionary Update.")
    params_dic, loss_du = updateDictionary(optimizer, optimizer_options, target, x, sparse_coefs, params_dic_init; epoch=epoch, shape=shape, vlevel=vlevel)
    vlevel >= 3 && println("end Dictionary Update Stage.")

    vlevel >= 2 && println("epoch #$itr, data #$k: loss(Sparse coding) = $(loss_sps[k]), loss(Dic. update) = $(loss_dus[k]).")

    setParamsDictionary!(target, params_dic)

    return (params_dictionary = params_dic,
            sparse_coefs = sparse_coefs,
            loss_sparse_coding = loss_sp,
            loss_dictionary_update = loss_du,)
end

gettrainingdata(filename::AbstractString) = gettrainingdata(FileIO.load(filename))
gettrainingdata(td::AbstractArray{T}) where {T<:AbstractFloat} = td
gettrainingdata(td::AbstractArray{Complex{T}}) where {T<:AbstractFloat} = td
function gettrainingdata(td::AbstractArray{T,D}, TP::Type=Float64) where {D,T<:Color}
    channelview(td) .|> TP
end

getvalidshape(shape::Shapes.AbstractShape, args...) = shape
getvalidshape(::Nothing, ::LearningTarget, sc, du, x::AbstractArray) = Shapes.Vec(size(x))
getvalidshape(::Nothing, ::Multiscale, ::SparseCoders.IHT{T}, args...) where {N,T<:NTuple{N}} = Shapes.Arrayed()

verbosenames() = Dict(:none => 0, :standard => 1, :specified => 2, :loquacious => 3)
verboselevel(sym::Symbol) = verbosenames()[sym]
verboselevel(lv::Integer) = lv

display_filters(::Val{false}, args...; kwargs...) = nothing

display_filters(::Val{true}, nsolt::AbstractNsolt) = display(plot(nsolt))

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
using CodecLz4
using CodecZstd

uncompress_cfst(path) = begin
	cfst_file = open(path, "r")

	# first few header bytes
	headerhash = read(cfst_file, UInt32)
	block_size = Int(read(cfst_file, UInt32))
	version = read(cfst_file, UInt32) |> Int
	compression_algorithm = read(cfst_file, UInt32) & 2147483647

	vec_length = read(cfst_file, UInt64)

	vec_length |> Int

	# calculate the number of blocks
	nr_of_blocks = Int(1 + div(vec_length - 1, block_size))

	block_hash = read(cfst_file, UInt64)

	block_offsets = Vector{UInt64}(undef, 1 + nr_of_blocks)
	# read block offsets
	read!(cfst_file, block_offsets)

	# the first number is the hash
	block_offsets_int = Int.(block_offsets)

	lo = block_offsets_int[1]
	hi = block_offsets_int[2]

	compressed_chunk = read(cfst_file, hi-lo)

	# LZ4
	if compression_algorithm == 1
		fnl_res =  transcode(LZ4Decompressor, compressed_chunk)
		# Blosc.set_compressor("lz4hc")
		# fnl_res = Blosc.decompress(UInt8, compressed_chunk)
		return compressed_chunk
	else
		fnl_res =  transcode(ZstdDecompressor, compressed_chunk)
	end

	for (lo, hi) in zip(@view(block_offsets_int[2:end-1]), @view(block_offsets_int[3:end]))
		compressed_chunk = read(cfst_file, hi-lo)
		if compression_algorithm == 1
			fnl_res = vcat(fnl_res, Blosc.decompress(UInt8, compressed_chunk))
		else
			fnl_res = vcat(fnl_res, transcode(ZstdDecompressor, compressed_chunk))
		end
	end

	close(cfst_file)
	fnl_res
end


text:
function grammar_simple_calculator(; erc=false)
    grammar = erc ?
        (@grammar begin
            R = R + R
            R = R - R
            R = R * R
            R = R / R
            R = _(rand(1:9))
        end) :
        (@grammar begin
            R = R + R
            R = R - R
            R = R * R
            R = R / R
            R = |(1:9)
        end)
end

function loss_approximate_pi(S::SymbolTable, tree::RuleNode, grammar::Grammar)
    ex = get_executable(tree, grammar)
    value = Core.eval(S, ex)
    if isinf(value) || isnan(value)
        return Inf
    end
    Δ = abs(value - π)
    return log(Δ) + length(tree) / 1e4
end

function run_approximate_pi(seed::Int=0, n_pop::Int=300)
    Random.seed!(seed)

    grammar = grammar_simple_calculator(erc=false)
    S = SymbolTable(grammar, ExprOptimization.Benchmarks)
    loss(tree::RuleNode, grammar::Grammar) = loss_approximate_pi(S, tree, grammar)

    init_method = RandomInit()
    select_method = TournamentSelection(3)
    p = GeneticProgram(n_pop, #pop_size
                       40,  #iterations
                       10,  #max_depth
                       0.4, #p_reproduction
                       0.5, #p_crossover
                       0.1, #p_mutation
                       init_method=init_method,
                       select_method=select_method)
    tstart = CPUtime_us()
    result = optimize(p, grammar, :R, loss)
    dt = CPUtime_us() - tstart
    result, dt
end


text:
###
### Functions for Epd based tests
###
function _L!(r::AbstractArray{T}, y::AbstractArray{T}, μ::T, σ::T, p::T) where {T <: Real}
    for i in 1:length(y)
        @inbounds r[i] = gamma(1 + 1/p) * abs(y[i] - μ) / σ
    end
    r
end

function logOrZero(x::Real)
    x == 0. ? 0 : log(x)
end

function effScore(y::AbstractArray{<:T}, μ::T, σ::T, p::T) where {T <: Real}
    r = Vector{T}(undef, length(y))
    _L!(r, y, μ, σ, p)
    Sₚ = r.^p .* (1/p * digamma(1 + 1/p) .- logOrZero.(r))
    Sₛ = p/σ .* r.^p .- 1/σ
    (Sₚ, Sₛ)
end

"""
    epdTest(y, μ, σ, p)

Computes the efficient score function wrt to the shape parameter, Eq. 8 of PAPER

```math
S_{n,p}^*
```

# Arguments
- `y::Array{<:Real, 1}`: data
- `μ::Real`: location parameter, μ ∈ ℝ
- `σ::Real`: scale parameter, σ > 0
- `p::Real`: shape parameter, p > 0
"""
function epdTest(y::AbstractArray{<:Real}, μ::Real, σ::Real, p::Real)
    n = length(y)
    S_p, S_σ = effScore(y, μ, σ, p)
    return sum(S_p .+ σ / 4 .* S_σ) / √(n*((p+1)/p^4 * trigamma((p+1)/p) - 1/p^3))
end

function loglikEPD(θ, p, x) where {T <: Real}
    μ, σ = θ
    -log.(pdf.(Epd(μ, exp(σ), p), x)) |> sum
end

function MLE(θ::AbstractVector{T}, p::T, x::Array{T, 1}) where {T <: Real}
    length(θ) === 2 || throw(ArgumentError("θ not of length 2"))
    func = TwiceDifferentiable(vars -> loglikEPD(vars, p, x), ones(2), autodiff =:forward)
    optimum = optimize(func, θ)
    Optim.converged(optimum) || throw(ConvergenceError("Optimizer did not converge"))
    mle = Optim.minimizer(optimum)
    mle[2] = exp(mle[2])
    mle
end

"""
    BivariateNormalTest(X, μ, Σ)

Computes the bivariate normal test based on the multivariate EPD

# Arguments
- `X::Array{<:Real, 2}`: data, n × p
- `μ::Array{<:Real, 1}`: location parameter, μ ∈ ℝᴾ
- `Σ::Array{<:Real, 2}`: Inverse of covariance matrix
"""
function BivariateNormalTest(X::Array{<:Real, 2}, μ::Array{<:Real, 1}, Σ::Array{<:Real, 2})
    n,p = size(X)
    p == length(μ) || throw(DomainError(μ, "μ and X dimension mismatch"))
    u = zeros(n)
    for i in 1:n
        u[i] = dot((X[i,:] - μ)' * S, (X[i,:] - μ))[1,1]
    end
    return √n * (1 + digamma(1) + log(2) - 0.5 * mean(u.* log.(u))) / √((π^2 - 9)/3)
end

"""
    BivariateNormalTest(X)

Computes the bivariate normal test based on the multivariate EPD

# Arguments
- `X::Array{<:Real, 2}`: data, n × p
"""
function BivariateNormalTest(X::Array{<:Real, 2})
    n,p = size(X)
    u = zeros(n)
    S = cov(X)^(-1)
    μ = reshape(sum(X, dims = 1)./n, 2)
    for i in 1:n
        u[i] = dot((X[i,:] - μ)' * S, (X[i,:] - μ))[1,1]
    end
    return √n * (1 + digamma(1) + log(2) - 0.5 * mean(u.* log.(u))) / √((π^2 - 9)/3)
end

function standard(X::AbstractMatrix{<:Real})
    n,d = size(X)
    S = inv(sqrt((n-1)/n * cov(X)))
    (S*(X .- mean(X, dims = 1))')'
end

"""
    DEHU(X)

Computes the test of multivariate normality based on the characteristic function

# Arguments
- `X::AbstractMatrix{<:Real}`: data, n × p
- `a::Real`: parameter > 0
"""
function DEHU(X::AbstractMatrix{<:Real}, a::Real=1)
    n,d = size(X)
    n > d || throw(ArgumentError("Rows of X must be larger than columns"))
    X = standard(X)
    Djk = X * X'
    Rquad = diag(Djk)
    # not sure how to do this the best way
    Dj = zeros(n,n)
    for i in 1:n
        Dj[:, i] = Rquad
    end
    Hilf = Dj - 2 .*Djk+Dj'

    Y=(Rquad*Rquad' - (Dj+Dj').*(Hilf.+2*d*a*(2*a-1))/(4*a^2)+
        (16*d^2*a^3*(a-1)+4*d*(d+2)*a^2 .+(8*d*a^2-(8+4*d)*a).*Hilf+Hilf.^2)/(16*a^4)).*exp.(-Hilf/(4*a))
    ret = (π/a)^(d/2) * sum(Y) / n
    ret*d^(-2) * (a/π)^(d/2)
end

"""
    JB(X)

Computes the Jarque-Bera test of Mardia
# Arguments
- `X::AbstractMatrix{<:Real}`: data, n × p
"""
function JB(X::AbstractMatrix{<:Real})
    n,p = size(X)
    μ = reshape(sum(X, dims = 1)./n, 2)
    S = cov(X)^(-1)
    u = zeros(n)
    for i in 1:n
        u[i] = dot((X[i,:] - μ)' * S, (X[i,:] - μ))[1,1]
    end
    √(n/64) * (mean(u.^2) - 8)
end


text:
@testset "Image tests" begin

end


text:
abstract type Metrics end
mutable struct Unit_e <: Metrics
end
mutable struct Diag_e <: Metrics
end

abstract type SamplingAlgorithm end
mutable struct Hmc <: SamplingAlgorithm
  # ...
  metric::Metrics
end
mutable struct Fixed_param <: SamplingAlgorithm
end

Hmc() = Hmc(Diag_e())
Hmc(m::DataType) = ( m <: Metrics ? Hmc(m()) : "Error: $(m) not a subtype of Metrics" )
Hmc(m::Symbol) = ( eval(m) <: Metrics ? Hmc(eval(m)()) : "Error: $(m) not a subtype of Metrics" )

h = Hmc()
l = Hmc(Diag_e)
k = Hmc(:Unit_e)

@assert isa(h.metric, Metrics)
@assert isa(l.metric, Metrics)
@assert isa(k.metric, Metrics)
@assert typeof(k.metric) <: Metrics


text:
function predict(net::Vector{Layer},input)
  net[1].rₒ=input
  for idx=1:length(net)-1
    forward!(net[idx])
  end
  ~,index=findmax(net[end-1].rₒ,1) #regression=>net[end-1].rₒ
  return map(x->mod1(x,size(net[end-1].rₒ,1)),index)
end


text:
dummy_state(stack) = Debugger.DebuggerState(stack, nothing)

# Steps through the whole expression using `s`
function step_through(frame)
    state = dummy_state(frame)
    while state.frame !== nothing
        execute_command(state, Val{:s}(), "s")
    end
    return state.overall_result
end


text:
module TinyAutoGrad

include("dualnumbers.jl")

end


text:
using Plots

# Model
include_model("planar_push_block_v2")

# Horizon
T = 26

# Time step
tf = 2.5
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[model.idx_u[1:2:8]] .= 1.0
_uu[model.idx_u[2:2:9]] .= 0.5 * π
_uu[model.idx_λ] .= 1.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[model.idx_u[1:2:8]] .= 0.0
_ul[model.idx_u[2:2:9]] .= -0.5 * π
_ul[model.idx_λ] .= 1.0

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.0, 0.0]
x1 = [q1; q1]
qT = [1.0; 0.0; 0.0]
xT = [qT; qT]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
include_objective(["velocity", "nonlinear_stage"])
obj_velocity = velocity_objective(
    [t > T / 2 ? Diagonal(10.0 * ones(model.nq)) : Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3]))

obj_tracking = quadratic_tracking_objective(
    [Diagonal(1.0 * ones(model.n)) for t = 1:T],
    # [Diagonal(0.1 * ones(model.m)) for t = 1:T-1],
	[Diagonal([0.1 * ones(8);
		zeros(4);
		zeros(model.nc);
		ones(model.nb);
		zeros(model.m - model.nu - model.nc - model.nb)]) for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj_penalty = PenaltyObjective(1.0e6, model.m)

function cost_l1(x, u, t)
	J = 0.0

	if t < T
		J += 1000.0 * sum(u[model.idx_u][9:12])
	end

	return J
end

cost_l1(x) = cost_l1(x, nothing, T)
obj_l1 = nonlinear_stage_objective(cost_l1, cost_l1)

obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity, obj_l1])

# Constraints
include_constraints(["contact", "stage"])
t_idx = vcat([t for t = 1:T-1])

function con_l1!(c, x, u, t)

	u_f = u[model.idx_u][1:2:8]
	u_l1 = u[model.idx_u][9:12]

	c[1:4] = u_l1 - u_f
	c[5:8] = u_f + u_l1

    nothing
end

n_l1 = 8
con_l1 = stage_constraints(con_l1!, n_l1, (1:8), t_idx)

# function control_limits_con!(c, x, u, t)
# 	q = x[nq .+ (1:nq)]
# 	θ = q[3]
# 	R = rotation_matrix(θ)
# 	u_ctrl = u[model.idx_u]
#
# 	c[1] = -1.0 * (R' * u_ctrl[1:2])[1]
# 	c[2] = (model.μ[end] * u_ctrl[1])^2.0 - u_ctrl[2]^2.0
#
# 	c[3] = -1.0 * (R' * u_ctrl[3:4])[2]
# 	c[4] = (model.μ[end] * u_ctrl[4])^2.0 - u_ctrl[3]^2.0
#
# 	c[5] = (R' * u_ctrl[5:6])[1]
# 	c[6] = (model.μ[end] * u_ctrl[5])^2.0 - u_ctrl[6]^2.0
#
# 	c[7] = (R' * u_ctrl[7:8])[2]
# 	c[8] = (model.μ[end] * u_ctrl[8])^2.0 - u_ctrl[7]^2.0
# end
#
# n_ctrl_lim = 8
# con_ctrl_lim = stage_constraints(control_limits_con!, n_ctrl_lim, (1:8), t_idx)

con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact, con_l1])#, con_ctrl_lim])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
#open(vis)
visualize!(vis, model,
    q̄, ū,
    Δt = h)

# settransform!(vis["/Cameras/default"],
# 	compose(Translation(0.0, 0.0, 50.0), LinearMap(RotZ(0.5 * pi) * RotY(-pi/2.5))))
# setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)

plot(hcat(q̄...)', color = :red, width = 1.0, labels = "")
plot(hcat([ū..., ū[end]]...)[model.idx_u[1:2:8], :]', linetype = :steppost, color = :red, width = 1.0, labels = "impulse magnitude")

[ū[t][model.idx_u[1:2:8]] for t = 1:T-1]

max_issues_repo_name: thowell/DirectMotionPlanning


text:
target = "Dtree/libdtree.$(Libdl.dlext)"
vers = "0.0.1"

if !isfile(target)
    @static if is_linux()
        LibGit2.clone("https://github.com/kpamnany/Dtree", "Dtree")
        println("Compiling libdtree...")
        run(`make -C Dtree`)
    else
	error("Dtree is Linux-only right now")
    end
end