text:
using Random, LinearAlgebra

import ..timeevolution: dschroedinger!

# TODO: Remove imports
import RecursiveArrayTools.copyat_or_push!

"""
    mcwf_h(tspan, rho0, Hnh, J; <keyword arguments>)

Calculate MCWF trajectory where the Hamiltonian is given in hermitian form.

For more information see: [`mcwf`](@ref)
"""
function mcwf_h(tspan, psi0::Ket, H::AbstractOperator, J;
        seed=rand(UInt), rates=nothing,
        fout=nothing, Jdagger=dagger.(J),
        tmp=copy(psi0),
        display_beforeevent=false, display_afterevent=false,
        kwargs...)
    check_mcwf(psi0, H, J, Jdagger, rates)
    f(t, psi, dpsi) = dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, tmp)
    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
    integrate_mcwf(f, j, tspan, psi0, seed, fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    mcwf_nh(tspan, rho0, Hnh, J; <keyword arguments>)

Calculate MCWF trajectory where the Hamiltonian is given in non-hermitian form.

```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```

For more information see: [`mcwf`](@ref)
"""
function mcwf_nh(tspan, psi0::Ket, Hnh::AbstractOperator, J;
        seed=rand(UInt), fout=nothing,
        display_beforeevent=false, display_afterevent=false,
        kwargs...)
    check_mcwf(psi0, Hnh, J, J, nothing)
    f(t, psi, dpsi) = dschroedinger!(dpsi, Hnh, psi)
    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, nothing)
    integrate_mcwf(f, j, tspan, psi0, seed, fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    mcwf(tspan, psi0, H, J; <keyword arguments>)

Integrate the master equation using the MCWF method.

There are two implementations for integrating the non-hermitian
schroedinger equation:

* [`mcwf_h`](@ref): Usual formulation with Hamiltonian + jump operators
separately.
* [`mcwf_nh`](@ref): Variant with non-hermitian Hamiltonian.

The `mcwf` function takes a normal Hamiltonian, calculates the
non-hermitian Hamiltonian and then calls [`mcwf_nh`](@ref) which is
slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
be displayed.
* `psi0`: Initial state vector.
* `H`: Arbitrary Operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
operator type.
* `seed=rand()`: Seed used for the random number generator.
* `rates=ones()`: Vector of decay rates.
* `fout`: If given, this function `fout(t, psi)` is called every time an
output should be displayed. ATTENTION: The state `psi` is neither
normalized nor permanent! It is still in use by the ode solve
and therefore must not be changed.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
operators. If they are not given they are calculated automatically.
* `display_beforeevent=false`: `fout` is called before every jump.
* `display_afterevent=false`: `fout` is called after every jump.
* `display_jumps=false`: If set to true, an additional list of times and indices
is returned. These correspond to the times at which a jump occured and the index
of the jump operators with which the jump occured, respectively.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function mcwf(tspan, psi0::Ket, H::AbstractOperator, J;
        seed=rand(UInt), rates=nothing,
        fout=nothing, Jdagger=dagger.(J),
        display_beforeevent=false, display_afterevent=false,
        kwargs...)
    isreducible = check_mcwf(psi0, H, J, Jdagger, rates)
    if !isreducible
        tmp = copy(psi0)
        dmcwf_h_(t, psi, dpsi) = dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, tmp)
        j_h(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
        integrate_mcwf(dmcwf_h_, j_h, tspan, psi0, seed,
            fout;
            display_beforeevent=display_beforeevent,
            display_afterevent=display_afterevent,
            kwargs...)
    else
        Hnh = copy(H)
        if isa(rates, Nothing)
            for i=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im)*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im*rates[i])*Jdagger[i]*J[i]
            end
        end
        dmcwf_nh_(t, psi, dpsi) = dschroedinger!(dpsi, Hnh, psi)
        j_nh(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates)
        integrate_mcwf(dmcwf_nh_, j_nh, tspan, psi0, seed,
            fout;
            display_beforeevent=display_beforeevent,
            display_afterevent=display_afterevent,
            kwargs...)
    end
end

"""
    mcwf_dynamic(tspan, psi0, f; <keyword arguments>)

Integrate the master equation using the MCWF method with dynamic
Hamiltonian and Jump operators.

The `mcwf` function takes a normal Hamiltonian, calculates the
non-hermitian Hamiltonian and then calls [`mcwf_nh`](@ref) which is
slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
be displayed.
* `psi0`: Initial state vector.
* `f`: Function `f(t, psi) -> (H, J, Jdagger)` or `f(t, psi) -> (H, J, Jdagger, rates)`
    that returns the time-dependent Hamiltonian and Jump operators.
* `seed=rand()`: Seed used for the random number generator.
* `rates=ones()`: Vector of decay rates.
* `fout`: If given, this function `fout(t, psi)` is called every time an
output should be displayed. ATTENTION: The state `psi` is neither
normalized nor permanent! It is still in use by the ode solve
and therefore must not be changed.
* `display_beforeevent=false`: `fout` is called before every jump.
* `display_afterevent=false`: `fout` is called after every jump.
* `display_jumps=false`: If set to true, an additional list of times and indices
is returned. These correspond to the times at which a jump occured and the index
of the jump operators with which the jump occured, respectively.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function mcwf_dynamic(tspan, psi0::Ket, f;
    seed=rand(UInt), rates=nothing,
    fout=nothing, display_beforeevent=false, display_afterevent=false,
    kwargs...)
    tmp = copy(psi0)
    dmcwf_(t, psi, dpsi) = dmcwf_h_dynamic!(dpsi, f, rates, psi, tmp, t)
    j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, f, psi_new, rates)
    integrate_mcwf(dmcwf_, j_, tspan, psi0, seed,
        fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    mcwf_nh_dynamic(tspan, rho0, f; <keyword arguments>)

Calculate MCWF trajectory where the dynamic Hamiltonian is given in non-hermitian form.

For more information see: [`mcwf_dynamic`](@ref)
"""
function mcwf_nh_dynamic(tspan, psi0::Ket, f;
    seed=rand(UInt), rates=nothing,
    fout=nothing, display_beforeevent=false, display_afterevent=false,
    kwargs...)
    dmcwf_(t, psi, dpsi) = dmcwf_nh_dynamic!(dpsi, f, psi, t)
    j_(rng, t, psi, psi_new) = jump_dynamic(rng, t, psi, f, psi_new, rates)
    integrate_mcwf(dmcwf_, j_, tspan, psi0, seed,
        fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    dmcwf_h_dynamic!(dpsi, f, rates, psi, dpsi_cache, t)

Compute the Hamiltonian and jump operators as `H,J,Jdagger=f(t,psi)` and
update `dpsi` according to a non-Hermitian Schrödinger equation.

See also: [`mcwf_dynamic`](@ref), [`dmcwf_h!`](@ref), [`dmcwf_nh_dynamic`](@ref)
"""
function dmcwf_h_dynamic!(dpsi, f, rates, psi, dpsi_cache, t)
    result = f(t, psi)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    if length(result) == 3
        H, J, Jdagger = result
        rates_ = rates
    else
        H, J, Jdagger, rates_ = result
    end
    QO_CHECKS[] && check_mcwf(psi, H, J, Jdagger, rates_)
    dmcwf_h!(dpsi, H, J, Jdagger, rates_, psi, dpsi_cache)
end

"""
    dmcwf_nh_dynamic!(dpsi, f, psi, t)

Compute the non-Hermitian Hamiltonian and jump operators as `H,J,Jdagger=f(t,psi)`
and update `dpsi` according to a Schrödinger equation.

See also: [`mcwf_nh_dynamic`](@ref), [`dmcwf_nh!`](@ref), [`dschroedinger!`](@ref)
"""
function dmcwf_nh_dynamic!(dpsi, f, psi, t)
    result = f(t, psi)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    H, J, Jdagger = result[1:3]
    QO_CHECKS[] && check_mcwf(psi, H, J, Jdagger, nothing)
    dschroedinger!(dpsi, H, psi)
end

function jump_dynamic(rng, t, psi, f, psi_new, rates)
    result = f(t, psi)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    J = result[2]
    if length(result) == 3
        rates_ = rates
    else
        rates_ = result[4]
    end
    jump(rng, t, psi, J, psi_new, rates_)
end

"""
    integrate_mcwf(dmcwf, jumpfun, tspan, psi0, seed; fout, kwargs...)

Integrate a single Monte Carlo wave function trajectory.

# Arguments
* `dmcwf`: A function `f(t, psi, dpsi)` that calculates the time-derivative of
        `psi` at time `t` and stores the result in `dpsi`.
* `jumpfun`: A function `f(rng, t, psi, dpsi)` that uses the random number
        generator `rng` to determine if a jump is performed and stores the
        result in `dpsi`.
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `psi0`: Initial state vector.
* `seed`: Seed used for the random number generator.
* `fout`: If given, this function `fout(t, psi)` is called every time an
        output should be displayed. ATTENTION: The state `psi` is neither
        normalized nor permanent! It is still in use by the ode solver
        and therefore must not be changed.
* `kwargs`: Further arguments are passed on to the ode solver.
"""
function integrate_mcwf(dmcwf, jumpfun, tspan,
                        psi0, seed, fout::Function;
                        display_beforeevent=false, display_afterevent=false,
                        display_jumps=false,
                        save_everystep=false, callback=nothing,
                        saveat=tspan,
                        alg=OrdinaryDiffEq.DP5(),
                        kwargs...)

    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    # Display before or after events
    function save_func!(affect!,integrator)
        affect!.saveiter += 1
        copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
        copyat_or_push!(affect!.saved_values.saveval, affect!.saveiter,
            affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
        return nothing
    end
    save_before! = display_beforeevent ? save_func! : (affect!,integrator)->nothing
    save_after! = display_afterevent ? save_func! : (affect!,integrator)->nothing

    # Display jump operator index and times
    jump_t = eltype(tspan_)[]
    jump_index = Int[]
    save_t_index = if display_jumps
        function(t,i)
            push!(jump_t,t)
            push!(jump_index,i)
            return nothing
        end
    else
        (t,i)->nothing
    end

    function fout_(x, t, integrator)
        recast!(state,x)
        fout(t, state)
    end

    state = copy(psi0)
    dstate = copy(psi0)
    out_type = pure_inference(fout, Tuple{eltype(tspan_),typeof(state)})
    out = DiffEqCallbacks.SavedValues(eltype(tspan_),out_type)
    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan_,
                                         save_everystep=save_everystep,
                                         save_start = false)

    cb = jump_callback(jumpfun, seed, scb, save_before!, save_after!, save_t_index, psi0)
    full_cb = OrdinaryDiffEq.CallbackSet(callback,cb,scb)

    function df_(dx, x, p, t)
        recast!(state,x)
        recast!(dstate,dx)
        dmcwf(t, state, dstate)
        recast!(dx,dstate)
    end

    prob = OrdinaryDiffEq.ODEProblem{true}(df_, as_vector(psi0), (tspan_[1],tspan_[end]))

    sol = OrdinaryDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-6,
                abstol = 1.0e-8,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=full_cb, kwargs...)

    if display_jumps
        return out.t, out.saveval, jump_t, jump_index
    else
        return out.t, out.saveval
    end
end

function integrate_mcwf(dmcwf, jumpfun, tspan,
                        psi0, seed, fout::Nothing;
                        kwargs...)
    function fout_(t, x)
        return normalize(x)
    end
    integrate_mcwf(dmcwf, jumpfun, tspan, psi0, seed, fout_; kwargs...)
end

function jump_callback(jumpfun, seed, scb, save_before!,
                        save_after!, save_t_index, psi0)

    tmp = copy(psi0)
    psi_tmp = copy(psi0)

    rng = MersenneTwister(convert(UInt, seed))
    jumpnorm = Ref(rand(rng))
    djumpnorm(x, t, integrator) = norm(x)^2 - (1-jumpnorm[])

    function dojump(integrator)
        x = integrator.u
        t = integrator.t

        affect! = scb.affect!
        save_before!(affect!,integrator)
        recast!(psi_tmp,x)
        i = jumpfun(rng, t, psi_tmp, tmp)
        x .= tmp.data
        save_after!(affect!,integrator)
        save_t_index(t,i)

        jumpnorm[] = rand(rng)
        return nothing
    end

    return OrdinaryDiffEq.ContinuousCallback(djumpnorm,dojump,
                     save_positions = (false,false))
end
as_vector(psi::StateVector) = psi.data

"""
    jump(rng, t, psi, J, psi_new)

Default jump function.

# Arguments
* `rng:` Random number generator
* `t`: Point of time where the jump is performed.
* `psi`: State vector before the jump.
* `J`: List of jump operators.
* `psi_new`: Result of jump.
"""
function jump(rng, t, psi, J, psi_new, rates::Nothing)
    if length(J)==1
        QuantumOpticsBase.mul!(psi_new,J[1],psi,true,false)
        psi_new.data ./= norm(psi_new)
        i=1
    else
        probs = zeros(real(eltype(psi)), length(J))
        for i=1:length(J)
            QuantumOpticsBase.mul!(psi_new,J[i],psi,true,false)
            probs[i] = dot(psi_new.data, psi_new.data)
        end
        cumprobs = cumsum(probs./sum(probs))
        r = rand(rng)
        i = findfirst(cumprobs.>r)
        QuantumOpticsBase.mul!(psi_new,J[i],psi,one(eltype(psi))/sqrt(probs[i]),zero(eltype(psi)))
    end
    return i
end

function jump(rng, t, psi, J, psi_new, rates::AbstractVector)
    if length(J)==1
        QuantumOpticsBase.mul!(psi_new,J[1],psi,eltype(psi)(sqrt(rates[1])),zero(eltype(psi)))
        psi_new.data ./= norm(psi_new)
        i=1
    else
        probs = zeros(real(eltype(psi)), length(J))
        for i=1:length(J)
            QuantumOpticsBase.mul!(psi_new,J[i],psi,eltype(psi)(sqrt(rates[i])),zero(eltype(psi)))
            probs[i] = dot(psi_new.data, psi_new.data)
        end
        cumprobs = cumsum(probs./sum(probs))
        r = rand(rng)
        i = findfirst(cumprobs.>r)
        QuantumOpticsBase.mul!(psi_new,J[i],psi,eltype(psi)(sqrt(rates[i]/probs[i])),zero(eltype(psi)))
    end
    return i
end

"""
    dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, dpsi_cache)

Update `dpsi` according to a non-hermitian Schrödinger equation. The
non-hermitian Hamiltonian is given in two parts - the hermitian part H and
the jump operators J.

See also: [`mcwf`](@ref)
"""
function dmcwf_h!(dpsi, H, J, Jdagger, rates::Nothing, psi, dpsi_cache)
    QuantumOpticsBase.mul!(dpsi,H,psi,eltype(psi)(-im),zero(eltype(psi)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(dpsi_cache,J[i],psi,true,false)
        QuantumOpticsBase.mul!(dpsi,Jdagger[i],dpsi_cache,eltype(psi)(-0.5),one(eltype(psi)))
    end
    return dpsi
end

function dmcwf_h!(dpsi, H, J, Jdagger, rates::AbstractVector, psi, dpsi_cache)
    QuantumOpticsBase.mul!(dpsi,H,psi,eltype(psi)(-im),zero(eltype(psi)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(dpsi_cache,J[i],psi,eltype(psi)(rates[i]),zero(eltype(psi)))
        QuantumOpticsBase.mul!(dpsi,Jdagger[i],dpsi_cache,eltype(psi)(-0.5),one(eltype(psi)))
    end
    return dpsi
end

"""
    check_mcwf(psi0, H, J, Jdagger, rates)

Check input of mcwf.
"""
function check_mcwf(psi0, H, J, Jdagger, rates)
    # TODO: replace type checks by dispatch; make types of J known
    check_psi0(psi0)
    isreducible = true
    if !(isa(H, DenseOpType) || isa(H, SparseOpType))
        isreducible = false
    end
    for j=J
        @assert isa(j, AbstractOperator)
        if !(isa(j, DenseOpType) || isa(j, SparseOpType))
            isreducible = false
        end
    end
    for j=Jdagger
        @assert isa(j, AbstractOperator)
        if !(isa(j, DenseOpType) || isa(j, SparseOpType))
            isreducible = false
        end
    end
    @assert length(J) == length(Jdagger)
    if isa(rates, Matrix)
        throw(ArgumentError("Matrix of decay rates not supported for MCWF!
            Use diagonaljumps(rates, J) to calculate new rates and jump operators."))
    elseif isa(rates, Vector)
        @assert length(rates) == length(J)
    end
    isreducible
end

"""
    diagonaljumps(rates, J)

Diagonalize jump operators.

The given matrix `rates` of decay rates is diagonalized and the
corresponding set of jump operators is calculated.

# Arguments
* `rates`: Matrix of decay rates.
* `J`: Vector of jump operators.
"""
function diagonaljumps(rates::AbstractMatrix, J)
    @assert length(J) == size(rates)[1] == size(rates)[2]
    d, v = eigen(rates)
    d, [sum([v[j, i]*J[j] for j=1:length(d)]) for i=1:length(d)]
end

function diagonaljumps(rates::AbstractMatrix, J::Vector{T}) where T<:Union{LazySum,LazyTensor,LazyProduct}
    @assert length(J) == size(rates)[1] == size(rates)[2]
    d, v = eigen(rates)
    d, [LazySum([v[j, i]*J[j] for j=1:length(d)]...) for i=1:length(d)]
end


text:
const m1 = ReTerms.ColMajorLowerTriangular(3)

@test size(m1) == (3,3)
@test size(m1,1) == 3
@test size(m1,2) == 3
@test size(m1,3) == 1
const m1c = copy(m1)
@test ReTerms.nθ(m1) == 6
@test m1c == m1
@test m1[:θ] == [1.,0,0,1,0,1]
@test lowerbd(m1) == [0.,-Inf,-Inf,0.,-Inf,0.]

m1[:θ] = collect(1.:6)                  # assignment of parameters
@test m1[:θ] == collect(1.:6)           # check assignment and extraction
@test m1c != m1
@test copy!(m1c,m1) == m1

@test full(m1) == reshape([1.,2,3,0,4,5,0,0,6],(3,3))

ext: jl


text:
## Canonical Form of Normal distribution

immutable NormalCanon <: ContinuousUnivariateDistribution
    h::Float64       # σ^(-2) * μ
    prec::Float64    # σ^(-2)
    μ::Float64       # μ

    function NormalCanon(h::Float64, prec::Float64)
    	prec > 0. || throw(ArgumentError("prec should be positive."))
    	new(h, prec, h / prec)
    end

    NormalCanon(h::Real, prec::Real) = NormalCanon(float64(h), float64(prec))
    NormalCanon() = new(0.0, 1.0, 0.0)
end

## conversion between Normal and NormalCanon

Base.convert(::Type{Normal}, cf::NormalCanon) = Normal(cf.μ, 1.0 / sqrt(cf.prec))
Base.convert(::Type{NormalCanon}, d::Normal) = (J = 1.0 / abs2(σ); NormalCanon(J * d.μ, J))
canonform(d::Normal) = convert(NormalCanon, d)

## Basic properties

@continuous_distr_support NormalCanon -Inf Inf

mean(cf::NormalCanon) = cf.μ
median(cf::NormalCanon) = mean(cf)
mode(cf::NormalCanon) = mean(cf)

skewness(cf::NormalCanon) = 0.0
kurtosis(cf::NormalCanon) = 0.0

var(cf::NormalCanon) = 1.0 / cf.prec
std(cf::NormalCanon) = sqrt(var(cf))

entropy(cf::NormalCanon) = 0.5 * (log2π + 1. - log(cf.prec))

# Evaluation

pdf(d::NormalCanon, x::Real) = (sqrt(d.prec) / sqrt2π) * exp(-0.5 * d.prec * abs2(x - d.μ))
logpdf(d::NormalCanon, x::Real) = 0.5 * (log(d.prec) - log2π - d.prec * abs2(x - d.μ))

zval(d::NormalCanon, x::Real) = (x - d.μ) * sqrt(d.prec)
xval(d::NormalCanon, z::Real) = d.μ + z / sqrt(d.prec)

cdf(d::NormalCanon, x::Real) = Φ(zval(d,x))
ccdf(d::NormalCanon, x::Real) = Φc(zval(d,x))
logcdf(d::NormalCanon, x::Real) = logΦ(zval(d,x))
logccdf(d::NormalCanon, x::Real) = logΦc(zval(d,x))

quantile(d::NormalCanon, p::Real) = xval(d, Φinv(p))
cquantile(d::NormalCanon, p::Real) = xval(d, -Φinv(p))
invlogcdf(d::NormalCanon, p::Real) = xval(d, logΦinv(p))
invlogccdf(d::NormalCanon, p::Real) = xval(d, -logΦinv(p))

## Sampling

rand(cf::NormalCanon) = cf.μ + randn() / sqrt(cf.prec)
rand!{T<:FloatingPoint}(cf::NormalCanon, r::Array{T}) = rand!(convert(Normal, cf), r)


text:
"""
  add_coord!(pt::Point, var::Variable, val::Number)

  Sparsely add `val` to the `var` coordinate of `pt`.
"""
function add_coord!(pt::Point, var::Variable, val::Number)
  return add_to_dict!(pt.coords, var, val)
end

## Overloads for transparent iteration on `coords` attribute
start(pt::Point) = start(pt.coords)
next(pt::Point, state) = next(pt.coords, state)
done(pt::Point, state) = done(pt.coords, state)
length(pt::Point) = length(pt.coords)
haskey(pt::Point, key) = haskey(pt.coords, key)
keys(pt::Point) = keys(pt.coords)
values(pt::Point) = values(pt.coords)
getindex(pt::Point, var::Variable) = pt.coords[var]
function setindex!(pt::Point, var::Variable, val::Number)
  if val != 0
    setindex!(pt.coords, var, val)
  end
  return
end

function add!(pt1::Point, pt2::Point)
  for (var, val) in pt2.coords
    add_coord!(pt1, var, val)
  end
end

function add(pt1::Point, pt2::Point)
  pt_coord = copy(pt1.coords)
  pt = Point(pt_coord)
  add!(pt, pt2)
  return pt
end

+(pt1::Point, pt2::Point) = add(pt1, pt2)
-(pt1::Point, pt2::Point) = pt1 + (-1)*pt2

function merge(x::Point ...)
  pt_res = Point()
  for pt in x
    add!(pt_res, pt)
  end
  return pt_res
end

function *(pt1::Point, λ::Number)
  pt = Point()
  if λ == 0
    return pt
  end
  for (var, val) in pt1
    pt[var] = λ*val
  end
  return pt
end
*(λ::Number, pt1::Point) = pt1*λ

function norm(pt::Point, p::Real = 2)
    p > 0 || error("norm(): p should be positive ($p ≤ 0).")
    if p == Inf
        maxi = -Inf
        for (var, val) in pt
            # iscomplex(var) && warn("norm(): var $(var) is not real, splitting real and imag parts.")
            if !iscomplex(var) && (val > maxi)
                maxi = val
            elseif iscomplex(var) && max(real(val), imag(val)) > maxi
                maxi = max(real(val), imag(val))
            end
        end
        return maxi
    elseif p == 1
      s = 0
      for (var, val) in pt
          if !iscomplex(var)
            s += abs(val)
          else
            # warn("norm(): var $(var) is not real, splitting real and imag parts.")
            s += abs(real(val)) + abs(imag(val))
          end
      end
      return s
    else
      s = 0
      for (var, val) in pt
          if !iscomplex(var)
            s += val^p
          else
            # warn("norm(): var $(var) is not real, splitting real and imag parts.")
            s += real(val)^p + imag(val)^p
          end
      end
      return s^(1/p)
  end
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Documenter, ActuaryUtilities

# Setup for doctests embedded in docstrings.
DocMeta.setdocmeta!(ActuaryUtilities, :DocTestSetup, :(import Pkg; Pkg.add("DayCounts");using ActuaryUtilities, Dates,DayCounts),recursive=true)

makedocs(;
    modules=[ActuaryUtilities],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaActuary/ActuaryUtilities.jl/blob/{commit}{path}#L{line}",
    sitename="ActuaryUtilities.jl",
    authors="Alec Loudenback",
)

deploydocs(;
    repo="github.com/JuliaActuary/ActuaryUtilities.jl",
)

max_stars_repo_name: alecloudenback/ActuaryUtilities

max_issues_repo_name: alecloudenback/ActuaryUtilities

max_forks_repo_name: alecloudenback/ActuaryUtilities


text:
struct ValidType <: Options end
struct InvalidType end

@testset "Options" begin
    @test Bioview.isavalidtype(Options) == true
    @test Bioview.isavalidtype(ValidType) == true
    @test_throws BioviewError Bioview.isavalidtype(InvalidType)

    opt = ValidType()
    oopt = ValidType()
    opts = (opt, oopt)
    @test Bioview.isavalidtype(opts) == true

    invalid_opt = InvalidType()
    opts = (opt, oopt, invalid_opt)
    @test_throws BioviewError Bioview.isavalidtype(opts)
end


text:
Base.@kwdef mutable struct Mesh
    type   ::Union{String, Missing}          = missing
    nc     ::Union{Int64,  Missing}          = missing
    nn     ::Union{Int64,  Missing}          = missing
    nn_el  ::Union{Int64,  Missing}          = missing
    nf_el  ::Union{Int64,  Missing}          = missing
    xn     ::Union{Vector{Float64}, Missing} = missing # node x coordinate
    yn     ::Union{Vector{Float64}, Missing} = missing # node y coordinate
    xf     ::Union{Vector{Float64}, Missing} = missing # face x coordinate
    yf     ::Union{Vector{Float64}, Missing} = missing # face y coordinate
    c2n    ::Union{Vector{Int64},   Missing} = missing # cell 2 node numbering
    c2f    ::Union{Vector{Int64},   Missing} = missing # cell 2 face numbering
    xyviz  ::Union{Vector{Float64}, Missing} = missing # table of coordinates for triangle visualisation with PyPlot
    numviz ::Union{Vector{Float64}, Missing} = missing # table of numbering for triangle visualisation with PyPlot
end

# Create structure with default missing fields
mesh = Mesh()

# Check if value is missing
ismissing(mesh.xn)

# Assign value to fields
mesh.type = "triangle"
mesh.nc   = 101
mesh.xf   = [1.1, 1.5, 1.9]

# Assigns value based on attribute value
if mesh.type == "triangle"
    mesh.nn_el = 3
    mesh.nf_el = 3
end

# Loop through field names and fields: standard
for fname in fieldnames(typeof(mesh))
    println("Field name: ", fname)
    println("Content   : ", getfield(mesh, fname))
end

# Loop through struct fields: meta-programming
@generated function test4(p::P) where P
    assignments = [
        :( println("Field name: ", p.$name, "\n Content: ")  ) for name in fieldnames(P) #println(p.$name)
    ]
    quote $(assignments...) end
end
test4(mesh)


text:
"""
    Mie(m, n, nl_only)

The Mie generalized interaction.
When `m` equals 6 and `n` equals 12 this is equivalent to the Lennard Jones interaction.
"""
struct Mie{S, C, T} <: GeneralInteraction
    m::T
    n::T
    cutoff::C
    nl_only::Bool
    mn_fac::T
end

Mie{S}(m, n, cutoff, nl_only, mn_fac) where S =
    Mie{S, typeof(cutoff), typeof(m)}(m, n, cutoff, nl_only, mn_fac)

Mie(m, n, cutoff, nl_only, mn_fac) =
    Mie{false, typeof(cutoff), typeof(m)}(m, n, cutoff, nl_only, mn_fac)

Mie(m, n, nl_only=false) = Mie(m, n, ShiftedPotentialCutoff(3.0), nl_only,
                                convert(typeof(m), (n / (n - m)) * (n / m) ^ (m / (n - m))))

@fastmath @inbounds function force(inter::Mie{S, C, T},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where {S, C, T}
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)
    r = √r2

    if !S && iszero(atom_i.σ) || iszero(atom_j.σ)
        return zero(coord_i)
    end

    σ = sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    m = inter.m
    n = inter.n
    # Derivative obtained via wolfram
    const_mn = inter.mn_fac * ϵ
    σ_r = σ / r
    σ2 = σ^2
    params = (m, n, σ_r, const_mn)

    if cutoff_points(C) == 0
        f = force_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        r2 > sqdist_cutoff && return zero(coord_i)

        f = force_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        activation_dist = cutoff.activation_dist * σ2

        r2 > sqdist_cutoff && return zero(coord_i)

        if r2 < activation_dist
            f = force_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_cutoff(cutoff, r2, inter, params)
        end
    end

    return f * dr
end

@fastmath function force_nocutoff(::Mie, r2, invr2, (m, n, σ_r, const_mn))
    return -const_mn / r2 * (m * σ_r ^ m - n * σ_r ^ n)
end

@inline @inbounds function potential_energy(inter::Mie{S, C, T},
                                            s::Simulation,
                                            i::Integer,
                                            j::Integer) where {S, C, T}
    U = eltype(s.coords[i]) # this is not Unitful compatible
    i == j && return zero(T)

    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)
    r = √r2

    if !S && iszero(s.atoms[i].σ) || iszero(s.atoms[j].σ)
        return zero(U)
    end

    σ = sqrt(s.atoms[i].σ * s.atoms[j].σ)
    ϵ = sqrt(s.atoms[i].ϵ * s.atoms[j].ϵ)

    cutoff = inter.cutoff
    m = inter.m
    n = inter.n
    const_mn = inter.mn_fac * ϵ
    σ_r = σ / r
    σ2 = σ^2
    params = (m, n, σ_r, const_mn)

    if cutoff_points(C) == 0
        potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        r2 > sqdist_cutoff && return zero(U)

        potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > sqdist_cutoff && return zero(U)

        if r2 < activation_dist
            potential(inter, r2, inv(r2), params)
        else
            potential_cutoff(cutoff, r2, inter, params)
        end
    end
end

@fastmath function potential(::Mie, r2, invr2, (m, n, σ_r, const_mn))
    return const_mn * (σ_r ^ m - σ_r ^ m)
end


text:
using MyExample
using Test
my_f(2,1)
MyExample.greet2()
MyExample.my_f2(1,1)

@testset "MyExample.jl" begin
    @test MyExample.my_f(2,1)==5
    @test MyExample.my_f(2,3)==7
    @test my_f(2,3)=9
end

max_stars_repo_name: giandopal/MyExample

max_issues_repo_name: giandopal/MyExample

max_forks_repo_name: giandopal/MyExample


text:
export parallel_replicates

"""
A function to be used in `pmap` in `parallel_replicates`. It runs the `step!` function, but has a `dummyvar` parameter that does nothing, but is required for the `pmap` function.
"""
function parallel_step_dummy!(model::ABM, agent_step!, model_step!, n::Int, properties, when::AbstractArray{V}, step0::Bool, dummyvar) where {V<:Integer}
  data = step!(deepcopy(model), agent_step!, model_step!, n, properties, when=when, step0=step0);
  return data
end

"""
    parallel_replicates(agent_step!, model::ABM, n::Integer, agent_properties::Array{Symbol}, when::AbstractArray{Integer}, replicates::Integer)

Runs `replicates` number of simulations in parallel and returns a `DataFrame`.
"""
function parallel_replicates(model::ABM, agent_step!, model_step!, n::T, properties;
  when::AbstractArray{T}, replicates::T, step0::Bool) where {T<:Integer}

  all_data = pmap(j-> parallel_step_dummy!(model, agent_step!, model_step!, n,
  properties, when, step0, j), 1:replicates)

  dd = DataFrame()
  for (rep, d) in enumerate(all_data)
    d[!, :replicate] = [rep for i in 1:size(d, 1)]
    dd = vcat(dd, d)
  end

  return dd

end


text:
using Plots
using Distributions

# https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988/2
ENV["GKSwstype"] = "100"

mu = 3.14;
sigma = 1.0;
x = [mu - 5 * sigma:0.1:mu + 5 * sigma;];
N = Normal(mu, sigma^2);
p = plot(x, pdf.(N, x));
display(p)


text:
using ParetoRecipes
using Test

fourpoints = [[1, 0], [0, 1], [1, 1], [2, 2]]

@testset "pareto front" begin
    @test pareto_front(fourpoints) == [[0, 1], [1, 0]]
end

@testset "pareto plot" begin
    m = ParetoRecipes.Pareto(fourpoints)
    @test m.args == fourpoints
end


text:
"""
    GlmSpeedCallback(; glm_scale=0.5, cfl)

Update the divergence cleaning wave speed `c_h` according to the time step
computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations.
The `cfl` number should be set to the same value as for the time step size calculation. The
`glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD
solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0`
deactivates the divergence cleaning.
"""
struct GlmSpeedCallback{RealT<:Real}
  glm_scale::RealT
  cfl::RealT
end

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})
  @nospecialize cb # reduce precompilation time

  glm_speed_callback = cb.affect!
  @unpack glm_scale, cfl = glm_speed_callback
  print(io, "GlmSpeedCallback(glm_scale=", glm_scale, ", cfl=", cfl, ")")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:GlmSpeedCallback})
  @nospecialize cb # reduce precompilation time

  if get(io, :compact, false)
    show(io, cb)
  else
    glm_speed_callback = cb.affect!

    setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
            ]
    summary_box(io, "GlmSpeedCallback", setup)
  end
end

function GlmSpeedCallback(; glm_scale=0.5, cfl)

  @assert 0 <= glm_scale <= 1 "glm_scale must be between 0 and 1"

  glm_speed_callback = GlmSpeedCallback(glm_scale, cfl)

  DiscreteCallback(glm_speed_callback, glm_speed_callback, # the first one is the condition, the second the affect!
                   save_positions=(false,false),
                   initialize=initialize!)
end

function initialize!(cb::DiscreteCallback{Condition,Affect!}, u, t, integrator) where {Condition, Affect!<:GlmSpeedCallback}
  cb.affect!(integrator)
end

# this method is called to determine whether the callback should be activated
function (glm_speed_callback::GlmSpeedCallback)(u, t, integrator)
  return true
end

# This method is called as callback after the StepsizeCallback during the time integration.
@inline function (glm_speed_callback::GlmSpeedCallback)(integrator)

  dt = get_proposed_dt(integrator)
  semi = integrator.p
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  @unpack glm_scale, cfl = glm_speed_callback

  # compute time step for GLM linear advection equation with c_h=1 (redone due to the possible AMR)
  c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)

  # c_h is proportional to its own time step divided by the complete MHD time step
  equations.c_h = glm_scale * c_h_deltat / dt

  return nothing
end

include("glm_speed_dg.jl")


text:
@testset "Indexing" begin
    @testset "Int" begin
        leaf = Leaf(Tape(), 5 * [1, 1, 1, 1, 1])
        y = getindex(leaf, 1)
        @test unbox(y) == 5
        @test ∇(y, one(unbox(y)))[leaf] == [1, 0, 0, 0, 0]
    end

    @testset "Vector" begin
        x = Leaf(Tape(), 10 * [1, 1, 1])
        y = x[2:3]
        @test unbox(y) == [10, 10]
        @test ∇(y, oneslike(unbox(y)))[x] == [0, 1, 1]
    end

    @testset "Overlapping indices (#139)" begin
        x = Leaf(Tape(), 10 * [1, 1, 1])
        y = x[[2, 3, 3]]
        @test unbox(y) == [10, 10, 10]
        @test ∇(y, oneslike(unbox(y)))[x] == [0, 1, 2]
    end

    @testset "Ref" begin
        @test ref_equal(∇(getindex)(Ref(4))[1], Ref(1))
    end
end


text:
"""
"""
module Server # begin submodule PredictMD.Server

############################################################################
# PredictMD.Server source files ############################################
############################################################################

# submodules/Server/cryptography/
# include("cryptography/mbedtls.jl")
# include("cryptography/transportlayersecurity.jl")

# submodules/Server/web/
# include("web/genie.jl")
include("web/http.jl")
# include("web/httpclient.jl")
# include("web/httpserver.jl")
# include("web/juliawebapi.jl")
# include("web/mux.jl")
# include("web/websockets.jl")

end # end submodule PredictMD.Server


text:
include("double_integrator.jl")

## Double Integrator Example
mutable struct DoubleIntegrator{D,NM} <: DI{D}
    cfg::ForwardDiff.JacobianConfig{
        Nothing,
        Float64,
        NM,
        Tuple{
            Vector{ForwardDiff.Dual{Nothing,Float64,NM}},
            Vector{ForwardDiff.Dual{Nothing,Float64,NM}},
        },
    }
    cache::FiniteDiff.JacobianCache{
        Vector{Float64},
        Vector{Float64},
        Vector{Float64},
        UnitRange{Int64},
        Nothing,
        Val{:forward}(),
        Float64,
    }
    function DoubleIntegrator{D}() where {D}
        n = 2D
        m = D
        y = zeros(n)
        z = zeros(n + m)
        cfg = ForwardDiff.JacobianConfig(nothing, y, z)
        cache = FiniteDiff.JacobianCache(z, y)
        new{D,length(cfg.seeds)}(cfg, cache)
    end
end

# These methods are automatically defined by @autodiff
function RobotDynamics.jacobian!(::StaticReturn, ::ForwardAD, model::DoubleIntegrator, J, y, z)
    f(_z) = RD.evaluate(model, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    J .= ForwardDiff.jacobian(f, getdata(z))
end

function RobotDynamics.jacobian!(::InPlace, ::ForwardAD, model::DoubleIntegrator, J, y, z)
    f!(_y, _z) = RD.evaluate!(model, _y, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    J .= ForwardDiff.jacobian!(J, f!, y, getdata(z), model.cfg)
end

function RobotDynamics.jacobian!(::StaticReturn, ::FiniteDifference, model::DoubleIntegrator, J, y, z)
    f!(_y,_z) = _y .= RD.evaluate(model, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, RD.getdata(z), model.cache)
end

function RobotDynamics.jacobian!(::InPlace, ::FiniteDifference, model::DoubleIntegrator, J, y, z)
    f!(_y,_z) = RD.evaluate!(model, _y, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, RD.getdata(z), model.cache)
end

dim = 3
model = DoubleIntegrator{dim}()
@test RobotDynamics.state_dim(model) == 2dim
@test RobotDynamics.control_dim(model) == dim
@test RobotDynamics.output_dim(model) == 2dim
test_model(model)

## Model with auto-generated jacobian methods
@autodiff struct DoubleIntegratorAuto{D} <: DI{D} end
model = DoubleIntegratorAuto{dim}()
test_model(model)


text:
"Return the indices of u,v pairs belonging to the same natural variable and harmonic from hvars."
function _get_uv_pairs(hvars::Vector{HarmonicVariable})
    u_idx = findall(x -> x.type == "u", hvars)
    pairs = []
    for i in u_idx
        j = findall(x -> isequal(x.ω, hvars[i].ω) && x.type == "v" && isequal(hvars[i].natural_variable, x.natural_variable), hvars)
        j = length(j) != 1 && error("no v coordinate found for ", hvars[i]) || j[1]
        push!(pairs, [i, j])
    end
    pairs
end


text:
"""
    Hamming(n)

Create a Hamming window with n samples, where n is type Int.

# Example
```julia
julia> w = Hamming(51); plot(w);
```
"""
function Hamming(n::Int)

    n==1 && return [1.]
    n -= 1
    k = collect(0:1:n)
    hamm = 0.54 .- 0.46 .* cos.(2pi.*k/n)

end


text:
include("./src/ShellModel.jl")
using .ShellModel

# similar to samplerun()
function run()
    ### for compilation
    println("For JIT compilation...")
    @time main("./snts/x_mass.snt","Be8",5,[0])
    @time main("./snts/x_mass.snt","Be8",5,[0];is_block=true,q=3)
    println("Done: JIT compilation")
    ### input
    tJ = 0 # target J
    tJs =[ tJ ]
    num_ev = 4; n_block = 4
    num_ev_target = 10
    sntf = "./snts/usdb.snt"; target_nuc = "Si28"
    ###

    println("\n### sample (a). 10-lowest states for $target_nuc")
    @time main(sntf,target_nuc,10,[])

    println("\n### sample (b). 10-lowest states with J=$tJs")
    @time main(sntf,target_nuc,10,tJs)

    println("\n### sample (c). calc. EC estimates for $num_ev_target J=$tJs states ")
    Hs = ["./snts/usdb.snt"]
    tJNs = [ [tJ,num_ev_target] ]
    solveEC(Hs,target_nuc,tJNs;
            verbose=true, exact_logf = "runlog_sigma1_valid.dat")

    println("\n### sample (d).  (b). with the preprocessing")
    println("Block Lanczos method with n_block=$n_block")
    @time main(sntf,target_nuc,num_ev,tJs;is_block=true,q=n_block)
    println("\n+ preprocessing")
    @time main(sntf,target_nuc,num_ev,tJs;is_block=true,q=n_block,
               in_wf="./appwavs/"*target_nuc*"_usdb_j"*string(tJ)*".appwav")

end
run()


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8efe4e7f-a63c-4a73-bea3-48c02b292d12
begin
    import Pkg
    # careful: this is _not_ a reproducible environment
    # activate the global environment
    Pkg.activate()

    using Plots, PlutoUI
	import HDF5
	import Pda
	using Statistics
end

# ╔═╡ e5f8844e-7a41-4ac6-a843-9904309f9fa8
file = "U10.000singlewallDT1.000.h5"

# ╔═╡ ad67db84-c555-438d-a4d1-aae9f81d3981
@bind frame PlutoUI.Slider(0:98, show_value=true)

# ╔═╡ f1c78cf0-a784-4273-a6c6-090f64cfed36
begin
	x = Pda.get_h5data(file, "config/$frame/x")
	theta = Pda.get_h5data(file, "config/$frame/theta")
end

# ╔═╡ be3457aa-a2a0-43e6-b240-fa0e076d396c
function compute_inf_linedensity(file; frame::Int=80, )
		right = Pda.get_h5attr(file, "right")
		threshold = right/2.0
		Lx = right-threshold
		xin = [xi for xi in x if xi > threshold]
		ninf = length(xin) / Lx
	return ninf
end

# ╔═╡ c6cd5b84-d00d-4898-a769-e9ebd819ac93
begin
	U0 = Pda.get_h5attr(file, "U0")
	tauR = Pda.get_h5attr(file, "tauR")
	ell = U0*tauR
	DT = Pda.get_h5attr(file, "DT")
	delta = sqrt(DT*tauR)
	# Lx = 20 # box length in x
	# # bulk line density, take only the latter half of the box
	# ninf = sum( x .>= Lx/2)/(Lx/2)
	# # number of particles on the wall
	# Nw = 2*ninf*ell
end

# ╔═╡ 09a21a89-e9a0-4236-88f8-d954d4165a55
sum( x .<= eps(Float64) ) /(compute_inf_linedensity(file)*ell)

# ╔═╡ ef496cd8-2a47-4c1e-9aac-c4fdbfd98ff6
begin
	qx_wall = Float64[]
	qy_wall = Float64[]
	theta_wall = Float64[]
	N = convert(Int, Pda.get_h5attr(file, "N") )
	for i in 1:N
		if x[i]  <= eps(Float64)
			push!(qx_wall, cos(theta[i]) )
			push!(qy_wall, sin(theta[i]) )
			push!(theta_wall, mod(theta[i], 2pi))
		end
	end
	h1 = histogram(qx_wall, label = "qx")
	h2 = histogram(qy_wall, label = "qy")
	plot(h1, h2, framestyle=:box, layout = (2,1) )
end

# ╔═╡ 65f836ea-1839-4c63-889c-9ad062ccf6b5
begin
	histogram(theta_wall, bins=50, framestyle=:box, xticks = pi*(0:1/2:2))
	let
		theta = LinRange(pi/2, 3pi/2, 100)
		lee = [
		0.47883923959186386 0.09804957410222426;
		0.5198535717065695 0.2995920825021341;
		0.5576166579070441 0.5216095918830788;
		0.5936994405206624 0.6869355109518869;
		0.6330218896210373 0.8176132098108629;
		0.6706910537986635 0.9246701632248963;
		0.7116140370966255 1.014401719930961;
		0.7476234833362387 1.0899639172173865;
		0.7476234833362387 1.0899639172173865;
		0.788528454191463 1.1576482440119886;
		0.8277930061544658 1.2174598467269784;
		0.8277930061544658 1.2174598467269784;
		0.8654081330038772 1.2583751104066245;
		0.8654081330038772 1.2583751104066245;
		0.9013789811519084 1.286693243597059;
		0.9406126546416456 1.308709595035256;
		0.9765719233622026 1.3228545089968935;
		1.0239549405878874 1.3196675932352986;
		1.0549903794256732 1.3070447306850326;
		1.0925540421529758 1.2849679089033574;
		1.1301151316741729 1.2597414828486164;
		1.174194152260465 1.212462680470202;
		1.2084589647605855 1.1525931806178398;
		1.2492545743563286 1.0864193258071344;
		1.2884084784567968 1.0107979447802835;
		1.2884084784567968 1.0107979447802835;
		1.325912957443674 0.9162802247180892;
		1.325912957443674 0.9162802247180892;
		1.3634058570030767 0.8075892854270975;
		1.3634058570030767 0.8075892854270975;
		1.4008846039288998 0.6815755226342425;
		1.4008846039288998 0.6815755226342425;
		1.4416017307384266 0.5193387374950225;
		1.487170637659772 0.29568080922185613;
		1.5196843834051366 0.092505601548041;
		1.5196843834051366 0.092505601548041];
		#scatter!(lee[:,1]*pi, lee[:,2]*210, )
		plot!(theta, 40*cos.(theta).^2 )
		# This simple theory doesnt work!???? why
	end
end

# ╔═╡ f9faa467-6f14-46e6-84b5-fb9762214e00
let

scatter(lee[:,1]*pi, lee[:,2], )
theta = LinRange(0.5*pi, 1.5*pi,200)
plot!(theta, cos.(theta).^2*1.35)
end

# ╔═╡ e47394f4-c336-4937-aa0c-aa174bf409b4
 sum( x .> sqrt(eps(Float64)))

# ╔═╡ cffe8f9b-b7ab-448b-bf6c-c2123bdc2616
begin
	dx = Pda.get_h5data(file, "dx/m")
	dxtotal =mean( dx*N) # average over frames
	dt = Pda.get_h5attr(file, "dt")
	ninf = compute_inf_linedensity(file)
	# 1st scale by ninf kbt
	# 2nd scale by ninf ksts, for DT=0 this should be unity
	dxtotal/(dt*ninf*DT), dxtotal/(dt*ninf*0.5*U0^2*tauR)
end

# ╔═╡ eaea4043-ad2e-4194-b2d7-52a21012ee4c
begin
 	1 + 0.5*(ell/delta)^2
end

# ╔═╡ 089cef24-5c87-4a59-b23e-cb7b952cd3d3
let
	y = [xi for xi in x if xi > 0.]

end

# ╔═╡ Cell order:
# ╠═8efe4e7f-a63c-4a73-bea3-48c02b292d12
# ╠═e5f8844e-7a41-4ac6-a843-9904309f9fa8
# ╠═ad67db84-c555-438d-a4d1-aae9f81d3981
# ╠═f1c78cf0-a784-4273-a6c6-090f64cfed36
# ╟─be3457aa-a2a0-43e6-b240-fa0e076d396c
# ╟─c6cd5b84-d00d-4898-a769-e9ebd819ac93
# ╠═09a21a89-e9a0-4236-88f8-d954d4165a55
# ╠═ef496cd8-2a47-4c1e-9aac-c4fdbfd98ff6
# ╠═65f836ea-1839-4c63-889c-9ad062ccf6b5
# ╠═f9faa467-6f14-46e6-84b5-fb9762214e00
# ╠═e47394f4-c336-4937-aa0c-aa174bf409b4
# ╠═cffe8f9b-b7ab-448b-bf6c-c2123bdc2616
# ╠═eaea4043-ad2e-4194-b2d7-52a21012ee4c
# ╠═089cef24-5c87-4a59-b23e-cb7b952cd3d3


text:
# Sample postprocessing script for customized figure.
#
# Reading 3D simulation data and plotting density in 3 cuts in the same figure.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot, PyCall
axes_grid1 = pyimport("mpl_toolkits.axes_grid1")
ImageGrid = axes_grid1.ImageGrid

file = "bulk1.0001000.vlsv"
nameρ = "proton/vg_rho"
colorscale = Log
addcolorbar = false

meta = load(file)

# normal cuts in the x,y,z directions
fig = plt.figure(figsize=(12, 4))
grid = ImageGrid(fig, 111,
                 nrows_ncols=(1, 3),
                 axes_pad=0.85,
                 cbar_mode="single",
                 cbar_location="right",
                 cbar_pad=0.1,
                 label_mode="all"
                 )

c1 = pcolormesh(meta, nameρ, grid[1]; normal=:x, addcolorbar, colorscale)
c2 = pcolormesh(meta, nameρ, grid[2]; normal=:y, addcolorbar, colorscale)
c3 = pcolormesh(meta, nameρ, grid[3]; normal=:z, addcolorbar, colorscale)

cb = fig.colorbar(c3, cax=grid.cbar_axes[1])
datainfo = readvariablemeta(meta, nameρ)

cb_title_str = datainfo.variableLaTeX
cb_title_str *= ",["*datainfo.unitLaTeX*"]"
cb_title = cb.ax.set_title(cb_title_str, fontsize=14, fontweight="bold")

plt.savefig("test.png",bbox_inches="tight")


text:
function Base.show(io::IO, wb::WorkBook)
    print(io, basename(wb.path), " (with ", length(wb), " WorkSheet)")
end

function Base.show(io::IO, ws::WorkSheet)
    print("WorkSheet[$(ws.name)] ")
    display(ws.data)
end

function Base.show(io::IO, x::Cell)
    print(io, x.v)
end


text:
#///////////////////////////////////////
#// File Name: crowd_nav_controller.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/04/02
#// Description: Controller based on CrowdNav RL Policies
#///////////////////////////////////////

struct CrowdNavControlParameter <: Parameter
    model_dir::String # Directory of the trained policy
    env_config::String # environment config file name
    policy_config::String # policy config file name
    policy_name::String # policy name
    goal_pos::Vector{Float64} # [x, y] goal position
    tcalc::Float64 # Allocated control computation time
    dtr::Float64 # Replanning time interval
    target_speed::Float64 # target speed for robot
end

mutable struct CrowdNavController
    sim_param::SimulationParameter
    cnt_param::CrowdNavControlParameter

    rl_robot::PyObject

    tcalc_actual::Union{Nothing, Float64}
    u_value::Union{Nothing, Vector{Float64}}

    control_update_task::Union{Nothing, Task}

    tcalc_actual_tmp::Union{Nothing, Float64}
    u_value_tmp::Union{Nothing, Vector{Float64}}
end

function CrowdNavController(sim_param::SimulationParameter,
                            cnt_param::CrowdNavControlParameter)
    rl_robot = py"configure_rl_robot"(cnt_param.model_dir, cnt_param.env_config,
                                      cnt_param.policy_config, cnt_param.policy_name);
    @assert rl_robot.visible == false "RL Robot has to be invisible."
    @assert rl_robot.time_step == sim_param.dto "RL Robot has to use the same timestep as sim_param.dto."
    rl_robot.v_pref = cnt_param.target_speed;
    rl_robot.gx, rl_robot.gy = cnt_param.goal_pos[1], cnt_param.goal_pos[2];
    rl_robot.theta = 0.0
    return CrowdNavController(sim_param, cnt_param, rl_robot, nothing, nothing,
                              nothing, nothing, nothing);
end

function get_action!(controller::CrowdNavController,
                     e_init::RobotState,
                     ado_state_dict::Dict{T, Vector{Float64}} where T <: Union{PyObject, String})
    # update rl_robot state
    controller.rl_robot.set_position(copy(get_position(e_init)));
    controller.rl_robot.set_velocity(copy(get_velocity(e_init)));

    # ado_state_dict has to have velocity information
    ado_state_array = [];
    radius::Float64 = controller.rl_robot.radius
    if length(ado_state_dict) < 2 # crowdnav seems to need >1 ado agent
        vec = controller.cnt_param.goal_pos - get_position(e_init);
        action = vec./norm(vec).*controller.rl_robot.v_pref;
        return action
    else
        for key in keys(ado_state_dict)
            @assert length(ado_state_dict[key]) >= 4 "ado_state_dict has to have velocity information."
            px, py = ado_state_dict[key][1], ado_state_dict[key][2];
            vx, vy = ado_state_dict[key][3], ado_state_dict[key][4];
            observable_state = py"ObservableState"(px, py, vx, vy, radius);
            push!(ado_state_array, observable_state)
        end
        action = controller.rl_robot.act(ado_state_array);
        return action
    end
end

@inline function crowdnav_control_update!(controller::CrowdNavController,
                                          e_init::RobotState,
                                          ado_state_dict::Dict{T, Vector{Float64}} where T <: Union{PyObject, String})
    tcalc_actual =
        @elapsed u_vel = get_action!(controller, e_init, ado_state_dict);
    tcalc_actual +=
        @elapsed u_acc = (u_vel .- get_velocity(e_init)) ./ controller.sim_param.dto;
        #@elapsed u_acc = 2 * (u_vel .- get_velocity(e_init)) ./ controller.sim_param.dto;

    if tcalc_actual >= controller.cnt_param.tcalc
        # tcalc_actual has exceeded allowable computation time
        time = @sprintf "Time %.2f" round(to_sec(e_init.t), digits=5)
        @warn "$(time) [sec]: CrowdNav computation took $(round(tcalc_actual, digits=3)) [sec], which exceeds the maximum computation time allowed."
    end

    return tcalc_actual, u_acc
end

function control!(controller::CrowdNavController,
                  current_time::Time;
                  log::Union{Nothing, Vector{Tuple{Time, String}}}=nothing)

    if !isnothing(controller.control_update_task) &&
        istaskdone(controller.control_update_task)
        if !isnothing(log)
            msg = "New CrowdNav control is available to the controller."
            push!(log, (current_time, msg))
        end
        controller.tcalc_actual = controller.tcalc_actual_tmp;
        controller.u_value = copy(controller.u_value_tmp);
        controller.control_update_task = nothing;
    end
    u = copy(controller.u_value);
    if !isnothing(log)
        msg = "Control: $(u) is applied to the system."
        push!(log, (current_time, msg))
    end
    return u
end

function schedule_control_update!(controller::CrowdNavController,
                                  e_init::RobotState,
                                  ado_state_dict::Dict{T, Vector{Float64}} where T <: Union{PyObject, String};
                                  log::Union{Nothing, Vector{Tuple{Time, String}}}=nothing)
    if !isnothing(log)
        msg = "New CrowdNav control computation is scheduled."
        push!(log, (e_init.t, msg))
    end

    controller.control_update_task = @task begin
        controller.tcalc_actual_tmp, controller.u_value_tmp =
            crowdnav_control_update!(controller, e_init, ado_state_dict);
    end
    schedule(controller.control_update_task);
end


text:
using Mocking
Mocking.activate()

using Dates: Dates, Hour
using Test

using Mocking: apply, splitdef, combinedef
using Mocking: anon_morespecific, anonymous_signature, dispatch, type_morespecific

@testset "Mocking" begin
    include("expr.jl")
    include("dispatch.jl")
    include("patch.jl")

    include("concept.jl")
    include("targets.jl")
    include("functions.jl")
    include("import.jl")
    include("real-open.jl")
    include("real-isfile.jl")
    include("real-nested.jl")
    include("mock-in-patch.jl")
    include("readme.jl")
    include("optional.jl")
    include("patch-gen.jl")
    include("anonymous-param.jl")
    include("reuse.jl")
    include("args.jl")
end


text:
function evaluate_decision(structure::HorizontalStructure, decision::AbstractVector)
    return _eval_subproblems(structure, decision)
end

function statistically_evaluate_decision(structure::HorizontalStructure, decision::AbstractVector)
    𝔼Q, σ² = _stat_eval_subproblems(structure, decision)
    return 𝔼Q, sqrt(σ²)
end

function _eval_subproblems(structure::HorizontalStructure{2,1,Tuple{SP}},
                           decision::AbstractVector) where SP <: ScenarioProblems
    # Update decisions
    take_decisions!(structure.decisions[1], decision)
    map(subprob -> update_decisions!(subprob), subproblems(structure))
    # Cache result
    result = outcome_mean(subproblems(structure), probability.(scenarios(structure)))
    # Revert back to untaken decisions
    untake_decisions!(structure)
    # Return evaluation result
    return result
end
function _eval_subproblems(structure::HorizontalStructure{2,1,Tuple{SP}},
                           decision::AbstractVector) where SP <: DistributedScenarioProblems
    Qs = Vector{Float64}(undef, nworkers())
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async Qs[i] = remotecall_fetch(
                w,
                sp[w-1],
                sp.decisions[w-1],
                decision) do sp, d, x
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0
                    take_decisions!(fetch(d), x)
                    map(subprob -> update_decisions!(subprob), subproblems(scenarioproblems))
                    return outcome_mean(subproblems(scenarioproblems),
                                        probability.(scenarios(scenarioproblems)))
                end
        end
    end
    return sum(Qs)
end

function _stat_eval_subproblems(structure::HorizontalStructure{2,1,Tuple{SP}},
                                decision::AbstractVector) where SP <: ScenarioProblems
    # Update decisions
    take_decisions!(structure.decisions[1], decision)
    map(subprob -> update_decisions!(subprob), subproblems(structure))
    return welford(subproblems(structure), probability.(scenarios(structure)))
end
function _stat_eval_subproblems(structure::HorizontalStructure{2,1,Tuple{SP}},
                                decision::AbstractVector) where SP <: DistributedScenarioProblems
    partial_welfords = Vector{Tuple{Float64,Float64,Float64,Int}}(undef, nworkers())
    sp = scenarioproblems(structure)
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch(
                w,
                sp[w-1],
                sp.decisions[w-1],
                x) do sp, d, x
                    scenarioproblems = fetch(sp)
                    num_scenarios(scenarioproblems) == 0 && return 0.0, 0.0, 0.0, 0
                    take_decisions!(fetch(d), x)
                    map(subprob -> update_decisions!(subprob), subproblems(scenarioproblems))
                    return welford(subproblems(scenarioproblems),
                                   probability.(scenarios(scenarioproblems)))
                end
        end
    end
    𝔼Q, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼Q, σ²
end


text:
"""
    (dist, r, z, d) = powerline_distances(catenary, point, [swing = 0])
    (dist, r, z, d) = powerline_distances(point, catenary, [swing = 0])

Distance between the given point and the closest point on the specified catenary. If
`swing` is specified, it returns the distance of closest approach as the catenary is swung
across an arc of `±swing` from vertical. The `swing` angle is specified in radians.

 * `dist` is 3D distance to closest point on catenary
 * `r` is a horizontal coordinate (positive or negative) transverse to the catenary
 * `z` is the vertical distance from the point on the catenary at given `d` (or nearest end)
 * `d` is a horizontal coordiante longitudanal to the catenary.
"""
function powerline_distances(cat::Catenary, p::AbstractVector{T}, swing = zero(T)) where {T <: Real}
    if length(p) !== 3
        error("Expected 3d point")
    end
    if swing != 0
        error("Swing is not implemented yet")
    end

    inv_a = inv(cat.a)

    # Point transformed into the frame of the catenary
    p_trans = cat.transform(p)

    # The beginning and end of the catenary
    p_begin = SVector(cat.lmin, 0, cat.a * (cosh(cat.lmin * inv_a) - 1))
    p_end = SVector(cat.lmax, 0, cat.a * (cosh(cat.lmax * inv_a) - 1))

    # Find coordinates relative to line
    @inbounds d = p_trans[1]
    @inbounds r = p_trans[2]
    @inbounds z = p_trans[3] - cat.a * (cosh(d * inv_a) - 1)

    # Determine whether we are closest in 3D to one of the ends, or the body of the catenary
    @inbounds if d > (p_end[3] - p_trans[3]) * sinh(cat.lmax * inv_a) + cat.lmax
        dist = norm(p_trans - p_end)
    elseif d < (p_begin[3] - p_trans[3]) * sinh(cat.lmin * inv_a) + cat.lmin
        dist = norm(p_trans - p_begin)
    else
        # With a local gradient approximation, find the distance to the closest point on the
        # line in the d-z plane. Effectively performs the first Newton-Raphson step.
        sinh_d_a = sinh(d * inv_a)
        d_perp = d + z * sinh_d_a / (1 + sinh_d_a*sinh_d_a)
        Δz = p_trans[3] - cat.a * (cosh(d_perp * inv_a) - 1) # recalculating z at new d helps accuracy signficantly
        Δd = d_perp - p_trans[1]

        dist = sqrt(r*r + Δz*Δz + Δd*Δd)
    end

    # output z depends on whether outside catenary region or not
    @inbounds if d > p_end[1] # beyond end
        z = p_trans[3] - p_end[3]
    elseif d < p_begin[1] # beyond begin
        z = p_trans[3] - p_begin[3]
    end

    return (dist, r, z, d)
end

powerline_distances(p::AbstractVector{<:Real}, cat::Catenary) = powerline_distances(cat, p)

"""
    p = closest_point(geom1, geom2)

Return the point on `geom1` closest to `geom2`.
"""
closest_point(a, b) = closest_points(a, b)[1]

"""
    (p1, p2) = closest_points(geom1, geom2)

Return the two closest points between two geometries.
"""
closest_points(p1::AbstractVector{<:Real}, p2::AbstractVector{<:Real}) = (p1, p2)

function closest_points(line::Line, p::AbstractVector{T}) where {T <: Real}
    d = line.p2 - line.p1
    sd = sum(abs2, d)
    closest = sd > 0 ? line.p1 + d * clamp(dot(p - line.p1, d) / sd, zero(T), one(T)) : line.p1
    return (closest, p)
end
closest_points(p::AbstractVector{<:Real}, line::Line) = reverse(closest_points(line, p))

function closest_points(l1::Line{2}, l2::Line{2})
    # Parameterize each line as pi + ti*di
    p1 = l1.p1
    p2 = l2.p1
    d1 = l1.p2 - l1.p1
    d2 = l2.p2 - l2.p1

    denom = d1 × d2 # 2D cross product from StaticArrays - returns scalar

    # Check that denominator is reasonable - otherwise assume parallel lines
    if abs(denom) > 10*eps(norm(d1)*norm(d2))
        r = p1 - p2
        t1 = (d2 × r) / denom
        t2 = (d1 × r) / denom

        # The segments cross iff the ti are in 0..1
        if t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1
            return (p1 + t1*d1, p2 + t2*d2)
        end
    end

    # Otherwise, one of the end points is closest - but there is no hint as to which
    (out1, out2) = closest_points(l1.p1, l2)
    dist = distance(out1, out2)

    (tmp1, tmp2) = closest_points(l1.p2, l2)
    dist_tmp = distance(tmp1, tmp2)
    if dist_tmp < dist
        dist = dist_tmp
        out1 = tmp1
        out2 = tmp2
    end

    (tmp1, tmp2) = closest_points(l1, l2.p1)
    dist_tmp = distance(tmp1, tmp2)
    if dist_tmp < dist
        dist = dist_tmp
        out1 = tmp1
        out2 = tmp2
    end

    (tmp1, tmp2) = closest_points(l1, l2.p2)
    dist_tmp = distance(tmp1, tmp2)
    if dist_tmp < dist
        dist = dist_tmp
        out1 = tmp1
        out2 = tmp2
    end

    return (out1, out2)
end

function closest_points(l1::Line{3}, l2::Line{3})
    # Parameterize each line as pi + ti*di
    p1 = l1.p1
    p2 = l2.p1
    d1 = l1.p2 - l1.p1
    d2 = l2.p2 - l2.p1

    denom = d1 × d2 # 2D cross product from StaticArrays - returns scalar

    # The line connecting points of closest approach is perendicular to both d1 and d2
    n = d1 × d2

    # Check that the cross product is reasonable - otherwise assume parallel lines
    if norm(n) > 10*eps(norm(d1)*norm(d2))
        r = p1 - p2
        t1 = -(d2 ⋅ (n × r)) / (d2 ⋅ (n × d1))
        t2 =  (d1 ⋅ (n × r)) / (d1 ⋅ (n × d2))

        # The segments reach closest approach (of infinite lines) iff the ti are in 0..1
        if t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1
            return (p1 + t1*d1, p2 + t2*d2)
        end
    end

    # Otherwise, one of the end points is closest - but there is no hint as to which
    (out1, out2) = closest_points(l1.p1, l2)
    dist = distance(out1, out2)

    (tmp1, tmp2) = closest_points(l1.p2, l2)
    dist_tmp = distance(tmp1, tmp2)
    if dist_tmp < dist
        dist = dist_tmp
        out1 = tmp1
        out2 = tmp2
    end

    (tmp1, tmp2) = closest_points(l1, l2.p1)
    dist_tmp = distance(tmp1, tmp2)
    if dist_tmp < dist
        dist = dist_tmp
        out1 = tmp1
        out2 = tmp2
    end

    (tmp1, tmp2) = closest_points(l1, l2.p2)
    dist_tmp = distance(tmp1, tmp2)
    if dist_tmp < dist
        dist = dist_tmp
        out1 = tmp1
        out2 = tmp2
    end

    return (out1, out2)
end

"""
    distance(geom1, geom2)

Calculate the minimum Euclidean distance between two geometries.
"""
@inline distance(p1::AbstractVector{<:Real}, p2::AbstractVector{<:Real}) = norm(p1 - p2)

@inline function distance(p1::StaticVector{n, <:Real}, p2::StaticVector{n, <:Real}) where {n}
    norm(p1 - p2)
end
@inline function distance(p1::StaticVector{2, <:Real}, p2::StaticVector{3, <:Real})
    norm(p1 - convert2d(p2))
end
@inline function distance(p1::StaticVector{3, <:Real}, p2::StaticVector{2, <:Real})
    norm(convert2d(p1) - p2)
end

# Distance between a point and the closest point on a Line
function distance(line::Line, p::AbstractVector{<:Real})
    return distance(p, closest_point(line, p))
end
distance(line::Line{2}, p::StaticVector{3, <:Real}) = distance(line, convert2d(p))
distance(line::Line{3}, p::StaticVector{2, <:Real}) = distance(convert2d(line), p)
distance(p::AbstractVector{<:Real}, line::Line) = distance(line, p)

# Shortest distance between two Lines
function distance(l1::Line{2}, l2::Line{2})
    # Parameterize each line as pi + ti*di
    p1 = l1.p1
    p2 = l2.p1
    d1 = l1.p2 - l1.p1
    d2 = l2.p2 - l2.p1

    denom = d1 × d2 # 2D cross product from StaticArrays - returns scalar

    # Check that denominator is reasonable - otherwise assume parallel lines
    if abs(denom) > 10*eps(norm(d1)*norm(d2))
        r = p1 - p2
        t1 = (d2 × r) / denom
        t2 = (d1 × r) / denom

        # The segments cross iff the ti are in 0..1
        if t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1
            return zero(t1)
        end
    end

    # Otherwise, one of the end points is closest - but there is no hint as to which
    dist = distance(l1.p1, l2)
    dist = min(dist, distance(l1.p2, l2))
    dist = min(dist, distance(l1, l2.p1))
    dist = min(dist, distance(l1, l2.p2))

    return dist
end

function distance(l1::Line{3}, l2::Line{3})
    # Parameterize each line as pi + ti*di
    p1 = l1.p1
    p2 = l2.p1
    d1 = l1.p2 - l1.p1
    d2 = l2.p2 - l2.p1

    # The line connecting points of closest approach is perendicular to both d1 and d2
    n = d1 × d2

    # Check that the cross product is reasonable - otherwise assume parallel lines
    if norm(n) > 10*eps(norm(d1)*norm(d2))
        r = p1 - p2
        t1 = -(d2 ⋅ (n × r)) / (d2 ⋅ (n × d1))
        t2 =  (d1 ⋅ (n × r)) / (d1 ⋅ (n × d2))

        # The segments reach closest approach (of infinite lines) iff the ti are in 0..1
        if t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1
            return n ⋅ r / norm(n)
        end
    end

    # Otherwise, one of the end points is closest - but there is no hint as to which
    dist = distance(l1.p1, l2)
    dist = min(dist, distance(l1.p2, l2))
    dist = min(dist, distance(l1, l2.p1))
    dist = min(dist, distance(l1, l2.p2))

    return dist
end

distance(l1::Line{2}, l2::Line{3}) = distance(l1, convert2d(l2))
distance(l1::Line{3}, l2::Line{2}) = distance(convert2d(l1), l2)

# Distance between a point and the closest point on a LineString
function distance(ls::LineString, p::AbstractVector{<:Real})
    mapreduce(l -> distance(l, p), min, ls)
end
distance(p::AbstractVector{<:Real}, ls::LineString) = distance(ls, p)

function distance(polygon::Polygon, p::AbstractVector{T}) where {T <: Real}
    if p in polygon
        return zero(T)
    end

    dist = distance(polygon.exterior, p)
    for ls in polygon.interiors
        dist = min(dist, distance(ls, p))
    end

    return dist
end
distance(p::AbstractVector{<:Real}, polygon::Polygon) = distance(polygon, p)

# Distance between a point and the closest point on a Catenary
distance(p::AbstractVector{<:Real}, cat::Catenary) = powerline_distances(cat, p)[1]
distance(cat::Catenary, p::AbstractVector{<:Real}) = powerline_distances(cat, p)[1]


text:
struct GhostConditions

    Proj :: SparseMatrixCSC
    Extr :: SparseMatrixCSC
    swaps :: Matrix{Float64}

end

function GhostConditions(mx, my, MinvT, bars::Vector{Barrier}, grd::Grid; w=0.5)

    ks = [Int64[] for _=1:length(bars)+1]
    js = [Int64[] for _=1:length(bars)+1]
    dats = [Float64[] for _=1:length(bars)+1]
    swaps = ones(Float64, (grd.Nk, length(bars)))
    ng = 0
    dr = 1.5*sqrt(grd.dx^2 + grd.dy^2)

    for k=1:grd.Nk
        us = [smoothstep(grd.points[:,k]..., 1.0, bar) for bar in bars]
        if !isempty(us[us .< 0.5])

            us[us .> 0.5] .= -Inf
            l = sortperm(us)[end]
            ng += 1

            # add permutation row
            ks[l] = [ks[l]; k]

            # compute reflection
            x, y = bars[l].rmap(grd.points[:,k]...)

            # for robustness, compute coefficients in terms of neighboring points.
            dist = Float64[x, y] - grd.points[:,k]
            er = dist / norm(dist)
            x1 = x + dr*er[1]
            x2 = x + 2*dr*er[1]
            y1 = y + dr*er[2]
            y2 = y + 2*dr*er[2]

            # add reflection rows
            row_dat1, row_j1 = interpolation_row(x1, y1, mx, my, MinvT, grd)
            row_dat2, row_j2 = interpolation_row(x2, y2, mx, my, MinvT, grd)
            js[l] = [js[l]; row_j1]
            js[l] = [js[l]; row_j2]
            dats[l] = [dats[l]; 2*row_dat1]
            dats[l] = [dats[l]; -1*row_dat2]

            # add identity rows
            row_dat, row_j = interpolation_row(grd.points[:,k]..., mx, my, MinvT, grd)
            js[l] = [js[l]; row_j]
            dats[l] = [dats[l]; row_dat]

            # change swaps
            swaps[k,l] = -1.0

        else

            # add permutation row
            ks[end] = [ks[end]; k]
        end
    end

    # compute permutation and reflection condition matrix
    Perm = sparse([1:grd.Nk...], vcat(ks...), ones(grd.Nk), grd.Nk, grd.Nk)
    M = sparse(vcat([k*ones(3*mx*my) for k=1:ng]...), vcat(js...), vcat(dats...), ng, grd._Nx*grd._Ny) * transpose(grd.Proj) * transpose(Perm)

    # invert to get extrapolation operation, delete small values
    N = inv(Matrix(M[1:ng, 1:ng]))
    N[abs.(N) .< 1e-8] .= 0
    R = -sparse(N) * M[:, ng+1:end]

    # make extrapolation matrix
    Extr_perm = sparse(Diagonal(ones(Float64, grd.Nk)))
    # Extr_perm = sparse([1:grd.Nk]..., [1:grd.Nk]..., ones(Float64, grd.Nk), grd.Nk, grd.Nk) # XXX
    Extr_perm[1:ng, ng+1:end] = R
    Extr_perm[1:ng, 1:ng] .= 0
    Extr = transpose(Perm) * Extr_perm * Perm

    # compute projection
    Proj = sparse([1:grd.Nk-ng...], ks[end], ones(grd.Nk-ng), grd.Nk-ng, grd.Nk)

    return GhostConditions(Proj, Extr, swaps)
end

function swap_sign(gc::GhostConditions, inds...)

    Extr = gc.Extr

    for i in inds
        Extr = Diagonal(gc.swaps[:,i]) * Extr
    end

    return GhostConditions(gc.Proj, Extr, gc.swaps)
end

function extrapolate_ghosts(x::Vector{Float64}, xb::Vector{Float64}, gc::GhostConditions)

    return gc.Extr*transpose(gc.Proj)*x + (I - gc.Extr)*xb
end

function require_boundary_conditions(A::SparseMatrixCSC, gc::GhostConditions)

    P = gc.Proj
    Pt = transpose(P)

    Anew = P*A*gc.Extr*Pt
    Bterm = P*A*(I - gc.Extr)

    return Anew, Bterm
end


text:
module EspyInsideFunction
export @request
export makekey,forloop,scalar
export @espy,@espydbg

# manipulation of data-structures representing Julia expressions (::Expr)
using Printf

## Generate unique symbol
newsym(name) = Symbol(name,"_",string(gensym())[3:end])

## Tests
issymbol(e::Symbol)         = true
issymbol(e)                 = false
isquote(e::QuoteNode)       = true
isquote(e)                  = false
isnumber(e::Number)         = true
isnumber(e)                 = false
isline(e::LineNumberNode)   = true
isline(e)                   = false
isexpr(e::Expr)             = true
isexpr(e,f...)              = false
isexpr(s,e::Expr)           = e.head == s

istuple(e)                  = isexpr(e) && e.head == :tuple  # (a,b)
isref(e)                    = isexpr(e) && e.head == :ref    # a[b]
isdot(e)                    = isexpr(e) && e.head == :(.)    # a.b
issub(e)                    = isexpr(e) && e.head == :(.)  && isquote(getright(e))  # a.b
isdeclare(e)                = isexpr(e) && e.head == :(::)   # a::b
isassign(e)                 = isexpr(e) && e.head == :(=)
iscall(e)                   = isexpr(e) && e.head == :call
isfunc(e)                   = isexpr(e) && e.head == :function
isfor(e)                    = isexpr(e) && e.head == :for
isescsymbol(e)              = isexpr(e) && e.head == (:escape) && issymbol(getleft(e))
issymbolish(e)              = issymbol(e)||isescsymbol(e)

## Analysis
getargs(e)                  = (e.args...,)
getright(e)                 = e.args[2]
getleft(e)                  = e.args[1]
getfor(e)                   = (var=e.args[1].args[1],from=e.args[1].args[2].args[2],to=e.args[1].args[2].args[3],body=e.args[2]) #(var,from,to,body)

## Synthesis
maketuple(e...)             = Expr(:tuple,          e...)
makeref(e...)               = Expr(:ref,            e...)
makecall(e...)              = Expr(:call,           e...)
makeblock(e...)             = Expr(:block,          e...)
makewhere(e...)             = Expr(:where,          e...)
makeassign(e...)            = Expr(:(=),            e...)
makefunction(f,e...)        = Expr(:function,f,makeblock(e...))
makefor(e...)               = Expr(:for,e...)
makemacro(m,e...)           = Expr(:macrocall,      Symbol("@",m),:LineNumberHere,e...)
makeespymacro(m,e...)       = Expr(:macrocall,      makesub(:EspyInsideFunction,Symbol("@",m)),:LineNumberHere,e...)
makestructure(e...)         = Expr(:struct,false,   e...)
makemutable(e...)           = Expr(:struct,true ,   e...)
makedeclare(e...)           = Expr(:(::),e[1],e[2])
makedot(e...)               = Expr(:(.) ,e[1],e[2])
makesub(e...)               = Expr(:(.) ,e[1],QuoteNode(e[2]))

## deline: delete LineNumberNode
deline(ex) = ex
function deline(ex::Expr)
    args = deline.(ex.args[.! isline.(ex.args)])
    return Expr(ex.head,args...)
end

## Take out QuoteNodes
function dequote(ex::Expr)
    if isexpr(:(.),ex)
        exo = makesub(dequote(getleft(ex)),dequote(getright(ex)))
    else
        exo = Expr(ex.head,[dequote(a) for a ∈ ex.args]...)
    end
    return exo
end
dequote(ex::QuoteNode) = ex.value
dequote(ex)            = ex

## For human readable code, use println(ex) or println(deline(ex))
# To understand the tree structure of an Expression, use pretty(ex)
function dent(i)
    for _ = 1:i
        @printf("   ")
    end
end
pretty(e::Union{Expr,Symbol,LineNumberNode,Symbol,Nothing})   = pretty(e,0)
function pretty(e::Expr,ind)
    dent(ind)
    @printf("%s\n",e.head)
    for iarg = 1:length( e.args)
        pretty(e.args[iarg],ind+1)
    end
end
function pretty(e::Symbol,ind)
    dent(ind)
    @printf("%s\n",e)
end
function pretty(::Nothing,ind)
    dent(ind)
    @printf("Nothing()\n")
end
function pretty(e::Number,ind)
    dent(ind)
    @printf("%g\n",e)
end
function pretty(e::LineNumberNode,ind)
    dent(ind)
    @printf("LineNumberNode\n")
end
function pretty(x,ind)
    dent(ind)
    display(x)
end
function pretty(m::String)
    @printf("%s\n",m)
end
macro pretty(ex)
    pretty(ex)
    return ex
end

## Definition by element of what can be requested
"""

    forloop

Component to build the `requestable` input to [`makekey`](@ref)
See also: [`makekey`](@ref), [`scalar`](@ref)
"""
struct forloop
    range :: Int
    body
end
"""

    scalar

Component to build the `requestable` input to [`makekey`](@ref)
See also: [`makekey`](@ref), [`forloop`](@ref)
"""
const scalar = ()#Int64[]
## Request definition

struct UnknownSymbolInOutputRequest <: Exception
    symbol::Symbol
    eltyp::DataType
end

# Julia parses a.(b::Tb,c.d).(e::Te,f::Tf) as (a.(b::Tb,c.d)).(e::Te,f::Tf)
# transform to a.((b::Tb,c.d).(e::Te,f::Tf))
function leftparse(e)
    if     issub(e)     && issub(getleft(e))    leftparse(makedot(getleft(getleft(e)),makedot(    getright(getleft(e)),getright(e))))
    elseif isref(e)     && issub(getleft(e))    leftparse(makedot(getleft(getleft(e)),makeref(    getright(getleft(e))            )))
    elseif isdeclare(e) && issub(getleft(e))    leftparse(makedot(getleft(getleft(e)),makedeclare(getright(getleft(e)),getright(e))))
    elseif isexpr(e)                            Expr(e.head,leftparse.(e.args)...)
    elseif isquote(e)                           e.value
    else                                        e
    end
end

# Given a request, recursively generate the "key" to access big output arrays
# The input is an Expr (or Symbol), the output a tree structure of arrays (over GP) and
# NamedTuple, with leaves that are Int64 or AbstractArray of Int64 - containing indices into the output array
namedtuple(sym,val)         = (;zip(sym,val)...)
namedtuple(sym::Symbol,val) = (;zip((sym,),(val,))...)
isloopreq(ele) = isdot(ele) && isref(getleft(ele))
iscallreq(ele) = isdot(ele) && ~isref(getleft(ele))
function makekey_tuple(cnt,ex,reqabl) # ex=(a,gp[].(...),foo.(...))   reqabl=(a=[...],gp=forloop(ngp,(...)),foo=(...))
    if ex == () return NamedTuple(),0 end
    if  isquote(ex) ex = maketuple(ex.value) end # allow user to type a.b for (a.(b,),)
    if ~istuple(ex) ex = maketuple(ex)       end
    len   = length(ex.args)
    elkey = Vector{Any   }(undef,len) # "Any" is OK here, will be used to build a tuple
    sym   = Vector{Symbol}(undef,len)
    for (i,ele) = enumerate(ex.args)
        if issymbol(ele)    # :ele
            sym[i]       = ele
            elkey[i],cnt = makekey_symbol(cnt,reqabl[ele])
        elseif isloopreq(ele) # sym[].(body)
            sym[i]       = getleft(getleft(ele))
            body         = getright(ele)
            elkey[i],cnt = makekey_loop(cnt,body,reqabl[sym[i]])
        elseif iscallreq(ele) # sym.(body)
            sym[i]       = getleft(ele)
            body         = getright(ele)
            elkey[i],cnt = makekey_tuple(cnt,body,reqabl[sym[i]])
        else
            error("Illegal request expression")
        end
    end
    return namedtuple(sym,elkey),cnt
end
function makekey_symbol(cnt,siz) # ex=:a  reqabl=[...]
    len     = prod(siz)
    ndim    = length(siz)
    key     = ndim>0 ? collect(reshape(cnt+1:cnt+len,siz)) : cnt+1
    cnt    += len
    return key,cnt
end
function makekey_loop(cnt,body,reqabl) # loop[].(body) reqabl=forloop(range,(rbody))
    range   = reqabl.range
    rbody   = reqabl.body
    tmp,cnt = makekey_tuple(cnt,body,rbody)
    key     = Vector{typeof(tmp)}(undef,range)
    key[1]  = tmp
    for ia  = 2:range
        key[ia],cnt = makekey_tuple(cnt,body,rbody)
    end
    return key,cnt
end

"""
    key = makekey(requested,requestable)

Create a "key" i.e. a data structure of indices into an array `out` of internal
results, returned by the code generated by `@espy`.

Inputs are
- `requested` a data structure defining a request. This input is provided
   by the user of the code to specify what results are to be extracted.
- `requestable` a named tuple defining the names and sizes of intermediate results
   that can be requested from a given function: this input is provided

# Example

    requestable  = (gp=forloop(2, (z=scalar,s=scalar, material=(a=scalar,b=scalar))),)
    requested    = @request gp[].(s,z,material.(a,b))
    key,nkey     = makekey(requested,requestable)

returns `key` such that

    key.gp[1] == (s=1, z=2, material = (a=3, b=4))
    key.gp[2] == (s=5, z=6, material = (a=7, b=8))
    key.gp[2].material.a == 7
    nkey      == 8

See also: [`@espy`](@ref), [`@espydbg`](@ref), [`@request`](@ref), [`forloop`](@ref), [`scalar`](@ref)
"""
makekey(requested,requestable) = makekey_tuple(0,leftparse(requested),requestable)
"""

    req = @request expr

Create a request of internal results wanted from a function. Considering the function
presented as example for [`@espy`](@ref), examples of possible syntax include

    req       = @request gp[].(s,z,material.(a,b))
    req       = @request gp[].(s)
    req       = @request gp[].(material.(a))

The first expression can be read as follows: "In the function, there is a `for` loop over variable `igp`,
and the results are wanted as a vector (one element for each cycle of the loop).  Each element of the vector
shall be a `type` (a structure) with a field `material`, because a function of that name is called in the
for loop.  Within that function, a variable `a` is to be retrieved.

Note the need to use parentheses also for single-element lists, as in `(s)`.

See also: [`@espy`](@ref), [`@espydbg`](@ref), [`makekey`](@ref)
"""
macro request(ex)
    return QuoteNode(ex)
end

## Data extraction macros
#@espy_loop key gp igp → key_gp=key.gp[igp]
macro espy_loop(key,loopname)                                                   # :key, :gp
    n       = Symbol("n"    ,loopname)                                          # :ngp
    i       = Symbol("i"    ,loopname)                                          # :igp
    key_sub = Symbol(key,"_",loopname)                                          # :key_gp
    return esc(quote
        if haskey($key,$(QuoteNode(loopname)))                                  # if haskey(key,:gp) QuoteNode to produce a Symbol in the code
            $key_sub   = $key.$loopname[$i]                                     # key_gp=key.gp[igp]
        else
            $key_sub =  Nothing
        end
    end)
end
#@espy_record out key var tmp → out[key.var] = var
macro espy_record(out,key,var)
    return esc(quote
        if haskey($key,$(QuoteNode(var)))                                       # if haskey(key,:x)
            if typeof($key.$var) == Int64
                $out[$key.$var] = $var                                          # out[key.x] = x
            else
                $out[$key.$var] .= $var                                         # out[key.x] = x
            end
        end
    end)
end
# @espy_call out key foo(args) → foo(out,key.foo,args...)
macro espy_call(out,key,f)                                                      # out,key,foo(args)
    if f.head ≠ :call error("EspyInsideFunction, @espy_call internal error") end
    foo     = f.args[1]                                                         # foo
    key_sub = makesub(key,foo)                                                  # key.foo
    fp      = makecall(foo,out,key_sub,f.args[2:end]...)                        # foo(out,key.foo,args)
    return esc(quote
        haskey($key,$(QuoteNode(foo))) ? $fp : $f                               # haskey(key,:foo) ? foo(out,key.foo,args) : foo(args)
    end)
end
## @espy
function newinput(ex,newargs...)
    if isexpr(:call,ex)                                                         # f(x)
        exo = makecall(ex.args[1],newargs...,ex.args[2:end]...)                 # f(newargs...,x)
    elseif isexpr(:where,ex)                                                    # f(x::R) where{R<:Number}
        exo = makewhere(newinput(ex.args[1],newargs...),ex.args[2:end]...)      # f(newargs...,x) where{R<:Number}
    else error("espy-newinput")
    end
end
function espy_(ex::Expr,out,key,trace=false)
    if     isexpr(:function,ex)                                                 # foo(a,b,c)...end
        trace && println("function")
        header,body = getleft(ex),getright(ex)                                  # foo(a,b,c) , ...
        header = newinput(header,out,key)                                       # foo(out,key,a,b,c)
        body   = espy_(body,out,key,trace)                                      # recursion
        exo    = makefunction(header,body.args...)                              # foo(out,key,a,b,c)...end
    elseif isexpr(:for,ex)
        trace && println("for")
        (var,lo,hi,body) = getfor(ex)                                           # igp,1,ngp, code-inside-the-block
        loopname = Symbol(string(var)[2:end])                                   # gp
        subkey   = Symbol(key,"_",loopname)                                     # key_gp
        body     = espy_(body,out,subkey,trace)                                 # recursion
        loopmac  = makeespymacro(:espy_loop,key,loopname)                       # @espy_loop key gp → key_gp = key.gp[igp]
        body.args = [loopmac,body.args...]                                      # prepend macro to body
        exo      = makefor(ex.args[1],body)                                     # rebuild the loop
    elseif isexpr(:(=),ex)
        trace && println("assign")
        right  = getright(ex)                                                   # rhs
        left   = getleft(ex)                                                    # lhs
        rec    = Vector{Expr}(undef,0)                                          # will contain the macros to insert
        if iscall(right)
            foo             = getleft(right)
            if isquote(foo)                                                     # ... = foo, or ... = :foo?
                right.args[1]   = foo.value                                     # right: foo(x,y,z)
                right           = makeespymacro(:espy_call,out,key,right)       # ... = @espy_call out key foo(x,y,z)
            end
        end
        if isexpr(:tuple,left)                                                  # (a,:b) = ...
            args = left.args
            left = ()                                                           # will contain (a,b)
            for o ∈ args
                if isquote(o)                                                   # ...,:b =
                    left = (left...,o.value)                                    # (a,b)
                    push!(rec    ,makeespymacro(:espy_record,out,key,o.value))  # @espy_record out key b
                else                                                            #  a,... =
                    left = (left...,o)                                          # (a,...)
                end
            end
            left = maketuple(left...)
        elseif isquote(left)                                                    #:a = ...
            left = left.value                                                   # a = ...
            push!(rec    ,makeespymacro(:espy_record,out,key,left))             # @espy_record out key a
        end
        exo = makeblock(makeassign(left,right),rec...)
    else
        trace && println("recursion")
        exo = Expr(ex.head,[espy_(a,out,key,trace) for a ∈ ex.args]...)
    end

    return exo
end
function espy_(ex,out,key,trace=false)
    trace && println("default")
    return ex
end
"""

    @espy function residual(x,y)
        ngp=2
        r = 0
        for igp=1:ngp
            :z = x[igp]+y[igp]
            :s,dum  = :material(z)
            r += s
        end
        return r
    end
    @espy function material(z)
        :a = z+1
        :b = a*z
        return b,3.
    end

Transform the code of a function in which variables and function calls have been annotated with `:`
in order to allow the extraction of intermediate results.

The above annotated code will result in the generation of "clean" code in which the `:` annotations
have been taken out

    function residual(x,y)
        ngp=2
        r = 0
        for igp=1:ngp
            z = x[igp]+y[igp]
            s,dum  = material(z)
            r += s
        end
        return r
    end
    function material(z)
        a = z+1
        b = a*z
        return b,3.
    end

The macro will also generate code with additional `out` and `key` arguments:

    function residual(out,key,x,y)
        ngp = 2
        r   = 0
        for igp = 1:ngp
            @espy_loop key gp                     # key_gp = key.gp[igp]
            z = x[igp]+y[igp]
            @espy_record out key_gp z             # out[key_gp.z] = z
            s = @espy_call out key_gp material(z) # s = material(out,key_gp.material,z)
            @espy_record out key_gp s             # out[key_gp.s] = s
            r += s
        end
        return r
    end
    function material(out,key,z)
        a = z+1
        @espy_record out key a                    # out[key.a] = a
        b = a*z
        @espy_record out key b                    # out[key.b] = b
        return b
    end

The above code contains more macros, which in turn evaluate as
shown in the comments.  More precisely,

    @espy_record out key a

evaluates to

    if haskey(key,a)
        out[key.a] = a
    end

`key` is a data structure generated by [`makekey`](@ref) based on a [`@request`](@ref).

When the version of `residual` with additional parameter `out` has been called, the content
of this output is accessed using `key`:

    requestable  = (gp=forloop(2, (z=scalar,s=scalar, material=(a=scalar,b=scalar))),)
    requested    = @request gp[].(s,z,material.(a,b))
    key,nkey     = makekey(requested,requestable)
    residual(out,key,x,y)
    igp          = 2
    z            = out[key.gp[igp].z]

See also: [`@espydbg`](@ref), [`@request`](@ref), [`makekey`](@ref)
"""
macro espy(ex)
    cleancode = dequote(ex)
    espycode  = espy_(ex,newsym(:espy_key),newsym(:espy_out),false)
    return makeblock(esc(espycode),esc(cleancode))
end
"""
    @espydbg function ...
    end

Run [`@espy`](@ref) and to generate code and print the output code (for debug purposes).

See also: [`@espy`](@ref), [`@request`](@ref), [`makekey`](@ref), [`forloop`](@ref), [`scalar`](@ref)"""
macro espydbg(ex)
    println(">>>>>> espy clean code")
    cleancode = dequote(ex)
    println(deline(cleancode))
    println("###### espy extractor code")
    espycode  = espy_(ex,newsym(:espy_key),newsym(:espy_out),false)
    println(deline(espycode))
    println("<<<<<<")
    return makeblock(esc(espycode),esc(cleancode))
end

end

max_stars_repo_name: PhilippeMaincon/Espy