text:
export result_quad_force,result_quad_displacement

"""
    result_quad_force(assembly,lc_id,elm_id)
读取已求解的assembly实例的单个quad单元内力
# Arguments
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `elm_id`: 单元id
# 返回
- 包含梁单元坐标系下的12个节点力向量
"""
function result_quad_force(assembly,lc_id,elm_id)
    path=assembly.working_path
    elm_id=string(elm_id)
    if !(id in keys(assembly.structure.quads))
        throw("quad id "*string(quad_id)*" doesn't exists!")
    end
    lc_id=string(lc_id)
    if !(lc_id in keys(assembly.lcset.cases))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    elm=assembly.structure.quads[quad_ids]
    i=elm.node1.hid
    j=elm.node2.hid
    k=elm.node3.hid
    l=elm.node4.hid
    iDOFs=6*i-5:6*i
    jDOFs=6*j-5:6*j
    kDOFs=6*k-5:6*k
    lDOFs=6*l-5:6*l
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    T=elm.T
    uᵉ=T*[u[iDOFs];u[jDOFs];u[kDOFs];u[lDOFs]]
    Kᵉ=elm.Kᵉ
    return Kᵉ*uᵉ
end

"""
    result_quad_displacement(assembly,lc_id,elm_id)
读取已求解的assembly实例的单个quad单元节点位移
# 参数
- `assembly`: Assembly类型实例
- `lc_id`: 工况id
- `elm_id`: 单元id
# 返回
- 包含quad单元坐标系下的12个节点位移向量
"""
function result_quad_displacement(assembly,lc_id,elm_id)
    path=assembly.working_path
    elm_id=string(elm_id)
    if !(id in keys(assembly.structure.quads))
        throw("quad id "*string(quad_id)*" doesn't exists!")
    end
    lc_id=string(lc_id)
    if !(lc_id in keys(assembly.lcset.cases))
        throw("loadcase with id "*lc_id*" doesn't exist!")
    end
    loadcase=assembly.lcset.cases[lc_id]
    u=read_vector(path*"/.analysis",lc_id*"_u.v")
    elm=assembly.structure.quads[quad_id]
    i=elm.node1.hid
    j=elm.node2.hid
    k=elm.node3.hid
    l=elm.node4.hid
    iDOFs=6*i-5:6*i
    jDOFs=6*j-5:6*j
    kDOFs=6*k-5:6*k
    lDOFs=6*l-5:6*l
    T=elm.T
    uᵉ=T*[u[iDOFs];u[jDOFs];u[kDOFs];u[lDOFs]]
    return ūᵉ
end

function result_element()
    Δdᵉ=(elm.T)'*[Δd[6*hid1-5:6*hid1];Δd[6*hid2-5:6*hid2]]
    Δϵᵉ=Bᵉ*Δdᵉ
    Δσᵉ=∫D*Δϵᵉ
    σᵉ+=Δσᵉ
end


text:
# test for gmi

using gmi
using PumiConfig

function test_gmi()

  if !HAVE_SIMMETRIX
    return nothing
  end

  @testset "GMI" begin

    fname = "meshes/airfoil.smd"
    gmi.sim_start()
    gmi.register_sim()

    g = gmi.load(fname)

    it = gmi._gmi_begin(g, 1)
    ge = gmi._gmi_next(g, it)
    edge_tags = Array{Cint}(0)
    while ( ge != gmi.NullModelEntity)
      @test gmi.dim(g, ge) == 1
      push!(edge_tags, gmi.tag(g, ge))

      ge = gmi._gmi_next(g, it)
    end

    gmi.gmi_end(g, it)
    sort!(edge_tags)
    @test edge_tags == Cint[5, 8]

    # test iteration on Gmi_iter
    edge_tags2 = Array{Cint}(0)
    edges = Array{gmi.ModelEntity}(0)
    for ge in gmi.Gmi_iter(g, 1)
      @test gmi.dim(g, ge) == 1
      push!(edge_tags2, gmi.tag(g, ge))
      push!(edges, ge)
    end

    _edge_tags = sort(edge_tags2)
    @test _edge_tags == Cint[5, 8]

    @test eltype(gmi.Gmi_iter) == gmi.ModelEntity

    for i=1:2
      ge = gmi.find(g, 1, edge_tags2[1])
      @test ge == edges[1]
    end

    # edge 8 is the free stream (circle)
    circle = gmi.find(g, 1, 8)
    r = zeros(Float64, 2)
    gmi.range(g, circle, 0, r)
    @test gmi.periodic(g, circle, 0)
    @test gmi.can_eval(g)
    # the circle is centered at the origin, so traversing half the parameter
    # range should rotate the point 180 degrees about the origin
    p1 = [0.25, 0.0]
    p2 = p1 + [ (r[2]-r[1])/2, 0]
    x1 = zeros(Float64, 3)
    x2 = zeros(Float64, 3)
    gmi.geval(g, circle, p1, x1)
    gmi.geval(g, circle, p2, x2)
    println("x1 = ", x1)
    println("x2 = ", x2)
    @test maximum(abs.(x1 + x2)) < 1e-13
    @test abs(norm(x1) - 10) < 1e-12
    @test abs(norm(x2) - 10) < 1e-12

    # lets see what happens when we go beyond the parametric range for a
    # periodic entity
    # It appears for p values greater than r[2], it returns the x values at
    # r[2].  It also appears the parameter is not linear
    # p = 0 -> (10.0, 0)
    # p = 0.25 -> (0, 10.0)
    # p = 0.125 -> (0.8, 0.6) ???
    #=
    p3 = copy(p1)
    p3[1] += 1
    x3 = zeros(Float64, 3)
    gmi.geval(g, circle, p3, x3)
    println("p1 = ", p1)
    println("p3 = ", p3)
    println("x1 = ", x1)
    println("x3 = ", x3)

    p0 = zeros(2)
    x0 = zeros(3)
    gmi.geval(g, circle, p0, x0)
    println("p0 = ", p0)
    println("x0 = ", x0)
    =#

    # test reparam
    p0 = [0.125, 0.0]
    x0 = zeros(3)
    gmi.geval(g, circle, p0, x0)
    faces = gmi.adjacent(g, circle, 2)
    @test length(faces.v) == 1
    face = faces.v[1]

    p1 = zeros(2)
    x1 = zeros(3)
    gmi.reparam(g, circle, p0, face, p1)
    gmi.geval(g, face, p1, x1)

    @test maximum(abs.(x1 - x0)) < 1e-10  # face evaluations are not so accurate

    # test closest point, normal
    x0 = zeros(3)
    x1 = zeros(3)
    p1 = zeros(2)
    gmi.closest_point(g, circle, x0, x1, p1)
    @test (norm(x1) - 10) < 1e-10

    # normal only works for geometric faces
    p = zeros(2)
    n = zeros(3)
    gmi.normal(g, face, p, n)
    println("n = ", n)
    @test abs(norm(n) - 1) < 1e-13
    @test abs(n[3] -1) < 1e-13

    # testing first_derivative is tricky, because the parameter is not
    # theta, do finite difference instead
    p = zeros(2)
    t0 = zeros(3)
    t1 = zeros(3)
    gmi.first_derivative(g, circle, p, t0, t1)

    h = 1e-6
    p1 = [h, 0.0]
    x0 = zeros(3)
    x1 = zeros(3)
    gmi.geval(g, circle, p, x0)
    gmi.geval(g, circle, p1, x1)
    t0_fd = (x1 - x0)/h
    @test maximum(abs.(t0_fd - t0)) < 1e-5*maximum(t0_fd)

    # is_in_closure_of
    @test gmi.is_in_closure_of(g, circle, face)

    gmi.destroy(g)
  end

end

ext: jl


text:
using MINLPTests, JuMP, Ipopt, Juniper, Test

const OPTIMIZER = MINLPTests.JuMP.with_optimizer(
    Juniper.Optimizer, nl_solver=with_optimizer(Ipopt.Optimizer, print_level=0), atol=1e-7
)

@testset "MINLPTests" begin
    ###
    ### src/nlp-mi tests.
    ###

    MINLPTests.test_nlp_mi(OPTIMIZER)
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, ParameterTests

makedocs(
    modules = [ParameterTests],
    clean = false,
    format = Documenter.HTML(),
    sitename = "ParameterTests.jl",
    authors = "Iain Skett",
    pages = [
        "Introduction" => "index.md"
        "Getting Started" => "getting_started.md"
    ],
)

deploydocs(
    repo = "github.com/InfiniteChai/ParameterTests.jl.git",
    devbranch = "main"
)


text:
using Documenter
using CorpusLoaders

makedocs(modules=[CorpusLoaders])

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
			repo = "github.com/oxinabox/CorpusLoaders.jl.git",
			osname = "linux")


text:
# -*- coding: utf-8 -*-
# $$
# \frac{d f}{dt} +  (v \frac{d f}{dx} - x \frac{d f}{dv}) = 0
# $$

# $$
# x \in [-\pi, \pi],\qquad y \in [-\pi, \pi] \qquad \mbox{ and } \qquad t \in [0, 200\pi]
# $$

using  FFTW
using  LinearAlgebra
using  Plots, ProgressMeter
using  BenchmarkTools

# ## Mesh parameters (matlab code)
#
# ```m
# Nx=128;Ny=256;
# xmin=-pi; xmax=pi;
#
# dx=(xmax-xmin)/Nx;
# x=xmin:dx:xmax-dx;
#
# ymin=-pi; ymax=pi;
# dy=(ymax-ymin)/Ny;
# y=ymin:dy:ymax-dy;
# ```

# ## Julia type for mesh information
#
# ```jl
# struct Mesh
#     nx   :: Int
#     ny   :: Int
#     xmin :: Float64
#     xmax :: Float64
#     ymin :: Float64
#     ymax :: Float64
#     dx   :: Float64
#     dy   :: Float64
#     x    :: Vector{Float64}
#     y    :: Vector{Float64}
# end
#
# mesh = Mesh( 128, 256, -π, π, -π, π, 2π/128, 2π/256, ...)
#
# ```
#

struct Mesh

    nx   :: Int
    ny   :: Int
    xmin :: Float64
    xmax :: Float64
    ymin :: Float64
    ymax :: Float64
    dx   :: Float64
    dy   :: Float64
    x    :: Vector{Float64}
    y    :: Vector{Float64}

    function Mesh( xmin, xmax, nx, ymin, ymax, ny)
        dx, dy = (xmax-xmin)/nx, (ymax-ymin)/ny
        x = range(xmin, stop=xmax, length=nx+1)[1:end-1]  ## we remove the end point
        y = range(ymin, stop=ymax, length=ny+1)[1:end-1]  ## periodic boundary condition
        new( nx, ny, xmin, xmax, ymin, ymax, dx, dy, x, y)
    end
end

mesh = Mesh(-π, π, 128, -π, π, 256)

@show mesh.xmin, mesh.xmax, mesh.nx, mesh.dx

# # Initialization of f : 2d array of double float
#
# ```m
# f=zeros(Nx,Ny);
# for i=1:Nx
#     xx=xmin+(i-1)*dx;
#     for j=1:Ny
#         yy=ymin+(j-1)*dy;
#         f(i,j)=exp(-(xx-1)*(xx-1)/0.1)*exp(-(yy-1)*(yy-1)/0.1);
#     end
# end
# ```

# ```jl
# f = zeros(Float64,(mesh.nx,mesh.ny))
#
# for (i, x) in enumerate(mesh.x), (j, y) in enumerate(mesh.y)
#
#     f[i,j] = exp(-(x-1)*(x-1)/0.1)*exp(-(y-1)*(y-1)/0.1)
#
# end
# ```

# ### Julia function to compute exact solution

function exact(tf, mesh)

    f = zeros(Float64,(mesh.nx, mesh.ny))
    for (i, x) in enumerate(mesh.x), (j, y) in enumerate(mesh.y)
        xn = cos(tf)*x - sin(tf)*y
        yn = sin(tf)*x + cos(tf)*y
        f[i,j] = exp(-(xn-1)*(xn-1)/0.1)*exp(-(yn-1)*(yn-1)/0.1)
    end

    f
end

f = exact(0.0, mesh)
surface(f)

# ## Create the gif to show what we are computing

function create_gif()
    x = y = range(-π, stop=π, length=40)
    n = 100

    prog = Progress(n,1) ## progress bar

    @gif for t in range(0, stop=2π, length=n)
        f(x,y) = exp(-((cos(t)*x-sin(t)*y)-1)^2/0.2)*exp(-((sin(t)*x+cos(t)*y)-1)^2/0.2)

        p = plot(x, y, f, st = [:surface])

        plot!(p[1])
        plot!(zlims=(-0.01,1.01))

        next!(prog) ## increment the progress bar
    end
end

create_gif();

# ![](tmp.gif)

# ## Advection and loop over time
# ```m
# tf=200*pi;Nt=1000;dt=tf/Nt;
# kx=2*pi/(xmax-xmin)*[0:Nx/2-1,Nx/2-Nx:Nx-1-Nx];
# ky=2*pi/(ymax-ymin)*[0:Ny/2-1,Ny/2-Ny:Ny-1-Ny];
#
# fnx=zeros(1,Nx);ffx=zeros(1,Nx);fny=zeros(1,Ny);ffy=zeros(1,Ny);
#
# for n=1:Nt
#     for i=1:Nx
#         xx=xmin+(i-1)*dx;
#         ffy=fft(f(i,:));
#         fny=real(ifft(exp(sqrt(-1)*xx*ky*tan(dt/2)).*ffy));
#         f(i,:)=fny;
#     end
#
#     for j=1:Ny
#         yy=ymin+(j-1)*dy;
#         ffx=fft(f(:,j));
#         fnx=real(ifft(exp(-sqrt(-1)*yy*kx*sin(dt)).*transpose(ffx)));
#         f(:,j)=fnx;
#     end
#
#     for i=1:Nx
#         xx=xmin+(i-1)*dx;
#         ffy=fft(f(i,:));
#         fny=real(ifft(exp(sqrt(-1)*xx*ky*tan(dt/2)).*ffy));
#         f(i,:)=fny;
#     end
# end
# ```

# ## Function to compute error
#
# ```m
#
# % compute errors in Linfty norm
# error1=max(max(f-f_exact))
#
# ```
#
# - In julia the max value of an array is `maximum`.

function error1(f, f_exact)
    maximum(abs.(f .- f_exact))
end
#----------------------------------------------------------------------------

# ## Naive translation of the matlab code

function naive_from_matlab(tf, nt, mesh::Mesh)

    dt = tf/nt

    kx = 2π/(mesh.xmax-mesh.xmin)*[0:mesh.nx÷2-1;mesh.nx÷2-mesh.nx:-1]
    ky = 2π/(mesh.ymax-mesh.ymin)*[0:mesh.ny÷2-1;mesh.ny÷2-mesh.ny:-1]

    f = exact(0.0, mesh)

    for n=1:nt

       for (i, x) in enumerate(mesh.x)
           f[i,:]=real(ifft(exp.(1im*x*ky*tan(dt/2)).*fft(f[i,:])))
       end

       for (j, y) in enumerate(mesh.y)
           f[:,j]=real(ifft(exp.(-1im*y*kx*sin(dt)).*fft(f[:,j])))
       end

       for (i, x) in enumerate(mesh.x)
           f[i,:]=real(ifft(exp.(1im*x*ky*tan(dt/2)).*fft(f[i,:])))
       end
   end

    f
end
#----------------------------------------------------------------------------

nt, tf = 1000, 200
println( " error = ", error1(naive_from_matlab(tf, nt, mesh), exact(tf, mesh)))
@btime naive_from_matlab(tf, nt, mesh);
#----------------------------------------------------------------------------

# ###  Vectorized version
#
# - We remove the for loops over direction x and y by creating the 2d arrays `exky` and `ekxy`.
# - We save cpu time by computing them before the loop over time

function vectorized(tf, nt, mesh::Mesh)

    dt = tf/nt

    kx = 2π/(mesh.xmax-mesh.xmin)*[0:mesh.nx÷2-1;mesh.nx÷2-mesh.nx:-1]
    ky = 2π/(mesh.ymax-mesh.ymin)*[0:mesh.ny÷2-1;mesh.ny÷2-mesh.ny:-1]

    f = exact(0.0, mesh)

    exky = exp.( 1im*tan(dt/2) .* mesh.x  .* ky')
    ekxy = exp.(-1im*sin(dt)   .* mesh.y' .* kx )

    for n = 1:nt
        f = real(ifft(exky .* fft(f, 2), 2))
        f = real(ifft(ekxy .* fft(f, 1), 1))
        f = real(ifft(exky .* fft(f, 2), 2))
    end

    f
end
#----------------------------------------------------------------------------

nt, tf = 1000, 200
println( " error = ", error1(vectorized(tf, nt, mesh), exact(tf, mesh)))
@btime vectorized(tf, nt, mesh);
#----------------------------------------------------------------------------

# ## Inplace computation
# - We remove the Float64-Complex128 conversion by allocating the distribution function `f` as a Complex array
# - Note that we need to use the inplace assignement operator ".="  to initialize the `f` array.
# - We use inplace computation for fft with the "bang" operator `!`

function inplace(tf, nt, mesh::Mesh)

    dt = tf/nt

    kx = 2π/(mesh.xmax-mesh.xmin)*[0:mesh.nx÷2-1;mesh.nx÷2-mesh.nx:-1]
    ky = 2π/(mesh.ymax-mesh.ymin)*[0:mesh.ny÷2-1;mesh.ny÷2-mesh.ny:-1]

    f  = zeros(Complex{Float64},(mesh.nx,mesh.ny))
    f .= exact(0.0, mesh)

    exky = exp.( 1im*tan(dt/2) .* mesh.x  .* ky')
    ekxy = exp.(-1im*sin(dt)   .* mesh.y' .* kx )

    for n = 1:nt
        fft!(f, 2)
        f .= exky .* f
        ifft!(f,2)
        fft!(f, 1)
        f .= ekxy .* f
        ifft!(f, 1)
        fft!(f, 2)
        f .= exky .* f
        ifft!(f,2)
    end

    real(f)
end
#----------------------------------------------------------------------------

nt, tf = 1000, 200
println( " error = ", error1(inplace(tf, nt, mesh), exact(tf, mesh)))
@btime inplace(tf, nt, mesh);
#----------------------------------------------------------------------------

# ### Use plans for fft

# - When you apply multiple fft on array with same shape and size, it is recommended to use fftw plan to improve computations.
# - Let's try to initialize our two fft along x and y with plans.

function with_fft_plans(tf, nt, mesh::Mesh)

    dt = tf/nt

    kx = 2π/(mesh.xmax-mesh.xmin)*[0:mesh.nx÷2-1;mesh.nx÷2-mesh.nx:-1]
    ky = 2π/(mesh.ymax-mesh.ymin)*[0:mesh.ny÷2-1;mesh.ny÷2-mesh.ny:-1]

    f  = zeros(Complex{Float64},(mesh.nx,mesh.ny))
    f .= exact(0.0, mesh)
    f̂  = similar(f)

    exky = exp.( 1im*tan(dt/2) .* mesh.x  .* ky')
    ekxy = exp.(-1im*sin(dt)   .* mesh.y' .* kx )

    Px = plan_fft(f, 1)
    Py = plan_fft(f, 2)

    for n = 1:nt

        f̂ .= Py * f
        f̂ .= f̂  .* exky
        f .= Py \ f̂

        f̂ .= Px * f
        f̂ .= f̂  .* ekxy
        f .= Px \ f̂

        f̂ .= Py * f
        f̂ .= f̂  .* exky
        f .= Py \ f̂

    end

    real(f)
end
#----------------------------------------------------------------------------

nt, tf = 1000, 200
println( " error = ", error1(with_fft_plans(tf, nt, mesh), exact(tf, mesh)))
@btime with_fft_plans(tf, nt, mesh);
#----------------------------------------------------------------------------

# ## Inplace computation and fft plans
#
# To apply fft plan to an array A, we use a preallocated output array Â by calling `mul!(Â, plan, A)`.
# The input array A must be a complex floating-point array like the output Â.
# The inverse-transform is computed inplace by applying `inv(P)` with `ldiv!(A, P, Â)`.

function with_fft_plans_inplace(tf, nt, mesh::Mesh)

    dt = tf/nt

    kx = 2π/(mesh.xmax-mesh.xmin)*[0:mesh.nx÷2-1;mesh.nx÷2-mesh.nx:-1]
    ky = 2π/(mesh.ymax-mesh.ymin)*[0:mesh.ny÷2-1;mesh.ny÷2-mesh.ny:-1]

    f  = zeros(Complex{Float64},(mesh.nx,mesh.ny))
    f .= exact(0.0, mesh)
    f̂  = similar(f)

    exky = exp.( 1im*tan(dt/2) .* mesh.x  .* ky')
    ekxy = exp.(-1im*sin(dt)   .* mesh.y' .* kx )

    Px = plan_fft(f, 1)
    Py = plan_fft(f, 2)

    for n = 1:nt

        mul!(f̂, Py, f)
        f̂ .= f̂ .* exky
        ldiv!(f, Py, f̂)

        mul!(f̂, Px, f)
        f̂ .= f̂ .* ekxy
        ldiv!(f, Px, f̂)

        mul!(f̂, Py, f)
        f̂ .= f̂ .* exky
        ldiv!(f, Py, f̂)

    end

    real(f)
end
#----------------------------------------------------------------------------

nt, tf = 1000, 200
println( " error = ", error1(with_fft_plans_inplace(tf, nt, mesh), exact(tf, mesh)))
@btime with_fft_plans_inplace(tf, nt, mesh);
#----------------------------------------------------------------------------

# ## Explicit transpose of `f`
#
# - Multidimensional arrays in Julia are stored in column-major order.
# - FFTs along y are slower than FFTs along x
# - We can speed-up the computation by allocating the transposed `f`
# and transpose f for each advection along y.

function with_fft_transposed(tf, nt, mesh::Mesh)

    dt = tf/nt

    kx = 2π/(mesh.xmax-mesh.xmin)*[0:mesh.nx÷2-1;mesh.nx÷2-mesh.nx:-1]
    ky = 2π/(mesh.ymax-mesh.ymin)*[0:mesh.ny÷2-1;mesh.ny÷2-mesh.ny:-1]

    f  = zeros(Complex{Float64},(mesh.nx,mesh.ny))
    f̂  = similar(f)
    fᵗ = zeros(Complex{Float64},(mesh.ny,mesh.nx))
    f̂ᵗ = similar(fᵗ)

    exky = exp.( 1im*tan(dt/2) .* mesh.x' .* ky )
    ekxy = exp.(-1im*sin(dt)   .* mesh.y' .* kx )

    FFTW.set_num_threads(4)
    Px = plan_fft(f,  1, flags=FFTW.PATIENT)
    Py = plan_fft(fᵗ, 1, flags=FFTW.PATIENT)

    f .= exact(0.0, mesh)

    for n = 1:nt
        transpose!(fᵗ,f)
        mul!(f̂ᵗ, Py, fᵗ)
        f̂ᵗ .= f̂ᵗ .* exky
        ldiv!(fᵗ, Py, f̂ᵗ)
        transpose!(f,fᵗ)

        mul!(f̂, Px, f)
        f̂ .= f̂ .* ekxy
        ldiv!(f, Px, f̂)

        transpose!(fᵗ,f)
        mul!(f̂ᵗ, Py, fᵗ)
        f̂ᵗ .= f̂ᵗ .* exky
        ldiv!(fᵗ, Py, f̂ᵗ)
        transpose!(f,fᵗ)
    end
    real(f)
end
#----------------------------------------------------------------------------

nt, tf = 1000, 200
println( " error = ", error1(with_fft_transposed(tf, nt, mesh), exact(tf, mesh)))
@btime with_fft_transposed(tf, nt, mesh);
#----------------------------------------------------------------------------

tf, nt = 400π, 1000
mesh = Mesh(-π, π, 512, -π, π, 256)

inplace_bench = @benchmark inplace(tf, nt, mesh)
vectorized_bench = @benchmark vectorized(tf, nt, mesh)
with_fft_plans_bench = @benchmark with_fft_plans(tf, nt, mesh)
with_fft_plans_inplace_bench = @benchmark with_fft_plans_inplace(tf, nt, mesh)
with_fft_transposed_bench = @benchmark with_fft_transposed(tf, nt, mesh)
#----------------------------------------------------------------------------

d = Dict()
d["vectorized"] = minimum(vectorized_bench.times) / 1e6
d["inplace"] = minimum(inplace_bench.times) / 1e6
d["with_fft_plans"] = minimum(with_fft_plans_bench.times) / 1e6
d["with_fft_plans_inplace"] = minimum(with_fft_plans_inplace_bench.times) / 1e6
d["with_fft_transposed"] = minimum(with_fft_transposed_bench.times) / 1e6;
#----------------------------------------------------------------------------

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value, digits=1), 6, "."))
end
#----------------------------------------------------------------------------

# ## Conclusion
# - Use pre-allocations of memory and inplace computation are very important
# - Try to always do computation on data contiguous in memory
# - In this notebook, use btime to not taking account of time consumed in
# compilation.

max_stars_repo_name: pnavaro/math-julia

max_issues_repo_name: pnavaro/math-julia

max_forks_repo_name: pnavaro/math-julia


text:
@testset "F" begin
    @testset "fovray" begin
        try
            furnsh(path(CORE, :lsk),
                   path(CORE, :pck),
                   path(CORE, :spk),
                   path(CASSINI, :pck),
                   path(CASSINI, :sat_spk),
                   path(CASSINI, :tour_spk),
                   path(CASSINI, :fk),
                   path(CASSINI, :ck),
                   path(CASSINI, :sclk),
                   path(CASSINI, :ik))
            camid = bodn2c("CASSINI_ISS_NAC")
            shape, frame, bsight, bounds = getfov(camid, 4)
            et = str2et("2013 FEB 25 11:50:00 UTC")
            visible = fovray("CASSINI_ISS_NAC", [0.0, 0.0, 1.0], frame, "S", "CASSINI", et)
            @test visible
            @test_throws ArgumentError fovray("CASSINI_ISS_NAC", [0.0, 0.0], frame, "S", "CASSINI", et)
            @test_throws SpiceError fovray("norbert", [0.0, 0.0, 1.0], frame, "S", "CASSINI", et)
            @test_throws SpiceError fovray("CASSINI_ISS_NAC", [0.0, 0.0, 1.0], frame, "S", "norbert", et)
            @test_throws SpiceError fovray("CASSINI_ISS_NAC", [0.0, 0.0, 1.0], "norbert", "S", "CASSINI", et)
        finally
            kclear()
        end
    end
    @testset "fovtrg" begin
        try
            furnsh(path(CORE, :lsk),
                   path(CORE, :pck),
                   path(CORE, :spk),
                   path(CASSINI, :pck),
                   path(CASSINI, :sat_spk),
                   path(CASSINI, :tour_spk),
                   path(CASSINI, :fk),
                   path(CASSINI, :ck),
                   path(CASSINI, :sclk),
                   path(CASSINI, :ik))
            et = str2et("2013 FEB 25 11:50:00 UTC")
            visible = fovtrg("CASSINI_ISS_NAC", "Enceladus", "Ellipsoid",
                             "IAU_ENCELADUS", "LT+S", "CASSINI", et)
            @test visible
            @test_throws SpiceError fovtrg("norbert", "Enceladus", "Ellipsoid",
                             "IAU_ENCELADUS", "LT+S", "CASSINI", et)
        finally
            kclear()
        end
    end
    @testset "frame" begin
        vec = [23.0, -3.0, 18.0]
        x, y, z = frame(vec)
        expected_x = [0.78338311, -0.10218041, 0.61308243]
        expected_y = [0.61630826, 0.0000000, -0.78750500]
        expected_z = [0.080467580, 0.99476588, 0.062974628]
        @test expected_x ≈ x
        @test expected_y ≈ y
        @test expected_z ≈ z
    end
    @testset "frinfo" begin
        @test frinfo(13000) == (399, 2, 3000)
        @test frinfo(42) === nothing
    end
    @testset "frmnam" begin
        @test frmnam(13000) == "ITRF93"
        @test frmnam(13000) == "ITRF93"
        @test isempty(frmnam(42))
    end
    @testset "furnsh" begin
        try
            furnsh(path(CORE, :pck),path(CORE, :spk),path(CORE, :gm_pck),path(CORE, :lsk))
            @test ktotal("ALL") == 4
        finally
            kclear()
        end
    end
end


text:
@testset "fsm_active_close.jl" begin
    base_seq = WrappingInt32(1 << 31)
    DEFAULT_CAPACITY = 64000
    TIMEOUT_DFLT = 1000
    @testset "START -> SYN_SENT -> SYN/ACK -> ACK" begin
        conn = TCPConnection()
        connect!(conn)
        tick!(conn, 1)
        seg1 = expect_one_seg(conn; ack=false, syn=true)
        expect_state(conn, JLSponge.SYN_SENT)

        isn = WrappingInt32(rand(UInt32))
        send_syn!(conn, isn, seg1.header.seqno + 1)
        tick!(conn, 1)
        expect_state(conn, JLSponge.ESTABLISHED)
        expect_one_seg(conn; ack=true, syn=false, ackno=isn+1)
        @test bytes_in_flight(conn) == 0
    end

    @testset "START -> SYN_SENT -> SYN/ACK -> ACK" begin
        conn = TCPConnection()
        connect!(conn)
        tick!(conn, 1)
        seg1 = expect_one_seg(conn; ack=false, syn=true)
        seg_hdr = seg1.header
        expect_state(conn, JLSponge.SYN_SENT)

        isn = WrappingInt32(rand(UInt32))
        send_syn!(conn, isn)
        tick!(conn, 1)
        expect_one_seg(conn, syn=false, ack=true, ackno=isn + 1)

        expect_state(conn, JLSponge.SYN_RCVD)
        send_ack!(conn, isn, seg_hdr.seqno + 1)
        tick!(conn, 1)
        expect_no_seg(conn)
        expect_state(conn, JLSponge.ESTABLISHED)
    end

    @testset "START -> SYN_SENT -> SYN/ACK -> ESTABLISHED" begin
        conn = TCPConnection()
        connect!(conn)
        tick!(conn, 1)
        seg1 = expect_one_seg(conn; ack=false, syn=true)
        seg_hdr = seg1.header
        expect_state(conn, JLSponge.SYN_SENT)

        isn = WrappingInt32(rand(UInt32))
        send_syn!(conn, isn, seg_hdr.seqno + 1)
        tick!(conn, 1)
        expect_one_seg(conn; ack=true, ackno=isn + 1, syn=false)
        expect_state(conn, JLSponge.ESTABLISHED)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


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

# Need to be refactored and improved

Base.@kwdef mutable struct TestPars
  refimgs::AbstractVector # [(img,prop1),(img,prop2)]
  smooth::AbstractVector{Number}
  variogram::AbstractVector # [(:Y,γ1),(:X,γ2)]

  ratio1::AbstractArray = [(0.2,1.0),(0.5,1.0)]
  ratio2::AbstractArray = [nothing]
  meths::AbstractVector{Symbol} = [:MovingWindows,:KernelConvolution]
  maxneighbors::AbstractArray{Int} = [15,25,40]

  #search::AbstractArray{Search}=[KNN(15),KNN(30),KNN(50)]

  folds::Int = 10
  forcegradients::Bool = false

end

macro name(arg)
   string(arg)
end

function localanisotropies(t::TestPars, problem::EstimationProblem)
    ip = Iterators.product
    vars = [v for (v,V) in variables(problem)]

    bestε = Dict(v => Inf for v in vars)
    optpars = Dict()
    listε = Dict(v => [] for v in vars)

    ids = grid2hd_ids(data(problem),domain(problem))

    for (img,w) in ip(t.refimgs, t.smooth)
        M = typeof(img[1]) <: CartesianGrid || t.forcegradients
        #lpars = M ? gradients(img[1],img[2],w) : geometry(img,w)
        lpars = localanisotropies(Gradients,img[1],img[2],w)

        for (r1,r2) in ip(t.ratio1, t.ratio2)
            lx = rescale_magnitude(lpars, r1, r2)

            for (v,n,e,γ) in ip(vars, t.maxneighbors, t.meths, t.variogram)
                p = (variogram=γ, localaniso=lx, method=e, maxneighbors=n)
                solver = LocalKriging(v => p)
                cv = CrossValidation(t.folds)
                ε = cverror(solver, problem, cv)
                push!(listε[v],ε[v])

                if bestε[v]>ε[v]
                    bestε[v] = ε[v]
                    optpars[v] = p
                end
            end
        end
    end

    optpars
end


text:
module Aqua

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest")
end Aqua

using Base: PkgId, UUID
using Pkg: Pkg, TOML
using Test

try
    findnext('a', "a", 1)
catch
    import Compat
end

const AQUA_VERSION = let path = joinpath(dirname(@__DIR__), "Project.toml")
    proj = TOML.parsefile(path)
    include_dependency(path)
    VersionNumber(proj["version"])
end

_lt05(y::AbstractString, n::AbstractString = "") = AQUA_VERSION < v"0.5-" ? y : n
_ge05(y::AbstractString, n::AbstractString = "") = AQUA_VERSION >= v"0.5-" ? y : n

include("pkg/Versions.jl")
using .Versions: VersionSpec, semver_spec

include("utils.jl")
include("ambiguities.jl")
include("unbound_args.jl")
include("exports.jl")
include("project_extras.jl")
include("stale_deps.jl")
include("deps_compat.jl")
include("project_toml_formatting.jl")

"""
    test_all(testtarget::Module)

Run following tests in isolated testset:

* [`test_ambiguities([testtarget, Base$(_ge05(", Core"))])`](@ref test_ambiguities)
* [`test_unbound_args(testtarget)`](@ref test_unbound_args)
* [`test_undefined_exports(testtarget)`](@ref test_undefined_exports)
* [`test_project_extras(testtarget)`](@ref test_project_extras) $(_lt05("(optional)"))
* [`test_stale_deps(testtarget)`](@ref test_stale_deps) $(_lt05("(optional)"))
* [`test_deps_compat(testtarget)`](@ref test_deps_compat) $(_lt05("(optional)"))
* [`test_project_toml_formatting(testtarget)`](@ref test_project_toml_formatting)
  $(_lt05("(optional)"))

!!! compat "Aqua.jl 0.5"

    Since Aqua.jl 0.5:

    * `test_all` runs [`test_ambiguities`](@ref) with `Core`.  This
       means method ambiguities of constructors can now be detected.
       In Aqua.jl 0.4, `test_ambiguities` was invoked with
       `[testtarget, Base]`.

    * `test_all` runs [`test_stale_deps`](@ref).  In Aqua.jl 0.4, this
      check was opt-in.

The keyword argument `\$x` (e.g., `ambiguities`) can be used to
control whether or not to run `test_\$x` (e.g., `test_ambiguities`).
If `test_\$x` supports keyword arguments, a `NamedTuple` can also be
passed to `\$x` to specify the keyword arguments for `test_\$x`.

# Keyword Arguments
- `ambiguities = true`
- `unbound_args = true`
- `undefined_exports = true`
- `project_extras = $(_lt05("false", "true"))`
- `stale_deps = $(_lt05("false", "true"))`
- `deps_compat = $(_lt05("false", "true"))`
- `project_toml_formatting = $(_lt05("false", "true"))`
"""
function test_all(
    testtarget::Module;
    ambiguities = true,
    unbound_args = true,
    undefined_exports = true,
    project_extras = AQUA_VERSION >= v"0.5-",
    stale_deps = AQUA_VERSION >= v"0.5-",
    deps_compat = AQUA_VERSION >= v"0.5-",
    project_toml_formatting = AQUA_VERSION >= v"0.5-",
)
    @testset "Method ambiguity" begin
        if ambiguities !== false
            if v"1.6.0-DEV.816" <= VERSION < v"1.6.0-DEV.875"
                # Maybe remove this branch?
                @warn "Ignoring ambiguities from `Base` to workaround JuliaLang/julia#36962"
                test_ambiguities([testtarget]; askwargs(ambiguities)...)
            elseif AQUA_VERSION >= v"0.5-"
                test_ambiguities([testtarget, Base, Core]; askwargs(ambiguities)...)
            else
                test_ambiguities([testtarget, Base]; askwargs(ambiguities)...)
            end
        end
    end
    @testset "Unbound type parameters" begin
        if unbound_args !== false
            test_unbound_args(testtarget; askwargs(unbound_args)...)
        end
    end
    @testset "Undefined exports" begin
        if undefined_exports !== false
            test_undefined_exports(testtarget; askwargs(undefined_exports)...)
        end
    end
    @testset "Compare Project.toml and test/Project.toml" begin
        if project_extras !== false
            test_project_extras(testtarget; askwargs(project_extras)...)
        end
    end
    @testset "Stale dependencies" begin
        if stale_deps !== false
            test_stale_deps(testtarget; askwargs(stale_deps)...)
        end
    end
    @testset "Compat bounds" begin
        if deps_compat !== false
            test_deps_compat(testtarget; askwargs(deps_compat)...)
        end
    end
    @testset "Project.toml formatting" begin
        if project_toml_formatting !== false
            test_project_toml_formatting(testtarget; askwargs(project_toml_formatting)...)
        end
    end
end

end # module


text:
#module roundPerRound

function σtot(σ::Float64,ρ::Float64)
    ρ*σ*1e-3
end

function SSHP(projectile::AbstractPenetrator, target::TargetCirc)
    D = abs(projectile.position-target.position)
    Stot=σtot(projectile.σ,D)
    y = (x)->x[1]/Stot^2*exp(-0.5*x[1]^2/Stot^2)
    g=hcubature(y,0,0.5*target.ρ)[1]
    return g
end
"""
    SSHP(target,error)
Returns the single shot hit probability (sshp) for a rectangular target.
"""
function SSHP(target::TargetRect,error::RectErrorB)
    y = x->1/(sqrt(2*pi)*error.σ_x)*exp(-((x[1]-error.μ_x)^2)/(2*error.σ_x^2))*1/(sqrt(2*pi)*error.σ_y)*exp(-((x[2]-error.μ_y)^2)/(2*error.σ_y^2))
    g=hcubature(y,[-target.a,-target.b],[target.a, target.b])[1]
    return g
end

"""
    pHit(target,error)

Returns the probability of hitting the target defined by'target' if the estimated error is provided by 'error'.
```math
\\int_{-a}^{a}\\int_{-b}^{b}\\frac{1}{2\\pi\\sigma_x\\sigma_y}\\text{exp}[\\frac{(x-\\mu_x)^2}{2\\sigma_x^2}\\frac{(y-\\mu_y)^2}{2\\sigma_y^2}]dxdy
```
"""
function pHit(target::AbstractTarget, error::RectErrorB)
    SSHP(target, error)
end
#end
#\\[\\frac{(x-\\mu_x)^2}{2\\sigma_x^2}
#+\\frac{(y-\\mu_y)^2}{2\\sigma_y^2}\\]dxdy\$


text:
# Use Julia's internal syntax trees as the representation for the SuperTux files
function write_scheme(io::IO, elem)
  write(io, elem)
end

function write_scheme(io::IO, expr::Expr)
  if(expr.head == :call || expr.head == :tuple)
    # tuples and calls are written the same here
    write(io, "(")
    for elem in expr.args
      write_scheme(io, elem)
      write(io, " ")
    end
    write(io, ")")
  elseif(expr.head == :block)
    # Scheme doesn't have blocks the way Julia does
    # Just recurse and write the first entry that is an expression
    for arg in expr.args
      if isa(arg, Expr) && arg.head != :line
        write_scheme(io, arg)
        break
      end
    end
  else
    # Unrecognized contructs show up here
    head = expr.head
    write(io, "<$head>")
  end
end

function write_scheme(io::IO, str::String)
  print(io, "\"$str\"")
end

function write_scheme(io::IO, num::Number)
  print(io, num)
end


text:
export Events

module Events

using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

export
    KeyboardEvent, MouseMotionEvent, MouseWheelEvent, MouseButtonEvent,
    get_key_code_sym, get_scancode, get_modifier, get_repeat,
    poll_event!, get_event_type, set!, print

using ..Ranger

import Base.print

# --------------------------------------------------------------------
# Event utilities
# --------------------------------------------------------------------
# returns a UInt32 in host-endian
function extract_4byte_host(start_index, arr::Array{UInt8})
    # First combine into a single 32 bit value
    word = UInt32(arr[start_index]) << 24 |
        UInt32(arr[start_index + 1]) << 16 |
        UInt32(arr[start_index + 2]) << 8 |
        arr[start_index + 3]
    # Convert to host-endian
    # Note: SDL has a note about endianess: https://wiki.libsdl.org/CategoryEndian
    ntoh(word)
end

function extract_2byte_host(start_index, arr::Array{UInt8})
    # First combine into a single 16 bit value
    word = UInt16(arr[start_index]) << 8 |
        arr[start_index + 1]
    ntoh(word)
end

# Grabs the first 4 bytes and places them into a UInt32 in host-ordering
function get_event_type(event::Array{UInt8})
    extract_4byte_host(1, event) # return UInt32
end

function get_event_type(e::SDL2.Event)
    e._Event[1] # UInt8
end

get_key_code_sym(event::Array{UInt8}) = extract_4byte_host(21, event)
get_scancode(event::Array{UInt8}) = extract_4byte_host(17, event)
get_modifier(event::Array{UInt8}) = extract_2byte_host(25, event)
get_state(event::Array{UInt8}) = event[13]
get_repeat(event::Array{UInt8}) = event[14]

SDL_Event() = Array{UInt8}(zeros(56))
event = SDL_Event()

function poll_event!()
    success = (SDL2.PollEvent(event) ≠ 0)
    event, success
end

# --------------------------------------------------------------
# Key maps
# --------------------------------------------------------------
const KEY_A = SDL2.SDLK_a
const KEY_E = SDL2.SDLK_e
const KEY_F = SDL2.SDLK_f
const KEY_S = SDL2.SDLK_s
const KEY_W = SDL2.SDLK_w
const KEY_UP = SDL2.SDLK_UP
const KEY_DOWN = SDL2.SDLK_DOWN
const KEY_LEFT = SDL2.SDLK_LEFT
const KEY_RIGHT = SDL2.SDLK_RIGHT
const KEY_PRESSED = 1
const KEY_RELEASED = 0
const KEY_REPEATING = 1

# --------------------------------------------------------------
# Debugging printing and misc.
# --------------------------------------------------------------
to_hex(num) = "0x" * string(num, base = 16)

function print_key_event_type(event::Array{UInt8})
    e32 = event[1]
    e24 = event[2]
    e16 = event[3]
    e08 = event[4]
    b32 = bitstring(e32)
    b24 = bitstring(e24)
    b16 = bitstring(e16)
    b08 = bitstring(e08)
    t = get_event_type(event)
    println("type: $b08 $b16 $b24 $b32 <= ", t, " = ", to_hex(t));
end

function print_key_event_state(event::Array{UInt8})
    # SDL_PRESSED or SDL_RELEASED
    e13 = event[13]
    println("state: ($e13) ", to_hex(e13));
end

function print_key_event_repeat(event::Array{UInt8})
    # non-zero if this is a key repeat
    e13 = event[14]
    println("repeat: ($e13) ", to_hex(e13));
end

function print_key_event_sym(event::Array{UInt8})
    scancode = get_scancode(event)
    println("scancode: ($scancode) ", to_hex(scancode))
    keycode = get_key_code_sym(event)
    println("keycode: ($keycode) ", to_hex(keycode))
    mod = get_modifier(event)
    println("mod: ($mod) ", to_hex(mod));
end

# --------------------------------------------------------------
# Keyboard
# --------------------------------------------------------------
function print_keyboard_event(event::Array{UInt8})
    println("*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--")
    println("event: ", event[1:31])
    println("f-type:     ", event[1:4])   # Ex: SDL2.KEYDOWN
    println("f-state:    ", event[13:13])
    println("f-repeat:   ", event[14:14])
    println("f-scancode: ", event[17:20]) # Ex: SDL2.SCANCODE_C
    println("f-keycode:  ", event[21:24]) # Ex: SDL2.SDLK_c
    println("f-mod:      ", event[25:25]) # Ex: SDL2.KMOD_LCTRL

    print_key_event_type(event)
    print_key_event_state(event)
    print_key_event_repeat(event)
    print_key_event_sym(event)
    println("--------------------------------------------------------------")
end

mutable struct KeyboardEvent <: Ranger.AbstractIOEvent
    keycode::UInt32
    scancode::UInt32
    modifier::UInt16
    repeat::UInt8
    state::UInt8

    function KeyboardEvent()
        new(0, 0, 0, 0)
    end
end

function set!(keyboard::KeyboardEvent, event::Array{UInt8,1})
    keyboard.keycode = get_key_code_sym(event)
    keyboard.scancode = get_scancode(event)
    keyboard.modifier = get_modifier(event)
    keyboard.repeat = get_repeat(event)
    keyboard.state = get_state(event)
end

function print(keyboard::KeyboardEvent)
    println("key: ", keyboard.keycode,
    ", scan: ", keyboard.scancode,
    ", mod: ", keyboard.modifier,
    ", rep: ", keyboard.repeat,
    ", state: ", keyboard.state)
end

# --------------------------------------------------------------
# Mouse
# --------------------------------------------------------------
function print_mouse_event(event::Array{UInt8})
    # mutable struct MouseMotionEvent <: AbstractEvent
    #     _type::Uint32        1:4
    #     timestamp::Uint32    5:8
    #     windowID::Uint32     9:12
    #     which::Uint32        13:16
    #     state::Uint32
    #     x::Sint32
    #     y::Sint32
    #     xrel::Sint32
    #     yrel::Sint32
    # end
    println("*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--")
    println("event: ", event[1:36])

    type = extract_4byte_host(1, event)
    which = extract_4byte_host(13, event)
    state = extract_4byte_host(17, event)
    x = extract_4byte_host(21, event)
    y = extract_4byte_host(25, event)
    xrel = extract_4byte_host(29, event)
    yrel = extract_4byte_host(33, event)

    println("m-type:     ", event[1:4], " <= ", type)   # Ex: SDL2.MOUSEMOTION
    println("m-which:    ", event[13:16], " <= ", which)
    println("m-state:    ", event[17:20], " <= ", state)
    println("m-x:        ", event[21:24], " <= ", signed(x))
    println("m-y:        ", event[25:28], " <= ", signed(y))
    println("m-xrel:     ", event[29:32], " <= ", signed(xrel))
    println("m-yrel:     ", event[33:36], " <= ", signed(yrel))

    println("--------------------------------------------------------------")
end

mutable struct MouseMotionEvent <: Ranger.AbstractIOEvent
    type::UInt32
    which::UInt32
    state::UInt32
    x::Int32
    y::Int32
    xrel::Int32
    yrel::Int32

    function MouseMotionEvent()
        new(0, 0, 0, 0, 0, 0)
    end
end

function set!(mouse::MouseMotionEvent, event::Array{UInt8,1})
    mouse.type = extract_4byte_host(1, event)
    mouse.which = extract_4byte_host(13, event)
    mouse.state = extract_4byte_host(17, event)
    mouse.x = signed(extract_4byte_host(21, event))
    mouse.y = signed(extract_4byte_host(25, event))
    mouse.xrel = signed(extract_4byte_host(29, event))
    mouse.yrel = signed(extract_4byte_host(33, event))
end

mutable struct MouseWheelEvent <: Ranger.AbstractIOEvent
    type::UInt32
    which::UInt32
    x::Int32
    y::Int32
    direction::UInt32

    function MouseWheelEvent()
        new(0, 0, 0, 0, 0)
    end
end

function set!(mouse::MouseWheelEvent, event::Array{UInt8,1})
    mouse.type = extract_4byte_host(1, event)
    mouse.which = extract_4byte_host(13, event)
    mouse.x = signed(extract_4byte_host(17, event))
    mouse.y = signed(extract_4byte_host(21, event))
    mouse.direction = signed(extract_4byte_host(25, event))
end

mutable struct MouseButtonEvent <: Ranger.AbstractIOEvent
    type::UInt32
    which::UInt32
    button::UInt8
    state::UInt8
    clicks::UInt8
    # padding1 UInt8
    x::Int32
    y::Int32

    function MouseButtonEvent()
        new(0, 0, 0, 0, 0, 0, 0)
    end
end

function set!(mouse::MouseButtonEvent, event::Array{UInt8,1})
    mouse.type = extract_4byte_host(1, event)
    mouse.which = extract_4byte_host(13, event)
    mouse.button = event[17]

    mouse.state = event[18]
    mouse.clicks = event[19]
    # padding at 20
    mouse.x = signed(extract_4byte_host(21, event))
    mouse.y = signed(extract_4byte_host(25, event))
end

function is_left_mouse_down(mouse::MouseButtonEvent)
    mouse.button == 1 && mouse.state == 1
end

function is_left_mouse_up(mouse::MouseButtonEvent)
    mouse.button == 1 && mouse.state == 0
end

end # Module ------------------------------------------------------------

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using DataStructures

"Results structure"
struct UninformedSearchResult{TState}
  status::Symbol
  path::Vector{TState}
  closedsetsize::Int64
  opensetsize::Int64
end

struct UninformedSearchNode{TState}
  data::TState
  depth::Int64
  parent::Union{UninformedSearchNode{TState}, Nothing}
end

mutable struct DepthFirstSearchState{TState, THash}
  openset::Stack{UninformedSearchNode{TState}}
  closedset::Set{THash}
  start_time::Float64
end

function _depthfirst!(
  search_state::DepthFirstSearchState{TState, THash},
  neighbours,
  start,
  goal,
  isgoal,
  hashfn,
  timeout,
  maxdepth,
  enable_closedset,
) where {TState, THash}
  while !isempty(search_state.openset)
    node = pop!(search_state.openset)
    node_data = node.data
    if isgoal(node_data, goal)
      return UninformedSearchResult(
        :success,
        reconstructpath(node),
        length(search_state.closedset),
        length(search_state.openset),
      )
    end

    if timeout < Inf && time() - search_state.start_time >= timeout
      return UninformedSearchResult(
        :timeout,
        reconstructpath(node),
        length(search_state.closedset),
        length(search_state.openset),
      )
    end

    neighbours_data = collect(neighbours(node_data))
    # reverse so that the neighbours are checked in order
    reverse!(neighbours_data)
    for neighbour in neighbours_data
      new_depth = node.depth + 1
      if (enable_closedset && hashfn(neighbour) in search_state.closedset) ||
         new_depth > maxdepth
        continue
      end

      neighbour_node = UninformedSearchNode(neighbour, new_depth, node)
      push!(search_state.openset, neighbour_node)
    end

    if enable_closedset
      push!(search_state.closedset, hashfn(node_data))
    end
  end

  return UninformedSearchResult(
    :nopath,
    [start],
    length(search_state.closedset),
    length(search_state.openset),
  )
end

"""
Execute a depth-first search given the neighbours function, a start and a goal.
It takes the same parameters as astar, but it ignores those about cost and heuristic.
Use the `maxdepth` parameter to limit the search to a certain depth.
"""
function depthfirst(
  neighbours,
  start,
  goal;
  isgoal = defaultisgoal,
  hashfn = hash,
  timeout = Inf,
  maxdepth = typemax(Int64),
  enable_closedset = true,
  kwargs...,
)
  start_time = time()
  start_hash = hashfn(start)

  start_node = UninformedSearchNode(start, 0, nothing)
  stack = Stack{typeof(start_node)}()
  push!(stack, start_node)
  search_state = DepthFirstSearchState(stack, Set{typeof(start_hash)}(), start_time)

  return _depthfirst!(
    search_state,
    neighbours,
    start,
    goal,
    isgoal,
    hashfn,
    timeout,
    maxdepth,
    enable_closedset,
  )
end

"""
Execute iteratively a depth first search starting from depth 0 to maxdepth.
"""
function iterative_deepening(
  neighbours,
  start,
  goal;
  isgoal = defaultisgoal,
  hashfn = hash,
  timeout = Inf,
  maxdepth = typemax(Int64),
  enable_closedset = true,
  kwargs...,
)
  start_time = time()
  end_time = start_time + timeout
  closedsetsize = 0
  opensetsize = 0

  for depth = 0:maxdepth
    depth_first_timeout = max(end_time - time(), 0.0)
    res = depthfirst(
      neighbours,
      start,
      goal;
      isgoal,
      hashfn,
      enable_closedset,
      maxdepth = depth,
      timeout = depth_first_timeout,
      kwargs...,
    )
    if res.status in (:success, :timeout)
      return res
    end
    closedsetsize = res.closedsetsize
    opensetsize = res.opensetsize
  end
  return UninformedSearchResult(:nopath, [start], closedsetsize, opensetsize)
end

mutable struct BreadthFirstSearchState{TState, THash}
  openset::Deque{UninformedSearchNode{TState}}
  closedset::Set{THash}
  start_time::Float64
end

function _breadthfirst!(
  search_state::BreadthFirstSearchState{TState, THash},
  neighbours,
  start,
  goal,
  isgoal,
  hashfn,
  timeout,
  maxdepth,
  enable_closedset,
) where {TState, THash}
  while !isempty(search_state.openset)
    node = popfirst!(search_state.openset)
    node_data = node.data
    if isgoal(node_data, goal)
      return UninformedSearchResult(
        :success,
        reconstructpath(node),
        length(search_state.closedset),
        length(search_state.openset),
      )
    end

    if timeout < Inf && time() - search_state.start_time >= timeout
      return UninformedSearchResult(
        :timeout,
        reconstructpath(node),
        length(search_state.closedset),
        length(search_state.openset),
      )
    end

    for neighbour in neighbours(node_data)
      new_depth = node.depth + 1
      if (enable_closedset && hashfn(neighbour) in search_state.closedset) ||
         new_depth > maxdepth
        continue
      end

      neighbour_node = UninformedSearchNode(neighbour, new_depth, node)
      push!(search_state.openset, neighbour_node)
    end

    if enable_closedset
      push!(search_state.closedset, hashfn(node_data))
    end
  end

  return UninformedSearchResult(
    :nopath,
    [start],
    length(search_state.closedset),
    length(search_state.openset),
  )
end

"""
Executes a breadth first search given the neightbours function, the start and goal state.
"""
function breadthfirst(
  neighbours,
  start,
  goal;
  isgoal = defaultisgoal,
  hashfn = hash,
  timeout = Inf,
  maxdepth = typemax(Int64),
  enable_closedset = true,
  kwargs...,
)
  start_time = time()
  start_hash = hashfn(start)

  start_node = UninformedSearchNode(start, 0, nothing)
  deque = Deque{typeof(start_node)}()
  push!(deque, start_node)
  search_state = BreadthFirstSearchState(deque, Set{typeof(start_hash)}(), start_time)

  return _breadthfirst!(
    search_state,
    neighbours,
    start,
    goal,
    isgoal,
    hashfn,
    timeout,
    maxdepth,
    enable_closedset,
  )
end


text:
export AbstractSize, time_space, time_space_src, time_src, rec_space, space, space_src

struct AbstractSize
    dims::Dict
end

AbstractSize(n::NTuple{N, Symbol}, sizes) where N = AbstractSize(Dict(zip(n, sizes)))

# Base
function getindex(S::AbstractSize, sinds::Union{AbstractRange, Vector{<:Integer}})
    k = keys(S.dims)
    dim_map = copy(S.dims)
    for (k, v) in S.dims
        new_v = typeof(v) <: Integer ? v : v[sinds]
        dim_map[k] = new_v
    end
    :src ∈ keys(dim_map) && (dim_map[:src] = length(sinds))
    AbstractSize(dim_map)
end
getindex(S::AbstractSize, i::Integer) = getindex(S, i:i)

getindex(S::AbstractSize, I...) = AbstractSize(keys(S.dims[I...]), values(S.dims[I...]))
getindex(S::AbstractSize, I::Symbol) = S.dims[I]
setindex!(S::AbstractSize, v, I::Symbol) = setindex!(S.dims, v, I)
iterate(S::AbstractSize) = iterate(S.dims)
iterate(S::AbstractSize, state) = iterate(S.dims, state)

convert(::Type{T}, S::AbstractSize) where T<:Number = convert(T, nsamples(S))
(::Type{T})(S::AbstractSize) where T<:Number = convert(T, nsamples(S))
(c::Colon)(i::T, S::AbstractSize) where T = c(i, T(S))

Base.keys(S::AbstractSize) = keys(S.dims)
Base.values(S::AbstractSize) = values(S.dims)
Base.merge!(S1::AbstractSize, S2::AbstractSize) = merge!(S1.dims, S2.dims)

==(S1::AbstractSize, S2::AbstractSize) = S1.dims == S2.dims
==(S1::Integer, S2::AbstractSize) = nsamples(S2) == S1
==(S1::AbstractSize, S2::Integer) = nsamples(S1) == S2

Base.repr(S::AbstractSize) = "($(join(keys(S.dims), " * ")))"

similar(x::judiMultiSourceVector, ::Type{ET}, dims::AbstractSize) where ET = similar(x, ET)
similar(x::AbstractVector, ::Type{ET}, dims::AbstractSize) where ET = similar(x, ET, Int(dims))

# Update size
set_space_size!(S::AbstractSize, vals::NTuple{2, Integer}) = begin pop!(S.dims, :y, 0); merge!(S.dims, Dict(zip((:x, :z), vals))) end
set_space_size!(S::AbstractSize, vals::NTuple{3, Integer}) = merge!(S.dims, Dict(zip((:x, :y, :z), vals)))

# Actual integer size
nsamples(dims::Dict{Symbol, Any}) = sum(.*(values(dims)...)) ÷ get(dims, :src, 1)
nsamples(dims::Dict{N, <:Integer}) where N = prod(values(dims))
nsamples(D::AbstractSize) = nsamples(D.dims)

# Constructors
space(N::NTuple{2, Integer}) = AbstractSize((:x, :z), (N[1], N[2]))
space(N::NTuple{3, Integer}) = AbstractSize((:x, :y, :z), (N[1], N[2], N[3]))

time_space(N::NTuple{2, Integer}) = AbstractSize((:time, :x, :z), ([0], N[1], N[2]))
time_space(N::NTuple{3, Integer}) = AbstractSize((:time, :x, :y, :z), ([0], N[1], N[2], N[3]))

time_space_src(nsrc::Integer, nt, N::Integer) = AbstractSize((:src, :time, :x, :y, :z), (nsrc, nt, zeros(Int, N)...))
time_space_src(nsrc::Integer, nt, N::NTuple{2, Integer}) = AbstractSize((:src, :time, :x, :z), (nsrc, nt, N...))
time_space_src(nsrc::Integer, nt, N::NTuple{3, Integer}) = AbstractSize((:src, :time, :x, :y, :z), (nsrc, nt, N...))
time_space_src(nsrc::Integer, nt) = AbstractSize((:src, :time, :x, :y, :z), (nsrc, nt, 0, 0, 0))

space_src(nsrc::Integer) = AbstractSize((:src, :x, :y, :z), (nsrc, 0, 0, 0))

time_src(nsrc::Integer, nt) = AbstractSize((:src, :time), (nsrc, nt))

rec_space(G::Geometry) = AbstractSize((:src, :time, :rec), (get_nsrc(G), G.nt, G.nrec))


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    InpaintAlgo

An inpainting algorithm with given parameters.
"""
abstract type InpaintAlgo end

"""
    inpaint(img, mask, algo)

Inpaint `img` on pixels marked as `true` in `mask` using
algorithm `algo`.
"""
function inpaint(img::AbstractArray, mask::BitArray, algo::InpaintAlgo)
  # sanity checks
  @assert size(img) == size(mask) "image and mask must have same size"

  # dispatch appropriate implementation
  inpaint_impl(img, mask, algo)
end

#------------------
# IMPLEMENTATIONS
#------------------
include("crimisini.jl")


text:
struct FDJONSWAP <: FDSpectralModel
    θ::Vector{Float64}
end

parameter(J::FDJONSWAP) = J.θ

function sdf(J::Union{FDJONSWAP,FDSwellJS}, ω)
    α = J.θ[1]
    ωₚ = J.θ[2]
    γ = J.θ[3]
    r = J.θ[4]
    ω = abs(ω)
    s = 4.0
    σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
    δ = exp(-1 / (2 * σ1²) * (ω / ωₚ - 1)^2)
    sdf = α * ω^(-r) * exp(-(r / s) * (ω / ωₚ)^(-s)) * γ^δ / 2
    sdf *= (!isnan(sdf))
    return sdf
end

function spreading(J::FDJONSWAP, ω, ϕ)
    ω = abs(ω)
    halfPeakSep = J.θ[6] / exp(J.θ[7] * ((ω > J.θ[2])*((J.θ[2] / ω) - 1.0) + 1.0)) / 2
    φₘ₁ = J.θ[5] + halfPeakSep
    φₘ₂ = J.θ[5] - halfPeakSep
    σ = J.θ[8] - (J.θ[9]/3)*(4*(J.θ[2]/ω)^2 - (J.θ[2]/ω)^8)
    return 0.5 * (WrappedGaussian(ϕ, φₘ₁, σ) + WrappedGaussian(ϕ, φₘ₂, σ))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f38a5752-3e3e-11eb-256c-eb1f8de56937
begin
	import Pda
	using Plots
	using PlutoUI
	using HDF5
	using Statistics
end

# ╔═╡ 09068b50-3e3f-11eb-330c-ab8cf9ce1600
file = "/media/HDD1/singlewall/U2.154singlewall.h5"

# ╔═╡ 14f7b2da-3e40-11eb-1aca-593125a5c40c
frame =99

# ╔═╡ a7c0b072-3e3f-11eb-1498-4343ce3bfccc
begin
	t = Pda.get_h5data(file, "config/$(frame)/t")
	x = Pda.get_h5data(file, "config/$(frame)/x")
end

# ╔═╡ be7aa564-3e41-11eb-343a-6f8d2a3367e9
function average_density()
	bins=100
	n = zeros(bins)
	loc = zeros(bins)
	for i in 0:99
		loc, n_ = Pda.density1d(Pda.get_h5data(file, "config/$(i)/x"),
			bins=bins,
			normalize=false)
		n += n_
	end
	return loc, n
end

# ╔═╡ 2bed370c-3e41-11eb-1fb8-c7f187b21521
loc,n = average_density()

# ╔═╡ 7e1e3a7a-3e42-11eb-06c8-679cebbbc9ec
function density_Qclosure(x, Pel::Real)
	laml = sqrt(Pel+0.5*Pel^2)
	n = 1 + 0.5*Pel*exp(-laml*x)
	return n
end

# ╔═╡ abfa5be6-3e41-11eb-0c1a-4d5d565887d3
begin
	scatter(loc[2:end-1],n[2:end-1]/(n[end-10:end-1] |>mean),
		framestyle=:box,
		label="BD")
	U0 = Pda.get_h5attr(file, "U0")
	tauR = Pda.get_h5attr(file, "tauR")
	DT = Pda.get_h5attr(file,"DT")
	Pel = (U0*tauR/sqrt(DT*tauR))^2
	plot!(loc, density_Qclosure.(loc, Pel),
		lw=2,
		label="Q closure")
end

# ╔═╡ Cell order:
# ╠═f38a5752-3e3e-11eb-256c-eb1f8de56937
# ╠═09068b50-3e3f-11eb-330c-ab8cf9ce1600
# ╠═14f7b2da-3e40-11eb-1aca-593125a5c40c
# ╠═a7c0b072-3e3f-11eb-1498-4343ce3bfccc
# ╠═be7aa564-3e41-11eb-343a-6f8d2a3367e9
# ╠═2bed370c-3e41-11eb-1fb8-c7f187b21521
# ╠═7e1e3a7a-3e42-11eb-06c8-679cebbbc9ec
# ╠═abfa5be6-3e41-11eb-0c1a-4d5d565887d3


text:
export genAmpl2dDataFEM
function
genAmpl2dDataFEM(M::RegularMesh,SM::SparseMatrixCSC{Float64},MC::SparseMatrixCSC{Float64},b::Vector,P::SparseMatrixCSC,dobs::Vector,nbx::Vector,nby::Vector,pNoise::Int,ns::Int)

#-------------------------------------------------------------------------------
# Generating lookup table for variables u and w
#-------------------------------------------------------------------------------
  nn1  = M.n.+1
  nn2 = nn1-1
  np1 = prod(nn1)
  np2 = prod(nn2)
  wx  = zeros(Int64,np2,1)
  wy  = zeros(Int64,np2,1)
  ux  = zeros(Int64,np1,1)
  uy  = zeros(Int64,np1,1)
  #bc=(:dir,:neu,:neu,:neu,:neu,:neu)
  #iddir = getBoundaryIndicesFEM(M, bc)  # Dirichlet nodes
  #(dirX, dirY) = ind2sub(zeros(nn1[1],nn1[2]),iddir)
  for ii=1:nn2[2]
  	for jj=1:nn2[1]
  	wx[(ii-1)*nn2[1]+jj] = jj;
  	wy[(ii-1)*nn2[1]+jj] = ii;
  	nodes =	[jj+(nn2[1]+1)*(ii-1);jj+(nn2[1]+1)*(ii-1)+1;jj+(nn2[1]+1)*(ii)+1;jj+(nn2[1]+1)*ii];
  	ux[nodes[1]] = jj;  uy[nodes[1]] = ii;
  	ux[nodes[2]] = jj+1; uy[nodes[2]] = ii;
  	ux[nodes[3]] = jj+1; uy[nodes[3]] = ii+1;
  	ux[nodes[4]] = jj;   uy[nodes[4]] = ii+1;
    end
  end

#-------------------------------------------------------------------------------
# Generating data file for AMPL
#-------------------------------------------------------------------------------

	printAMPL = true
	fname = string(M.n[1], "x",	M.n[2],"S",ns,"P",length(dobs),"N",pNoise,".dat");
	fname1 = string(M.n[1], "x",M.n[2],"S",ns,"P",length(dobs),"N",pNoise);
	if printAMPL
		open(fname,"w") do f
			println(f,"data;")
			println(f, "param outfile:=",fname1, ";")
			println(f,"param Ly  := $(M.domain[2]-M.domain[1]) ;")
			println(f,"param Lx  := $(M.domain[4]-M.domain[3]) ;")
			println(f,"param Ny  := $((M.n[1])+1) ;")
			println(f,"param Nx  := $((M.n[2])+1) ;")
			println(f,"param NBH := $(length(dobs)) ;")
			@printf(f,"param nbx := ")
			for k=1:length(dobs)
				@printf(f, " %d %.8f",k, nbx[k])
			end
			@printf(f, ";\n")
			@printf(f,"param nby := ")
			for k=1:length(dobs)
				@printf(f, " %d %.8f",k, nby[k])
			end
			@printf(f, ";\n")
			for k=1:length(dobs)
				irow = find(P[k,:])
				for i=1:length(irow)
					@printf(f, "let Uhat[%d,%d]:=%.8f;\n",ux[irow[i]],uy[irow[i]],dobs[k])
				end
			end
			#for k=1:length(dirX)
			  #@printf(f, "fix U[%d,%d]:=0;\n",dirX[k], dirY[k])
			#end
			for k=1:length(dobs)
				irow = find(P[k,:])
				for i=1:length(irow)
					@printf(f, "let P[%d,%d]:=%d;\n",ux[irow[i]],uy[irow[i]],1)
				end
			end

			println(f, "model;")
			println(f, "subject to")
			for k=1:nn1[1]*nn1[2]
			  #if mod(k,nn1[1]) != 1
			    irows = find(SM[k,:])
			    print(f, "FEM_",k,": ")
			    for i=1:length(irows)
				  #@printf(f,"+ U[%d,%d]*SM[%d,%d]",ux[irows[i]],uy[irows[i]],k,irows[i])
				  if i==1
					@printf(f,"%.8f* U[%d,%d]", SM[k,irows[i]],ux[irows[i]],uy[irows[i]])

				  else
					@printf(f,"+ %.8f* U[%d,%d]", SM[k,irows[i]],ux[irows[i]],uy[irows[i]])

				  end
			    end
			    print(f,"=")
			    irows = find(MC[k,:])
			    for i=1:length(irows)
				  #@printf(f,"+ W[%d,%d]*MC[%d,%d]",wx[irows[i]],wy[irows[i]],k,irows[i])
				  if i==1
					@printf(f,"%.8f* W[%d,%d] + %.8f", MC[k,irows[i]],wx[irows[i]],wy[irows[i]], b[k])
				  else
					@printf(f,"+ %.8f* W[%d,%d] + %.8f", MC[k,irows[i]],wx[irows[i]],wy[irows[i]], b[k])
				  end
			    end
				if length(irows) == 0
			      @printf(f,"0")
				end
			    println(f,";")
			  #end
			end
		end
	end
	return
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module ShaleDrillingLikelihood_OptimizeStatic

DOBTIME = false
DOPROFILE = false

using ShaleDrillingLikelihood
using Test
using SparseArrays
using BenchmarkTools
using Calculus
using Random
using InteractiveUtils
using Distributions
using Dates
using StatsBase
using Optim
using CountPlus
using LinearAlgebra
using Profile
using Distributed
using CountPlus
# using PProf
using Juno
# if DOPROFILE
#     using ProfileView
# end

using Optim: minimizer, minimum, maximize, maximizer, maximum

using Base.Iterators: product, OneTo

using ShaleDrillingLikelihood: ObservationDrill,
    SimulationDraw,
    idx_ρ,
    actionspace,
    ichars_sample,
    j1chars,
    tptr,
    _x, _y,
    zchars,
    ztrange,
    InitialDrilling,
    _i,
    j1_sample,
    j1_range,
    jtstart,
    exploratory_terminal,
    ssprime,
    total_leases,
    zero_out_small_probs,
    lrvar,
    lrstd,
    lrmean,
    simulate,
    theta_royalty_ρ,
    revenue,
    cost,
    extend,
    theta_ρ,
    vw_revenue,
    vw_cost,
    vw_extend,
    DataDynamicDrill,
    max_states,
    _D,
    simloglik!,
    simloglik_drill_data!,
    SimulationDraws,
    DCDPTmpVars,
    _model,
    theta_drill_ρ,
    update!,
    value_function,
    _x, _y,
    end_ex0, end_ex1, end_lrn, end_inf,
    state_idx,
    j2ptr, tptr,
    idx_produce_ψ, idx_drill_ψ,
    idx_produce_g, idx_drill_g,
    merge_thetas,
    split_thetas,
    theta_produce,
    serial_simloglik!,
    parallel_simloglik!,
    TestDrillModel,
    drill,
    DataFull,
    coeftable

println("print to keep from blowing up")

@testset "simulate all data" begin

    Random.seed!(7)

    discount = ((0x1.006b55c832502p+0)^12 / 1.125) ^ (1/4)  # real discount rate

    # set up coefs
    θρ =  0.0
    αψ = 0.33
    αg = 0.56

    model_drill = TestDrillModel()
    θ_drill   = [αψ, 2.0, -2.0, -0.75, θρ]
    θ_royalty = [θρ, 1.0,    0.5, 0.5,    -0.6, 0.6]  # dψdρ, ψ, β, κ
    θ_produce = vcat(αψ, αg, 0.2, 0.3, 0.4)

    θ = vcat(θ_drill, θ_royalty[2:end], θ_produce[2:end])
    println(sprintf_binary(θ))

    M = 250
    num_i = 300
    num_zt = 200          # drilling

    # observations
    obs_per_well = 10:20  # pdxn

    # roylaty
    royalty_rates = [1/8, 3/16, 1/4,]
    num_royalty_rates = length(royalty_rates)
    nk_royalty = length(θ_royalty) - (length(royalty_rates)-1)-2

    # geology
    ogip_dist = Normal(4.68,0.31)

    # price parameters
    zrho = 0.8                # rho
    zmean = 1.33*(1-zrho)     # mu
    zvar = 0.265^2*(1-zrho^2) # sigsq

    # simulate price process
    zprocess = AR1process(zmean, zrho, zvar)
    zvec = simulate(zprocess, num_zt)
    ztrng = range(Date(2003,10); step=Dates.Quarter(1), length=num_zt)
    _zchars = ExogTimeVars(tuple.(zvec), ztrng)

    # random vars
    log_ogip = rand(ogip_dist, num_i)
    Xroyalty = vcat(log_ogip', randn(nk_royalty-1, num_i))
    u,v = randn(num_i), randn(num_i)

    # construct royalty data
    data_roy = DataRoyalty(u, v, Xroyalty, θ_royalty, num_royalty_rates)
    royrates_endog = royalty_rates[_y(data_roy)]
    royrates_exog = sample(royalty_rates, num_i)

    # construct ichars for drilling
    _ichars = [gr for gr in zip(log_ogip, royrates_exog)]
    @test _ichars isa Vector{<:NTuple{2,Number}}

    data_drill_opt = (minmaxleases=1:5, nper_initial=40:40, nper_development=40:40, tstart=1:50,)
    data_drill   = DataDrill(u, v, _zchars, _ichars, TestDrillModel(), θ_drill; data_drill_opt...)
    data_produce = DataProduce(u, 10, obs_per_well, θ_produce, log_ogip)

    coefnames(data_produce)
    coefnames(data_roy)
    coefnames(data_drill)

    coef_links = [(idx_produce_ψ, idx_drill_ψ,),]
    data_full = DataSetofSets(data_drill, data_roy, data_produce, coef_links)
    data_small = DataSetofSets(EmptyDataSet(), EmptyDataSet(), data_produce)
    data_drill = DataSetofSets(data_drill, EmptyDataSet(), EmptyDataSet())

    coefnames(data_full)
    coefnames(data_small)
    coefnames(data_drill)

    @test data_full isa DataFull

    DOPAR = true

    if DOPAR
        rmprocs(workers())
        pids = addprocs()
    else
        pids = [1,]
    end
    println_time_flush("Putting pkg on workers")
    @everywhere using ShaleDrillingLikelihood
    println_time_flush("Package on workers")

    datas  = [ data_small,   data_drill,   data_full,]
    datanm = ["data_small", "data_drill", "data_full",]
    thetas = [θ_produce, θ_drill, θ, ]

    theta_tuple = (θ_drill, θ_royalty, θ_produce)
    theta_full = merge_thetas(theta_tuple, data_full)
    @test theta_full == θ
    @test all(theta_tuple .== split_thetas(data_full, theta_full))

    let leo = LocalEstObj(data_small, θ_produce)
        coeftable(leo)
        coeftable(leo, 0.1)
    end

    for (d,t,nm) in zip(datas, 0.9 .* thetas, datanm)

        @test _nparm(d) .== length(t)
        println("rocking out dataset $nm")
        resetcount!()

        leo = LocalEstObj(d,t)
        s = SimulationDraws(M, ShaleDrillingLikelihood.data(leo))
        reo = RemoteEstObj(leo, M)
        ew = EstimationWrapper(leo, reo)

        leograd = ShaleDrillingLikelihood.grad(leo)

        stopcount!()
        for dograd in false:true
            @eval @everywhere set_g_RemoteEstObj($reo)

            simloglik!(1, t, dograd, reo)
            ShaleDrillingLikelihood.reset!(leo)
            sll = serial_simloglik!(ew, t, dograd)
            sgrad = copy(leograd)

            pll = parallel_simloglik!(ew, t, dograd)
            pgrad = copy(leograd)

            @test sll == pll
            @test sgrad == pgrad
        end

        fill!(leograd, 0)
        serial_simloglik!(ew, t, true)
        gs = copy(leograd)
        @test any(gs .!= 0)

        fill!(leograd, 0)
        parallel_simloglik!(ew, t, true)
        gp = copy(leograd)
        @test any(gp .!= 0)
        @test gs == gp

        fds = Calculus.gradient(x -> serial_simloglik!(  ew, x, false), t, :central)
        fdp = Calculus.gradient(x -> parallel_simloglik!(ew, x, false), t, :central)
        @test fds == fdp
        @test gs ≈ fds
        @test gp ≈ fdp

        if true
            # startcount!([1, 100000,], [100, 100,])
            resetcount!()
            startcount!([100, 500, 100000,], [10, 10, 100,])
            opts = Optim.Options(show_trace=true, time_limit=20, allow_f_increases=true)
            res = solve_model(ew, 0.9*t; OptimOpts=opts)
            @show res
            @test minimizer(res) == theta1(leo)
            println(sprintf_binary(minimizer(res)))

            if drill(d) != EmptyDataSet()
                @test theta_ρ(d, t) == θρ
            end
            println(coeftable(leo))
        end
    end
end

end # module


text:
"""
A *finite polynomial* is a polynomial endofunctor on FinSet or equivalently a
polynomial functor whose sets of positions and directions are finite. These
objects are simple but perhaps not very useful on their own. For a richer class
of polynomials, see the module [`SumProdPolynomials`](@ref).
"""
module FinPolynomials
export AbstractFinPoly, FinPoly, FinSet, FinFunction, positions, directions,
  proj, plus, times, otimes, ⊗

using Catlab, Catlab.CategoricalAlgebra, Catlab.CategoricalAlgebra.FinSets
using Catlab.Theories
import Catlab.Theories: proj, plus, otimes, ⊗

@present FinPolySchema(FreeSchema) begin
  (Pos, Dir)::Ob
  pos::Hom(Dir, Pos)
end

""" Abstract type for finite polynomial. See [`FinPolynomial`](@ref).
"""
const AbstractFinPoly = AbstractACSetType(FinPolySchema)

""" Finite polynomial, i.e., polynomial endofunctor on FinSet.

Finite polynomials can be constructed from a list of exponents or from the
position-direction projection. For example, both `FinPoly([3,2,0])` and
`FinPoly(FinFunction([1,1,1,2,2],3))` yield the polynomial ``y³+y²+1``.
"""
const FinPoly = ACSetType(FinPolySchema, index=[:pos])

function (::Type{P})(ndirs::AbstractVector{Int}) where P <: AbstractFinPoly
  p = P()
  add_parts!(p, :Pos, length(ndirs))
  for (pos, ndir) in enumerate(ndirs)
    add_parts!(p, :Dir, ndir, pos=pos)
  end
  p
end

function (::Type{P})(f::FinFunction) where P <: AbstractFinPoly
  p = P()
  add_parts!(p, :Pos, length(codom(f)))
  add_parts!(p, :Dir, length(dom(f)), pos=collect(f))
  p
end

""" Positions of a polynomial.
"""
positions(p::AbstractFinPoly) = FinSet(p, :Pos)

""" All directions of a polynomial, or directions at a given position.
"""
directions(p::AbstractFinPoly) = FinSet(p, :Dir)
directions(p::AbstractFinPoly, pos) = incident(p, :pos, pos)

""" Projection from directions to positions in a polynomial.
"""
proj(p::AbstractFinPoly) = FinFunction(p, :pos)

plus(p::AbstractFinPoly, q::AbstractFinPoly) = ob(coproduct(p, q))
Base.:+(p::AbstractFinPoly, q::AbstractFinPoly) = plus(p, q)

times(p::P, q::P) where P <: AbstractFinPoly = P(fiber_sum(proj(p), proj(q)))
Base.:*(p::AbstractFinPoly, q::AbstractFinPoly) = times(p, q)

otimes(p::AbstractFinPoly, q::AbstractFinPoly) = ob(product(p, q))
⊗(p, q) = otimes(p, q)

# Set constructions
###################

""" Fiber sum of indexed sets.

Given indexed sets ``Xᵢ``, ``i ∈ I``, and ``Yⱼ``, ``j ∈ J``, represented by
functions ``f: X → I`` and ``g: Y → J``, compute the "fiber sum" ``Xᵢ+Yⱼ``,
``(i,j) ∈ I×J``, as a function ``X×J + I×Y → I×J``.
"""
function fiber_sum(f, g)
  X, Y, I, J = dom(f), dom(g), codom(f), codom(g)
  XJ, IY, IJ = product(X,J), product(I,Y), product(I,J)
  fJ = pair(IJ, proj1(XJ)⋅f, proj2(XJ)) # f×id(J) : X×J → I×J
  Ig = pair(IJ, proj1(IY), proj2(IY)⋅g) # id(I)×g : I×Y → I×J
  copair(coproduct(ob(XJ), ob(IY)), fJ, Ig)
end

end


text:
using Dates

mutable struct ForwardRateAgreement{DC <: DayCount, BC <: BusinessCalendar, C <: BusinessDayConvention, Y <: YieldTermStructure, P <: PositionType, TP <: TenorPeriod, CUR <: AbstractCurrency, IB <: BusinessCalendar, IC <: BusinessDayConvention, IDC <: DayCount, IT <: TermStructure} <: AbstractForward
  lazyMixin::LazyMixin
  underlyingIncome::Float64
  underlyingSpotValue::Float64
  dc::DC
  calendar::BC
  convention::C
  settlementDays::Int
  payoff::ForwardTypePayoff
  valueDate::Date
  maturityDate::Date
  discountCurve::Y
  fraType::P
  forwardRate::InterestRate{IDC, SimpleCompounding, Once}
  strikeForwardRate::InterestRate{IDC, SimpleCompounding, Once}
  notionalAmount::Float64
  iborIndex::IborIndex{TP, CUR, IB, IC, IDC, IT}
end

function ForwardRateAgreement(valueDate::Date,
                              maturityDate::Date,
                              position::P,
                              strikeForward::Float64,
                              notionalAmount::Float64,
                              iborIndex::IborIndex{TP, CUR, IB, IC, IDC, IT},
                              discountCurve::Y) where {P <: PositionType, TP <: TenorPeriod, CUR <: AbstractCurrency, IB <: BusinessCalendar, IC <: BusinessDayConvention, IDC <: DayCount, IT <: TermStructure, Y <: YieldTermStructure}
  calendar = iborIndex.fixingCalendar
  convention = iborIndex.convention
  settlementDays = iborIndex.fixingDays
  maturityDate = adjust(calendar, convention, maturityDate)
  fixingDate = advance(Dates.Day(-settlementDays), calendar, valueDate)
  forwardRate = InterestRate(fixing(iborIndex, iborIndex.ts, fixingDate), iborIndex.dc, SimpleCompounding(), Once())
  strikeForwardRate = InterestRate(strikeForward, iborIndex.dc, SimpleCompounding(), Once())

  strike = notionalAmount * compound_factor(strikeForwardRate, valueDate, maturityDate)
  payoff = ForwardTypePayoff(position, strike)

  return ForwardRateAgreement{IDC, IB, IC, Y, P, TP, CUR, IB, IC, IDC, IT}(LazyMixin(), 0.0, 0.0, iborIndex.dc, calendar, convention, settlementDays, payoff, valueDate,
                              maturityDate, discountCurve, position, forwardRate, strikeForwardRate, notionalAmount, iborIndex)
end

function spot_value(fra::ForwardRateAgreement)
  calculate!(fra)
  return fra.notionalAmount * compound_factor(forward_rate(fra), fra.valueDate, fra.maturityDate) * discount(fra.discountCurve, fra.maturityDate)
end

function forward_rate(fra::ForwardRateAgreement)
  calculate!(fra)
  return fra.forwardRate
end

function forward_value(fra::ForwardRateAgreement)
  calculate!(fra)
  return (fra.underlyingSpotValue - fra.underlyingIncome) / discount(fra.discountCurve, fra.maturityDate)
end

spot_income(::ForwardRateAgreement, ::YieldTermStructure) = 0.0

function implied_yield(fra::ForwardRateAgreement, underlyingSpotValue::Float64, forwardValue::Float64, settlementDate::Date, comp::CompoundingType, dc::DayCount)
  t = year_fraction(dc, settlementDate, fra.maturityDate)
  compoundingFactor = forwardValue / (underlyingSpotValue - spot_income(fra, fra.discountCurve))

  return implied_rate(compoundingFactor, dc, comp, t, Annual())
end

function perform_calculations!(fra::ForwardRateAgreement)
  fixingDate = advance(Dates.Day(-fra.settlementDays), fra.calendar, fra.valueDate)
  fra.forwardRate = InterestRate(fixing(fra.iborIndex, fra.iborIndex.ts, fixingDate), fra.iborIndex.dc, SimpleCompounding(), Once())
  fra.underlyingSpotValue = spot_value(fra)
  fra.underlyingIncome = 0.0
  # NPV calc

  return fra
end


text:
# This does simple Poisson estimation using the meps1996.data

# The MEPS data

# Define dep and expl vbls
#
# 	The dep. vbls, in corresponding column
# 	Office based doctor visits	1
# 	Outpatient doctor visits	2
# 	Emergency room visits		3
# 	Inpatient visits		4
# 	Dental visits			5
# 	Prescriptions			6
#
using Econometrics, DelimitedFiles, Statistics
which_dep = 1
if (which_dep == 1) dep = "OBDV"   end
if (which_dep == 2) dep = "OPV"    end
if (which_dep == 3) dep = "IPV"    end
if (which_dep == 4) dep = "ERV"    end
if (which_dep == 5) dep = "DV"     end
if (which_dep == 6) dep = "PRESCR" end
cd(@__DIR__)
data = readdlm("../Data/meps1996.data")
y = data[:,which_dep]
x = data[:,7:12]
n = size(x,1)
x = [ones(n,1) x]
# to check effects of poor scaling, try commenting the next line
x[:,end] = x[:,end]/1000.0 # put income in thousands
names = ["constant", "pub. ins.","priv. ins.", "sex", "age","edu","inc"]
title = "Poisson model, "*dep* ",  MEPS 1996 full data set"
model = theta -> poisson(theta, y, x)
θstart = zeros(size(x,2)) # start values for estimation
# try adding the option vc=2 for "Hessian" or vc=3 for OPG
θhat, objvalue, V, converged = mleresults(model, θstart, title, names, vc=1)

λhat = exp.(x*θhat)
mv = round(mean(λhat),digits=2)
println("Estimated marginal variance: $mv")
sv = round(var(y),digits=2)
println("Sample marginal variance: $sv")

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
using Test: get_test_counts
using CancellationTokens
using Test

@testset "CancellationTokens" begin

    src = CancellationTokenSource()
    cancel(src)
    wait(get_token(src))

    src = CancellationTokenSource()
    @async begin
        sleep(0.1)
        cancel(src)
    end
    wait(get_token(src))

    src = CancellationTokenSource(0.1)
    wait(get_token(src))

    src = CancellationTokenSource()
    @async begin
        sleep(0.1)
        cancel(src)
    end
    wait(get_token(src))

    src = CancellationTokenSource()
    sleep(0.1, get_token(src))

    src = CancellationTokenSource()
    @async begin
        sleep(0.1)
        cancel(src)
    end
    @test_throws OperationCanceledException sleep(20.0, get_token(src))
end


text:
"""
    solve_objbound_edge_cases(model::MOI.ModelLike, config::TestConfig)

Test a variety of edge cases related to the ObjectiveBound attribute.
"""
function solve_objbound_edge_cases(model::MOI.ModelLike, config::TestConfig)
    @testset "Min IP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(model, """
            variables: x
            minobjective: 2.0x + -1.0
            c1: x >= 1.5
            c2: x in Integer()
        """)
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(model, config;
            objective_value = 3.0,
            variable_primal = [(x, 2.0)]
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0
        end
    end

    @testset "Max IP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(model, """
            variables: x
            maxobjective: 2.0x + 1.0
            c1: x <= 1.5
            c2: x in Integer()
        """)
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(model, config;
            objective_value = 3.0,
            variable_primal = [(x, 1.0)]
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) >= 3.0
        end
    end

    @testset "Min LP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(model, """
            variables: x
            minobjective: 2.0x + -1.0
            c1: x >= 1.5
        """)
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(model, config;
            objective_value = 2.0,
            variable_primal = [(x, 1.5)]
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0
        end
    end

    @testset "Max LP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(model, """
            variables: x
            maxobjective: 2.0x + 1.0
            c1: x <= 1.5
        """)
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(model, config;
            objective_value = 4.0,
            variable_primal = [(x, 1.5)]
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) >= 4.0
        end
    end
end
unittests["solve_objbound_edge_cases"] = solve_objbound_edge_cases


text:
struct InvalidTopology <: Exception
    msg::String
end

struct SpeciesNameCollision <: Exception
    sppnames::Vector{String}
end

SpeciesNameCollision() = SpeciesNameCollision(String[])

function Base.showerror(io::IO, ex::SpeciesNameCollision)
    if isempty(ex.sppnames)
        print(io, "All the species names in the directory must be unique.")
    else
        msg = "The following names are already present in the species name list: "
        print(io, msg, ex.sppnames)
    end

    return nothing
end


text:
###################################################
## SchoolBusRouting.jl
##      Module root
## Authors: Arthur Delarue, Sébastien Martin, 2018
###################################################

module SchoolBusRouting

using DataFrames, CSV
using LightGraphs, NearestNeighbors, ProgressMeter
using Geodesy, Colors
using JuMP, Gurobi
import JLD
import Base: convert
using SFML, Plots

include("problem.jl")
include("load.jl")
include("scenarios.jl")
include("selection.jl")
include("tests.jl")
include("output.jl")
include("lbh.jl")
include("viz.jl")
include("synthetic.jl")

end

max_stars_repo_name: jihoyeo/SchoolBusRouting

max_issues_repo_name: jihoyeo/SchoolBusRouting

max_forks_repo_name: jihoyeo/SchoolBusRouting


text:
function freqresp(sys::DelayLtiSystem, ω::AbstractVector{T}) where {T <: Real}
    ny = noutputs(sys)
    nu = ninputs(sys)

    P_fr = freqresp(sys.P.P, ω);

    G_fr = zeros(eltype(P_fr), length(ω), ny, nu)

    for ω_idx=1:length(ω)
        P11_fr = P_fr[ω_idx, 1:ny, 1:nu]
        P12_fr = P_fr[ω_idx, 1:ny, nu+1:end]
        P21_fr = P_fr[ω_idx, ny+1:end, 1:nu]
        P22_fr = P_fr[ω_idx, ny+1:end, nu+1:end]

        delay_matrix_inv_fr = Diagonal(exp.(im*ω[ω_idx]*sys.Tau)) # Frequency response of the diagonal matrix with delays
        # Inverse of the delay matrix, so there should not be any minus signs in the exponents

        G_fr[ω_idx,:,:] .= P11_fr + P12_fr/(delay_matrix_inv_fr - P22_fr)*P21_fr # The matrix is invertible (?!)
    end

    return G_fr
end

function evalfr(sys::DelayLtiSystem, s)
    (ny, nu) = size(sys)

    P_fr = evalfr(sys.P.P, s)

    P11_fr = P_fr[1:ny, 1:nu]
    P12_fr = P_fr[1:ny, nu+1:end]
    P21_fr = P_fr[ny+1:end, 1:nu]
    P22_fr = P_fr[ny+1:end, nu+1:end]

    delay_matrix_inv_fr = Diagonal(exp.(s*sys.Tau))

    return P11_fr + P12_fr/(delay_matrix_inv_fr - P22_fr)*P21_fr
end

"""
    delayd_ss(τ, Ts)
Discrete-time statespace realization of a delay τ sampled with period Ts,
i.e. of z^-N where N = τ / Ts.

τ must be a multiple of Ts.
"""
function delayd_ss(τ::Number, Ts::Number)
    n = Int(round(τ / Ts))
    if !(τ - n*Ts ≈ 0)
        error("The delay τ must be a multiple of the sample time Ts")
    end
    ss(diagm(1 => ones(n-1)), [zeros(n-1,1); 1], [1 zeros(1,n-1)], 0, Ts)
end

"""
    c2d(G::DelayLtiSystem, Ts, method=:zoh)
"""
function c2d(G::DelayLtiSystem, Ts::Real, method=:zoh)
    if !(method === :zoh)
        error("c2d for DelayLtiSystems only supports zero-order hold")
    end
    X = append([delayd_ss(τ, Ts) for τ in G.Tau]...)
    Pd = c2d(G.P.P, Ts)
    return lft(Pd, X)
end

"""
    `y, t, x = lsim(sys::DelayLtiSystem, u, t::AbstractArray{<:Real}; x0=fill(0.0, nstates(sys)), alg=MethodOfSteps(Tsit5()), abstol=1e-6, reltol=1e-6, kwargs...)`

    Simulate system `sys`, over time `t`, using input signal `u`, with initial state `x0`, using method `alg` .

    Arguments:

    `t`: Has to be an `AbstractVector` with equidistant time samples (`t[i] - t[i-1]` constant)
    `u`: Function to determine control signal `ut` at a time `t`, on any of the following forms:
        Can be a constant `Number` or `Vector`, interpreted as `ut .= u` , or
        Function `ut .= u(t)`, or
        In-place function `u(ut, t)`. (Slightly more effienct)
    `alg, abstol, reltol` and `kwargs...`: are sent to `DelayDiffEq.solve`.

    Returns: times `t`, and `y` and `x` at those times.
"""
function lsim(sys::DelayLtiSystem{T,S}, u, t::AbstractArray{<:Real};
        x0=fill(zero(T), nstates(sys)),
        alg=DelayDiffEq.MethodOfSteps(Tsit5()),
        abstol=1e-6, reltol=1e-6,
        kwargs...) where {T,S}

    # Make u! in-place function of u
    u! = if isa(u, Number) || isa(u,AbstractVector) # Allow for u to be a constant number or vector
        (uout, t) -> uout .= u
    elseif DiffEqBase.isinplace(u, 2)               # If u is an inplace (more than 1 argument function)
        u
    else                                            # If u is a regular u(t) function
        (out, t) -> (out .= u(t))
    end

    _lsim(sys, u!, t, x0, alg; abstol=abstol, reltol=reltol, kwargs...)
end

# Generic parametrized dde used for simulating DelayLtiSystem
# We simulate the integral of the states
# u, du are from the notation of variables in DifferentialEquations.jl
# The state from the control system is x
# u!, uout is the control law and its output
function dde_param(du, u, h, p, t)
    A, B1, B2, C1, C2, D11, D12, D21, D22, y, Tau, u!, uout, hout, tmpy, tsave = p

    nx = size(A,1)
    nd = length(Tau)
    ny = size(C1,1)

    dx = view(du, 1:nx)
    dY = view(du, (nx+1):(nx+ny))
    dD = view(du, (nx+ny+1):(nx+ny+nd))
    x = view(u, 1:nx)

    # uout = u(t)
    u!(uout, t)

    # hout = d(t-Tau)
    for k=1:length(Tau)
        # Get the derivative of the history for d at indices corresponding to Tau
        hout[k] = h(p, t-Tau[k], Val{1}, idxs=(nx+ny+k))
    end

    #dx(t) .= A*x(t) + B1*u(t) +B2*d(t-tau)
    mul!(dx, A, x)
    mul!(dx, B1, uout, true, true)
    mul!(dx, B2, hout, true, true)

    # dY = y(t) = C1*x + D11*u(t) + D12*d(t-Tau)
    mul!(dY, C1, x)
    mul!(dY, D11, uout, true, true)
    mul!(dY, D12, hout, true, true)

    # dD = d(t) = C2*x + D21*u(t) + D22*d(t-Tau)
    mul!(dD, C2, x)
    mul!(dD, D21, uout, true, true)
    mul!(dD, D22, hout, true, true)

    # Save y value in tmpy to be used by dde_saver
    if t in tsave
        # The value of y at this time is given by the derivative
        tmpy .= dY
    end
    return
end

# Save x(t) and y(t) to output
function dde_saver(u,t,integrator)
    A, B1, B2, C1, C2, D11, D12, D21, D22, y, Tau, u!, uout, hout, tmpy, tsave = integrator.p
    nx = size(A,1)
    # y is already saved in tmpy
    u[1:nx], copy(tmpy)
end

function _lsim(sys::DelayLtiSystem{T,S}, Base.@nospecialize(u!), t::AbstractArray{<:Real}, x0::Vector{T}, alg; kwargs...) where {T,S}

    P = sys.P

    t0 = first(t)
    dt = t[2] - t[1]

    # Get all matrices to save on allocations
    A, B1, B2, C1, C2, D11, D12, D21, D22 = P.A, P.B1, P.B2, P.C1, P.C2, P.D11, P.D12, P.D21, P.D22
    Tau = sys.Tau
    nx = size(A,1)
    nd = length(Tau)
    ny = size(C1,1)
    nt = length(t)

    hout = fill(zero(T), nd)            # in place storage for delays
    uout = fill(zero(T), ninputs(sys))  # in place storage for u
    tmpy = similar(x0, ny)              # in place storage for output
    y = Matrix{T}(undef, ny, nt)        # Output matrix
    x = Matrix{T}(undef, nx, nt)        # State matrix

    p = (A, B1, B2, C1, C2, D11, D12, D21, D22, y, Tau, u!, uout, hout, tmpy, t)

    # This callback computes and stores the delay term
    sv = SavedValues(Float64, Tuple{Vector{Float64},Vector{Float64}})
    cb = SavingCallback(dde_saver, sv, saveat = t)

    # History function, only used for d
    # The true d(t) is stored as a derivative
    h!(p, t, deriv::Type{Val{1}}; idxs=0) = zero(T)

    # The states are x(t), Y(t), D(t), where Y, D are integrals of y(t), d(t)
    u0 = [x0;zeros(T,ny);zeros(T,nd)]
    prob = DDEProblem{true}(dde_param, u0, h!,
                (T(t[1]), T(t[end])),
                p,
                constant_lags=sort(Tau),# Not sure if sort needed
                neutral=true,           # We have derivatives on RHS (d(t))
                callback=cb)
    # Important to stop at t since we can not access derivatives in SavingCallback
    sol = DelayDiffEq.solve(prob, alg; tstops=t, saveat=t, kwargs...)

    # Retrive the saved values
    for k = 1:nt
        x[:,k] .= sv.saveval[k][1]
        y[:,k] .= sv.saveval[k][2]
    end

    return y', t, x'
end

# We have to default to something, look at the sys.P.P and delays
function _bounds_and_features(sys::DelayLtiSystem, plot::Val)
    ws, pz =  _bounds_and_features(sys.P.P, plot)
    logtau = log10.(abs.(sys.Tau))
    logtau = filter(x->x>4, logtau) # Ignore low frequency
    if isempty(logtau)
        return ws, pz
    end
    extreme = extrema(logtau)
    return [min(ws[1], floor(extreme[1]-0.2)), max(ws[2], ceil(extreme[2]+0.2))], pz
end

# Againm we have to do something for default vectors, more or less a copy from timeresp.jl
function _default_dt(sys::DelayLtiSystem)
    if !isstable(sys.P.P)
        return 0.05   # Something small
    else
        ps = pole(sys.P.P)
        r = minimum([abs.(real.(ps));0]) # Find the fastest pole of sys.P.P
        r = min(r, minimum([sys.Tau;0])) # Find the fastest delay
        if r == 0.0
            r = 1.0
        end
        return 0.07/r
    end
end

function Base.step(sys::DelayLtiSystem{T}, t::AbstractVector; kwargs...) where T
    nu = ninputs(sys)
    if t[1] != 0
        throw(ArgumentError("First time point must be 0 in step"))
    end
    u = (out, t) -> (t < 0 ? out .= 0 : out .= 1)
    x0=fill(zero(T), nstates(sys))
    if nu == 1
        y, tout, x = lsim(sys, u, t; x0=x0, kwargs...)
    else
        x = Array{T}(undef, length(t), nstates(sys), nu)
        y = Array{T}(undef, length(t), noutputs(sys), nu)
        for i=1:nu
            y[:,:,i], tout, x[:,:,i] = lsim(sys[:,i], u, t; x0=x0, kwargs...)
        end
    end
    return y, tout, x
end

function impulse(sys::DelayLtiSystem{T}, t::AbstractVector; alg=MethodOfSteps(BS3()), kwargs...) where T
    nu = ninputs(sys)
    iszero(sys.P.D12) || @warn("Impulse with a direct term from input to delay vector leads to poor accuracy.")
    if t[1] != 0
        throw(ArgumentError("First time point must be 0 in impulse"))
    end
    u = (out, t) -> (out .= 0)
    if nu == 1
        y, tout, x = lsim(sys, u, t; alg=alg, x0=sys.P.B[:,1], kwargs...)
    else
        x = Array{T}(undef, length(t), nstates(sys), nu)
        y = Array{T}(undef, length(t), noutputs(sys), nu)
        for i=1:nu
            y[:,:,i], tout, x[:,:,i] = lsim(sys[:,i], u, t; alg=alg, x0=sys.P.B[:,i], kwargs...)
        end
    end
    return y, tout, x
end

# Used for pade approximation
"""
`p2 = _linscale(p::Polynomial, a)`

Given a polynomial `p` and a number `a, returns the polynomials `p2` such that
`p2(s) == p(a*s)`.
"""
function _linscale(p::Polynomial, a)
    # This function should perhaps be implemented in Polynomials.jl
    coeffs_scaled = similar(p.coeffs, promote_type(eltype(p), typeof(a)))
    a_pow = 1
    coeffs_scaled[1] = p.coeffs[1]
    for k=2:length(p.coeffs)
        a_pow *= a
        coeffs_scaled[k] = p.coeffs[k]*a_pow
    end
    return Polynomial(coeffs_scaled)
end

# Coefficeints for Padé approximations
# PADE_Q_COEFFS = [Polynomial([binomial(N,i)*prod(N+1:2*N-i) for i=0:N]) for N=1:10]
const PADE_Q_COEFFS =  [[2, 1],
 [12, 6, 1],
 [120, 60, 12, 1],
 [1680, 840, 180, 20, 1],
 [30240, 15120, 3360, 420, 30, 1],
 [665280, 332640, 75600, 10080, 840, 42, 1],
 [17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1],
 [518918400, 259459200, 60540480, 8648640, 831600, 55440, 2520, 72, 1],
 [17643225600, 8821612800, 2075673600, 302702400, 30270240, 2162160, 110880, 3960, 90, 1],
 [670442572800, 335221286400, 79394515200, 11762150400, 1210809600, 90810720, 5045040, 205920, 5940, 110, 1]]

"""
    pade(τ::Real, N::Int)

Compute the `N`th order Padé approximation of a time-delay of length `τ`.
"""
function pade(τ::Real, N::Int)
    if !(1 <= N <= 10); error("Order of Padé approximation must be between 1 and 10. Got $N."); end

    Q = Polynomials.Polynomial(PADE_Q_COEFFS[N])

    return tf(_linscale(Q, -τ), _linscale(Q, τ)) # return Q(-τs)/Q(τs)
end

"""
    pade(G::DelayLtiSystem, N)

Approximate all time-delays in `G` by Padé approximations of degree `N`.
"""
function pade(G::DelayLtiSystem, N)
    ny, nu = size(G)
    nTau = length(G.Tau)
    X = append([ss(pade(τ,N)) for τ in G.Tau]...) # Perhaps append should be renamed blockdiag
    return lft(G.P.P, X)
end


text:
using Router

route("/", "todos#TodosController.index", named = :root)
route("/todos", "todos#TodosController.index", named = :todos)
route("/todos/:id::Int", "todos#TodosController.show", named = :todo_item)
route("/todos/:id::Int/edit", "todos#TodosController.edit")
route("/todos/:id::Int", "todos#TodosController.update", method = POST, named = :update_todo_item)
route("/todos/new", "todos#TodosController.new")
route("/todos", "todos#TodosController.create", method = POST, named = :create_todo_item)
route("/todos/:id::Int/delete", "todos#TodosController.delete", method = POST, named = :delete_todo_item)

# Authentication
route("/login", "user_sessions#UserSessionsController.show_login", named = :show_login)
route("/login", "user_sessions#UserSessionsController.login", method = POST, named = :login)
route("/logout", "user_sessions#UserSessionsController.logout", named = :logout)

route("/error") do
  error("foo!")
end

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
using CalibrationAnalysis
using Documenter

DocMeta.setdocmeta!(
    CalibrationAnalysis, :DocTestSetup, :(using CalibrationAnalysis); recursive=true
)

makedocs(;
    modules=[CalibrationAnalysis],
    authors="David Widmann",
    repo="https://github.com/devmotion/CalibrationAnalysis.jl/blob/{commit}{path}#{line}",
    sitename="CalibrationAnalysis.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/CalibrationAnalysis.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
    checkdocs=:exports,
    strict=true,
)

deploydocs(;
    repo="github.com/devmotion/CalibrationAnalysis.jl", push_preview=true, devbranch="main"
)