text:
include("TypeSynonyms.jl")
include("Design.jl")
include("General.jl")
include("Input.jl")
include("Ellipticals.jl")
include("Symbolic.jl")
include("Lshear-main.jl")
include("PlaneProblems.jl")
include("Assemble.jl")


text:
using Documenter, OneHotArrays

DocMeta.setdocmeta!(OneHotArrays, :DocTestSetup, :(using OneHotArrays); recursive = true)
makedocs(sitename = "OneHotArrays", doctest = false,
         pages = ["Overview" => "index.md",
                  "Reference" => "reference.md"])

deploydocs(repo = "github.com/FluxML/OneHotArrays.jl.git",
           target = "build",
           push_preview = true)

ext: jl


text:
line_pattern = r"(\d+)-(\d+) (\w): (\w+)"

input = map(eachline("a2.txt")) do line
    m = match(line_pattern, line)
    first, last, required, word = m.captures
    parse(Int16, first), parse(Int16, last), required[1], collect(word)
end

part1 = count(input) do (min, max, required, word)
  min <= count(==(required), word) <= max
end
print("$part1\n")

part2 = count(input) do (first, last, required, word)
  (word[first] == required) ⊻ (word[last] == required)
end
print("$part2\n")


text:
using CUDAnative

function kernel_vadd(a, b, c, i)
    c[i] = a[i] + b[i]
    return nothing
end

Base.code_llvm(stdout, kernel_vadd, Tuple{CuDeviceArray{Float32,1,AS.Global},
                                          CuDeviceArray{Float32,1,AS.Global},
                                          CuDeviceArray{Float32,1,AS.Global},
                                          Int})

using TinyCassette
struct Ctx end

Base.code_llvm(stdout, TinyCassette.execute,
               Tuple{Ctx, typeof(kernel_vadd),
                     CuDeviceArray{Float32,1,AS.Global},
                     CuDeviceArray{Float32,1,AS.Global},
                     CuDeviceArray{Float32,1,AS.Global},
                     Int})


text:
"""
    add_disturbance(sys::StateSpace, Ad::Matrix, Cd::Matrix)

See CCS pp. 144

# Arguments:
- `sys`: System to augment
- `Ad`: The dynamics of the disturbance
- `Cd`: How the disturbance states affect the states of `sys`. This matrix has the shape (sys.nx, size(Ad, 1))

See also [`add_low_frequency_disturbance`](@ref), [`add_resonant_disturbance`](@ref)
"""
function add_disturbance(sys::AbstractStateSpace, Ad::AbstractMatrix, Cd::AbstractMatrix)
    A,B,C,D = ControlSystems.ssdata(sys)
    T = eltype(A)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    Ae = [A Cd; zeros(T, size(Ad, 1), nx) Ad]
    Be = [B; zeros(T, size(Ad, 1), nu)]
    Ce = [C zeros(T, ny, size(Ad, 1))]
    De = D
    ss(Ae,Be,Ce,De,sys.timeevol)
end

"""
    add_measurement_disturbance(sys::StateSpace{Continuous}, Ad::Matrix, Cd::Matrix)
"""
function add_measurement_disturbance(sys::AbstractStateSpace{Continuous}, Ad::AbstractMatrix, Cd::AbstractMatrix)
    A,B,C,D = ControlSystems.ssdata(sys)
    T = eltype(A)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    Ae = [A zeros(T, nx, size(Ad, 1)); zeros(T, size(Ad, 1), nx) Ad]
    Be = [B; zeros(T, size(Ad, 1), nu)]
    Ce = [C Cd]
    De = D
    ss(Ae,Be,Ce,De)
end

"""
    add_low_frequency_disturbance(sys::StateSpace, Ai::Integer; ϵ = 0)

A disturbance affecting only state `Ai`.
"""
function add_low_frequency_disturbance(sys::AbstractStateSpace, Ai::Integer; ϵ=0)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    1 ≤ Ai ≤ nx || throw(ArgumentError("Ai must be a valid state index"))
    Cd = zeros(nx, 1)
    Cd[Ai] = 1
    Ad = -ϵ*I(nu)
    isdiscrete(sys) && (Ad += I)
    add_disturbance(sys, Ad, Cd)
end

"""
    add_low_frequency_disturbance(sys::StateSpace; ϵ = 0, measurement = false)

Augment `sys` with a low-frequency (integrating if `ϵ=0`) disturbance model.
If an integrating input disturbance is used together with an observer, the controller will have integral action.

# Arguments:
- `ϵ`: Move the integrator pole `ϵ` into the stable region.
- `measurement`: If true, the disturbance is a measurement disturbance, otherwise it's an input diturbance.
"""
function add_low_frequency_disturbance(sys::AbstractStateSpace; ϵ=0, measurement=false)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    Ad = -ϵ*I(nu)
    isdiscrete(sys) && (Ad += I)
    if measurement
        Cd = I(nu)
        add_measurement_disturbance(sys, Ad, Cd)
    else
        Cd = sys.B
        add_disturbance(sys, Ad, Cd)
    end
end

"""
    add_resonant_disturbance(sys::StateSpace{Continuous}, ω, ζ, Ai::Int; measurement = false)

Augment `sys` with a resonant disturbance model.

# Arguments:
- `ω`: Frequency
- `ζ`: Relative damping.
- `Ai`: The affected state
- `measurement`: If true, the disturbace is acting on the output, this will cause the controller to have zeros at ω (roots of poly s² + 2ζωs + ω²). If false, the disturbance is acting on the input, this will cause the controller to have poles at ω (roots of poly s² + 2ζωs + ω²).
"""
function add_resonant_disturbance(sys::AbstractStateSpace{Continuous}, ω, ζ, Ai::Integer; measurement=false)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    if measurement
        1 ≤ Ai ≤ sys.ny || throw(ArgumentError("Ai must be a valid output index"))
        Cd = zeros(ny, 2)
        Cd[Ai, 1] = 1
    else
        1 ≤ Ai ≤ sys.nx || throw(ArgumentError("Ai must be a valid state index"))
        Cd = zeros(nx, 2)
        Cd[Ai, 1] = 1
    end
    Ad = [-ζ -ω; ω -ζ]
    measurement ? add_measurement_disturbance(sys, Ad, Cd) : add_disturbance(sys, Ad, Cd)
end

"""
    add_differentiator(sys::StateSpace{<:Discrete})

Augment the output of `sys` with the numerical difference (discrete-time derivative) of output, i.e.,
`y_aug = [y; (y-y_prev)/sys.Ts]`
To add both an integrator and a differentiator to a SISO system, use
```julia
Gd = add_output_integrator(add_output_differentiator(G), 1)
```
"""
function add_output_differentiator(sys::AbstractStateSpace{<: Discrete}, diffsys=sys)
    A,B,C,D = ssdata(diffsys)
    all(iszero, D) || throw(ArgumentError("Can't add a differentiator to a system with non-zero D matrix. The system would not be proper."))
    C = C ./ sys.Ts
    Cd = C*(A-I)
    Dd = C*B
    A,B,C,D = ssdata(sys)
    ss(A, B, [C; Cd], [D; Dd], sys.timeevol)
end

function ControlSystems.tf(M::AbstractArray{TransferFunction{TE,ControlSystems.SisoRational{T}}}) where {TE, T<:Number}
    all(ControlSystems.issiso, M) || throw(ArgumentError("To make a MIMO system out of several MIMO systems is not yet supported"))
    matrix = first.(getproperty.(M, :matrix))
    TransferFunction{TE,ControlSystems.SisoRational{T}}(matrix, M[1].timeevol)
end

"""
add_output_integrator(sys::StateSpace{<:Discrete}, ind = 1; ϵ = 0)

Augment the output of `sys` with the integral of output at index `ind`, i.e.,
`y_aug = [y; ∫y[ind]]`
To add both an integrator and a differentiator to a SISO system, use
```julia
Gd = add_output_integrator(add_output_differentiator(G), 1)
```

Note: numerical integration is subject to numerical drift. If the output of the system corresponds to, e.g., a velocity reference and the integral to position reference, consider methods for mitigating this drift.
"""
function add_output_integrator(sys::AbstractStateSpace{<: Discrete}, ind=1; ϵ=0)
    int = tf(1.0*sys.Ts, [1, -(1-ϵ)], sys.Ts)
    𝟏 = tf(1.0,sys.Ts)
    𝟎 = tf(0.0,sys.Ts)
    M = [i==j ? 𝟏 : 𝟎 for i = 1:sys.ny, j = 1:sys.ny]
    M = [M; permutedims([i==ind ? int : 𝟎 for i = 1:sys.ny])]
    tf(M)*sys
end

"""
    add_input_integrator(sys::StateSpace, ui = 1, ϵ = 0)

Augment the output of `sys` with the integral of input at index `ui`, i.e.,
`y_aug = [y; ∫u[ui]]`
See also [`add_low_frequency_disturbance`](@ref)
"""
function add_input_integrator(sys::AbstractStateSpace, ui=1; ϵ=0)
    A,B,C,D = ControlSystems.ssdata(sys)
    T = eltype(A)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    1 ≤ ui ≤ nu || throw(ArgumentError("ui must be a valid input index"))
    Cd = zeros(T, 1, nx+1)
    Cd[end] = 1
    Bd = zeros(T, 1, nu)
    Bd[ui] = ControlSystems.isdiscrete(sys) ? sys.Ts : 1
    Ad = -ϵ*I(1)
    isdiscrete(sys) && (Ad += I)

    Ae = [A zeros(T, nx, 1); zeros(T, size(Ad, 1), nx) Ad]
    Be = [B; Bd]
    Ce = [[C zeros(T, ny, 1)]; Cd]
    De = [D; zeros(T, 1, nu)]
    ss(Ae,Be,Ce,De,sys.timeevol)

end

"""
    add_input_differentiator(sys::StateSpace, ui = 1:sys.nu; goodwin=false)

Augment the output of `sys` with the difference `u(k+1)-u(k)`

# Arguments:
- `ui`: An index or vector of indices indicating which inputs to differentiate.
- `goodwin`: If true, the difference operator will use the Goodwin δ operator, i.e., `(u(k+1)-u(k)) / sys.Ts`.

The augmented system will have the matrices
```
[A 0; 0 0]  [B; I]  [C 0; 0 -I]  [D; I]
```
with `length(ui)` added states and outputs.
"""
function add_input_differentiator(sys::AbstractStateSpace{<:Discrete}, ui=1:sys.nu; goodwin=false)
    A,B,C,D = ControlSystems.ssdata(sys)
    T = eltype(A)
    nx,nu,ny = sys.nx,sys.nu,sys.ny
    all(1 .≤ ui .≤ nu) || throw(ArgumentError("ui must be a valid input index"))
    nnu = length(ui) # number of new states and outputs

    den = goodwin ? 1/sys.Ts : 1

    Cd = zeros(T, nnu, nx+nnu)
    Cd[:, nx+1:end] .= -den*I(nnu)
    Bd = zeros(T, nnu, nu)
    for (i, ui) in enumerate(ui)
        Bd[i, ui] = 1
    end
    Ad = zeros(nnu, nnu)
    Dd = den*I(nnu)

    Ae = [A zeros(T, nx, nnu); zeros(T, size(Ad, 1), nx) Ad]
    Be = [B; Bd]
    Ce = [[C zeros(T, ny, nnu)]; Cd]
    De = [D; Dd]
    ss(Ae,Be,Ce,De,sys.timeevol)

end

# using ControlSystems.DemoSystems
# sys = DemoSystems.resonant()
# sys2 = add_low_frequency_disturbance(sys, 2)
# sys25 = add_low_frequency_disturbance(sys)
# sys3 = add_resonant_disturbance(sys, 1, 0.5, 1)
# ss([0], [0], [1], 1)*sys
# sys + ss([0.0], [0], [1], 1)


text:
@test sprint(show, h1) == "{5, 0} undirected graph"
@test sprint(show, h3) == "empty undirected graph"

@test Graph(DiGraph(g3)) == g3

@test degree(g3, 1) == 1
# @test all_neighbors(g3, 3) == [2, 4]
@test density(g3) == 0.4

g = Graph(5)
@test add_edge!(g, 1, 2)

e2 = Edge(1,3)
e3 = Edge(1,4)
e4 = Edge(2,5)
e5 = Edge(3,5)

@test add_edge!(g, e2)
@test add_edge!(g, e3)
@test add_edge!(g, e4)
@test add_edge!(g, e5)

h = DiGraph(5)
@test add_edge!(h, 1, 2)
@test add_edge!(h, e2)
@test add_edge!(h, e3)
@test add_edge!(h, e4)
@test add_edge!(h, e5)

@test LightGraphs.fadj(g)[1] == LightGraphs.fadj(g,1) ==
    LightGraphs.badj(g)[1] == LightGraphs.badj(g,1) ==
    LightGraphs.adj(g)[1] == LightGraphs.adj(g,1) ==
    [2,3,4]

@test sprint(show, h4) == "{7, 0} directed graph"
@test sprint(show, h5) == "empty directed graph"
@test has_edge(g, e1)
@test has_edge(h, e1)
@test !has_edge(g, e0)
@test !has_edge(h, e0)

@test degree(g4, 1) == 1
# @test all_neighbors(g4, 3) == [4]
@test density(g4) == 0.2

@test nv(a1) == 3
@test ne(a1) == 2
@test nv(a2) == 3
@test ne(a2) == 5

badadjmx = [ 0 1 0; 1 0 1]
@test_throws ErrorException Graph(badadjmx)
@test_throws ErrorException Graph(sparse(badadjmx))
@test_throws ErrorException DiGraph(badadjmx)
@test_throws ErrorException DiGraph(sparse(badadjmx))
@test_throws ErrorException Graph([1 0; 1 1])

@test !add_edge!(g, 100, 100)
@test !add_edge!(h, 100, 100)

@test_throws ErrorException Graph(sparse(adjmx2))

g = Graph(sparse(adjmx1))
h = DiGraph(sparse(adjmx1))

@test (nv(g), ne(g)) == (3, 2)
@test (nv(h), ne(h)) == (3, 4)
@test Graph(h) == g

@test sort(all_neighbors(WheelDiGraph(10),2)) == [1, 3, 10]


text:
using Setfield
replacein(x, oldnode, newnode) = x
replacein(x::Tuple, oldnode, newnode) = tuple([replacein(m, oldnode, newnode) for m in x]...)
replacein(x::NamedTuple, oldnode, newnode) = (;[k => replacein(x[k], oldnode, newnode) for k in keys(x)]...)

function replacein(x::T, oldnode, newnode) where {T<:Union{AbstractNode, AbstractMillModel}}
    x === oldnode && return(newnode)
    fields = map(f -> replacein(getproperty(x, f), oldnode, newnode), fieldnames(T))
    n = nameof(T)
    p = parentmodule(T)
    eval(:($p.$n))(fields...)
end

function replacein(x::LazyNode{N}, oldnode, newnode) where {N}
    x === oldnode && return(newnode)
    LazyNode{N}(replacein(x.data, oldnode, newnode))
end

function replacein(x::LazyModel{N}, oldnode, newnode) where {N}
    x === oldnode && return(newnode)
    LazyModel{N}(replacein(x.m, oldnode, newnode))
end

function findin(x, node)
    x === node && return(@lens _)
    return(nothing)
end

function findin(x::T, node) where {T<:Union{AbstractNode, AbstractMillModel}}
    x === node && return(@lens _)
    for k in fieldnames(T)
        l = findin(getproperty(x, k), node)
        if l != nothing
            lo = Setfield.PropertyLens{k}() ∘ l
            return(lo)
        end
    end
    return(nothing)
end

function findin(x::NamedTuple, node)
	x === node && return(@lens _)
    for k in keys(x)
    	l = findin(x[k], node)
    	if l != nothing
    		lo = Setfield.PropertyLens{k}() ∘ l
    		return(lo)
    	end
    end
    return(nothing)
end

function findin(x::Tuple, node)
    error("findin does not support Tuples due to restrinctions of Lens from Setfield.")
end

# function findin(x::T, node) where {T<:ArrayNode}
# 	x === node && return(@lens _)
# 	x.data === node && return(@lens _.data)
#     return(nothing)
# end

# function findin(x::LazyNode{N}, node) where {N}
# 	x === node && return(@lens _)
# 	x.data === node && return(@lens _.data)
#     return(nothing)
# end

# function findin(x::AbstractBagNode, node)
# 	x === node && return(@lens _)
# 	l = findin(x.data, node)
# 	if l != nothing
# 		return((@lens _.data) ∘  l)
# 	end
#     return(nothing)
# end


text:
import Interpolations
import Interpolations: interpolate, extrapolate, scale, BSpline, Quadratic, OnCell
using SparseArrays

"""
Interpolate a function expressed in a basis `basis_in` to a basis `basis_out`
This interpolation uses a very basic real-space algorithm, and makes
a DWIM-y attempt to take into account the fact that basis_out can be a supercell of basis_in
"""
function interpolate_density(ρ_in, basis_in::PlaneWaveBasis, basis_out::PlaneWaveBasis)
    ρ_out = interpolate_density(ρ_in, basis_in.fft_size, basis_out.fft_size,
                                basis_in.model.lattice, basis_out.model.lattice)
end

# TODO Specialization for the common case lattice_out = lattice_in
function interpolate_density(ρ_in::AbstractArray, grid_in, grid_out, lattice_in, lattice_out=lattice_in)
    T = real(eltype(ρ_in))
    @assert size(ρ_in) == grid_in

    # First, build supercell, array of 3 ints
    supercell = zeros(Int, 3)
    for i = 1:3
        if norm(lattice_in[:, i]) == 0
            @assert norm(lattice_out[:, i]) == 0
            supercell[i] = 1
        else
            supercell[i] = round(Int, norm(lattice_out[:, i]) / norm(lattice_in[:, i]))
        end
        if norm(lattice_out[:, i] - supercell[i]*lattice_in[:, i]) > .3*norm(lattice_out[:, i])
            @warn "In direction $i, the output lattice is very different from the input lattice"
        end
    end

    # ρ_in represents a periodic function, on a grid 0, 1/N, ... (N-1)/N
    grid_supercell = grid_in .* supercell
    ρ_in_supercell = similar(ρ_in, (grid_supercell...))
    for i = 1:supercell[1]
        for j = 1:supercell[2]
            for k = 1:supercell[3]
                ρ_in_supercell[
                    1 + (i-1)*grid_in[1] : i*grid_in[1],
                    1 + (j-1)*grid_in[2] : j*grid_in[2],
                    1 + (k-1)*grid_in[3] : k*grid_in[3]] = ρ_in
            end
        end
    end

    # interpolate ρ_in_supercell from grid grid_supercell to grid_out
    axes_in = (range(0, 1, length=grid_supercell[i]+1)[1:end-1] for i=1:3)
    itp = interpolate(ρ_in_supercell, BSpline(Quadratic(Interpolations.Periodic(OnCell()))))
    sitp = scale(itp, axes_in...)
    ρ_interp = extrapolate(sitp, Periodic())
    ρ_out = similar(ρ_in, grid_out)
    for i = 1:grid_out[1]
        for j = 1:grid_out[2]
            for k = 1:grid_out[3]
                ρ_out[i, j, k] = ρ_interp((i-1)/grid_out[1],
                                          (j-1)/grid_out[2],
                                          (k-1)/grid_out[3])
            end
        end
    end

    ρ_out
end

"""
Interpolate some data from one ``k``-point to another. The interpolation is fast, but not
necessarily exact or even normalized. Intended only to construct guesses for iterative
solvers
"""
function interpolate_kpoint(data_in::AbstractVecOrMat,
                            basis_in::PlaneWaveBasis,  kpoint_in::Kpoint,
                            basis_out::PlaneWaveBasis, kpoint_out::Kpoint)
    # TODO merge with transfer_blochwave_kpt
    if kpoint_in == kpoint_out
        return copy(data_in)
    end
    @assert length(G_vectors(basis_in, kpoint_in)) == size(data_in, 1)

    n_bands  = size(data_in, 2)
    n_Gk_out = length(G_vectors(basis_out, kpoint_out))
    data_out = similar(data_in, n_Gk_out, n_bands) .= 0
    for iin in 1:size(data_in, 1)
        idx_fft = kpoint_in.mapping[iin]
        idx_fft in keys(kpoint_out.mapping_inv) || continue
        iout = kpoint_out.mapping_inv[idx_fft]
        data_out[iout, :] = data_in[iin, :]
    end
    data_out
end

"""
Compute the index mapping between two bases. Returns two arrays
`idcs_in` and `idcs_out` such that `ψkout[idcs_out] = ψkin[idcs_in]` does
the transfer from `ψkin` (defined on `basis_in` and `kpt_in`) to `ψkout`
(defined on `basis_out` and `kpt_out`).
"""
function transfer_mapping(basis_in::PlaneWaveBasis{T},  kpt_in::Kpoint,
                          basis_out::PlaneWaveBasis{T}, kpt_out::Kpoint) where T
    idcs_in  = 1:length(G_vectors(basis_in, kpt_in))  # All entries from idcs_in
    kpt_in == kpt_out && return idcs_in, idcs_in

    # Get indices of the G vectors of the old basis inside the new basis.
    idcs_out = index_G_vectors.(Ref(basis_out), G_vectors(basis_in, kpt_in))

    # In the case where G_vectors(basis_in.kpoints[ik]) are bigger than vectors
    # in the fft_size box of basis_out, we need to filter out the "nothings" to
    # make sure that the index linearization works. It is not an issue to
    # filter these vectors as this can only happen if Ecut_in > Ecut_out.
    if any(isnothing, idcs_out)
        idcs_in  = idcs_in[idcs_out .!= nothing]
        idcs_out = idcs_out[idcs_out .!= nothing]
    end
    idcs_out = getindex.(Ref(LinearIndices(basis_out.fft_size)), idcs_out)

    # Map to the indices of the corresponding G-vectors in
    # G_vectors(basis_out, kpt_out) this array might contains some nothings if
    # basis_out has less G_vectors than basis_in at this k-point
    idcs_out = indexin(idcs_out, kpt_out.mapping)
    if any(isnothing, idcs_out)
        idcs_in  = idcs_in[idcs_out .!= nothing]
        idcs_out = idcs_out[idcs_out .!= nothing]
    end

    idcs_in, idcs_out
end

"""
Return a sparse matrix that maps quantities given on `basis_in` and `kpt_in`
to quantities on `basis_out` and `kpt_out`.
"""
function compute_transfer_matrix(basis_in::PlaneWaveBasis{T}, kpt_in::Kpoint,
                                 basis_out::PlaneWaveBasis{T}, kpt_out::Kpoint) where T
    idcs_in, idcs_out = transfer_mapping(basis_in, kpt_in, basis_out, kpt_out)
    sparse(idcs_out, idcs_in, true)
end

"""
Return a list of sparse matrices (one per ``k``-point) that map quantities given in the
`basis_in` basis to quantities given in the `basis_out` basis.
"""
function compute_transfer_matrix(basis_in::PlaneWaveBasis{T}, basis_out::PlaneWaveBasis{T}) where T
    @assert basis_in.model.lattice == basis_out.model.lattice
    @assert length(basis_in.kpoints) == length(basis_out.kpoints)
    @assert all(basis_in.kpoints[ik].coordinate == basis_out.kpoints[ik].coordinate
                for ik in 1:length(basis_in.kpoints))
    [compute_transfer_matrix(basis_in, kpt_in, basis_out, kpt_out)
     for (kpt_in, kpt_out) in zip(basis_in.kpoints, basis_out.kpoints)]
end

"""
Transfer an array ψk defined on basis_in ``k``-point kpt_in to basis_out ``k``-point kpt_out.
"""
function transfer_blochwave_kpt(ψk_in, basis_in::PlaneWaveBasis{T}, kpt_in::Kpoint,
                                basis_out::PlaneWaveBasis{T}, kpt_out::Kpoint) where T
    kpt_in == kpt_out && return copy(ψk_in)
    @assert length(G_vectors(basis_in, kpt_in)) == size(ψk_in, 1)
    idcsk_in, idcsk_out = transfer_mapping(basis_in, kpt_in, basis_out, kpt_out)

    n_bands = size(ψk_in, 2)
    ψk_out  = similar(ψk_in, length(G_vectors(basis_out, kpt_out)), n_bands)
    ψk_out .= 0
    ψk_out[idcsk_out, :] .= ψk_in[idcsk_in, :]

    ψk_out
end

"""
Transfer Bloch wave between two basis sets. Limited feature set.
"""
function transfer_blochwave(ψ_in, basis_in::PlaneWaveBasis{T},
                            basis_out::PlaneWaveBasis{T}) where T
    @assert basis_in.model.lattice == basis_out.model.lattice
    @assert length(basis_in.kpoints) == length(basis_out.kpoints)
    @assert all(basis_in.kpoints[ik].coordinate == basis_out.kpoints[ik].coordinate
                for ik in 1:length(basis_in.kpoints))

    # If, for some kpt ik, basis_in has less vectors than basis_out, then idcs_out[ik] is
    # the array of the indices of the G_vectors from basis_in in basis_out.
    # It is then of size G_vectors(basis_in.kpoints[ik]) and the transfer can be done with
    # ψ_out[ik] .= 0
    # ψ_out[ik][idcs_out[ik], :] .= ψ_in[ik]

    # Otherwise, if, for some kpt ik, basis_in has more vectors than basis_out, then
    # idcs_out[ik] just keep the indices of the G_vectors from basis_in that are in basis_out.
    # It is then of size G_vectors(basis_out.kpoints[ik]) and the transfer can be done with
    # ψ_out[ik] .= ψ_in[ik][idcs_in[ik], :]

    map(enumerate(basis_out.kpoints)) do (ik, kpt_out)
        transfer_blochwave_kpt(ψ_in[ik], basis_in, basis_in.kpoints[ik], basis_out, kpt_out)
    end
end


text:
module BasicTestTargetTests
using Test
using Markdown
using BasicTestTarget

end


text:
include("repl.jl")

listing(code, output=showoutput(code)) = begin
    output = output == nothing ? showoutput(code) : output
    input = Signal(code)
    cell = hbox(
        code_io(code, input) |> width(30em),
        hskip(0.5em),
        vline(),
        hskip(0.5em),
        vbox(
            map(showoutput, input, typ=Any, init=output)
        ) |> wrap |> width(30em) |> Escher.pad(0.5em) |> fillcolor("white")
    ) |> Escher.pad(1em) |> paper(1)
    hbox(cell, flex())
end


text:
using Glob
using JLD

function searchAndSave(FolderPath,Name,SaveFilename,saveKey)
    Filename=glob(Name,FolderPath)
    n=sizeof(FolderPath)
    for k=1:size(Filename,1)
        println(Filename[k][n+1:end-4]) #
        Filename[k]=Filename[k][n+1:end-4]
    end

    save(SaveFilename, saveKey, Filename)

end

function searchAndSaveMTiXsens(FolderPath,NameID,filenameMTiXsensID,SaveFilename,saveKey)
    MTiXsensFilename=glob(string("*",NameID),FolderPath)
    n=sizeof(FolderPath)
    m=sizeof(NameID)
    for k=1:size(MTiXsensFilename,1)
        println(MTiXsensFilename[k][n+1:end-m]) #
        MTiXsensFilename[k]=MTiXsensFilename[k][n+1:end-m]
    end

    sensorXsensNames=collect(keys(filenameMTiXsensID));

    save(SaveFilename, saveKey[1], MTiXsensFilename,
                        saveKey[2], filenameMTiXsensID,
                        saveKey[3], sensorXsensNames)

end


text:
function fib(i::Int64)::Int64
    if i == 0 || i == 1
        return 1
    end
    return fib(i - 1) + fib(i - 2)
end

if abspath(PROGRAM_FILE) == @__FILE__
    @assert(fib(0) == 1)
    @assert(fib(1) == 1)
    @assert(fib(5) == 8)
    @assert(fib(30) == 1346269)
    println("OK")
end


text:
@testset "Test fft_view_tuple"  begin
    N = 100
    x = randn((N, N))

    yi, yr, ya = fft_view_tuple(x)
    @test ya ≈ abs.(yi * 1im + yr)

    yi, yr, ya = fft_view_tuple(x, false)
    @test ya ≈ abs.(yi * 1im + yr)

    yi, yr, ya = fft_view_tuple(x, false, [1])
    @test ya ≈ abs.(yi * 1im + yr)

    yi, yr, ya = fft_view_tuple(x, false, logarithmic=true)
    x_fft = fft(x)
    @test yi ≈ log.(1 .+ abs.(imag(x_fft)))

end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
using Serialization

function deref(pointer::Int, program::Array{Int,1}, mode::Int, base::Int)::Int
    if mode == 0
        return program[pointer+1]
    elseif mode == 1
        return pointer
    elseif mode == 2
        return program[pointer+base+1]
    else
        throw(ErrorException("Invalid mode: " * mode))
    end
end

function write_pointer(pointer::Int, mode::Int, base::Int)::Int
    if mode == 0
        return pointer+1
    elseif mode == 2
        return pointer+base+1
    else
        throw(ErrorException("Invalid mode for writing: " * mode))
    end
end

function run_intcode(program::Vector{Int}, input::Channel, output::Channel)::Int
    program = vcat(program, zeros(Int, 1000))
    pointer = 1
    ret = nothing
    base = 0
    while pointer <= length(program)
        # println(program)
        instruction = reverse(digits(program[pointer], pad=5))
        # println(instruction)
        opcode = instruction[end] + instruction[end-1]*10
        modes = reverse(instruction[1:end-2])
        # println(join([opcode,": ",modes]))
        if opcode == 1 # add
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 + val2
            pointer = pointer + 4
        elseif opcode == 2 # multiply
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 * val2
            pointer = pointer + 4
        elseif opcode == 3 # store input
            pt1 = program[pointer+1]
            program[write_pointer(pt1, modes[1], base)] = take!(input)
            pointer = pointer + 2
        elseif opcode == 4 # output
            pt1 = program[pointer+1]
            val1 = deref(pt1, program, modes[1], base)
            pointer = pointer + 2
            ret = val1
            put!(output, val1)
        elseif opcode == 5 # jump if true
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            if val1 != 0
                pointer = val2 + 1
            else
                pointer = pointer + 3
            end
        elseif opcode == 6 # jump if false
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            if val1 == 0
                pointer = val2 + 1
            else
                pointer = pointer + 3
            end
        elseif opcode == 7 # less than
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 < val2 ? 1 : 0
            pointer = pointer + 4
        elseif opcode == 8 # equal
            pt1 = program[pointer+1]
            pt2 = program[pointer+2]
            pt3 = program[pointer+3]
            val1 = deref(pt1, program, modes[1], base)
            val2 = deref(pt2, program, modes[2], base)
            program[write_pointer(pt3, modes[3], base)] = val1 == val2 ? 1 : 0
            pointer = pointer + 4
        elseif opcode == 9 # modify base
            pt1 = program[pointer+1]
            val1 = deref(pt1, program, modes[1], base)
            base = base + val1
            pointer = pointer + 2
        elseif opcode == 99
            close(output)
            return ret
        else
            throw(ErrorException("Invalid opcode: " * string(opcode)))
        end
    end
end

function draw_map(map::Array{Int, 2}, min::Vector{Int}, max::Vector{Int})
    map = string.(map)
    for row in eachrow(map[min[2]-1:max[2]+1,:])
        println(replace(replace(replace(replace(replace(replace(join(row[min[1]-1:max[1]+1]), "0" => " "), "1" => "."), "2" => "#"), "4" => "S"), "3" => "F"), "7" => "O"))
    end
end

function draw_map(map::Array{Int, 2}, robot::Vector{Int}, target::Vector{Int}, min::Vector{Int}, max::Vector{Int})
    map = string.(map)
    map[robot[2], robot[1]] = "^"
    map[target[2], target[1]] = "*"
    for row in eachrow(map[min[2]-1:max[2]+1,:])
        println(replace(replace(replace(replace(replace(join(row[min[1]-1:max[1]+1]), "0" => " "), "1" => "."), "2" => "#"), "4" => "S"), "3" => "F"))
    end
end

function draw_map(map::Array{Int, 2})
    min = findfirst(x -> x != 0, map)
    max = findlast(x -> x != 0, map)
    draw_map(map, [min[2], min[1]], [max[2], max[1]])
end

function draw_map(map::Array{Int, 2}, robot::Vector{Int}, direction::Int, queue::Vector{Vector{Int}})
    map = deepcopy(map)
    for q in queue
        map[q[2], q[1]] = 5 # queue
    end
    map[queue[1][2], queue[1][1]] = 6 # next objective
    map = string.(map)
    map[robot[2], robot[1]] = if direction == 0
        "^"
    elseif direction == 90
        "<"
    elseif direction == 180
        "v"
    elseif direction == 270
        ">"
    end
    for row in eachrow(map[robot[2]-2:robot[2]+2,:])
        println(replace(replace(replace(replace(replace(replace(replace(join(row[robot[1]-3:robot[1]+3]), "0" => " "), "1" => "."), "2" => "#"), "4" => "S"), "5" => "?"), "6" => "*"), "3" => "F"))
    end
end

function update_bounding_box(point, min, max)
    if point[1] < min[1]
        min[1] = point[1]
    end
    if point[1] > max[1]
        max[1] = point[1]
    end
    if point[2] < min[2]
        min[2] = point[2]
    end
    if point[2] > max[2]
        max[2] = point[2]
    end

    return min, max
end

function route(start::Vector{Int}, target::Vector{Int}, arr)
    i = 0
    arr = map(x -> x == 0 || x == 2 ? -99 : -1, arr) # Only consider visited points walkable
    # arr = map(x -> x == 2 ? -99 : -1, arr) # Consider unvisited points walkable, -1 => walkable

    # 1. init
    arr[start[2], start[1]] = i
    arr[target[2], target[1]] = -1 # NOTE: target has to be reachable

    # 2. wave expansion
    target_reached = false
    while !target_reached # && i < 10
        # println("Wave expansion #", i)
        res = findall(x -> x == i, arr)
        if length(res) == 0
            throw(ErrorException("Could not find path to "*string(target)*" from "*string(start)))
        end
        for p in res
            # println(p)
            neighbors = [
                [p[1]-1, p[2]],
                [p[1]+1, p[2]],
                [p[1], p[2]+1],
                [p[1], p[2]-1],
            ]
            for n in neighbors
                # println(n," ",arr[n[1], n[2]])
                if arr[n[1], n[2]] == -1
                    # println(n," ",arr[n[2], n[1]])
                    arr[n[1], n[2]] = i + 1
                    if n == reverse(target)
                        # println("Target found at ", n, " on iteration #",i)
                        target_reached = true
                        break
                    end
                end
            end
        end
        i = i + 1
    end

    # 3. backtrace
    path = []
    p = reverse(target)
    while p != reverse(start)
        neighbors = [
            [p[1]-1, p[2]],
            [p[1]+1, p[2]],
            [p[1], p[2]+1],
            [p[1], p[2]-1],
        ]
        for n in neighbors
            if arr[n[1],n[2]] == arr[p[1],p[2]] - 1
                # println(n)
                push!(path, n)
                p = n
                break
            end
        end
    end

    return path
end

# https://en.wikipedia.org/wiki/Lee_algorithm
function find_path_to(target::Vector, start::Vector, arr)::Int
    global cache
    # 0 = unvisited
    # 1 = visited
    # 2 = wall
    # 3 = finish
    # 4 = start
    distance = sqrt((start[1]-target[1])^2 + (start[2]-target[2])^2)
    # println("Start: ", start)
    # println("Target: ", target)
    # println("Distance: ", round(distance, digits=1))

    if arr[target[2], target[1]] == 2
        throw(ErrorException("Target is a wall"))
    end

    min = deepcopy(start)
    max = deepcopy(start)
    min, max = update_bounding_box(target, min, max)

    if distance > 1 && haskey(cache, target)
        # println("Possible cached paths")
        for path in cache[target]
            res = findlast(x -> x == reverse(start), path)
            if res != nothing
                # println("Cache hit!")
                target = reverse(path[res-1])
                break
            end
        end
    end

    if target == [start[1]-1, start[2]]
        # println("shortcircuit to left")
        return 90
    elseif target == [start[1]+1, start[2]]
        # println("shortcircuit to right")
        return 270
    elseif target == [start[1], start[2]+1]
        # println("shortcircuit to down")
        return 180
    elseif target == [start[1], start[2]-1]
        # println("shortcircuit to up")
        return 0
    end

    path = route(start, target, arr)

    # println("Write cache from ", start, " to ", target)
    if !haskey(cache, target)
        cache[target] = [path]
    else
        push!(cache[target], path)
    end

    res = reverse(path[end-1])

    if  res == [start[1]-1, start[2]]
        # println("to reach target go left")
        return 90
    elseif res == [start[1]+1, start[2]]
        # println("to reach target go right")
        return 270
    elseif res == [start[1], start[2]+1]
        # println("to reach target go down")
        return 180
    elseif res == [start[1], start[2]-1]
        # println("to reach target go up")
        return 0
    else
        throw(ErrorException("Node not next to start"))
    end

    throw(ErrorException("You shoudln't be here"))
end

function process_output(in::Channel, out::Channel)
    room = zeros(Int, 1000, 1000) # unvisited
    room[500, 500] = 4 # start
    robot = [500, 500]
    max = [500, 500] # drawing boundaries
    min = [500, 500]
    direction = 0 # up
    queue = [
        [robot[1]-1, robot[2]],
        [robot[1]+1, robot[2]],
        [robot[1], robot[2]+1],
        [robot[1], robot[2]-1],
        ]
    counter = 0
    status = -1
    # println(counter, ": robot at ",robot,":")
    # draw_map(room, robot, direction, queue)
    # println("Queue: ", length(queue))
    while true
        counter = counter + 1

        target = deepcopy(robot)
        if direction == 0
            put!(out, 1) # north/up
            target[2] = target[2] - 1
        elseif direction == 90
            put!(out, 3) # west/left
            target[1] = target[1] - 1
        elseif direction == 180
            put!(out, 2) # south/down
            target[2] = target[2] + 1
        elseif direction == 270
            put!(out, 4) # east/right
            target[1] = target[1] + 1
        else
            throw(ErrorException("invalid direction: "*direction ))
        end

        min, max = update_bounding_box(target, min, max)

        try
            status = take!(in)
        catch e
            if isa(e, InvalidStateException)
                break
            else
                rethrow
            end
        end

        if status == 0
            if room[target[2], target[1]] == 0
                room[target[2], target[1]] = 2 # wall
            end
        elseif status == 1
            robot = target
            if room[target[2], target[1]] == 0
                room[target[2], target[1]] = 1 # visited
            end
        elseif status == 2
            robot = target
            room[target[2], target[1]] = 3 # oxygen
            # break # FIXME: we use robot to exhaust the map, not to find target
        else
            throw(ErrorException("invalid status: "*direction ))
        end

        neighbors = [
            [robot[1]-1, robot[2]],
            [robot[1]+1, robot[2]],
            [robot[1], robot[2]+1],
            [robot[1], robot[2]-1],
        ]

        # println("Queue (", length(queue), "): ", queue)
        # println("Neighbors: ", neighbors)
        queue = filter(x -> room[x[2], x[1]] == 0, queue)
        neighbors = filter(x -> room[x[2], x[1]] == 0, neighbors)
        append!(queue, neighbors)
        if length(queue) == 0
            break # map exhausted
        end
        prepend!(queue, filter(x -> sqrt((x[1]-robot[1])^2 + (x[2]-robot[2])^2) == 1, queue))
        unique!(queue)
        # sort!(queue, by = x -> sqrt((x[1]-robot[1])^2 + (x[2]-robot[2])^2))

        # if counter % 200 == 0
            # println(counter, ": robot at ",robot,":")
            #draw_map(room, robot, direction, queue)
            #  draw_map(room, robot, queue[1], min, max)
            # println("Queue: ", length(queue))
            #    println(queue)
#        end

       direction = find_path_to(queue[1], robot, room)

    end
    # println("Mapped in ", counter) # 2444
    return room
end

function draw_oxygen(arr)
    arr = map(x -> x >= 0 ? 7 : x, arr) # oxygen
    arr = map(x -> x == -1 ? 1 : x, arr) # visited (empty of oxygen)
    arr = map(x -> x == -99 ? 2 : x, arr) # walls back
    arr = map(x -> x == -98 ? 0 : x, arr) # emptiness back

    draw_map(arr)
end

function fill(arr, start)
    i = 0
    arr = map(x -> x == 2 ? -99 : x == 0 ? -98 : -1, arr) # block walls

    # 1. init
    arr[start[2], start[1]] = i

    # 2. wave expansion
    while true
        # println("Wave expansion #", i)
        # draw_oxygen(arr)
        res = findall(x -> x == i, arr)
        if length(res) == 0
            break
        end
        for p in res
            # println(p)
            neighbors = [
                [p[1]-1, p[2]],
                [p[1]+1, p[2]],
                [p[1], p[2]+1],
                [p[1], p[2]-1],
            ]
            for n in neighbors
                # println(n," ",arr[n[1], n[2]])
                if arr[n[1], n[2]] == -1
                    # println(n," ",arr[n[2], n[1]])
                    arr[n[1], n[2]] = i + 1
                end
            end
        end
        i = i + 1
    end

    # draw_oxygen(arr)

    return i-1
end

cache = Dict()

function map_room()
    program = open("droid.csv") do f
        parse.(Int, split(readlines(f)[1], ","))
    end

    in = Channel(1)
    out = Channel(2)

    @async run_intcode(deepcopy(program), in, out)

    room = process_output(out, in)
    serialize("room.jls", room)

    return room
end

room =  if isfile("room.jls")
            deserialize("room.jls")
        else
            map_room()
        end
# draw_map(room)
target = findfirst(x -> x == 3, room)
steps = length(route([500,500], [target[2], target[1]], room))
println("Shortest path to oxygen in ", steps, " steps")
steps = fill(room, [target[2], target[1]])
println("It takes oxygen ", steps, " minutes to fill the space")


text:
using DSP
using FFTW

"Make window function in frequenVV domain with 3x3 size for crop=1."
function make_freq_window(sz, crop=1; option1=3.0, window_type="kaiser")
    # Ref: https://github.com/miaosiSari/Regridding/blob/master/get_window.m
    if window_type == "hanning"
        vec_ = DSP.hanning(sz)
    elseif window_type == "kaiser"
        vec_ = DSP.kaiser(sz, option1)
    end

    ws = vec_ * vec_'
    wf = fftshift(fft(ifftshift(ws)))

    if Bool(sz & 1)
        wf_crop = wf[div(sz+1,2)-crop:div(sz+1,2)+crop, div(sz+1,2)-crop:div(sz+1,2)+crop]
    else
        wf_crop = wf[div(sz,2)-crop:div(sz+1,2)+crop, div(sz+1,2)-crop:div(sz+1,2)+crop]
    end

    return wf_crop, ws
end

"""
    recon2d_gridrec(p::Array{T, 3}, angles::Array{T, 1}) where {T<:AbstractFloat}

Reconstruct based on [1] gridrec from p

# Args
- p [nangles x nslices x detcount] : (for the moment) detcount should be a power of 2

[1] : Marone, F., Stampanoni, M., 2012. Regridding reconstruction algorithm
    for real-time tomographic imaging. Journal of Synchrotron Radiation
"""
function recon2d_gridrec(p::Array{T, 3}, angles::Array{T, 1}) where {T<:AbstractFloat}
    nangles, nslices, detcount = size(p)

    if detcount % 2 != 0
        println("!TODO: Odd number of detector size can have problems for the moment.")
    end

    #------------- step1 padding
    # we assume that detcount is a power of 2
    # pad 2 times
    sz_pad = Int(2^ceil(log(detcount+1)/log(2)))
    sz_half = sz_pad >> 1 # even number
    padding = sz_pad - detcount
    padding_half = padding >> 1

    UU = zeros(nangles, sz_half)
    VV = zeros(nangles, sz_half)

    tau = cat(Array(1:2:Int(sz_pad/2-1)), Array(Int(sz_pad/2-1):-2:0), dims=1)
    spatial_filter = zeros(sz_pad)
    spatial_filter[1] = 0.25
    spatial_filter[2:2:end] .= -1.0 ./ (pi * tau).^2
    ramlak = fftshift( fft(spatial_filter) )
    # ramlak = fftshift(fft(ifftshift(spatial_filter))) # even function

    UU = zeros(Int, sz_pad)
    VV = zeros(Int, sz_pad)

    sz_window = 9
    sz_window_half = div(sz_window-1, 2)
    kernel, kernel_spatial = make_freq_window(sz_pad, sz_window_half)
    # ws ./= sum(ws)
    kernel_idx = CartesianIndices(kernel)

    # assume that H=detcount, W=detcount
    H, W = detcount, detcount
    rec_img = zeros(H,W,nslices)

    # 1:sz_pad - div(detcount, 2) + 0.5 - 1.0
    offset = div(sz_pad, 2) + 0.5
    dtheta = angles[2]-angles[1]
    halfsz = Int(floor(sz_window/2))

    UU0 = collect(1:sz_pad) .- offset
    one_index = CartesianIndex(1, 1)

    Threads.@threads for slice=1:nslices
        println("slice no. $slice")
        # (optional) add the current slice and the next slice
        # if slice < nslices-1
        #     p_pad[:,1:detcount] .= view(p,:,slice,:) + view(p,:,slice+1,:)
        # else
        #     p_pad[:,1:detcount] .= view(p,:,slice,:)
        # end

        # for multithreading, give up saving memory allocations
        Q = zeros(ComplexF64, sz_pad, sz_pad)
        q = zeros(sz_pad, sz_pad)
        p_pad = zeros(nangles, sz_pad)
        # fill!(Q, 0.0)
        # fill!(p_pad, 0.0)

        for i=1:detcount
            p_pad[:, padding_half+1:end-padding_half] .= view(p,:,slice,:)
        end

        P = fftshift(fft(ifftshift(p_pad), 2));
        P .*= reshape(ramlak, (1,length(ramlak)))
        out_conv = zeros(ComplexF64, sz_pad+sz_window-1, sz_pad+sz_window-1)

        for iang=1:nangles
            cos_ = cos(angles[iang])
            sin_ = sin(angles[iang])

            #---------------------------
            # convolution of P[nangles x detcount] and C
            # continuous H = F * W (4)
            # discrete   Q = P * C (6)
            #---------------------------
            UU .= Int.(round.( UU0 .* cos_ .+ offset ) )
            # UU .-= sz_window_half
            UU[UU.>sz_pad] .-= sz_pad
            UU[UU.<1] .+= sz_pad

            VV .= Int.(round.( UU0 .* sin_ .+ offset ) )
            # VV .-= sz_window_half
            VV[VV.>sz_pad] .-= sz_pad
            VV[VV.<1] .+= sz_pad

            P_idx = CartesianIndex.(VV,UU)
            fill!(out_conv, 0.0)

            for (idx_detcount, ci_P) in enumerate(P_idx)
                Pf = P[iang, idx_detcount]
                # println(ci_P)
                @simd for ci_C in kernel_idx
                     @inbounds out_conv[ci_P+ci_C-one_index] += Pf * kernel[ci_C]
                end
            end

            Q .+= out_conv[halfsz+1:end-halfsz, halfsz+1:end-halfsz]
        end

        Q .*= dtheta / (detcount^2)
        q .= real( fftshift(ifft!(ifftshift(Q))) ./ kernel_spatial ) #./ kernel_weights )
        # q .*= size(Q, 2)
        # print(size(Q))
        # q[q .< 0] .= 0
        # q[q .> 1] .= 1

        rec_img[:,:,slice] .= reverse( q[padding_half+1:end-padding_half, padding_half+1:end-padding_half], dims=1 )
    end

    return rec_img
end

"Gridrec for 2D image"
function recon2d_gridrec(p::Array{T, 2}, angles::Array{T, 1}) where {T<:AbstractFloat}
    p_3d_ = reshape(p, size(p)..., 1)
    p_3d = permutedims(p_3d_, 1, 3, 2)
    rec = recon2d_gridrec(p_3d, angles)
    return rec[:,:,1]
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
module JuliaTutor
    using Crayons

    greet() = print("Hi my name is JuliaTutor, \n
                    Here you can learn interactively through pre-made lessons.")

    exit_attempt(ui) = lowercase(ui[1:min(4,length(ui))]) == "exit"

    function inform(query)
        println( Crayon(foreground = :blue, bold = true),
                 "> ",
                 Crayon(foreground = :white, bold = false),
                 query)
        println( Crayon(foreground = :blue, italics = true),
                 "...Press any key to continue...")
        userinput = readline()
        exit_attempt(userinput) && exit()
    end

    function evaluate(ui, ans)::Tuple{ Bool, Bool }
        #try to run the users code
        syntax_error,incorrect_answer = false, false;
        try
            eval(ui)
        catch
            syntax_error = true;
        end
        incorrect_answer = !(ui == ans)
        return (syntax_error, incorrect_answer)
    end

    function request_read_evaluate(request, answer)
        print( Crayon( foreground = :red, bold = true ),
                 "> ")
        println( Crayon( foreground = :white, bold = true ),
                 request )
        println(" ")
        correct = false
        while !correct
            print( Crayon( foreground = :green ), "> ")
            userinput = readline()
            exit_attempt(userinput) && exit()
            (syntax_error, incorrect_answer) = evaluate(userinput, answer)
            if syntax_error
                println(Crayon( foreground = :white, italics = true ),
                        "Sorry it appears there is a syntax error in your code.")
            elseif incorrect_answer
                println(Crayon( foreground = :white, italics = true ),
                        "Good try but, this isn't the answer we are looking for. Please try again.")
            end
            correct = !syntax_error && !incorrect_answer
        end
    end

    export inform, request_read_evaluate, evaluate

end # module


text:
module QuadraticProgram

using LinearAlgebra, SparseArrays

using MathOptInterface
const MOI = MathOptInterface

import LazyArrays
import IterativeSolvers

import DiffOpt

struct Cache
    lhs::SparseMatrixCSC{Float64, Int}
end

Base.@kwdef struct ForwardReverseCache
    dz::Vector{Float64}
    dλ::Vector{Float64}
    dν::Vector{Float64}
end

MOI.Utilities.@product_of_sets(
    Equalities,
    MOI.EqualTo{T},
)

MOI.Utilities.@product_of_sets(
    Inequalities,
    MOI.LessThan{T},
)

MOI.Utilities.@struct_of_constraints_by_set_types(
    EqualitiesOrInequalities,
    MOI.EqualTo{T},
    MOI.LessThan{T},
)

const Form{T} = MOI.Utilities.GenericModel{
    T,
    MOI.Utilities.ObjectiveContainer{T},
    MOI.Utilities.FreeVariables,
    EqualitiesOrInequalities{T}{
        MOI.Utilities.MatrixOfConstraints{
            T,
            MOI.Utilities.MutableSparseMatrixCSC{
                T,
                Int,
                MOI.Utilities.OneBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{T},
            Equalities{T},
        },
        MOI.Utilities.MatrixOfConstraints{
            T,
            MOI.Utilities.MutableSparseMatrixCSC{
                T,
                Int,
                MOI.Utilities.OneBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{T},
            Inequalities{T},
        },
    },
}

"""
    DiffOpt.QuadraticProgram.Model <: DiffOpt.AbstractModel

Model to differentiate quadratic programs.

For the reverse differentiation, it differentiates the optimal solution `z` and return
product of jacobian matrices (`dz / dQ`, `dz / dq`, etc) with
the backward pass vector `dl / dz`

The method computes the product of
1. jacobian of problem solution `z*` with respect to
    problem parameters set with the [`DiffOpt.ReverseVariablePrimal`](@ref)
2. a backward pass vector `dl / dz`, where `l` can be a loss function

Note that this method *does not returns* the actual jacobians.

For more info refer eqn(7) and eqn(8) of https://arxiv.org/pdf/1703.00443.pdf
"""
mutable struct Model <: DiffOpt.AbstractModel
    # storage for problem data in matrix form
    model::Form{Float64}
    # includes maps from matrix indices to problem data held in `optimizer`
    # also includes KKT matrices
    # also includes the solution
    gradient_cache::Union{Nothing,Cache}

    # caches for sensitivity output
    # result from solving KKT/residualmap linear systems
    # this allows keeping the same `gradient_cache`
    # if only sensitivy input changes
    forw_grad_cache::Union{Nothing,ForwardReverseCache}
    back_grad_cache::Union{Nothing,ForwardReverseCache}

    # sensitivity input cache using MOI like sparse format
    input_cache::DiffOpt.InputCache

    x::Vector{Float64} # Primal
    λ::Vector{Float64} # Dual of inequalities
    ν::Vector{Float64} # Dual of equalities
end
function Model()
    return Model(Form{Float64}(), nothing, nothing, nothing, DiffOpt.InputCache(), Float64[], Float64[], Float64[])
end

function MOI.is_empty(model::Model)
    return MOI.is_empty(model.model)
end

function MOI.empty!(model::Model)
    MOI.empty!(model.model)
    model.gradient_cache = nothing
    model.forw_grad_cache = nothing
    model.back_grad_cache = nothing
    empty!(model.input_cache)
    empty!(model.x)
    empty!(model.λ)
    empty!(model.ν)
    return
end

const EQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}}
const LE = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}

function MOI.set(model::Model, ::MOI.ConstraintPrimalStart, ci::MOI.ConstraintIndex, value)
    MOI.throw_if_not_valid(model, ci)
    return # Ignored
end

# We want to stay consistent with the variable `ν` defined in (3) of
# Left hand side of eq. (6) in https://arxiv.org/pdf/1703.00443.pdf
# However, in eq. (6), they put it in the lagrangian as
# `+ ν ⋅ (Az - b)`
# while in MOI, we put it as
# `- ν ⋅ (Az - b)`
# so the we should reverse the sign if we want to use the same equations
# as in the paper.
function MOI.set(model::Model, ::MOI.ConstraintDualStart, ci::EQ, value)
    MOI.throw_if_not_valid(model, ci)
    DiffOpt._enlarge_set(model.ν, MOI.Utilities.rows(_equalities(model), ci), -value)
end

function MOI.set(model::Model, ::MOI.ConstraintDualStart, ci::LE, value)
    MOI.throw_if_not_valid(model, ci)
    DiffOpt._enlarge_set(model.λ, MOI.Utilities.rows(_inequalities(model), ci), -value)
end

function _gradient_cache(model::Model)
    if model.gradient_cache !== nothing
        return model.gradient_cache
    end
    A = convert(SparseMatrixCSC{Float64, Int}, _equalities(model).coefficients)
    b = _equalities(model).constants.upper
    G = convert(SparseMatrixCSC{Float64, Int}, _inequalities(model).coefficients)
    h = _inequalities(model).constants.upper

    nz = size(A, 2)
    # TODO ideally, the objective should be stored in matrix form in `model.model`
    objective_function = MOI.get(model.model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    sparse_array_obj = DiffOpt.sparse_array_representation(objective_function, nz)
    Q = sparse_array_obj.quadratic_terms
    q = sparse_array_obj.affine_terms

    LHS = create_LHS_matrix(model.x, model.λ, Q, G, h, A)

    model.gradient_cache = Cache(LHS)

    return model.gradient_cache
end

function _equalities(model::Model)
    return MOI.Utilities.constraints(
        model.model.constraints,
        MOI.ScalarAffineFunction{Float64},
        MOI.EqualTo{Float64},
    )
end
function _inequalities(model::Model)
    return MOI.Utilities.constraints(
        model.model.constraints,
        MOI.ScalarAffineFunction{Float64},
        MOI.LessThan{Float64},
    )
end

# TODO: create test functions for the methods

# """
#     Left hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_LHS_matrix(z, λ, Q, G, h, A=nothing)
#     if A == nothing || size(A)[1] == 0
#         return [Q                G';
#                 Diagonal(λ) * G    Diagonal(G * z - h)]
#     else
#         @assert size(A)[2] == size(G)[2]
#         p, n = size(A)
#         m    = size(G)[1]
#         return [Q                  G'                    A';
#                 Diagonal(λ) * G    Diagonal(G * z - h)   zeros(m, p);
#                 A                  zeros(p, m)           zeros(p, p)]
#     end
# end

"""
    create_LHS_matrix(z, λ, Q, G, h, A=nothing)

Inverse matrix specified on RHS of eqn(7) in https://arxiv.org/pdf/1703.00443.pdf

Helper method while calling `reverse_differentiate!`
"""
function create_LHS_matrix(z, λ, Q, G, h, A=nothing)::AbstractMatrix{Float64}
    if (A === nothing || size(A)[1] == 0) && (G === nothing || size(G)[1] == 0)
        return Q
    elseif A === nothing || size(A)[1] == 0
        return [Q         G' * Diagonal(λ);
                G         Diagonal(G * z - h)]
    elseif G === nothing || size(G)[1] == 0
        p, n = size(A)
        return [Q         A';
                A         spzeros(p, p)]
    else
        p, n  = size(A)
        m, n2 = size(G)
        if n != n2
            throw(DimensionError("Sizes of $A and $G do not match"))
        end
        return [Q         G' * Diagonal(λ)       A';
                G         Diagonal(G * z - h)    spzeros(m, p);
                A         spzeros(p, m)          spzeros(p, p)]
    end
end
# TODO: this is the transpose, check back for usage

# """
#     Right hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_RHS_matrix(z, dQ, dq, λ, dG, dh, ν=nothing, dA=nothing, db=nothing)
#     if dA == nothing || size(dA)[1] == 0
#         return -[dQ * z + dq + dG' * λ      ;
#                  Diagonal(λ) * (dG * z - dh)]
#     else
#         return -[dQ * z + dq + dG' * λ + dA' * ν;
#                  Diagonal(λ) * (dG * z - dh)    ;
#                  dA * z - db                    ]
#     end
# end

function MOI.get(model::Model, ::DiffOpt.ForwardVariablePrimal, vi::MOI.VariableIndex)
    return model.forw_grad_cache.dz[vi.value]
end

function DiffOpt._get_db(model::Model, ci::LE)
    i = ci.value
    # dh = -Diagonal(λ) * dλ
    return model.λ[i] * model.back_grad_cache.dλ[i]
end
function DiffOpt._get_db(model::Model, ci::EQ)
    return model.back_grad_cache.dν[ci.value]
end

function DiffOpt.reverse_differentiate!(model::Model)
    gradient_cache = _gradient_cache(model)
    LHS = gradient_cache.lhs

    neq = length(model.ν)
    nineq = length(model.λ)

    dl_dz = zeros(length(model.x))
    for (vi, value) in model.input_cache.dx
        dl_dz[vi.value] = value
    end

    RHS = [dl_dz; zeros(nineq + neq)]

    nv = length(model.x)
    Q = view(LHS, 1:nv, 1:nv)
    partial_grads = if norm(Q) ≈ 0
        -IterativeSolvers.lsqr(LHS, RHS)
    else
        -LHS \ RHS
    end

    dz = partial_grads[1:nv]
    dλ = partial_grads[nv+1:nv+nineq]
    dν = partial_grads[nv+nineq+1:end]

    model.back_grad_cache = ForwardReverseCache(dz, dλ, dν)
    return
    # dQ = 0.5 * (dz * z' + z * dz')
    # dq = dz
    # dG = Diagonal(λ) * (dλ * z' + λ * dz') # was: Diagonal(λ) * dλ * z' - λ * dz')
    # dh = -Diagonal(λ) * dλ
    # dA = dν * z'+ ν * dz' # was: dν * z' - ν * dz'
    # db = -dν
    # todo, check MOI signs for dA and dG
end

_linsolve(A, b) = A \ b
# See https://github.com/JuliaLang/julia/issues/32668
_linsolve(A, b::SparseVector) = A \ Vector(b)

# Just a hack that will be removed once we use `MOI.Utilities.MatrixOfConstraints`
struct _QPSets end
MOI.Utilities.rows(::_QPSets, ci::MOI.ConstraintIndex) = ci.value

function DiffOpt.forward_differentiate!(model::Model)
    gradient_cache = _gradient_cache(model)
    LHS = gradient_cache.lhs

    z = model.x
    nv = length(z)

    objective_function = DiffOpt._convert(MOI.ScalarQuadraticFunction{Float64}, model.input_cache.objective)
    sparse_array_obj = DiffOpt.sparse_array_representation(objective_function, nv)
    dQ = sparse_array_obj.quadratic_terms
    dq = sparse_array_obj.affine_terms

    # The user sets the constraint function in the sense `func`-in-`set` while
    # `db` and `dh` corresponds to the tangents of the set constants. Therefore,
    # we should multiply the constant by `-1`. For `GreaterThan`, we needed to
    # multiply by `-1` to transform it to `LessThan` so it cancels out.
    db = zero(model.ν)
    DiffOpt._fill(isequal(MOI.EqualTo{Float64}), (::Type{MOI.EqualTo{Float64}}) -> true, gradient_cache, model.input_cache, _QPSets(), db)
    dh = zero(model.λ)
    DiffOpt._fill(!isequal(MOI.EqualTo{Float64}), !isequal(MOI.GreaterThan{Float64}), gradient_cache, model.input_cache, _QPSets(), dh)

    dAi = zeros(Int, 0)
    dAj = zeros(Int, 0)
    dAv = zeros(Float64, 0)
    DiffOpt._fill(isequal(MOI.EqualTo{Float64}), isequal(MOI.GreaterThan{Float64}), gradient_cache, model.input_cache, _QPSets(), dAi, dAj, dAv)
    dA = sparse(dAi, dAj, dAv, length(model.ν), nv)

    dGi = zeros(Int, 0)
    dGj = zeros(Int, 0)
    dGv = zeros(Float64, 0)
    DiffOpt._fill(!isequal(MOI.EqualTo{Float64}), isequal(MOI.GreaterThan{Float64}), gradient_cache, model.input_cache, _QPSets(), dGi, dGj, dGv)
    dG = sparse(dGi, dGj, dGv, length(model.λ), nv)

    λ = model.λ
    ν = model.ν
    RHS = [
        dQ * z + dq + dG' * λ + dA' * ν
        λ .* (dG * z) - λ .* dh
        dA * z - db
    ]

    Q = view(LHS, 1:nv, 1:nv)
    partial_grads = if norm(Q) ≈ 0
        -IterativeSolvers.lsqr(LHS', RHS)
    else
        -_linsolve(LHS', RHS)
    end

    dz = partial_grads[1:nv]
    dλ = partial_grads[nv+1:nv+length(λ)]
    dν = partial_grads[nv+length(λ)+1:end]

    model.forw_grad_cache = ForwardReverseCache(dz, dλ, dν)
    return
end

function MOI.get(model::Model, ::DiffOpt.ReverseObjectiveFunction)
    ∇z = model.back_grad_cache.dz
    z = model.x
    # `∇z * z' + z * ∇z'` doesn't work, see
    # https://github.com/JuliaArrays/LazyArrays.jl/issues/178
    dQ = LazyArrays.@~ (∇z .* z' + z .* ∇z') / 2.0
    return DiffOpt.MatrixScalarQuadraticFunction(
        DiffOpt.VectorScalarAffineFunction(∇z, 0.0),
        dQ,
    )
end

# quadratic matrix indexes are split by type either == or <=
function DiffOpt._get_dA(model::Model, ci::EQ)
    i = ci.value
    dz = model.back_grad_cache.dz
    dν = model.back_grad_cache.dν
    return DiffOpt.lazy_combination(+, dν[i], model.x, model.ν[i], dz)
end
function DiffOpt._get_dA(model::Model, ci::LE)
    i = ci.value
    dz = model.back_grad_cache.dz
    dλ = model.back_grad_cache.dλ
    l = model.λ[i]
    return DiffOpt.lazy_combination(+, l * dλ[i], model.x, l, dz)
end

end


text:
using Test, LinearAlgebra, OrderedCollections
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Greens interpolant test" begin
    # construct interior solution
    @testset "Helmholtz" begin
        pde  = Helmholtz(dim=3,k=1+im)
        xout = SVector(3,3,3)
        u    = (x)   -> SingleLayerKernel(pde)(xout,x)
        dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
        Geometry.clear!()
        Ω,M  = WaveProp.IO.gmsh_sphere(dim=2,h=0.2)
        Γ    = boundary(Ω)
        mesh = NystromMesh(view(M,Γ);order=2)
        γ₀u   = γ₀(u,mesh)
        γ₁u   = γ₁(dudn,mesh)
        S     = SingleLayerOperator(pde,mesh)
        D     = DoubleLayerOperator(pde,mesh)
        e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
        Sdim  = Nystrom.assemble_dim(S)
        Ddim  = Nystrom.assemble_dim(D)
        e1    = WaveProp.Nystrom.error_interior_green_identity(Sdim,Ddim,γ₀u,γ₁u)
        @test 10*norm(e1,Inf) < norm(e0,Inf)
        S,D   = Nystrom.single_doublelayer_dim(pde,mesh)
        e2    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
        @test e1 ≈ e2
        K     = AdjointDoubleLayerOperator(pde,mesh)
        H     = HyperSingularOperator(pde,mesh)
        e0    = WaveProp.Nystrom.error_interior_derivative_green_identity(K,H,γ₀u,γ₁u)
        Kdim  = Nystrom.assemble_dim(K)
        Hdim  = Nystrom.assemble_dim(H)
        e1    = WaveProp.Nystrom.error_interior_derivative_green_identity(Kdim,Hdim,γ₀u,γ₁u)
        @test 10*norm(e1,Inf) < norm(e0,Inf)
        @test norm(e1,Inf) < 1e-4
    end
    @testset "Elastostatic" begin
        pde  = Elastostatic(dim=3,μ=1,λ=2.)
        xout = SVector(3,3,3)
        c    = ones(Nystrom.default_density_eltype(pde))
        u    = (x)   -> SingleLayerKernel(pde)(xout,x)*c
        dudn = (x,n) -> transpose(DoubleLayerKernel(pde)(xout,x,n))*c
        Geometry.clear!()
        Ω,M  = WaveProp.IO.gmsh_sphere(dim=2,h=0.2)
        Γ    = boundary(Ω)
        mesh = NystromMesh(view(M,Γ);order=2)
        γ₀u   = γ₀(u,mesh)
        γ₁u   = γ₁(dudn,mesh)
        S     = SingleLayerOperator(pde,mesh)
        D     = DoubleLayerOperator(pde,mesh)
        e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
        Sdim  = Nystrom.assemble_dim(S)
        Ddim  = Nystrom.assemble_dim(D)
        e1    = WaveProp.Nystrom.error_interior_green_identity(Sdim,Ddim,γ₀u,γ₁u)
        @test 10*norm(norm.(e1,Inf),Inf) < norm(norm.(e0,Inf),Inf)
        S,D   = Nystrom.single_doublelayer_dim(pde,mesh)
        e2    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
        @test e1 ≈ e2
        K     = AdjointDoubleLayerOperator(pde,mesh)
        H     = HyperSingularOperator(pde,mesh)
        e0    = WaveProp.Nystrom.error_interior_derivative_green_identity(K,H,γ₀u,γ₁u)
        Kdim  = Nystrom.assemble_dim(K)
        Hdim  = Nystrom.assemble_dim(H)
        e1    = WaveProp.Nystrom.error_interior_derivative_green_identity(Kdim,Hdim,γ₀u,γ₁u)
        @test 10*norm(norm(e1,Inf),Inf) < norm(norm(e0,Inf),Inf)
        @test norm(norm(e1,Inf),Inf) < 1e-4
    end
    @testset "Maxwell" begin
        pde  = Maxwell(dim=3,k=1.0)
        xout = SVector(3,3,3)
        c    = ones(Nystrom.default_density_eltype(pde))
        u    = (x)   -> SingleLayerKernel(pde)(xout,x)*c
        dudn = (x,n) -> transpose(DoubleLayerKernel(pde)(xout,x,n))*c
        Geometry.clear!()
        Ω,M  = WaveProp.IO.gmsh_sphere(dim=2,h=0.2)
        Γ    = boundary(Ω)
        mesh = NystromMesh(view(M,Γ);order=2)
        γ₀u   = γ₀(u,mesh)
        normals = mesh.qnormals
        γ₁u   = γ₁(dudn,mesh)
        S     = SingleLayerOperator(pde,mesh)
        D     = DoubleLayerOperator(pde,mesh)
        e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
        Sdim  = Nystrom.assemble_dim(S)
        Ddim  = Nystrom.assemble_dim(D)
        e1    = WaveProp.Nystrom.error_interior_green_identity(Sdim,Ddim,γ₀u,γ₁u)
        @test 10*norm(norm.(e1,Inf),Inf) < norm(norm.(e0,Inf),Inf)
        S,D   = Nystrom.single_doublelayer_dim(pde,mesh)
        e2    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
        @test e1 ≈ e2
    end
end

@testset "Multiple obstacles" begin
    # construct interior solution
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    Geometry.clear!()
    geo1  = Geometry.Kite()
    geo2  = Circle(center=(10,0))
    Ω     = Domain([geo1,geo2])
    M     = meshgen(Ω,h=0.1)
    Γ     = boundary(Ω)
    # generate a Nystrom mesh with Gauss-Legendre quadrature
    qrule   = GaussLegendre(5)
    e2qrule = OrderedDict(E=>qrule for E in etypes(M))
    mesh    = NystromMesh(M,Γ,e2qrule)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    Smk   = Nystrom.assemble_dim(S)
    Dmk   = Nystrom.assemble_dim(D)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    e1 = WaveProp.Nystrom.error_interior_green_identity(Smk,Dmk,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
    @test sum(qweights(mesh[Γ[2]])) ≈ 2π
end

@testset "Extract domain mesh" begin
    # construct interior solution
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    Geometry.clear!()
    geo1  = Geometry.Kite()
    geo2  = Circle(center=(10,0))
    Ω     = Domain([geo1,geo2])
    M   = meshgen(Ω,h=0.1)
    Γ     = boundary(Ω)
    # generate a Nystrom mesh with Gauss-Legendre quadrature
    qrule   = GaussLegendre(5)
    e2qrule = OrderedDict(E=>qrule for E in etypes(M))
    Γ_mesh    = NystromMesh(M,Γ,e2qrule)
    Γ1_mesh   = Γ_mesh[Γ[1]]
    Γ2_mesh   = Γ_mesh[Γ[2]]
    γ₀u   = γ₀(u,Γ2_mesh)
    γ₁u   = γ₁(dudn,Γ2_mesh)
    S     = SingleLayerOperator(pde,Γ2_mesh)
    D     = DoubleLayerOperator(pde,Γ2_mesh)
    Smk   = Nystrom.assemble_dim(S)
    Dmk   = Nystrom.assemble_dim(D)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    e1 = WaveProp.Nystrom.error_interior_green_identity(Smk,Dmk,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module Display
using Requires, Crayons, Dates, ANSIColoredPrinters
using ..Structures
using ..Calculations
using ..Jobs
using ..Client
using ..DFControl: Band, TimingData
using ..Utils

include("printing.jl")

const dfprintln = println
const dfprint = print

function __init__()
    @require Glimpse = "f6e19d58-12a4-5927-8606-ac30a9ce9b69" include("glimpse.jl")
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plotting.jl")
    if isdefined(Base, :active_repl)
        Display.include(joinpath(@__DIR__, "base.jl"))
    elseif isdefined(Main, :PlutoRunner)
        Display.include(joinpath(@__DIR__, "pluto.jl"))

    else
        @require IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" include("base.jl")
    end
end
end


text:
# discovering binary CUDA dependencies

using Pkg, Pkg.Artifacts
using Libdl

## global state

const __version = Ref{VersionNumber}()

"""
    version()

Returns the version of the CUDA toolkit in use.
"""
version() = @after_init(__version[])

"""
    release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
release() = @after_init(VersionNumber(__version[].major, __version[].minor))

const __libcublas = Ref{String}()
const __libcusparse = Ref{String}()
const __libcusolver = Ref{String}()
const __libcufft = Ref{String}()
const __libcurand = Ref{String}()
const __libcudnn = Ref{Union{Nothing,String}}(nothing)
const __libcutensor = Ref{Union{Nothing,String}}(nothing)

libcublas() = @after_init(__libcublas[])
libcusparse() = @after_init(__libcusparse[])
libcusolver() = @after_init(__libcusolver[])
libcufft() = @after_init(__libcufft[])
libcurand() = @after_init(__libcurand[])
function libcudnn()
    @after_init begin
        @assert has_cudnn() "This functionality is unavailabe as CUDNN is missing."
        __libcudnn[]
    end
end
function libcutensor()
    @after_init begin
        @assert has_cutensor() "This functionality is unavailabe as CUTENSOR is missing."
        __libcutensor[]
    end
end

export has_cudnn, has_cutensor
has_cudnn() = @after_init(__libcudnn[]) !== nothing
has_cutensor() = @after_init(__libcutensor[]) !== nothing

## discovery

# CUDA

# NOTE: we don't use autogenerated JLLs, because we have multiple artifacts and need to
#       decide at run time (i.e. not via package dependencies) which one to use.
const cuda_artifacts = Dict(
    v"10.2" => ()->artifact"CUDA10.2",
    v"10.1" => ()->artifact"CUDA10.1",
    v"10.0" => ()->artifact"CUDA10.0",
    v"9.2"  => ()->artifact"CUDA9.2",
    v"9.0"  => ()->artifact"CUDA9.0",
)

# utilities to look up stuff in the artifact (at known locations, so not using CUDAapi)
get_binary(artifact_dir, name) = joinpath(artifact_dir, "bin", Sys.iswindows() ? "$name.exe" : name)
function get_library(artifact_dir, name)
    filename = if Sys.iswindows()
        "$name.dll"
    elseif Sys.isapple()
        "lib$name.dylib"
    else
        "lib$name.so"
    end
    joinpath(artifact_dir, Sys.iswindows() ? "bin" : "lib", filename)
end

function use_artifact_cuda()
    @debug "Trying to use artifacts..."

    # select compatible artifacts
    if haskey(ENV, "JULIA_CUDA_VERSION")
        wanted_version = VersionNumber(ENV["JULIA_CUDA_VERSION"])
        filter!(((version,artifact),) -> version == wanted_version, cuda_artifacts)
    else
        driver_version = CUDAdrv.release()
        filter!(((version,artifact),) -> version <= driver_version, cuda_artifacts)
    end

    # download and install
    artifact = nothing
    for release in sort(collect(keys(cuda_artifacts)); rev=true)
        try
            artifact = (release=release, dir=cuda_artifacts[release]())
            break
        catch
        end
    end
    if artifact == nothing
        @debug "Could not find a compatible artifact."
        return false
    end

    nvdisasm = get_binary(artifact.dir, "nvdisasm")
    @assert isfile(nvdisasm)
    __version[] = parse_toolkit_version(nvdisasm)

    # discover libraries
    for name in  ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CuArrays, Symbol("__lib$name"))

        # on Windows, the library name is version dependent
        if Sys.iswindows()
            name *= if artifact.release >= v"10.1"
                "64_$(artifact.release.major)"
            else
                "64_$(artifact.release.major)$(artifact.release.minor)"
            end
        end

        handle[] = get_library(artifact.dir, name)
        Libdl.dlopen(handle[])
    end

    @debug "Using CUDA $(__version[]) from an artifact at $(artifact.dir)"
    use_artifact_cudnn(artifact.release)
    use_artifact_cutensor(artifact.release)
    return true
end

function use_local_cuda()
    @debug "Trying to use local installation..."

    cuda_dirs = find_toolkit()

    tool = find_cuda_binary("nvdisasm", cuda_dirs)
    if tool === nothing
        @debug "Could not find nvdisasm"
        return false
    end

    cuda_version = parse_toolkit_version(tool)
    __version[] = cuda_version

    # discover libraries
    for name in  ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CuArrays, Symbol("__lib$name"))

        path = find_cuda_library(name, cuda_dirs, [cuda_version])
        if path === nothing
            @debug "Could not find $name"
            return false
        end
        handle[] = path
    end

    @debug "Found local CUDA $(cuda_version) at $(join(cuda_dirs, ", "))"
    use_local_cudnn(cuda_dirs)
    use_local_cutensor(cuda_dirs)
    return true
end

# CUDNN

const cudnn_artifacts = Dict(
    v"10.2" => ()->artifact"CUDNN+CUDA10.2",
    v"10.1" => ()->artifact"CUDNN+CUDA10.1",
    v"10.0" => ()->artifact"CUDNN+CUDA10.0",
    v"9.2"  => ()->artifact"CUDNN+CUDA9.2",
    v"9.0"  => ()->artifact"CUDNN+CUDA9.0",
)

function use_artifact_cudnn(release)
    artifact_dir = try
        cudnn_artifacts[release]()
    catch ex
        @debug "Could not use CUDNN from artifacts" exception=(ex, catch_backtrace())
        return false
    end

    __libcudnn[] = get_library(artifact_dir, Sys.iswindows() ? "cudnn64_7" : "cudnn")
    Libdl.dlopen(__libcudnn[])
    @debug "Using CUDNN from an artifact at $(artifact_dir)"
    return true
end

function use_local_cudnn(cuda_dirs)
    path = find_cuda_library("cudnn", cuda_dirs, [v"7"])
    path === nothing && return false

    __libcudnn[] = path
    @debug "Using local CUDNN at $(path)"
    return true
end

# CUTENSOR

const cutensor_artifacts = Dict(
    v"10.2" => ()->artifact"CUTENSOR+CUDA10.2",
    v"10.1" => ()->artifact"CUTENSOR+CUDA10.1",
)

function use_artifact_cutensor(release)
    artifact_dir = try
        cutensor_artifacts[release]()
    catch ex
        @debug "Could not use CUTENSOR from artifacts" exception=(ex, catch_backtrace())
        return false
    end

    __libcutensor[] = get_library(artifact_dir, "cutensor")
    Libdl.dlopen(__libcutensor[])
    @debug "Using CUTENSOR from an artifact at $(artifact_dir)"
    return true
end

function use_local_cutensor(cuda_dirs)
    path = find_cuda_library("cutensor", cuda_dirs, [v"1"])
    path === nothing && return false

    __libcutensor[] = path
    @debug "Using local CUTENSOR at $(path)"
    return true
end

function __configure_dependencies__(show_reason::Bool)
    found = false

    if parse(Bool, get(ENV, "JULIA_CUDA_USE_BINARYBUILDER", "true"))
        found = use_artifact_cuda()
    end

    if !found
        found = use_local_cuda()
    end

    if !found
        show_reason && @error "Could not find a suitable CUDA installation"
    end

    return found
end


text:
using ..Inputs: QuantumESPRESSOInput
using ..Inputs.PWscf: PWInput

struct PhInput <: QuantumESPRESSOInput
    title_line::String
    inputph::PhNamelist
    q_points::Union{Nothing,QPointsCard}
end # struct PhInput
PhInput(inputph::PhNamelist, qpts::QPointsCard) = PhInput(inputph.prefix, inputph, qpts)
PhInput(inputph::PhNamelist) = PhInput(inputph.prefix, inputph, nothing)
PhInput() = PhInput(PhNamelist().prefix, PhNamelist(), nothing)

struct Q2rInput <: QuantumESPRESSOInput
    input::Q2rNamelist
end # struct Q2rInput
Q2rInput() = Q2rInput(Q2rNamelist())

struct MatdynInput <: QuantumESPRESSOInput
    input::MatdynNamelist
    q_points::Union{Nothing,QPointsCard}
end # struct MatdynInput
MatdynInput(input) = MatdynInput(input, nothing)
MatdynInput() = MatdynInput(MatdynNamelist(), nothing)

struct DynmatInput <: QuantumESPRESSOInput
    input::DynmatNamelist
end # struct DynmatInput
DynmatInput() = DynmatInput(DynmatNamelist())

function (x::VerbositySetter)(template::PhInput)
    @set! template.inputph.verbosity = x.v
    return template
end

"""
    relayinfo(from::PWInput, to::PhInput)

Relay shared information from a `PWInput` to a `PhInput`.

A `PWInput` before a `PhInput` has the information of `outdir` and `prefix`. They must keep the same in a
phonon calculation.
"""
function relayinfo(pw::PWInput, ph::PhInput)
    @set! ph.inputph.outdir = pw.control.outdir
    @set! ph.inputph.prefix = pw.control.prefix
    return ph
end # function relayinfo
"""
    relayinfo(from::PhInput, to::Q2rInput)

Relay shared information from a `PhInput` to a `Q2rInput`.

A `PhInput` before a `Q2rInput` has the information of `fildyn`. It must keep the same in a q2r calculation.
"""
function relayinfo(ph::PhInput, q2r::Q2rInput)
    @set! q2r.input.fildyn = ph.inputph.fildyn
    return q2r
end # function relayinfo
"""
    relayinfo(from::Q2rInput, to::MatdynInput)

Relay shared information from a `Q2rInput` to a `MatdynInput`.

A `Q2rInput` before a `MatdynInput` has the information of `fildyn`, `flfrc` and `loto_2d`. They must keep the same
in a matdyn calculation.
"""
function relayinfo(q2r::Q2rInput, matdyn::MatdynInput)
    @set! matdyn.input.flfrc = q2r.input.flfrc
    @set! matdyn.input.loto_2d = q2r.input.loto_2d
    return matdyn
end # function relayinfo
function relayinfo(ph::PhInput, matdyn::MatdynInput)
    @set! matdyn.input.amass = ph.inputph.amass
    @set! matdyn.input.q_in_band_form = ph.inputph.q_in_band_form
    return matdyn
end # function relayinfo
"""
    relayinfo(from::PhInput, to::DynmatInput)

Relay shared information from a `PhInput` to a `DynmatInput`.

A `PhInput` before a `DynmatInput` has the information of `asr`, `fildyn` and `amass`. They must keep the same
in a dynmat calculation.
"""
function relayinfo(ph::PhInput, dynmat::DynmatInput)
    # @set! dynmat.input.asr = ph.inputph.asr  # TODO
    @set! dynmat.input.fildyn = ph.inputph.fildyn
    @set! dynmat.input.amass = ph.inputph.amass
    return dynmat
end # function relayinfo


text:
immutable IDMMOBILBehavior <: BehaviorModel
	p_idm::IDMParam
	p_mobil::MOBILParam
	idx::Int
end
==(a::IDMMOBILBehavior,b::IDMMOBILBehavior) = (a.p_idm==b.p_idm) && (a.p_mobil==b.p_mobil)
Base.hash(a::IDMMOBILBehavior,h::UInt64=zero(UInt64)) = hash(a.p_idm,hash(a.p_mobil,h))

+(a::IDMMOBILBehavior, b::IDMMOBILBehavior) = IDMMOBILBehavior(a.p_idm+b.p_idm, a.p_mobil+b.p_mobil, 0)
-(a::IDMMOBILBehavior, b::IDMMOBILBehavior) = a+(-1.0*b)
*(a::Real, b::IDMMOBILBehavior) = IDMMOBILBehavior(a*b.p_idm, a*b.p_mobil, 0)
.*(a::Real, b::IDMMOBILBehavior) = a*b
.*(b::IDMMOBILBehavior, v::Vector{Float64}) = IDMMOBILBehavior(b.p_idm.*v[1:6], b.p_mobil.*v[7:9], 0)
^(b::IDMMOBILBehavior, p::Integer) = IDMMOBILBehavior(b.p_idm.^p, b.p_mobil.^p, 0)
.^(b::IDMMOBILBehavior, p::Integer) = IDMMOBILBehavior(b.p_idm.^p, b.p_mobil.^p, 0)
.-{B<:BehaviorModel}(v::Vector{B}, b::IDMMOBILBehavior) = B[v[i]-b for i in 1:length(v)]
/(b::IDMMOBILBehavior, f::Real) = 1/f*b
sqrt(b::IDMMOBILBehavior) = IDMMOBILBehavior(sqrt(b.p_idm), sqrt(b.p_mobil), 0)
zero(::Type{IDMMOBILBehavior}) = IDMMOBILBehavior(zero(IDMParam), zero(MOBILParam), 0)
nan(::Type{IDMMOBILBehavior}) = IDMMOBILBehavior(nan(IDMParam), nan(MOBILParam), 0)
abs(b::IDMMOBILBehavior) = IDMMOBILBehavior(abs(b.p_idm), abs(b.p_mobil), 0)
max(a::IDMMOBILBehavior, b::IDMMOBILBehavior) = IDMMOBILBehavior(max(a.p_idm, b.p_idm), max(a.p_mobil, b.p_mobil), 0)

function IDMMOBILBehavior(s::AbstractString,v0::Float64,s0::Float64,idx::Int)
	return IDMMOBILBehavior(IDMParam(s,v0,s0), MOBILParam(s), idx)
end

typical_velocity(b::IDMMOBILBehavior) = b.p_idm.v0

gen_accel(bmodel::BehaviorModel, dmodel::AbstractMLDynamicsModel, s::MLState, neighborhood::Array{Int,1}, idx::Int, rng::AbstractRNG) = error("Uninstantiated Behavior Model")

gen_lane_change(bmodel::BehaviorModel, dmodel::AbstractMLDynamicsModel, s::MLState, neighborhood::Array{Int,1}, idx::Int, rng::AbstractRNG) = error("Uninstantiated Behavior Model")

function gen_accel(bmodel::IDMMOBILBehavior, dmodel::AbstractMLDynamicsModel, s::MLState, neighborhood::Array{Int,1}, idx::Int, rng::AbstractRNG)

    d = accel_dist(bmodel, dmodel, s, neighborhood, idx)
    acc = rand(rng, d)

    return max(acc, -dmodel.phys_param.brake_limit)
end

function accel_dist(bmodel::IDMMOBILBehavior, dmodel::AbstractMLDynamicsModel, s::MLState, neighborhood::Array{Int,1}, idx::Int)
	pp = dmodel.phys_param
	dt = pp.dt
	car = s.cars[idx]
	vel = car.vel

	dv, ds = get_dv_ds(pp,s,neighborhood,idx,2)

	dvel = get_idm_dv(bmodel.p_idm,dt,vel,dv,ds) #call idm model
    acc = dvel/dt
#=
    if acc > (1.01*bmodel.p_idm.a)
      print(acc)
	  print(" || ")
	  println(bmodel.p_idm.a)
	  warn("acc is bigger than the settle value")
	end
=#
    @assert acc <= 1.01*bmodel.p_idm.a

    @if_debug if ds < 0.0
        Gallium.@enter accel_dist(bmodel, dmodel, s, neighborhood, idx)
    end
    #@assert ds >= 0.0 # can get rid of this
    if neighborhood[2] > 0
        @assert abs(s.cars[neighborhood[2]].vel - (vel-dv)) < 0.0001
    end

    max_safe = max_safe_acc(ds, vel, vel-dv, pp.brake_limit, dt)
	#max_safe = nullable_max_safe_acc(ds, vel, vel-dv, pp.brake_limit, dt)

    # if T, a, and b are small the idm may command something faster than the max_safe
    # @assert max_safe >= acc "max_safe=$max_safe; acc=$acc"
    if acc > max_safe
        acc = max_safe
    end

    lower_bound = min(-1e-5, max(-bmodel.p_idm.a/2, -pp.brake_limit-acc))
    upper_bound = min(bmodel.p_idm.a/2, max(max_safe-acc, 1e-5))

	if lower_bound > upper_bound
		lower_bound = -0.1
		upper_bound = 0.1
	end

	if acc > 8
		acc = 8
	end
	if acc<-8
		acc = -8
	end

    return TriangularDist(acc+lower_bound, acc+upper_bound, acc)
end

function max_accel(bmodel::IDMMOBILBehavior)
    return 1.5*bmodel.p_idm.a
end

function gen_lane_change(bmodel::IDMMOBILBehavior, dmodel::AbstractMLDynamicsModel, s::MLState, neighborhood::Array{Int,1}, idx::Int, rng::AbstractRNG)
    #print("LC2")
	pp = dmodel.phys_param
	dt = pp.dt
	car = s.cars[idx]
	lane_change = car.lane_change #this is a velocity in the y direction in LANES PER SECOND
	lane_ = car.y

	if mod(lane_-0.5,1) == 0. #in between lanes
        @assert lane_change != 0
		#println("in")
		return lane_change
	end

	#sample normally
	lanechange_::Int = get_mobil_lane_change(bmodel, pp, s, neighborhood, idx, rng)
	#gives +1, -1 or 0

    lanechange = lanechange_
    #println(lanechange * dmodel.lane_change_rate)
	return lanechange * dmodel.lane_change_rate
end

#############################################################################
#TODO How to make compatible with MOBIL?

type AvoidModel <: BehaviorModel
	jerk::Bool
end
AvoidModel() = AvoidModel(false)
# NOTE: JerkModel might be better modeled by a /very/ aggressive IDM car
JerkModel() = AvoidModel(true)
==(a::AvoidModel, b::AvoidModel) = (a.jerk == b.jerk)
Base.hash(a::AvoidModel, h::UInt64=zero(UInt64)) = hash(a.jerk,h)

function closest_car(dmodel::IDMMOBILModel, s::MLState, nbhd::Array{Int,1}, idx::Int, lookahead_only::Bool)
	x = s.cars[idx].x
    y = s.cars[idx].y
	dy = dmodel.phys_param.w_lane

	closest = 0
	min_dist = Inf

	for i = 1:3 #front cars
		if nbhd[i] == 0
			continue
		end
		dist = norm([x - car.x; (y - car.y)*dy])
		if dist < min_dist
			min_dist = dist
			closest = i
		end
	end

	if lookahead_only
		return closest
	end

	for i = 4:6 #rear cars
		if nbhd[i] == 0
			continue
		end
		car = s.cars[nbhd[i]]
		dists = norm([x - car.x; (y - car.y)*dy])
		if dist < min_dist
			min_dist = dist
			closest = i
		end
	end

	return closest

end

function gen_accel(bmodel::AvoidModel, dmodel::IDMMOBILModel, s::MLState, neighborhood::Array{Int,1}, idx::Int, rng::AbstractRNG)
	closest_car = closest_car(dmodel,s,neighborhood,idx,bmodel.jerk)
	if closest_car == 0
		return 0.
	end

	x = s.cars[idx].x
	cc = s.cars[closest_car]

	if p.jerk #TODO fix this so it doesn't always accelerate if there's no space
        accel = cs.x - x > 2*dmodel.phys_param.l_car ? 1 : -3. #slam brakes
    else
        accel = -sign(cs.x - x)
    end

	return accel*dmodel.phys_param.dt #XXX times accel interval?

end

function gen_lane_change(bmodel::AvoidModel, dmodel::IDMMOBILModel, s::MLState, neighborhood::Array{Int,1}, idx::Int, rng::AbstractRNG)

	closest_car = closest_car(dmodel,s,neighborhood,idx,bmodel.jerk)
    println("L-C1")
	if closest_car == 0
		return 0.
	end

	pos = s.cars[idx].y
	cc = s.cars[closest_car]

	lanechange = -sign(cc.y - pos)
	#if same lane, go in direction with more room
	if lanechange == 0
		if pos >= 2*dmodel.phys_param.nb_lanes - 1 - pos >= pos #more room to left >= biases to left
			lanechange = 1
		else
			lanechange = -1
		end
	end
	if lanechange > 0 && pos >= 2*dmodel.phys_param.nb_lanes - 1
		lanechange = 0
	elseif lanechange < 0 && pos <= 1
		lanechange = 0
	end

	return lanechange
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type flash_steady
	flash_steady()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed Stream",
				:PosX=>0,
				:PosY=>0.48,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.43,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0.43,
				:PosY=>0,
				:Symbol=>"_{out}^{Vapour}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat Duty",
				:PosX=>1,
				:PosY=>0.81,
				:Protected=>true,
				:Symbol=>"Q_{in}"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Vaporization fraction",
				:Symbol=>"\\phi",
				:Protected=>true
			)),
			enth_mol(Dict{Symbol,Any}(
				:Brief=>"Mixture enthalpy",
				:Hidden=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Symbol=>"P_{ratio}",
				:Protected=>true
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P",
				:Protected=>true
			)),
			[
				:([vfrac, OutletLiquid.z, OutletVapour.z] = PP.Flash(OutletVapour.T, OutletVapour.P, Inlet.z)),
				:([vfrac,OutletLiquid.z,OutletVapour.z]=PP.FlashPH(OutletLiquid.P,h,Inlet.z)),
				:(Inlet.F = OutletVapour.F + OutletLiquid.F),
				:(OutletVapour.F = Inlet.F * vfrac),
				:(Inlet.F*(h - Inlet.h) = InletQ),
				:(Inlet.F*h = Inlet.F*(1-vfrac)*OutletLiquid.h + Inlet.F*vfrac*OutletVapour.h),
				:(OutletVapour.T = OutletLiquid.T),
				:(OutletVapour.P = OutletLiquid.P),
				:(OutletLiquid.P = Inlet.P - Pdrop),
				:(OutletLiquid.P = Inlet.P * Pratio),
			],
			[
				"The flash calculation","Chemical equilibrium","Global Molar Balance","Vapour Fraction","Energy Balance","","Thermal Equilibrium","Mechanical Equilibrium","Pressure Drop","Pressure Ratio",
			],
			[:PP,:NComp,],
			[:Inlet,:OutletLiquid,:OutletVapour,:InletQ,:vfrac,:h,:Pratio,:Pdrop,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	InletQ::power
	vfrac::fraction
	h::enth_mol
	Pratio::positive
	Pdrop::press_delta
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export flash_steady
function setEquationFlow(in::flash_steady)
	if vfrac > 0 && vfrac <1
		addEquation(1)
	else
		addEquation(2)
	end
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
end
function atributes(in::flash_steady,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Flash"
	fields[:Brief]="Model of a static PH flash."
	fields[:Info]="This model is for using the flashPH routine available on VRTherm.

== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.
* static model.

== SPECIFY ==
* The Inlet stream;
* The heat duty;
* The outlet pressure.
"
	drive!(fields,_)
	return fields
end
flash_steady(_::Dict{Symbol,Any})=begin
	newModel=flash_steady()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module AstroImageView

using AstroImages, GtkReactive, Gtk.ShortNames, Graphics, Gtk, Colors, IntervalSets, Cairo
using AstroImages: render
using WCS, FITSIO

export ui_basic, ui_advanced

include("utils.jl")
include("basic.jl")
include("advanced.jl")

end # module


text:
using DashBootstrapComponents, DashHtmlComponents
using HTTP

LOREM = rstrip(String(HTTP.get("https://raw.githubusercontent.com/facultyai/dash-bootstrap-components/main/docs/components_page/components/modal/lorem.txt").body), ['\n'])

modal = html_div([
    dbc_button("Scrolling modal", id="open-scroll", className="mr-1", n_clicks=0),
    dbc_button("Modal with scrollable body", id="open-body-scroll", n_clicks=0),
    dbc_modal(
        [
            dbc_modalheader("Header"),
            dbc_modalbody(LOREM),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id="close-scroll",
                    className="ml-auto",
                    n_clicks=0,
                ),
            ),
        ],
        id="modal-scroll",
        is_open=false,
    ),
    dbc_modal(
        [
            dbc_modalheader("Header"),
            dbc_modalbody(LOREM),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id="close-body-scroll",
                    className="ml-auto",
                    n_clicks=0,
                ),
            ),
        ],
        id="modal-body-scroll",
        scrollable=true,
        is_open=false,
    ),
]);

function toggle_modal(n1, n2, is_open)
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

callback!(
    app,
    Output("modal-scroll", "is_open"),
    Input("open-scroll", "n_clicks"),
    Input("close-scroll", "n_clicks"),
    State("modal-scroll", "is_open"),
) do n1, n2, is_open
    return toggle_modal(n1, n2, is_open)
end;

callback!(
    app,
    Output("modal-body-scroll", "is_open"),
    Input("open-body-scroll", "n_clicks"),
    Input("close-body-scroll", "n_clicks"),
    State("modal-body-scroll", "is_open"),
) do n1, n2, is_open
    return toggle_modal(n1, n2, is_open)
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module mt_elasticity_t3

using Test
using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities
using MeshSteward: T3block
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs
using Elfel.FEIterators: jacjac
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!

E = 1.0;
nu = 1.0/3;
D = SMatrix{3, 3}(E / (1 - nu^2) * [1 nu 0
                                    nu 1 0
                                    0 0 (1 - nu) / 2])
A = 1.0 # length of the side of the square
N = 100;# number of subdivisions along the sides of the square domain

function genmesh()
    conn = T3block(A, A, N, N)
    mesh = Mesh()
    attach!(mesh, conn)
    return mesh
end

function assembleK(fesp, D)
    function integrateK!(ass, geom, elit, qpit, D)
        B = (g, k) -> k == 1 ? SVector{3}((g[1], 0, g[2])) : SVector{3}((0, g[2], g[1]))
        c = edofcompnt(elit.fesp)
        nedof = ndofsperel(elit)
        ke = LocalMatrixAssembler(nedof, nedof, 0.0)
        for el in elit
            init!(ke, eldofs(el), eldofs(el))
            for qp in qpit
                Jac, J = jacjac(el, qp)
                gradN = bfungrad(qp, Jac)
                JxW = J * weight(qp)
                for j in 1:nedof
                    DBj = D * B(gradN[j], c[j])
                    for i in 1:nedof
                        Bi = B(gradN[i], c[i])
                        ke[i, j] += dot(DBj, Bi) * JxW
                    end
                end
            end
            assemble!(ass, ke)
        end
        return ass
    end

    elit = FEIterator(fesp)
    qpit = QPIterator(fesp, (kind = :default,))
    geom = geometry(fesp.mesh)
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, ndofs(fesp), ndofs(fesp))
    integrateK!(ass, geom, elit, qpit, D)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function test()
    mesh = genmesh()
    fesp = FESpace(Float64, mesh, FEH1_T3(), 2)
    locs = geometry(mesh)
    inflate = A / N / 100
    box = [0.0 0.0 0.0 A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(fesp, 0, i, 1, 0.0)
        setebc!(fesp, 0, i, 2, 0.0)
    end
    box = [A A 0.0 A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(fesp, 0, i, 1, A / 10)
        setebc!(fesp, 0, i, 2, 0.0)
    end
    numberfreedofs!(fesp)
    numberdatadofs!(fesp)
    # @show nunknowns(fesp)
    K = assembleK(fesp, D)
    U = fill(0.0, ndofs(fesp))
    gathersysvec!(U, fesp)
    F = fill(0.0, ndofs(fesp))
    solve!(U, K, F, nunknowns(fesp))
    scattersysvec!(fesp, U)
    makeattribute(fesp, "U", 1:2)
    vtkwrite("elast_stretch_t3-U", baseincrel(mesh), [(name = "U", allxyz = true)])
    try rm("elast_stretch_t3-U.vtu"); catch end
end

end
using .mt_elasticity_t3
mt_elasticity_t3.test()

module mt_elast_stretch_t6

using Test
using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities
using MeshSteward: T6block
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs
using Elfel.FEIterators: jacjac
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!

E = 1.0;
nu = 1.0/3;
D = SMatrix{3, 3}(E / (1 - nu^2) * [1 nu 0
                                    nu 1 0
                                    0 0 (1 - nu) / 2])
A = 1.0 # length of the side of the square
N = 100;# number of subdivisions along the sides of the square domain

function genmesh()
    conn = T6block(A, A, N, N)
    mesh = Mesh()
    attach!(mesh, conn)
    return mesh
end

function assembleK(fesp, D)
    function integrateK!(ass, geom, elit, qpit, D)
        B = (g, k) -> k == 1 ? SVector{3}((g[1], 0, g[2])) : SVector{3}((0, g[2], g[1]))
        c = edofcompnt(elit.fesp)
        nedof = ndofsperel(elit)
        ke = LocalMatrixAssembler(nedof, nedof, 0.0)
        for el in elit
            init!(ke, eldofs(el), eldofs(el))
            for qp in qpit
                Jac, J = jacjac(el, qp)
                gradN = bfungrad(qp, Jac)
                JxW = J * weight(qp)
                for j in 1:nedof
                    DBj = D * B(gradN[j], c[j])
                    for i in 1:nedof
                        Bi = B(gradN[i], c[i])
                        ke[i, j] += dot(DBj, Bi) * JxW
                    end
                end
            end
            assemble!(ass, ke)
        end
        return ass
    end

    elit = FEIterator(fesp)
    qpit = QPIterator(fesp, (kind = :default, npts = 3,))
    geom = geometry(fesp.mesh)
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, ndofs(fesp), ndofs(fesp))
    integrateK!(ass, geom, elit, qpit, D)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function test()
    mesh = genmesh()
    fesp = FESpace(Float64, mesh, FEH1_T6(), 2)
    locs = geometry(mesh)
    inflate = A / N / 100
    box = [0.0 0.0 0.0 A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(fesp, 0, i, 1, 0.0)
        setebc!(fesp, 0, i, 2, 0.0)
    end
    box = [A A 0.0 A]
    vl = vselect(locs; box = box, inflate = inflate)
    for i in vl
        setebc!(fesp, 0, i, 1, A / 10)
        setebc!(fesp, 0, i, 2, 0.0)
    end
    numberfreedofs!(fesp)
    numberdatadofs!(fesp)
    # @show nunknowns(fesp)
    K = assembleK(fesp, D)
    U = fill(0.0, ndofs(fesp))
    gathersysvec!(U, fesp)
    F = fill(0.0, ndofs(fesp))
    solve!(U, K, F, nunknowns(fesp))
    scattersysvec!(fesp, U)
    makeattribute(fesp, "U", 1:2)
    vtkwrite("elast_stretch_t6-U", baseincrel(mesh), [(name = "U", allxyz = true)])
    try rm("elast_stretch_t6-U.vtu"); catch end
end

end
using .mt_elast_stretch_t6
mt_elast_stretch_t6.test()


text:
module test_bukdu_system

module NAStuff

using Bukdu
struct NA <: ApplicationController; conn::Conn end
index(c::NA) = nothing
hello(c::NA) = TEST_INTERNAL_ERROR

end # module NAStuff

using Test
using Bukdu
using .Bukdu.System

function index
end

routes() do
    get("/na", NAStuff.NA, index)
    post("/hello", NAStuff.NA, NAStuff.hello)
end

result = Router.call(get, "/na")
@test occursin("Bukdu.System.NotApplicableError", result.got)
@test result.resp.status == 500
@test result.route.action === System.not_applicable

Plug.Loggers.config[:error_stackframes_range] = 1:2
result = Router.call(post, "/hello")
@test result.route.action === System.internal_error

Routing.reset!()

@test Plug.Loggers._regularize_text("가1", 1) == "가"

struct Controller <: ApplicationController
    conn::Conn
end
function index(::Controller)
end

struct VeryLongNamedController <: ApplicationController
    conn::Conn
end
function index(::VeryLongNamedController)
end

routes() do
    get("/just", Controller, index)
    get("/long", VeryLongNamedController, index)
end
Router.call(get, "/just")
Router.call(get, "/long")

Routing.reset!()

end # module test_bukdu_system

module test_bukdu_system_proc_time

using Bukdu

function Bukdu.System.catch_request(route::Bukdu.Route, conn)
    conn.private[:req_time_ns] = time_ns()
end

using Logging: AbstractLogger

struct MyLogger <: AbstractLogger
    stream
end

function Plug.Loggers.info_response(logger::MyLogger, conn::Conn, route::Bukdu.RouteAction)
    io = logger.stream
    proc_time = (time_ns() - conn.private[:req_time_ns]) / 1e9
    print(io, proc_time, ' ')
    Plug.Loggers.default_info_response(io, conn, route)
end

plug(MyLogger, IOContext(Core.stdout, :color => Plug.Loggers.have_color()))

get("/") do conn::Conn
    42
end

Router.call(get, "/")

plug(Plug.Loggers.DefaultLogger)
Routing.reset!()

end # module test_bukdu_system_proc_time