text:
println("===============================================")
println("Testing zener.jl")
println("===============================================")

function Fract_Zener_G_reduce_SLS()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fraczen = relaxmod(Fract_Zener, t, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = relaxmod(SLS_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanSLS[i]), eachindex(t))
end
@test Fract_Zener_G_reduce_SLS()

function FractSLS_Zener_G_reduce_SLS()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracSLS = relaxmod(FractSLS_Zener, t, [2.0, 1.0, 1.0, 0.5])
    stanSLS = relaxmod(SLS_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(t))
end
@test FractSLS_Zener_G_reduce_SLS()

function FractSolid_G_reduce_SLS()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracSpec = relaxmod(FractSolid, t, [2.0, 1.0, 0.0, 0.5])
    stanSLS = relaxmod(SLS_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSpec[i], stanSLS[i]), eachindex(t))
end
@test FractSolid_G_reduce_SLS()

function Fract_Zener_J_reduce_SLS()
    dt = 0.001
    # NaN at 0.0 due to InverseLaplace so start at dt
    t = collect(dt:dt:10.0)

    fraczen = creepcomp(Fract_Zener, t, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = creepcomp(SLS_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanSLS[i]), eachindex(t))
end
@test Fract_Zener_J_reduce_SLS()

function FractSLS_Zener_J_reduce_SLS()
    dt = 0.001
    # NaN at 0.0 due to InverseLaplace so start at dt
    t = collect(dt:dt:10.0)

    fracSLS = creepcomp(FractSLS_Zener, t, [2.0, 1.0, 1.0, 0.5])
    stanSLS = creepcomp(SLS_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(t))
end
@test FractSLS_Zener_J_reduce_SLS()

function FractSolid_J_reduce_SLS()
    dt = 0.001
    # NaN at 0.0 due to InverseLaplace so start at dt
    t = collect(dt:dt:10.0)

    fracSpec = creepcomp(FractSolid, t, [2.0, 1.0, 0.0, 0.5])
    stanSLS = creepcomp(SLS_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSpec[i], stanSLS[i]), eachindex(t))
end
@test FractSolid_J_reduce_SLS()

function Fract_Zener_Gp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fraczen = storagemod(Fract_Zener, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = storagemod(SLS_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanSLS[i]), eachindex(chirp))
end
@test Fract_Zener_Gp_reduce_SLS()

function FractSLS_Zener_Gp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracSLS = storagemod(FractSLS_Zener, chirp, [2.0, 1.0, 1.0, 0.5])
    stanSLS = storagemod(SLS_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(chirp))
end
@test FractSLS_Zener_Gp_reduce_SLS()

function FractSolid_Gp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracSpec = storagemod(FractSolid, chirp, [2.0, 1.0, 0.0, 0.5])
    stanSLS = storagemod(SLS_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSpec[i], stanSLS[i]), eachindex(chirp))
end
@test FractSolid_Gp_reduce_SLS()

function Fract_Zener_Gpp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fraczen = lossmod(Fract_Zener, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = lossmod(SLS_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanSLS[i]), eachindex(chirp))
end
@test Fract_Zener_Gpp_reduce_SLS()

function FractSLS_Zener_Gpp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracSLS = lossmod(FractSLS_Zener, chirp, [2.0, 1.0, 1.0, 0.5])
    stanSLS = lossmod(SLS_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(chirp))
end
@test FractSLS_Zener_Gpp_reduce_SLS()

function FractSolid_Gpp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracSpec = lossmod(FractSolid, chirp, [2.0, 1.0, 0.0, 0.5])
    stanSLS = lossmod(SLS_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSpec[i], stanSLS[i]), eachindex(chirp))
end
@test FractSolid_Gpp_reduce_SLS()

function Fract_Zener_G_reduce_Jeffreys()
    dt = 0.001
    # NaN at 0.0 due to Gamma function so start at dt
    t = collect(dt:dt:10.0)

    fraczen = relaxmod(Fract_Zener, t, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = relaxmod(Jeffreys_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanJef[i]), eachindex(t))
end
@test Fract_Zener_G_reduce_Jeffreys()

function FractJeffreys_Zener_G_reduce_Jeffreys()
    dt = 0.001
    # NaN at 0.0 due to Gamma function so start at dt
    t = collect(dt:dt:10.0)

    fracJef = relaxmod(FractJeffreys_Zener, t, [2.0, 1.0, 0.0, 0.5])
    stanJef = relaxmod(Jeffreys_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(t))
end
@test FractJeffreys_Zener_G_reduce_Jeffreys()

function Fract_Zener_J_reduce_Jeffreys()
    dt = 0.001
    # NaN at 0.0 due to InverseLaplace function so start at dt
    t = collect(dt:dt:10.0)

    fraczen = creepcomp(Fract_Zener, t, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = creepcomp(Jeffreys_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanJef[i]), eachindex(t))
end
@test Fract_Zener_J_reduce_Jeffreys()

function FractJeffreys_Zener_J_reduce_Jeffreys()
    dt = 0.001
    # NaN at 0.0 due to InverseLaplace function so start at dt
    t = collect(dt:dt:10.0)

    fracJef = creepcomp(FractJeffreys_Zener, t, [2.0, 1.0, 0.0, 0.5])
    stanJef = creepcomp(Jeffreys_Zener, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(t))
end
@test FractJeffreys_Zener_J_reduce_Jeffreys()

function Fract_Zener_Gp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fraczen = storagemod(Fract_Zener, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = storagemod(Jeffreys_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanJef[i]), eachindex(chirp))
end
@test Fract_Zener_Gp_reduce_Jeffreys()

function FractJeffreys_Zener_Gp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracJef = storagemod(FractJeffreys_Zener, chirp, [2.0, 1.0, 0.0, 0.5])
    stanJef = storagemod(Jeffreys_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(chirp))
end
@test FractJeffreys_Zener_Gp_reduce_Jeffreys()

function Fract_Zener_Gpp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fraczen = lossmod(Fract_Zener, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = lossmod(Jeffreys_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fraczen[i], stanJef[i]), eachindex(chirp))
end
@test Fract_Zener_Gpp_reduce_Jeffreys()

function FractJeffreys_Zener_Gpp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracJef = lossmod(FractJeffreys_Zener, chirp, [2.0, 1.0, 0.0, 0.5])
    stanJef = lossmod(Jeffreys_Zener, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(chirp))
end
@test FractJeffreys_Zener_Gpp_reduce_Jeffreys()


text:
"""
Compute p(x) where p is the unique polynomial of degree length(xi),
such that p(x[i]) = y[i]) for all i.

    ti: interpolation nodes
    xi: interpolation values
    t:  evaluation point
    x:  evaluation value
"""
function aitken_neville(ti::Vector{TT}, xi::Matrix, t::TT, x::Vector) where {TT}
    @assert length(ti) == size(xi,2)
    @assert length(x)  == size(xi,1)

    for j in axes(t)
        for i in 1:(length(t)-j)
            for k in axes(x,1)
                xi[k,i] = xi[k,i+1] + (xi[k,i] - xi[k,i+1]) * (ti[i+j] - t) / (ti[i+j] - ti[i])
            end
        end
    end
    for k in eachindex(x)
        x[k] = xi[k,1]
    end
end

"""
Euler extrapolation method with arbitrary order p.

    v:  function to compute vector field
    t₀: initial time
    t₁: final   time
    x₀: initial value
    x₁: final   value
    s:  number of interpolations (order p=s+1)

# TODO This is probably broken!
"""
function euler_extrapolation(v::Function, t₀::TT, t₁::TT, x₀::Vector, x₁::Vector, s::Int) where {TT}
    @assert size(x₀) == size(x₁)

    F   = collect(1:(s+1))
    Δt  = t₁ - t₀
    σ   = Δt ./ F
    pts = repmat(x₀, 1, s+1)

    local xᵢ = zero(x₀)
    local vᵢ = zero(x₀)

    for i in 1:(s+1)
        for j in 1:(F[i]-1)
            tᵢ  = t₀ + σ[i]
            for k in axes(pts,1)
                xᵢ[k] = pts[k,i]
            end
            v(tᵢ, xᵢ, vᵢ)
            for k in axes(pts,1)
                pts[k,i] += σ[i] * vᵢ[k]
            end
        end
    end

    aitken_neville(σ, pts, zero(TT), x₁)
end

"""
Midpoint extrapolation method with arbitrary order p.

    v:  function to compute vector field
    t₀: initial time
    t₁: final   time
    x₀: initial value
    x₁: final   value
    s:  number of interpolations (order p=2s+2)
"""
function midpoint_extrapolation(v::Function, t₀::TT, t₁::TT, x₀::Vector, x₁::Vector, s::Int) where {TT}
    @assert size(x₀) == size(x₁)

    local F   = 2*collect(1:(s+1))
    local Δt  = t₁ - t₀
    local σ   = Δt ./ F
    local σ²  = σ.^2
    local pts = zeros(eltype(x₀), length(x₀), s+1)

    local xᵢ₁= zero(x₀)
    local xᵢ₂= zero(x₀)
    local xᵢₜ= zero(x₀)
    local vᵢ = zero(x₀)
    local v₀ = zero(x₀)

    v(t₀, x₀, v₀)

    for i in 1:s+1
        tᵢ   = t₀ + σ[i]
        xᵢ₁ .= x₀
        xᵢ₂ .= x₀ + σ[i] .* v₀
        for j in 1:(F[i]-1)
            v(tᵢ, xᵢ₂, vᵢ)
            xᵢₜ .= xᵢ₁ + 2σ[i] .* vᵢ
            xᵢ₁ .= xᵢ₂
            xᵢ₂ .= xᵢₜ
        end
        for k in axes(pts,1)
            pts[k,i] += xᵢ₂[k]
        end
    end

    aitken_neville(σ², pts, zero(TT), x₁)
end

"""
Midpoint extrapolation method with arbitrary order p.

    v:  function to compute vector field
    f:  function to compute force  field
    t₀: initial time
    t₁: final   time
    q₀: initial positions
    p₀: initial momenta
    q₁: final   positions
    p₁: final   momenta
    s:  number of interpolations (order p=2s+2)
"""
function midpoint_extrapolation(v::Function, f::Function, t₀::TT, t₁::TT, q₀::Vector, q₁::Vector, p₀::Vector, p₁::Vector, s::Int) where {TT}
    @assert size(q₀) == size(q₁) == size(p₀) == size(p₁)

    local F   = 2*collect(1:(s+1))
    local Δt  = t₁ - t₀
    local σ   = Δt ./ F
    local σ2  = σ.^2

    local qts = zeros(eltype(q₀), length(q₀), s+1)
    local pts = zeros(eltype(p₀), length(p₀), s+1)

    local qᵢ₁= zero(q₀)
    local qᵢ₂= zero(q₀)
    local qᵢₜ= zero(q₀)

    local pᵢ₁= zero(p₀)
    local pᵢ₂= zero(p₀)
    local pᵢₜ= zero(p₀)

    local v₀ = zero(q₀)
    local vᵢ = zero(q₀)

    local f₀ = zero(p₀)
    local fᵢ = zero(p₀)

    v(t₀, q₀, p₀, v₀)
    f(t₀, q₀, v₀, f₀)

    for i in 1:(s+1)
        tᵢ   = t₀ + σ[i]
        qᵢ₁ .= q₀
        qᵢ₂ .= q₀ + σ[i] .* v₀
        pᵢ₁ .= p₀
        pᵢ₂ .= p₀ + σ[i] .* f₀
        for j in 1:(F[i]-1)
            v(tᵢ, qᵢ₂, pᵢ₂, vᵢ)
            f(tᵢ, qᵢ₂, vᵢ,  fᵢ)
            qᵢₜ .= qᵢ₁ + 2σ[i] .* vᵢ
            qᵢ₁ .= qᵢ₂
            qᵢ₂ .= qᵢₜ
            pᵢₜ .= pᵢ₁ + 2σ[i] .* fᵢ
            pᵢ₁ .= pᵢ₂
            pᵢ₂ .= pᵢₜ
        end
        for k in axes(qts,1)
            qts[k,i] += qᵢ₂[k]
        end
        for k in axes(pts,1)
            pts[k,i] += pᵢ₂[k]
        end
    end

    aitken_neville(σ2, qts, zero(TT), q₁)
    aitken_neville(σ2, pts, zero(TT), p₁)
end

ext: jl


text:
"""
    boostimage(img, param; [iterations = 1, res_x = 500, res_y = 500])

Turns an image into an abstract figure by predicting its features with a gradient booster.

# Arguments
* `img::Matrix{T<:AbstractRGB}`: the image to boost.
* `param::Dict`: a dictionary with the XGBoost parameters to use in the prediction.
* `iterations::Int`: keyword argument that sets the number of boosting iterations. Defaults to `1`.
* `res_x::Int`: keyword argument that sets the horizontal resolution of the rendered abstract
    figure. Defaults to `500`.
* `res_y::Int`: keyword argument that sets the vertical resolution of the rendered abstract figure.
    Defaults to `500`.
"""
function boostimage{T<:AbstractRGB}(img::Matrix{T}, param::Dict;
                                    iterations::Int = 1, res_x::Int = 500, res_y::Int = 500)
    println("Preparing data structures...")
    bst_img = BoostingImage(img)
    pred_img = BoostingImage(res_x, res_y, bst_img)
    gr_img = GRImage(res_x, res_y)

    X = hcat(bst_img.x, bst_img.y)
    data_r = DMatrix(X, label = bst_img.r; transposed = false)
    data_g = DMatrix(X, label = bst_img.g; transposed = false)
    data_b = DMatrix(X, label = bst_img.b; transposed = false)
    pred_coords = DMatrix(hcat(pred_img.x, pred_img.y); transposed = false)

    bst_r = xgboost(data_r, 0, param = param, metrics = ["rmse"])
    bst_g = xgboost(data_g, 0, param = param, metrics = ["rmse"])
    bst_b = xgboost(data_b, 0, param = param, metrics = ["rmse"])

    initws(res_x, res_y)
    for iter in 1:iterations
        println("[", iter, "]: training Red channel...")
        XGBoost.update(bst_r, data_r, 1)
        pred_img.r = predict(bst_r, pred_coords)

        println("[", iter, "]: training Green channel...")
        XGBoost.update(bst_g, data_g, 1)
        pred_img.g = predict(bst_g, pred_coords)

        println("[", iter, "]: training Blue channel...")
        XGBoost.update(bst_b, data_b, 1)
        pred_img.b = predict(bst_b, pred_coords)

        println("[", iter, "]: Rendering...")
        update!(gr_img, pred_img)
        render(gr_img)
    end

    return nothing
end


text:
function splint(n, x, f )
    # arguments
    # integer, intent(in) :: n
    # real(8), intent(in) :: x(n),f(n)
    res = 0.0
    if n <= 4
        res = polynm(-1, n, x, f, x[n])
        return res
    end

    # fit piecewise cubic spline to data and integrate
    x0 = x[1]
    x1 = x[2] - x0
    x2 = x[3] - x0
    x3 = x[4] - x0
    t4 = x1 - x2
    t5 = x1 - x3
    t6 = x2 - x3

    y0 = f[1]
    y1 = f[2] - y0
    y2 = f[3] - y0
    y3 = f[4] - y0
    t1 = x1*x2*y3
    t2 = x2*x3*y1
    t3 = x1*x3
    t0 = 0.5/(t3*t4*t5*t6)
    t3 = t3*y2
    t7 = t1*t4 + t2*t6 - t3*t5
    t4 = x1^2
    t5 = x2^2
    t6 = x3^2
    y1 = t3*t6 - t1*t5
    y3 = t2*t5 - t3*t4
    y2 = t1*t4 - t2*t6
    t1 = x1*y1 + x2*y2 + x3*y3
    t2 = y1 + y2 + y3
    res = x2*( y0 + t0*(t1 + x2*(0.5*t7*x2 - 2*t2/3) ) )

    for i in 3:n-3
        x0 = x[i]
        x1 = x[i-1] - x0
        x2 = x[i+1] - x0
        x3 = x[i+2] - x0
        t4 = x1 - x2
        t5 = x1 - x3
        t6 = x2 - x3
        t3 = x1*x3
        y0 = f[i]
        y1 = f[i-1] - y0
        y2 = f[i+1] - y0
        y3 = f[i+2] - y0
        t1 = x1*x2*y3
        t2 = x2*x3*y1
        t0 = 0.5/(t3*t4*t5*t6)
        t3 = t3*y2
        t7 = t1*t4 + t2*t6 - t3*t5
        t4 = x1^2
        t5 = x2^2
        t6 = x3^2
        y1 = t3*t6 - t1*t5
        y2 = t1*t4 - t2*t6
        y3 = t2*t5 - t3*t4
        t1 = x1*y1 + x2*y2 + x3*y3
        t2 = y1 + y2 + y3
        res = res + x2*( y0 + t0*( t1 + x2*(0.5*t7*x2 - 2*t2/3) ) )
    end
    x0 = x[n-2]
    x1 = x[n-3] - x0
    x2 = x[n-1] - x0
    x3 = x[n] - x0
    t4 = x1 - x2
    t5 = x1 - x3
    t6 = x2 - x3
    y0 = f[n-2]
    y1 = f[n-3] - y0
    y2 = f[n-1] - y0
    y3 = f[n] - y0
    t1 = x1*x2
    t2 = x2*x3*y1
    t3 = x1*x3*y2
    t0 = 0.5/(t1*t4*t5*t6)
    t1 = t1*y3
    t7 = t1*t4 + t2*t6 - t3*t5
    t4 = x1^2
    t5 = x2^2
    t6 = x3^2
    y1 = t3*t6 - t1*t5
    y2 = t1*t4 - t2*t6
    y3 = t2*t5 - t3*t4
    t1 = x1*y1 + x2*y2+x3*y3
    t2 = y1 + y2 + y3
    res = res + x3*(y0 + t0*(t1 + x3*(0.5*t7*x3 - 2*t2/3) ) )

    return res
end


text:
#devsafety.jl - options and related macros that deal with developer mode safety
#checks.

const __LOADTIME_DEVMODE = options[:devmode]

doc"""
  `@devmode_on` initiates a section of code where developer checking is guaranteed to be "on".
  This code is not executed if developer mode was off at loadtime.
"""
macro devmode_on(expr)
  if __LOADTIME_DEVMODE

    #add this to the front
    unshift!(expr.args, quote
      cache = Unums.options[:devmode]
      Unums.options[:devmode] = true
    end)

    #add this to the end
    push!(expr.args, quote
      Unums.options[:devmode] = cache
    end)

    esc(expr)
  else
    :()
  end
end

doc"""
  `@devmode_off` initiates a section of code where developer checking is guaranteed to
  "off".
"""
macro devmode_off(expr)
  if __LOADTIME_DEVMODE

    #add this to the front
    unshift!(expr.args, quote
      cache = Unums.options[:devmode]
      Unums.options[:devmode] = false
    end)

    #add this to the end
    push!(expr.args, quote
      Unums.options[:devmode] = cache
    end)

    esc(expr)
  else
    :()
  end
end

export @devmode_on, @devmode_off

################################################################################
# @dev_check macro which uses reflection to automatically execute a function
# parameters integrity check before certain functions.

function __extract_function_string(fvar)
  #for variadic functions, the function name parameter appears as part of an
  #implicit untyped array.  This takes care of this.
  isa(fvar, Symbol) && return string(fvar)
  if isa(fvar, Expr) && (fvar.head == :curly)
    return string(fvar.args[1])
  end
  error("malformed function head.")
end

#strips the variable name off the expression, if it's just a symbol we return it
#if it's more, then we
function getvariable(e)
  isa(e, Symbol) && return e
  if isa(e, Expr)
    e.head == :(::) && return e.args[1]
  end
end

doc"""
the `@dev_check` macro injects a conditional check that sees if
`options[:devmode]` is set.  If it is, `@dev_check` inserts a line of code
that expects a companion function with symbol `__check_[functionname]`.  This
function is passed the same parameters as the parent function and is expected
to throw an error if parameter preconditions are violated.

For types it may be necessary to check the parameters passed to the constructor.
The type-variadic parameters will be passed to the check function immediately in
front of the type parameters.
"""
macro dev_check(expr)
  #bail with an unaltered expression if we haven't set the developer mode.
  options[:devmode] || return expr
  #make sure this macro precedes a function definition.
  isa(expr, Expr) || ArgumentError("checkable macro must precede a function or type definition")
  if (expr.head == :function)
    #prepend a __check_ prefix to the function name.
    #do a little bit of syntactic sugar to fix functions that already start with __
    fnstring = string("__check_", __extract_function_string(expr.args[1].args[1]))
    fparams = join(map(getvariable, expr.args[1].args[2:end]), ", ")
    #generate the line of code to be executed, pass it to parse.
    b2 = parse("$fnstring($fparams)")
    #inject the faux line number and the generated code into the ast.
    unshift!(expr.args[2].args, b2)
    #return the escaped function to the parser so that it generates the new function.

    return esc(expr)
  elseif (expr.head == :type)
    #what to do if dev_check is appied to a type statement.
    typesig = expr.args[2]    #extract the type signature
    typemembers = filter((e) -> (e.head != :line), expr.args[3].args)  #extract the type members, with typing info
    typevars = map((e) -> (e.args[1]), typemembers)                    #a list of type variables, without typing info

    #type signature might be a complex expression or it might be a single symbol.
    if isa(typesig, Expr)
      if typesig.head == :curly
        typename = typesig.args[1]
        varparams = typesig.args[2:end]
      elseif typesig.head == :<:
        typename = typesig.args[1].args[1]
        varparams = typesig.args[1].args[2:end]
      else
        throw(ArgumentError("malformed type expression passed to @dev_check"))
      end
    else
      typename = typesig
      varparams = []
    end

    constructor = Expr(:function)
    callexp = Expr(:call)
    push!(callexp.args, typename)
    append!(callexp.args, typemembers)

    #assemble the call to the check function
    checkname = Symbol("__check_", typename)
    checkcall = Expr(:call)
    push!(checkcall.args, checkname)
    append!(checkcall.args, varparams)  #first add the variadic parameters
    append!(checkcall.args, typevars)   #then add the the type parameters

    #create the check expression which is dependent on the global options variable
    check_expr = :(options[:devmode] && $checkcall)

    #assemble the call to the stealth "new" constructor.
    newcall = Expr(:call)
    push!(newcall.args, :new)
    append!(newcall.args, typevars)

    #assemble fnblock
    fnblock = Expr(:block)
    push!(fnblock.args, check_expr)
    push!(fnblock.args, newcall)

    #assemble the substitute constructor.
    push!(constructor.args, callexp)
    push!(constructor.args, fnblock)

    #append the substitute constructor onto the end of the type definition
    push!(expr.args[3].args, constructor)

    return esc(expr)
  else
    ArgumentError("checkable macro must precede a function or type definition")
  end
end

doc"""
  `Unums.@devmode(expr)` makes the expression disappear if we're not in developer mode.
"""
macro devmode(expr)
  if options[:devmode]
    expr
  else
    :()
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Random
using IntervalRootFinding

# px = p
# py = P
# x = q
# y = Q

""" Energy of the Dicke model """
function Energy(x, parameters)
    P, p, Q, q = x
    A, B, C = parameters

    Σ2 = p*p + P*P + q*q + Q*Q
    s2 = 2.0 - Σ2

    if s2 < 0.0
        return NaN
#        print("s negative!")
#        s2 = 0.0
    end

    return A * Σ2 + B * ((P*P + p*p) * s2 + (P*q - p*Q) ^ 2) + C * P * sqrt(s2)
end

""" Filling initial conditions (at a given energy calculates randomly all missing points) """
function InitialConditions(x0, e, parameters, coordinate)
    x = deepcopy(x0)
    x[coordinate] = 0

    P, p, Q, q = x
    A, B, C  = parameters

    result = []

    s2 = 2.0 - (P*P + p*p + Q*Q + q*q)
    if s2 >= 0
        try
            if coordinate == 1
                result = roots(x -> Energy((x, p, Q, q), parameters) - e, -sqrt(2)..sqrt(2))
            elseif coordinate == 2
                result = roots(x -> Energy((P, x, Q, q), parameters) - e, -sqrt(2)..sqrt(2))
            elseif coordinate == 3
                result = roots(x -> Energy((P, p, x, q), parameters) - e, -sqrt(2)..sqrt(2))
            elseif coordinate == 4
                result = roots(x -> Energy((P, p, Q, x), parameters) - e, -sqrt(2)..sqrt(2))
            end
        catch e
        end
    end

    r = zeros(Float64, length(result))
    for i = 1:length(result)
        r[i] = mid(result[i].interval)
    end

    return sort(r, rev=true)
end

""" True if given trajectory is within the kinematically accessible domain """
function CheckDomain(x, parameters, t)
    P, p, Q, q = x[1:4]
    s2 = 2.0 - (p*p + P*P + q*q + Q*Q)
    return s2 < 0
end

function EquationOfMotion!(dx, x, parameters, t)
    P, p, Q, q = x[1:4]             # Equations of motion for trajectories

    s2 = 2.0 - (p*p + P*P + q*q + Q*Q)
    if s2 < 0
        s2 = 0.0                    # With isoutofdomain=CheckDomain the calculation shouldn't enter here, but in enters anyway
        print("s negative!")
    end

    Φ = reshape(x[5:end], 4, 4)     # Tangent dynamics
    A, B, C = parameters.modelParameters

    s = sqrt(s2)
    a = p*p + P*P
    b = P*Q + p*q

    A2 = 2.0 * A
    B2 = 2.0 * B
    Cs = C / s
    Cs3 = Cs / s2

    dx[1] = -A2 * Q + B2 * P * b + Cs * P * Q
    dx[2] = -A2 * q + B2 * p * b + Cs * P * q
    dx[3] = A2 * P - B2 * (2.0 * P * (a - 1.0) + Q * b) + Cs * (s2 - P*P)
    dx[4] = A2 * p - B2 * (2.0 * p * (a - 1.0) + q * b) - Cs * p * P

    HPQ = -B2 * (p * q + 2.0 * P * Q) - Cs3 * Q * (s2 + P*P)
    HpQ = -B2 * P * q - Cs3 * p * P * Q
    HQQ = A2 - B2 * P*P - Cs3 * P * (s2 + Q*Q)
    HQq = -B2 * p * P - Cs3 * P * q * Q

    HPq = -B2 * p * Q - Cs3 * q * (s2 + P*P)
    Hpq = -B2 * (2.0 * p * q + P * Q) - Cs3 * p * P * q
    Hqq = A2 - B2 * p*p - Cs3 * P * (s2 + q*q)

    HPP = A2 - B2 * (6.0 * P*P + 2.0 * p*p + Q*Q - 2.0) - Cs3 * P * (P*P - 3.0 * a + 3.0 * (s2 + a))
    HPp = -B2 * (4.0 * p * P + q * Q) - Cs3 * p * (s2 + P*P)

    Hpp = A2 - B2 * (2.0 * P*P + 6.0 * p*p + q*q - 2.0) - Cs3 * P * (s2 + p*p)

    j = [-HPQ -HpQ -HQQ -HQq;
         -HPq -Hpq -HQq -Hqq;
         HPP HPp HPQ HPq;
         HPp Hpp HpQ Hpq]

    dx[5:end] = j * Φ
end

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# Define a special stride where x and z are updated.
function get_stride(model, traj)
    stride = zeros(model.nq)
    stride[1:2] = traj.q[end-1][1:2] - traj.q[1][1:2]
    return stride
end

# ## Simulation
s = get_simulation("hopper_2D", "flat_2D_lc", "flat");
s_sim = get_simulation("hopper_2D", "stairs3_2D_lc", "stairs");

# ## Reference Trajectory (stairs)
ref_traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_stair_ref.jld2"),
    load_type=:split_traj_alt));

H = ref_traj.H
h = ref_traj.h
h_sim = h / N_sample
H_sim = 240*N_sample

# ## MPC parameters (stairs)
N_sample = 10
H_mpc = 10
κ_mpc = 2.0e-4
n_opts = NewtonOptions(
    r_tol = 3e-4,
    max_iter = 5)
mpc_opts = CIMPCOptions(
    altitude_update = true,
    altitude_impact_threshold = 0.1,
    altitude_verbose = true)

obj = TrackingVelocityObjective(s.model, s.env, H_mpc,
    v = [Diagonal(1e-3 * [1e-2,1,1,10]) for t = 1:H_mpc],
    q = [[Diagonal(1e-0 * [1e1,1e-1,1,1])   for t = 1:H_mpc-5]; [Diagonal(1e-1 * [1,1e-1,1e1,0.1])   for t = 1:5]],
    u = [Diagonal(1e-0 * [1e0, 1e0]) for t = 1:H_mpc],
    γ = [Diagonal(1e-100 * ones(s.model.nc)) for t = 1:H_mpc],
    b = [Diagonal(1e-100 * ones(s.model.nc * friction_dim(s.env))) for t = 1:H_mpc])

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = n_opts,
    mpc_opts = mpc_opts);

# ## Initial conditions (stairs)
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator (stairs)
sim_stair = simulator(s, H_sim, h=h_sim, policy=p)

# ## Simulate (stairs)
@time status = simulate!(sim_stair, q1_sim, v1_sim)

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize (stairs)
anim = visualize_robot!(vis, model, sim_stair.traj, sample=10, name=:Sim, α=1.0)
stairs!(vis)

# ## Reference trajectory (flip)
ref_traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_tall_flip_ref.jld2"),
    load_type=:split_traj_alt));

for t = 1:ref_traj.H+2
    ref_traj.q[t][1] += sim_stair.traj.q[end][1] # shift flip to top of stairs
end

H = ref_traj.H
h = ref_traj.h
h_sim = h / N_sample
H_sim = 64*N_sample

# ## MPC setup (flip)
obj = TrackingVelocityObjective(s.model, s.env, H_mpc,
    v = [Diagonal(1e-13 * [1e-2,1,1,10]) for t = 1:H_mpc],
    q = [[Diagonal(1e-10 * [1e1,1e1,1,1])   for t = 1:H_mpc-5]; [Diagonal(1e-11 * [1,1e1,1e1,0.1])   for t = 1:5]],
    u = [Diagonal(1e-0 * [1e0, 1e0]) for t = 1:H_mpc],
    γ = [Diagonal(1e-100 * ones(s.model.nc)) for t = 1:H_mpc],
    b = [Diagonal(1e-100 * ones(s.model.nc * friction_dim(s.env))) for t = 1:H_mpc])

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = n_opts,
    mpc_opts = mpc_opts,
    );

# ## Set initial configurations to simulation result
q0_sim = deepcopy(ContactImplicitMPC.SVector{model.nq}(sim_stair.traj.q[end-1]))
q1_sim = deepcopy(ContactImplicitMPC.SVector{model.nq}(sim_stair.traj.q[end]))
v1_sim = (q1_sim - q0_sim) ./ sim_stair.traj.h

# ## Simulator (flip)
sim_flip = simulator(s_sim, H_sim, h=h_sim, policy=p)

# ## Simulate (flip)
@time status = ContactImplicitMPC.simulate!(sim_flip, q1_sim, v1_sim)

# ## Visualize (flip)
anim = visualize_robot!(vis, model, sim_flip.traj, sample=10, name=:Sim, α=1.0)
anim = visualize_robot!(vis, model, ref_traj, anim=anim, name=:Ref, α=0.3)

# ## Visualize (parkour)
ref_traj_full = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/hopper_stairs_flip_ref.jld2"),
    load_type=:split_traj_alt));

N_sample = 10
sim_traj_full = [sim_stair.traj.q[1:end-2]; sim_flip.traj.q]

anim = visualize_robot!(vis, model,
    [sim_traj_full[1:N_sample:end]..., [sim_traj_full[end] for i = 1:50]...],
    name = :Sim, α = 1.0, h = h_sim*N_sample);

anim = visualize_robot!(vis, model,
    [ref_traj_full.q..., [ref_traj_full.q[end] for i = 1:50]...],
    anim = anim, name = :Ref, α = 0.15, h = h_sim*N_sample);

ContactImplicitMPC.stairs!(vis)


text:
# We want to make sure that the first argument of change_t_via_interpolation!
# is specialized, yet, it needs to take both ODEIntegrator and DDEIntegrator.
# Hence, we need to have two separate functions.
function _change_t_via_interpolation!(integrator,t,modify_save_endpoint::Type{Val{T}}) where T
  # Can get rid of an allocation here with a function
  # get_tmp_arr(integrator.cache) which gives a pointer to some
  # cache array which can be modified.
  if integrator.tdir*t < integrator.tdir*integrator.tprev
    error("Current interpolant only works between tprev and t")
  elseif t != integrator.t
    if is_constant_cache(integrator.cache)
      integrator.u = integrator(t)
    else
      integrator(integrator.u,t)
    end
    integrator.t = t
    integrator.dt = integrator.t - integrator.tprev
    DiffEqBase.reeval_internals_due_to_modification!(integrator)
    if T
      solution_endpoint_match_cur_integrator!(integrator)
    end
  end
end
DiffEqBase.change_t_via_interpolation!(integrator::ODEIntegrator,
                                        t,modify_save_endpoint::Type{Val{T}}=Val{false}) where T =
                                          _change_t_via_interpolation!(integrator,t,modify_save_endpoint)

function DiffEqBase.reeval_internals_due_to_modification!(integrator::ODEIntegrator)
  if integrator.opts.calck
    resize!(integrator.k,integrator.kshortsize) # Reset k for next step!
    alg = unwrap_alg(integrator, false)
    if typeof(alg) <: BS5 || typeof(alg) <: Vern6 || typeof(alg) <: Vern7 ||
       typeof(alg) <: Vern8 || typeof(alg) <: Vern9
       DiffEqBase.addsteps!(integrator,integrator.f,true,false,!alg.lazy)
    else
      DiffEqBase.addsteps!(integrator,integrator.f,true,false)
    end
  end

  if integrator.isdae
    DiffEqBase.initialize_dae!(integrator)
  end

  integrator.u_modified = false
end

function u_modified!(integrator::ODEIntegrator,bool::Bool)
  integrator.u_modified = bool
end

get_proposed_dt(integrator::ODEIntegrator) = ifelse(integrator.opts.adaptive, integrator.dtpropose, integrator.dtcache)
set_proposed_dt!(integrator::ODEIntegrator,dt::Number) = (integrator.dtpropose = dt; integrator.dtcache = dt)

function set_proposed_dt!(integrator::ODEIntegrator,integrator2::ODEIntegrator)
  integrator.dtpropose = integrator2.dtpropose
  integrator.dtcache = integrator2.dtcache
  integrator.qold = integrator2.qold
  integrator.erracc = integrator2.erracc
  integrator.dtacc = integrator2.dtacc
end

@inline function DiffEqBase.get_du(integrator::ODEIntegrator)
  integrator.cache isa FunctionMapCache || integrator.cache isa FunctionMapConstantCache &&
                          error("Derivatives are not defined for this stepper.")
  return if isdefined(integrator, :fsallast)
    integrator.fsallast
  else
    integrator(integrator.t, Val{1})
  end
end

@inline function DiffEqBase.get_du!(out,integrator::ODEIntegrator)
  integrator.cache isa FunctionMapCache || integrator.cache isa FunctionMapConstantCache &&
                          error("Derivatives are not defined for this stepper.")
  if typeof(integrator.cache) <: FunctionMapCache
    out .= integrator.cache.tmp
  else
    return if isdefined(integrator, :fsallast) &&
      !(typeof(integrator.alg) <: Union{Rosenbrock23,Rosenbrock32,Rodas4,Rodas4P,Rodas4P2,Rodas5})
      # Special stiff interpolations do not store the right value in fsallast
      out .= integrator.fsallast
    else
      integrator(out, integrator.t, Val{1})
    end
  end
end

#TODO: Bigger caches for most algorithms
@inline DiffEqBase.get_tmp_cache(integrator::ODEIntegrator) =
          get_tmp_cache(integrator::ODEIntegrator,integrator.alg,integrator.cache)
# avoid method ambiguity
for typ in (OrdinaryDiffEqAlgorithm,Union{RadauIIA3,RadauIIA5},OrdinaryDiffEqNewtonAdaptiveAlgorithm,OrdinaryDiffEqRosenbrockAdaptiveAlgorithm,Union{SSPRK22,SSPRK33,SSPRK53_2N1,SSPRK53_2N2,SSPRK43,SSPRK432,SSPRK932})
  @eval @inline DiffEqBase.get_tmp_cache(integrator,alg::$typ,cache::OrdinaryDiffEqConstantCache) = nothing
end

# the ordering of the cache arrays is important!!!
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqAlgorithm,cache) = (cache.tmp,)
@inline DiffEqBase.get_tmp_cache(integrator,alg::Union{RadauIIA3,RadauIIA5},cache) = (cache.tmp,cache.atmp)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqNewtonAdaptiveAlgorithm,cache) = (cache.nlsolver.tmp,cache.atmp)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqNewtonAlgorithm,cache) = (cache.nlsolver.tmp,cache.nlsolver.z)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqRosenbrockAdaptiveAlgorithm,cache) = (cache.tmp,cache.linsolve_tmp)
@inline DiffEqBase.get_tmp_cache(integrator,alg::Union{SSPRK22,SSPRK33,SSPRK53_2N1,SSPRK53_2N2,SSPRK43,SSPRK432,SSPRK932},cache) = (cache.k,)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqImplicitExtrapolationAlgorithm,cache) = (cache.tmp,cache.utilde)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqAdaptiveExponentialAlgorithm,cache) = (cache.tmp,cache.utilde)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqExponentialAlgorithm,cache) = (cache.tmp,cache.dz)
@inline DiffEqBase.get_tmp_cache(integrator,alg::OrdinaryDiffEqLinearExponentialAlgorithm,cache) = (cache.tmp,)
@inline DiffEqBase.get_tmp_cache(integrator,alg::CompositeAlgorithm, cache) = get_tmp_cache(integrator, integrator.alg.algs[1], cache.caches[1])
@inline DiffEqBase.get_tmp_cache(integrator,alg::DAEAlgorithm,cache) = (cache.nlsolver.cache.dz,cache.atmp)

full_cache(integrator::ODEIntegrator) = full_cache(integrator.cache)
full_cache(integrator::CompositeCache) = Iterators.flatten(full_cache(c) for c in integrator.caches)

function add_tstop!(integrator::ODEIntegrator,t)
  integrator.tdir * (t - integrator.t) < zero(integrator.t) && error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
  push!(integrator.opts.tstops, integrator.tdir * t)
end

function DiffEqBase.add_saveat!(integrator::ODEIntegrator,t)
  integrator.tdir * (t - integrator.t) < zero(integrator.t) && error("Tried to add a saveat that is behind the current time. This is strictly forbidden")
  push!(integrator.opts.saveat, integrator.tdir * t)
end

function resize!(integrator::ODEIntegrator, i::Int)
  @unpack cache = integrator

  for c in full_cache(cache)
    # Skip nothings which may exist in the cache since extra variables
    # may be required for things like units
    c !== nothing && resize!(c,i)
  end
  resize_nlsolver!(integrator, i)
  resize_J_W!(cache, integrator, i)
  resize_non_user_cache!(integrator, cache, i)
end
# we can't use resize!(..., i::Union{Int, NTuple{N,Int}}) where {N} because of method ambiguities with DiffEqBase
function resize!(integrator::ODEIntegrator, i::NTuple{N,Int}) where {N}
  @unpack cache = integrator

  for c in full_cache(cache)
    resize!(c,i)
  end
  # TODO the parts below need to be adapted for implicit methods
  isdefined(integrator.cache, :nlsolver) && resize_nlsolver!(integrator, i)
  resize_J_W!(cache, integrator, i)
  resize_non_user_cache!(integrator, cache, i)
end

function resize_J_W!(cache, integrator, i)
  (isdefined(cache, :J) && isdefined(cache, :W)) || return

  @unpack f = integrator

  if cache.W isa WOperator
    nf = nlsolve_f(f, integrator.alg)
    islin = f isa Union{ODEFunction,SplitFunction} && islinear(nf.f)
    if !islin
      if isa(cache.J, DiffEqBase.AbstractDiffEqLinearOperator)
        resize!(cache.J,i)
      elseif f.jac_prototype !== nothing
        J = similar(f.jac_prototype, i, i)
        J = DiffEqArrayOperator(J; update_func=f.jac)
      elseif cache.J isa SparseDiffTools.JacVec
        resize!(cache.J.cache1,i)
        resize!(cache.J.cache2,i)
        resize!(cache.J.x,i)
      end
      if cache.W.jacvec !== nothing
        resize!(cache.W.jacvec.cache1,i)
        resize!(cache.W.jacvec.cache2,i)
        resize!(cache.W.jacvec.x,i)
      end
      cache.W = WOperator{DiffEqBase.isinplace(integrator.sol.prob)}(
                          f.mass_matrix, integrator.dt, cache.J, integrator.u,
                          cache.W.jacvec;
                          transform = cache.W.transform)
      cache.J = cache.W.J
    end
  else
    if cache.J !== nothing
      cache.J = similar(cache.J, i, i)
    end
    cache.W = similar(cache.W, i, i)
  end

  nothing
end

Base.resize!(p::LinearSolve.LinearCache,i) = p
resize_non_user_cache!(integrator::ODEIntegrator,i::Int) = resize_non_user_cache!(integrator,integrator.cache,i)
deleteat_non_user_cache!(integrator::ODEIntegrator,i) = deleteat_non_user_cache!(integrator,integrator.cache,i)
addat_non_user_cache!(integrator::ODEIntegrator,i) = addat_non_user_cache!(integrator,integrator.cache,i)

resize_non_user_cache!(integrator::ODEIntegrator,cache,i) = nothing

function resize_non_user_cache!(integrator::ODEIntegrator,cache::CompositeCache,i)
  for _cache in cache.caches
    resize_non_user_cache!(integrator,_cache,i)
  end
end

function deleteat_non_user_cache!(integrator::ODEIntegrator,cache::CompositeCache,i)
  for _cache in cache.caches
    deleteat_non_user_cache!(integrator,_cache,i)
  end
end

function addat_non_user_cache!(integrator::ODEIntegrator,cache::CompositeCache,i)
  for _cache in cache.caches
    addat_non_user_cache!(integrator,_cache,i)
  end
end

function resize_non_user_cache!(integrator::ODEIntegrator,
                      cache::RosenbrockMutableCache,i)
  cache.J = similar(cache.J,i,i)
  cache.W = similar(cache.W,i,i)
  resize_jac_config!(cache.jac_config, i)
  resize_grad_config!(cache.grad_config, i)
  nothing
end

function deleteat_non_user_cache!(integrator::ODEIntegrator,cache,idxs)
  # ordering doesn't matter in deterministic cache, so just resize
  # to match the size of u
  i = length(integrator.u)
  resize_non_user_cache!(integrator,cache,i)
end

function addat_non_user_cache!(integrator::ODEIntegrator,cache,idxs)
  # ordering doesn't matter in deterministic cache, so just resize
  # to match the size of u
  i = length(integrator.u)
  resize_non_user_cache!(integrator,cache,i)
end

function deleteat!(integrator::ODEIntegrator,idxs)
  for c in full_cache(integrator)
    deleteat!(c,idxs)
  end
  deleteat_non_user_cache!(integrator,integrator.cache,idxs)
end

function addat!(integrator::ODEIntegrator,idxs)
  for c in full_cache(integrator)
    addat!(c,idxs)
  end
  addat_non_user_cache!(integrator,integrator.cache,idxs)
end

function terminate!(integrator::ODEIntegrator, retcode = :Terminated)
  integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, retcode)
  integrator.opts.tstops.valtree = typeof(integrator.opts.tstops.valtree)()
end

DiffEqBase.has_reinit(integrator::ODEIntegrator) = true
function DiffEqBase.reinit!(integrator::ODEIntegrator,u0 = integrator.sol.prob.u0;
  t0 = integrator.sol.prob.tspan[1], tf = integrator.sol.prob.tspan[2],
  erase_sol = true,
  tstops = integrator.opts.tstops_cache,
  saveat = integrator.opts.saveat_cache,
  d_discontinuities = integrator.opts.d_discontinuities_cache,
  reset_dt = (integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive,
  reinit_callbacks = true, initialize_save = true,
  reinit_cache = true)

  if isinplace(integrator.sol.prob)
    recursivecopy!(integrator.u,u0)
    recursivecopy!(integrator.uprev,integrator.u)
  else
    integrator.u = u0
    integrator.uprev = integrator.u
  end

  if alg_extrapolates(integrator.alg)
    if isinplace(integrator.sol.prob)
      recursivecopy!(integrator.uprev2,integrator.uprev)
    else
      integrator.uprev2 = integrator.uprev
    end
  end

  integrator.t = t0
  integrator.tprev = t0

  tType = typeof(integrator.t)
  tspan = (tType(t0), tType(tf))
  integrator.opts.tstops = initialize_tstops(tType, tstops, d_discontinuities, tspan)
  integrator.opts.saveat = initialize_saveat(tType, saveat, tspan)
  integrator.opts.d_discontinuities =
    initialize_d_discontinuities(tType, d_discontinuities, tspan)

  if erase_sol
    if integrator.opts.save_start
      resize_start = 1
    else
      resize_start = 0
    end
    resize!(integrator.sol.u,resize_start)
    resize!(integrator.sol.t,resize_start)
    resize!(integrator.sol.k,resize_start)
    if integrator.opts.save_start
      copyat_or_push!(integrator.sol.t,1,t0)
      if integrator.opts.save_idxs === nothing
        copyat_or_push!(integrator.sol.u,1,u0)
      else
        u_initial = u0[integrator.opts.save_idxs]
        copyat_or_push!(integrator.sol.u,1,u_initial,Val{false})
      end
    end
    if integrator.sol.u_analytic !== nothing
      resize!(integrator.sol.u_analytic,0)
    end
    if typeof(integrator.alg) <: OrdinaryDiffEqCompositeAlgorithm
      resize!(integrator.sol.alg_choice,resize_start)
    end
    integrator.saveiter = resize_start
    if integrator.opts.dense
      integrator.saveiter_dense = resize_start
    end
  end
  integrator.iter = 0
  integrator.success_iter = 0
  integrator.u_modified = false

  # full re-initialize the PI in timestepping
  reinit!(integrator, integrator.opts.controller)
  integrator.qold = integrator.opts.qoldinit
  integrator.q11 = typeof(integrator.q11)(1)
  integrator.erracc = typeof(integrator.erracc)(1)
  integrator.dtacc = typeof(integrator.dtacc)(1)

  if reset_dt
    auto_dt_reset!(integrator)
  end

  if reinit_callbacks
    initialize_callbacks!(integrator, initialize_save)
  end

  if reinit_cache
    initialize!(integrator,integrator.cache)
  end
end

function DiffEqBase.auto_dt_reset!(integrator::ODEIntegrator)
  integrator.dt = ode_determine_initdt(integrator.u,integrator.t,
  integrator.tdir,integrator.opts.dtmax,integrator.opts.abstol,integrator.opts.reltol,
  integrator.opts.internalnorm,integrator.sol.prob,integrator)
  integrator.dtpropose = integrator.dt
  integrator.destats.nf += 2
end

function DiffEqBase.set_t!(integrator::ODEIntegrator, t::Real)
  if integrator.opts.save_everystep
    error("Integrator time cannot be reset unless it is initialized",
          " with save_everystep=false")
  end
  if alg_extrapolates(integrator.alg) || !isdtchangeable(integrator.alg)
    reinit!(integrator, integrator.u;
            t0 = t,
            reset_dt = false,
            reinit_callbacks = false,
            reinit_cache = false)
  else
    integrator.t = t
  end
end

function DiffEqBase.set_u!(integrator::ODEIntegrator, u)
  if integrator.opts.save_everystep
    error("Integrator state cannot be reset unless it is initialized",
          " with save_everystep=false")
  end
  integrator.u = u
  u_modified!(integrator, true)
end

DiffEqBase.has_destats(i::ODEIntegrator) = true


text:
export LT

"""
    LT{T <: Sampleable{Univariate, Discrete}} <: BinaryCode

Binary Luby Transform erasure correction code. Create a code object
with LT(K, dd).

Arguments

* K::Int: Number of source symbols.

* dd::Sampleable{Univariate, Discrete}: Degree distribution.

"""
struct LT{Tv,Td<:Sampleable{Univariate,Discrete}} <: AbstractErasureCode
    K::Int # number of source symbols
    dd::Td # degree distribution
    function LT{Tv,Td}(K, dd) where {Tv,Td<:Sampleable{Univariate,Discrete}}
        0 < minimum(support(dd)) || throw(ArgumentError("Degree distribution support must be positive"))
        maximum(support(dd)) <= K || throw(ArgumentError("Degree distribution support must be at most K"))
        new(K, dd)
    end
end
function LT{Tv}(K::Integer, dd::Sampleable{Univariate,Discrete}) where Tv
    LT{Tv,typeof(dd)}(K, dd)
end
function Base.show(io::IO, code::LT{Tv}) where Tv
    print(io, "LT{$Tv}($(code.K), $(code.dd)")
end
dimension(code::LT) = code.K

"""

Return a `SparseVector` corresponding to the LT symbol with ESI `X`.
"""
function lt_constraint(code::LT{Tv}, X::Integer) where Tv
    rng = MersenneTwister(X)
    d = rand(rng, code.dd)
    @assert d > 0
    Is = sample(rng, 1:code.K, d, replace=false)
    Vs = rand_nonzero(rng, Tv, d)
    sparsevec(Is, Vs, code.K)
end

generator_matrix(code::LT, Xs::AbstractVector{<:Integer}) = hcat([lt_constraint(code, X) for X in Xs]...)

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
function test_modeconversion_segmented(scenario)
    @unpack distances, tx, bfield, species, ground = scenario()
    params = LMPParams()

    waveguide = SegmentedWaveguide([HomogeneousWaveguide(bfield[i], species[i],
                                        ground[i], distances[i]) for i in 1:2])

    heighttype = typeof(params.wavefieldheights)
    wavefields_vec = Vector{LMP.Wavefields{heighttype}}(undef, 2)
    adjwavefields_vec = Vector{LMP.Wavefields{heighttype}}(undef, 2)

    mesh = LMP.defaultmesh(tx.frequency)
    for i = 1:2
        wvg = waveguide[i]
        adjwvg = LMP.adjoint(wvg)
        modeequation = PhysicalModeEquation(tx.frequency, wvg)
        modes = findmodes(modeequation, mesh; params=params)

        wavefields = LMP.Wavefields(params.wavefieldheights, modes)
        adjwavefields = LMP.Wavefields(params.wavefieldheights, modes)
        LMP.calculate_wavefields!(wavefields, adjwavefields, tx.frequency, wvg, adjwvg)

        wavefields_vec[i] = wavefields
        adjwavefields_vec[i] = adjwavefields
    end

    a = LMP.modeconversion(wavefields_vec[1], wavefields_vec[1], adjwavefields_vec[1])
    @test all(diag(a) .≈ complex(1))
    # a = modeconversion(wavefields_vec[1], wavefields_vec[2], adjwavefields_vec[2])
end


text:
using RunTests
using Base.Test

@testmodule ThirdSetOfTests begin
  function test_hello()
    @test true
  end
end

@testmodule FourthSetOfTests begin
  function test_hello()
    @test true
  end
end


text:
module TestStepSize
    using Base.Test
    import FiniteDiff

    @testset "Per-mode step size functions" begin
        modes = (
            FiniteDiff.ForwardMode,
            FiniteDiff.BackwardMode,
            FiniteDiff.CentralMode,
            FiniteDiff.ComplexMode,
            FiniteDiff.HessianMode,
        )

        for T in (Float32, Float64, BigFloat)
            for mode in modes
                ϵ0 = FiniteDiff.step_size(mode(), zero(T))
                @test isa(ϵ0, T)
                @test ϵ0 <= one(T)

                ϵ1 = FiniteDiff.step_size(mode(), one(T))
                @test isa(ϵ1, T)
                @test ϵ1 <= one(T)

                @test ϵ0 <= ϵ1
            end
        end
    end
end


text:
module MPFirst

# Write your package code here.

end


text:
@testset "Finds the eigenvector for the largest eigenvalue when that eigenvalue is much larger than the others" begin

    matrix = [5 0;0 1]
    guess = LinearAlgebra.normalize([1;0.5])
    stoppingcriteria = Eigen.stoppingcriteria(0.1, 30)

    actual = Eigen.rayleighiteration(matrix, guess, stoppingcriteria)
    expected = [1;0]

    @test LinearAlgebra.norm(actual.eigenvector - expected) <= 0.1
    @test actual.isverified == true
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9af04825-842b-4794-a317-6f1b78a64fb5
begin
	using Statistics, CSV, DataFrames
	using ForwardDiff
end

# ╔═╡ be7baa4f-6871-4d74-b4d2-1c2534cdfa17
# Row Orientation
begin
	using Random
	Random.seed!(0)

	struct Layer_dense
		w::Matrix{Float64}
		b::Matrix{Float64}
		σ::Function

		function Layer_dense(n_in::Int, n_out::Int; σ::Function = identity)
			w = 0.01 .* randn(n_out, n_in)
			b = reshape(zeros(n_out), (1, n_out))
			return new(w, b, σ)
		end
	end

	function (L::Layer_dense)(input)
		W, b, σ = L.w, L.b, L.σ
		σ( (input * W') .+ b )
	end
end

# ╔═╡ c4ea1432-abb4-11eb-1cf5-edac0735d67d
md"# LayerDense"

# ╔═╡ 45d31083-883b-4e84-913e-c040c6fdfd5f
md"""
>No futuro, podemos pensar em novas formas de fazer esse objeto. Por exemplo, passando a matriz W como input construtor do layer.
"""

# ╔═╡ dc76be1b-292e-4545-887e-1e62faf8e1b6
md"""
### Struct by Column Orientation
"""

# ╔═╡ fe237f5f-8f90-429a-8ddc-e5933b8fd808
begin
	Random.seed!(0)

	struct LayerDense{T<:Real, F<:Function}
		W::Matrix{T}
		b::Vector{T}
		σ::F

		function LayerDense(W::Matrix{T},
							b::Vector{T},
							σ::Function = identity) where {T<:Real}
			@assert size(W, 1) == length(b)

			new{T, typeof(σ)}(W, b, σ)
		end
	end

	function LayerDense(W::Matrix{T},
						b::Vector{S},
						σ::Function = identity) where {T<:Real, S<:Real}
		R = promote_type(T, S)
		return LayerDense(Matrix{R}(W), Vector{R}(b), σ)
	end

	function LayerDense(W::Matrix{T}, σ::Function = identity) where {T<:Real}
		b = zeros(T, size(W, 1))
		return LayerDense(W, b, σ)
	end

	function LayerDense(n_in::Int, n_out::Int, σ::Function = identity)
		W = 0.01 * randn(n_out, n_in)
		b = zeros(n_out)
		return LayerDense(W, b, σ)
	end

	function (L::LayerDense)(input::AbstractVector{<:Real})
		W, b, σ = L.W, L.b, L.σ
		return σ(W * input + b)
	end

	function (L::LayerDense)(input::AbstractMatrix{<:Real})
		W, b, σ = L.W, L.b, L.σ
		return σ(W * input .+ b)
	end
end

# ╔═╡ 4e6d165a-9707-4029-af4a-2356b0fe9c41
md"#### Example using the struct LayerDense"

# ╔═╡ 5d3ae0e6-7c5d-48b7-9845-d204a1df9934
W = rand(3, 4)

# ╔═╡ e85bc0f7-65ab-4a41-a4ba-e63c8b7c8184
B = LayerDense(W) # Layer de 4 inputs e 3 outputs

# ╔═╡ 142bb8f9-4188-46d2-94a2-8407d633c496
begin
	obs01 = [2, 3, 4, 5]
	obs02 = [2, 3, 4, 5]
	input = [obs01 obs02] # Each observation must be on a colunm
end

# ╔═╡ 2377f014-0bd0-457b-b9e9-3ace1ad7c331
B(input) # Each column is an output

# ╔═╡ c99c7f4f-7e62-4036-ba08-4478e5480d19
md"#### Adding sigmoid functions"

# ╔═╡ 96ce3aa7-9261-4e60-93ce-4d9321657d85
md"""
Função **Logística**: \
$$\;\;\;\;\;$$ $$\sigma(x) = {1 \over 1 + e^{-x}}$$

Na implementação via código é necessário fazer uma adaptação que evita problemas numéricos de *Overflow* resultando na seguinte função: \
$$\sigma(x) = \left\{
  \begin{array}{lr}
    {1 \over 1 + e^{-|x|}}, & x \ge 0\\
    {e^{-|x|} \over 1 + e^{-|x|}}, & x < 0
  \end{array}
\right.$$ \

Para isso, utilizamos o módulo de $$x$$ na expressão $$e^{-|x|}$$ que possui a seguinte propriedade: \
$$e^{-|x|} = \left\{
  \begin{array}{lr}
    e^{-x}, & x \ge 0\\
    e^x, & x < 0
  \end{array}
\right.$$

Com o modulo, quando x é positivo ou zero ($$x ≥ 0$$) a adaptação resulta na própria função logística, mas quando x é negativo são necessários os seguintes passos: \
$$Se \space x \ge 0:$$  $$\;\;\;\;\;$$  $$\sigma(x) = {1 \over 1 + e^{-|x|}}$$ \
$$Se \space x < 0:$$ $$\;\;\;\;\;$$ $$1 = e^{-x} ⋅ e^x → {1 \over 1 + e^{-x}} = {e^{-x} ⋅ e^x \over e^{-x} ⋅ e^x + e^{-x}} = {e^{-x} ⋅ e^x \over e^{-x}(e^x + 1)} = {e^x \over 1 + e^x} = {e^{-|x|} \over 1 + e^{-|x|}}$$
"""

# ╔═╡ 58265283-a36a-4c83-9101-3387165da021
begin
	function σ(x::Real)
		t = exp(-abs(x))
		return ifelse(x ≥ 0, inv(1 + t), t / (1 + t))
	end

	σ(x::AbstractVecOrMat{T}) where {T<:Real} = σ.(x)

	const sigmoid = σ
end

# ╔═╡ fc1ddc84-62dc-4af1-9b9e-9d37f581a252
md"""
#### ReLU function
$$ReLU(x) = \left\{
  \begin{array}{lr}
    x, & x \ge 0\\
    0, & x < 0
  \end{array}
\right.$$
$$Or$$
$$ReLU(x) = max(0, x)$$
"""

# ╔═╡ 45bafada-34ec-4c41-b0a2-000f35ebffef
begin
	relu(x::Real) = max(zero(x), x) # zero function keep the type of x
	relu(x::AbstractVecOrMat{<:Real}) = relu.(x)
end

# ╔═╡ c85cd202-7223-4c46-80ef-afd750f28ead
md"""
$$softplus(x) = ln(1 + e^x)$$
$$For \space x > 0:$$
$$1 = e^{-x} ⋅ e^x → ln(1 + e^x) = ln(e^{-x} ⋅ e^x + e^x) = ln(e^x(e^{-x} + 1)) = ln(e^x) + ln(1 + e^{-x})$$
$$softplus(x) = x + ln(1 + e^{-x})$$
"""

# ╔═╡ d6fce490-1498-469d-8ecc-55d0a71d1565
begin
	#log1p(x) is the seme as log(1 + x)
	softplus(x::Real) = ifelse(x > 0, x + log(1 + exp(-x)), log(1 + exp(x)))
	softplus(x::AbstractVecOrMat{<:Real}) = softplus.(x)
end

# ╔═╡ 04cadf6a-4fd1-4b5f-aa4d-140ea6e5cff2
function softmax(x::Vector{T}) where {T<:Real}
	m = maximum(x)
	exp_val = exp.(x .- m)
	s = sum(exp_val)
	return exp_val ./ s
end

# ╔═╡ a7898ffd-1acb-4220-bbc8-3b428b95db85
function softmax(x::Matrix{T}) where {T<:Real}
	m = maximum(x, dims = 1)
	exp_val = exp.(x .- m)
	s = sum(exp_val, dims = 1)
	return exp_val ./ s
end

# ╔═╡ c116a1e4-fd4d-446f-bdea-02301f43b2c0
md"##### Example with sigmoid functions"

# ╔═╡ d0c91a1c-8798-4985-978b-4c47d1abc212
C = LayerDense(2, 3, σ)

# ╔═╡ 3de0841e-f0f3-436c-a891-dd4af32d2af6
# C(randn(2, 3))
C([[1, 4] [2, 3] [3, 5]])

# ╔═╡ 963f6276-2440-4bfe-ae0f-539d8bfae0a2
D = LayerDense(2, 3, softmax)

# ╔═╡ 56191b0c-fe49-4ed4-9ca7-2e22e7b028ba
md"""
#### Cost Function
"""

# ╔═╡ 7b5a6859-a9e1-48b1-a0fd-1cde06249d03
begin
	function xlogy(x::Real, y::Real)
		result = x * log(y)
		ifelse(iszero(x) && !isnan(y), zero(result), result)
	end

	function crossentropy(ŷ::AbstractVecOrMat{<:Real},
						  y::AbstractVecOrMat{<:Real};
						  dims::Int = 1,
						  agg::Function = mean)
    	agg(.-sum(xlogy.(y, ŷ), dims = dims))
	end
end

# ╔═╡ 4b070348-fdbc-46cd-85f0-29b64432ae21
begin
	outp2 = D([1 4; 2 3; 3 5]')
	target2 = [0, 1, 0]
	crossentropy(outp2, target2)
end

# ╔═╡ e3ec0333-14ad-42a8-9cd3-4e3937dc5ff0
begin
	outp1 = [0.7 0.1 0.2
			 0.1 0.5 0.4
	   		 0.02 0.9 0.08]'
	target1 = [1 0 0
			   0 1 0
		  	   0 1 0]'
	crossentropy(outp1, target1)
end

# ╔═╡ 8f3cbc01-53ba-4ff0-8245-564f8fee53f2
md"##### Importing data from nnfs python package"

# ╔═╡ 099f5fa8-6687-4c77-9f61-119d75377bab
data = CSV.read("spiral_data.csv", DataFrame)

# ╔═╡ b5981a68-16e7-451f-9e13-d74b2b78a958
md"#### Chain type, params(), gradient()"

# ╔═╡ d2fa0422-7f39-4365-8f3b-bfa6997ae0ae
begin
	Random.seed!(1998)

	L1 = LayerDense(10, 5, σ)
	L2 = LayerDense(5, 2, softmax)

	output = L2(L1(rand(10)))
end

# ╔═╡ 72f68243-8ae1-465c-854d-a76c07c5e346
function params(layers)
	p = []
	for layer in layers
		push!(p, layer.W[:], layer.b )
	end
	return vcat(p...)
end

# ╔═╡ e84cff98-1bea-42cb-8e8d-6f5ef08226f4
begin
	Θ = params([L1, L2])

	U = (
		size(L1.W, 1)*size(L1.W, 2), length(L1.b),
		size(L2.W, 1)*size(L2.W, 2), length(L2.b),
	)
	U = cumsum(U)

	L1w = reshape( Θ[ 1:U[1] ], size(L1.W, 1), size(L1.W, 2))

	L1b = Θ[ U[1]+1:U[2] ]

	L2w = reshape( Θ[ U[2]+1:U[3] ], size(L2.W, 1), size(L2.W, 2))

	L2b = Θ[ U[3]+1:U[4] ]

	( length(Θ), U, (L1w, L1b, L2w, L2b) )
end

# ╔═╡ a5a511f7-38f0-4074-9192-7f11d611bf42
begin

	struct Chain2
		Layers::Vector{LayerDense{Float64}}
		Wl::Vector{Int}
		Wc::Vector{Int}
		b::Vector{Int}

		function Chain(L...)
			Wl = Int[]
			Wc = Int[]
			b = Int[]
			for (i, l) in enumerate(L)
				Wl[i], Wc[i] = size(l.W)
				b[i] = length(l.b)
			end
			return new{}(L, Wl, Wc, b)
		end

	end

end

# ╔═╡ 154e4560-0769-429e-89bd-3b0392209d6d
typeof( (1.3,2.12) )

# ╔═╡ 74ccd8c6-f513-4178-a41f-7b70835516e7
begin
	struct Chain{T<:Tuple}
		layers::T
		Chain(layers...) = new{typeof(layers)}(layers)
	end

	function (L::Chain)(input::AbstractVecOrMat{<:Real})
		return ∘(reverse(L.layers)...)(input)
	end
end

# ╔═╡ 4504eaff-0c86-4606-a410-c94a8e8a653e
begin
	#L1 = LayerDense(10, 5, σ)
	#L2 = LayerDense(5, 2, softmax)

	C1 = Chain(L1, L2)
	typeof(C1.layers)
end

# ╔═╡ ac678c79-a804-44e0-80d8-15a40dff327d
begin
	p = []
	for layer in C1.layers
		push!(p, layer.W[:], layer.b )
	end
	vcat(p...)
end

# ╔═╡ 4b019add-7eb4-496a-b0ef-3c3eee9d85bd
begin
	C1(rand(10,12))
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
ForwardDiff = "~0.10.19"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.0"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.0"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "b5e930ac60b613ef3406da6d4f42c35d8dc51419"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.19"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "3d682c07e6dd250ed082f883dc88aee7996bf2cc"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.0"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "2ca267b08821e86c5ef4376cffed98a46c2cb205"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─c4ea1432-abb4-11eb-1cf5-edac0735d67d
# ╟─45d31083-883b-4e84-913e-c040c6fdfd5f
# ╟─dc76be1b-292e-4545-887e-1e62faf8e1b6
# ╠═fe237f5f-8f90-429a-8ddc-e5933b8fd808
# ╟─4e6d165a-9707-4029-af4a-2356b0fe9c41
# ╠═5d3ae0e6-7c5d-48b7-9845-d204a1df9934
# ╠═e85bc0f7-65ab-4a41-a4ba-e63c8b7c8184
# ╠═142bb8f9-4188-46d2-94a2-8407d633c496
# ╠═2377f014-0bd0-457b-b9e9-3ace1ad7c331
# ╟─c99c7f4f-7e62-4036-ba08-4478e5480d19
# ╟─96ce3aa7-9261-4e60-93ce-4d9321657d85
# ╠═58265283-a36a-4c83-9101-3387165da021
# ╟─fc1ddc84-62dc-4af1-9b9e-9d37f581a252
# ╠═45bafada-34ec-4c41-b0a2-000f35ebffef
# ╟─c85cd202-7223-4c46-80ef-afd750f28ead
# ╠═d6fce490-1498-469d-8ecc-55d0a71d1565
# ╠═04cadf6a-4fd1-4b5f-aa4d-140ea6e5cff2
# ╠═a7898ffd-1acb-4220-bbc8-3b428b95db85
# ╟─c116a1e4-fd4d-446f-bdea-02301f43b2c0
# ╠═d0c91a1c-8798-4985-978b-4c47d1abc212
# ╠═3de0841e-f0f3-436c-a891-dd4af32d2af6
# ╠═963f6276-2440-4bfe-ae0f-539d8bfae0a2
# ╟─56191b0c-fe49-4ed4-9ca7-2e22e7b028ba
# ╠═7b5a6859-a9e1-48b1-a0fd-1cde06249d03
# ╠═4b070348-fdbc-46cd-85f0-29b64432ae21
# ╠═e3ec0333-14ad-42a8-9cd3-4e3937dc5ff0
# ╟─8f3cbc01-53ba-4ff0-8245-564f8fee53f2
# ╟─099f5fa8-6687-4c77-9f61-119d75377bab
# ╠═9af04825-842b-4794-a317-6f1b78a64fb5
# ╟─b5981a68-16e7-451f-9e13-d74b2b78a958
# ╠═d2fa0422-7f39-4365-8f3b-bfa6997ae0ae
# ╠═ac678c79-a804-44e0-80d8-15a40dff327d
# ╠═72f68243-8ae1-465c-854d-a76c07c5e346
# ╠═e84cff98-1bea-42cb-8e8d-6f5ef08226f4
# ╠═a5a511f7-38f0-4074-9192-7f11d611bf42
# ╠═154e4560-0769-429e-89bd-3b0392209d6d
# ╠═74ccd8c6-f513-4178-a41f-7b70835516e7
# ╠═4504eaff-0c86-4606-a410-c94a8e8a653e
# ╠═4b019add-7eb4-496a-b0ef-3c3eee9d85bd
# ╟─be7baa4f-6871-4d74-b4d2-1c2534cdfa17
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Newbe
using Test

@testset "Newbe.jl" begin
    # Write your tests here.
end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
module Tabulate

using   NamedArrays,
        VariableTypes,
        IterTools.Product,
        IterTools.product

export  tabulate

include("./util.jl")
include("./tab.jl")

end # module


text:
"""
  Three-dimensional panel method solver. All these solvers asumes that every
  panel is a planar polygon of any number of vertices higher or equal to three.

  # AUTHORSHIP
    * Author    : Eduardo J. Alvarez
    * Email     : Edo.AlvarezR@gmail.com
    * Created   : Jul 2018
    * License   : MIT License
"""
module PanelSolver

import ForwardDiff
using LinearAlgebra

# GeometricTools from https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
const gt = GeometricTools

const RType = Union{Float64,                    # Concrete real types
                    Int64,
                    ForwardDiff.Dual{Nothing,Float64,3},
                    ForwardDiff.Dual{Nothing,Int64,3}
                    }

const SMOOTH = 1e-2               # Smoothing radius of source panel
const SMOOTH2 = 2e-3              # Smoothing radius for vortex ring
# const SMOOTH2 = 1e-8              # Smoothing radius for vortex ring
const SMOOTH3 = SMOOTH2           # Smoothing radius for semi-infinite vortex
const SMOOTH4 = 1e-5              # Smoothing radius for doublet panel

"""
Calculates and returns the geometric matrix of a collection of constant-source
panels on the given control points with their associated normals.

**ARGUMENTS**
  * `nodes::Array{T,2}`                 : All nodes in the collection of
                                          panels.
  * `panels::Array{Array{Int64,1},1}`   : Node connectivity data defining each
                                          panel, where `panels[i][j]` is the
                                          index in `nodes` of the j-th node in
                                          the i-th panel.
  * `CPs::Array{Array{Float64,1},1}`    : Control points.
  * `normals::Array{Array{Float64,1},1}`: Normal associated to every CP.
"""
function G_constant_source(nodes::Array{T1,2},
                                panels::Array{Array{Int64,1},1},
                                CPs::Array{Array{T2,1},1},
                                normals::Array{Array{T3,1},1}
                            ) where{T1<:RType, T2<:RType, T3<:RType}
  N = size(panels, 1)
  G = Array{Float64}(undef, N, N)

  # Builds geometric matrix
  for j in 1:N # Iterates over columns (panels)
      Vconstant_source(
                        [nodes[:,ind] for ind in panels[j]], # Nodes in j-th panel
                        1.0,                               # Unitary strength,
                        CPs,                               # Targets
                        view(G, :, j);                     # Velocity of j-th
                                                           # panel on every CP
                        dot_with=normals                   # Normal of every CP
                      )
  end

  return G
end

"""
Returns the velocity induced by a panel of vertices `nodes` and constant
strength source `strength` on the targets `targets`. It adds the velocity at the
i-th target to out[i].
"""
function Vconstant_source(nodes::Array{Array{T1,1},1}, strength::RType,
                          targets::Array{Array{T2,1},1},
                          # out::Array{Array{T,1},1}
                          out;
                          dot_with=nothing
                          ) where{T1<:RType, T2<:RType}
  if size(out)!=size(targets)
    error("Invalid `out` argument."*
          " Expected size $(size(targets)), got $(size(out)).")
  end

  nn = size(nodes, 1)                      # Number of nodes

  # Tangent, oblique, and normal vectors
  t, o, n = gt._calc_unitvectors(nodes)

  # Coordinate system defined by Hess & Smith
  unitxi, uniteta, unitz = o, t, -n        # Unit vectors
  O = nodes[1]                             # Origin
  Oaxis = hcat(unitxi, uniteta, unitz)'    # Transformation matrix

  # Converts nodes to H&S coordinate system
  HSnodes = [Oaxis*(node-O) for node in nodes]

  # Iterates over targets
  for ti in 1:size(targets, 1)
    HSX = Oaxis*(targets[ti]-O)
    V = zeros(T2, 3)
    dtheta = 2*pi

    nR0 = 0

    for i in 1:nn
      xi, xj = HSnodes[i], HSnodes[i%nn + 1]

      dij = norm(xj-xi)
      ri = norm(HSX-xi)
      rj = norm(HSX-xj)

      #   println("ri,rj,dij=$ri,$rj,$dij")

      Qij = log( (ri+rj+dij)/(ri+rj-dij + SMOOTH) )

      Sij = (xj[2]-xi[2])/dij
      Cij = (xj[1]-xi[1])/dij

      siji = (xi[1]-HSX[1])*Cij + (xi[2]-HSX[2])*Sij
      sijj = (xj[1]-HSX[1])*Cij + (xj[2]-HSX[2])*Sij
      Rij = (HSX[1]-xi[1])*Sij - (HSX[2]-xi[2])*Cij

      Jij = atan( Rij*abs(HSX[3])*( ri*sijj - rj*siji ) ,
                   ri*rj*Rij^2 + HSX[3]^2*sijj*siji)

      #  println("Rij,ri,rj,siji,sijj,HSX=$Rij,$ri,$rj,$siji,$sijj,$HSX")
      #  println("Rij,ri,rj,siji,sijj,HSX,Sij,Cij=$Rij,$ri,$rj,$siji,$sijj,$HSX,$Sij,$Cij")

      V[1] -= Sij*Qij
      V[2] += Cij*Qij
      V[3] -= Jij

      dtheta *= Rij>=0
      nR0 += Rij==0
    end

    V[3] += dtheta
    V[3] *= sign(HSX[3])   # Isn't this sign already accounted for in atan2?
    V[3] *= !(nR0>1)       # Singularity fix of any z position aligned with node

    if dot_with!=nothing
      out[ti] += dot( strength*(V[1]*unitxi + V[2]*uniteta + V[3]*unitz),
                                                                  dot_with[ti])
    else
      out[ti] += strength*(V[1]*unitxi + V[2]*uniteta + V[3]*unitz)
    end

  end
end

"""
Calculates and returns the geometric matrix of a collection of constant-doublet
panels on the given control points with their associated normals.

**ARGUMENTS**
  * `nodes::Array{T,2}`                 : All nodes in the collection of
                                          panels.
  * `panels::Array{Array{Int64,1},1}`   : Node connectivity data defining each
                                          panel, where `panels[i][j]` is the
                                          index in `nodes` of the j-th node in
                                          the i-th panel.
  * `CPs::Array{Array{Float64,1},1}`    : Control points.
  * `normals::Array{Array{Float64,1},1}`: Normal associated to every CP.
"""
function G_constant_doublet(nodes::Array{T1,2},
                                panels::Array{Array{Int64,1},1},
                                CPs::Array{Array{T2,1},1},
                                normals::Array{Array{T3,1},1}
                            ) where{T1<:RType, T2<:RType, T3<:RType}
  N = size(panels, 1)
  G = Array{Float64}(undef, N, N)

  # Builds geometric matrix
  for j in 1:N # Iterates over columns (panels)
      Vconstant_doublet(
                        [nodes[:,ind] for ind in panels[j]], # Nodes in j-th panel
                        1.0,                               # Unitary strength,
                        CPs,                               # Targets
                        view(G, :, j);                     # Velocity of j-th
                                                           # panel on every CP
                        dot_with=normals                   # Normal of every CP
                      )
  end

  return G
end

"""
Returns the velocity induced by a panel of vertices `nodes` and constant
strength doublet `strength` on the targets `targets`. It adds the velocity at
the i-th target to out[i].
"""
function Vconstant_doublet(nodes::Array{Array{T1,1},1}, strength::RType,
                          targets::Array{Array{T2,1},1},
                          out;
                          dot_with::Union{Array{Array{T3,1},1}, Nothing}=nothing
                          ) where{T1<:RType, T2<:RType, T3<:RType}
  if size(out)!=size(targets)
    error("Invalid `out` argument."*
          " Expected size $(size(targets)), got $(size(out)).")
  end

  nn = size(nodes, 1)                      # Number of nodes

  # Tangent, oblique, and normal vectors
  t, o, n = gt._calc_unitvectors(nodes)

  # Coordinate system defined by Hess & Smith
  unitxi, uniteta, unitz = o, t, -n        # Unit vectors
  O = nodes[1]                             # Origin
  Oaxis = hcat(unitxi, uniteta, unitz)'    # Transformation matrix

  # Converts nodes to H&S coordinate system
  HSnodes = [Oaxis*(node-O) for node in nodes]

  # Iterates over targets
  for ti in 1:size(targets, 1)
    HSX = Oaxis*(targets[ti]-O)
    V = Array{T2}(undef, 3)

    for i in 1:nn
      xi, xj = HSnodes[i], HSnodes[i%nn + 1]

      ri = norm(HSX-xi)
      rj = norm(HSX-xj)

      if (ri>SMOOTH4 && rj>SMOOTH4 &&
                        abs( ri*rj - dot(HSX-xi, HSX-xj) )>SMOOTH4*SMOOTH4)

        aux = (ri+rj)/( ri*rj*( ri*rj - dot(HSX-xi, HSX-xj) ) )
        V[1] -= HSX[3]*(xj[2]-xi[2])*aux
        V[2] += HSX[3]*(xj[1]-xi[1])*aux
        V[3] += ( (HSX[1]-xj[1])*(HSX[2]-xi[2]) -
                                        (HSX[1]-xi[1])*(HSX[2]-xj[2]) )*aux
      end

    end

    if dot_with!=nothing
      out[ti] += strength/(4*pi)*dot( V[1]*unitxi + V[2]*uniteta + V[3]*unitz,
                                                                  dot_with[ti])
    else
      out[ti] += strength/(4*pi)*(V[1]*unitxi + V[2]*uniteta + V[3]*unitz)
    end

  end
end

"""
Calculates and returns the geometric matrix of a collection of vortex-ring
panels of a lifting body on the given control points with their associated
normals.

**ARGUMENTS**
  * `nodes::Array{T,2}`                 : All nodes in the collection of
                                          panels.
  * `panels::Array{Array{Int64,1},1}`   : Node connectivity data defining each
                                          panel, where `panels[i][j]` is the
                                          index in `nodes` of the j-th node in
                                          the i-th panel.
  * `U::Array{Int64,1}`                 : Indices of all panels along the upper
                                          side of the trailing edge.
  * `L::Array{Int64,1}`                 : Indices of all panels along the lower
                                          side of the trailing edge.
  * `CPs::Array{Array{Float64,1},1}`    : Control points.
  * `normals::Array{Array{Float64,1},1}`: Normal associated to every CP.
"""
function G_lifting_vortexring(nodes::Array{T1,2},
                                panels::Array{Array{Int64,1},1},
                                U::Array{Int64,1},
                                L::Array{Int64,1},
                                CPs::Array{Array{T2,1},1},
                                normals::Array{Array{T3,1},1}
                            ) where{T1<:RType, T2<:RType, T3<:RType}

  # Sorts trailing edge indices
  _U = sort(U)
  _L = sort(L)

  cur_u = 1                 # Index of current upper trailing edge cell
  cur_l = 1                 # Index of current lower trailing edge cell

  N = size(panels, 1)
  G = Array{Float64}(undef, N, N)

  # Builds geometric matrix --- Vortex rings
  for j in 1:N # Iterates over columns (panels)

    Vvortexring(# Nodes in j-th panel: puts first node last to match TE
                # [nodes[:,ind] for ind in vcat(panels[j][2:end], panels[j][1])],
                # Nodes in j-th panel: puts last node first to match TE
                [nodes[:,ind] for ind in vcat(panels[j][end], panels[j][1:end-1])],
                # Unitary strength,
                1.0,
                # Targets
                CPs,
                # Velocity of j-th panel on every CP
                view(G, :, j);
                # Normal of every CP
                dot_with=normals,
                # Checks for TE
                closed_ring= !( (cur_u<=size(_U,1) && j==_U[cur_u]) ||
                                (cur_l<=size(_L,1) && j==_L[cur_l]))
                )

      if cur_u<=size(_U,1) && j==_U[cur_u]; cur_u+=1; end;
      if cur_l<=size(_L,1) && j==_L[cur_l]; cur_l+=1; end;
  end

  return G
end

"""
Calculates and returns the geometric matrix of a collection of vortex-ring
panels on the given control points with their associated normals.

**ARGUMENTS**
  * `nodes::Array{T,2}`                 : All nodes in the collection of
                                          panels.
  * `panels::Array{Array{Int64,1},1}`   : Node connectivity data defining each
                                          panel, where `panels[i][j]` is the
                                          index in `nodes` of the j-th node in
                                          the i-th panel.
  * `CPs::Array{Array{Float64,1},1}`    : Control points.
  * `normals::Array{Array{Float64,1},1}`: Normal associated to every CP.
"""
function G_vortexring(nodes::Array{T1,2},
                                panels::Array{Array{Int64,1},1},
                                CPs::Array{Array{T2,1},1},
                                normals::Array{Array{T3,1},1}
                        ) where{T1<:RType, T2<:RType, T3<:RType}
  N = size(panels, 1)
  G = Array{Float64}(undef, N, N)

  # Builds geometric matrix
  for j in 1:N # Iterates over columns (panels)
      Vvortexring(
                        [nodes[:,ind] for ind in panels[j]], # Nodes in j-th panel
                        1.0,                               # Unitary strength,
                        CPs,                               # Targets
                        view(G, :, j);                     # Velocity of j-th
                                                           # panel on every CP
                        dot_with=normals                   # Normal of every CP
                      )
  end

  return G
end

"""
Returns the velocity induced by a vortex ring panel of vertices `nodes` and
vortex strength `strength` on the targets `targets`. It adds the velocity at the
i-th target to out[i].
"""
function Vvortexring(nodes::Array{Arr1,1}, strength::RType,
                          targets::Array{Array{T2,1},1},
                          out;
                          dot_with::Union{Array{Array{T3,1},1}, Nothing}=nothing,
                          closed_ring::Bool=true
                          ) where{Arr1<:AbstractArray, T2<:RType, T3<:RType}
  if size(out)!=size(targets)
    error("Invalid `out` argument."*
          " Expected size $(size(targets)), got $(size(out)).")
  end

  nn = size(nodes, 1)                      # Number of nodes

  # Iterates over targets
  for ti in 1:size(targets, 1)
    V = zeros(3)

    for i in 1:(nn - 1*!closed_ring)
      p1, p2 = nodes[i], nodes[i%nn + 1]
      r1 = targets[ti] - p1
      r2 = targets[ti] - p2
      crossr1r2 = cross(r1,r2)
      # This if statement avoids the singularity at the vortex line
      if dot(crossr1r2,crossr1r2) > SMOOTH2*SMOOTH2
        V += crossr1r2/dot(crossr1r2,crossr1r2) * dot(
                                              (p1-p2), r1/norm(r1) - r2/norm(r2) )
      end
    end

    if dot_with!=nothing
      out[ti] -= strength/(4*pi)*dot(V, dot_with[ti])
    else
      out[ti] -= strength/(4*pi)*V
    end

  end
end

"""
Returns the velocity induced by a semi-infinite vortex starting at point `p` in
the unitary direction `D` and vortex strength `strength` on the targets
`targets`. It adds the velocity at the i-th target to out[i].
"""
function Vsemiinfinitevortex(p::Array{T1,1}, D::Array{T2,1}, strength::RType,
                              targets::Array{Array{T3,1},1},
                              out;
                              dot_with::Union{Array{Array{T3,1},1}, Nothing}=nothing,
                              check::Bool=true
                              ) where{T1<:RType, T2<:RType, T3<:RType}
  # ERROR CASES
  if size(out)!=size(targets)
    error("Invalid `out` argument."*
          " Expected size $(size(targets)), got $(size(out)).")
  elseif check && abs(norm(D)-1.0)>1e-8
    error("Received non-unitary infinite direction D! (norm(D) = $(norm(D)))")
  end

  # Iterates over targets
  for ti in 1:size(targets, 1)

    p2 = p + dot(targets[ti]-p, D)*D
    h = norm(targets[ti]-p2)

    if h>SMOOTH3

      # Adds semi-infinite section
      if dot_with!=nothing
        out[ti] += dot(
                   strength / (4*pi*h) * cross(D, (targets[ti]-p2)/h ),
                                                                  dot_with[ti])
      else
        out[ti] += strength / (4*pi*h) * cross(D, (targets[ti]-p2)/h )
      end

      # Adds bound vortex section
      Vvortexring([p, p2], strength, targets[ti:ti], view(out, ti:ti);
                                          dot_with=dot_with!=nothing ? dot_with[ti:ti] : nothing,
                                          closed_ring=false)
    end

  end
end

end # END OF MODULE


text:
const IM_INTERRUPTED = -1
const IM_ALL = -1
const IM_COUNTER_CB = 0
const IM_RESOLUTION_CB = 1
const IM_GIF_TRANSPARENT_COLOR_CB = 0
const IM_TIF_IMAGE_DESCRIPTION_CB = 0
# begin enum imDataType
typealias imDataType UInt32
const IM_BYTE = 0
const IM_SHORT = 1
const IM_USHORT = 2
const IM_INT = 3
const IM_FLOAT = 4
const IM_CFLOAT = 5
# end enum imDataType
# begin enum imColorSpace
typealias imColorSpace UInt32
const IM_RGB = 0
const IM_MAP = 1
const IM_GRAY = 2
const IM_BINARY = 3
const IM_CMYK = 4
const IM_YCBCR = 5
const IM_LAB = 6
const IM_LUV = 7
const IM_XYZ = 8
# end enum imColorSpace
# begin enum imColorModeConfig
typealias imColorModeConfig UInt32
const IM_ALPHA = 256
const IM_PACKED = 512
const IM_TOPDOWN = 1024
# end enum imColorModeConfig
# begin enum imErrorCodes
typealias imErrorCodes UInt32
const IM_ERR_NONE = 0
const IM_ERR_OPEN = 1
const IM_ERR_ACCESS = 2
const IM_ERR_FORMAT = 3
const IM_ERR_DATA = 4
const IM_ERR_COMPRESS = 5
const IM_ERR_MEM = 6
const IM_ERR_COUNTER = 7
const IM_ERR_READ = IM_ERR_ACCESS
const IM_ERR_WRITE = IM_ERR_ACCESS
const IM_ERR_TYPE = IM_ERR_DATA
const IM_ERR_COMP = IM_ERR_COMPRESS
# end enum imErrorCodes
typealias imFile Void
# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const IM_BMP = 0
const IM_PCX = 1
const IM_GIF = 2
const IM_TIF = 3
const IM_RAS = 4
const IM_SGI = 5
const IM_JPG = 6
const IM_LED = 7
const IM_TGA = 8
# end enum ANONYMOUS_1
# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const IM_NONE = 0
const IM_DEFAULT = 256
const IM_COMPRESSED = 512
# end enum ANONYMOUS_2
typealias imCallback Ptr{Void}
typealias imFileCounterCallback Ptr{Void}
typealias imResolutionCallback Ptr{Void}
# begin enum ANONYMOUS_3
typealias ANONYMOUS_3 UInt32
const IM_RES_NONE = 0
const IM_RES_DPI = 1
const IM_RES_DPC = 2
# end enum ANONYMOUS_3
typealias imGifTranspIndex Ptr{Void}
typealias imTiffImageDesc Ptr{Void}


text:
import Base: size, show, summary, similar, parent, parentindices

const NDIM = 2

"""
    ScalarData <: PointData

A wrapper for a one-dimensional array of scalar-valued data. The resulting wrapper
can be indexed in the same way as the array itself.

# Constructors
- `ScalarData(d::AbstractVector[,dtype=Float64])` constructs a wrapper for the one-dimensional array of data `d`
- `ScalarData(n::Int)` constructs a wrapper for an array of zeros of length `n`.
- `ScalarData(x::PointData)` constructs a wrapper for an array of zeros of the
   same length as that wrapped by `x`.
- `ScalarData(n::Int,dtype=ComplexF64)` constructs a wrapper for complex-valued data.
- `ScalarData(x::PointData,dtype=ComplexF64)` constructs a wrapper for an array of
   complex zeros of the same length as that wrapped by `x`.

# Example

```jldoctest
julia> f = ScalarData(10);

julia> f[5] = 1.0;

julia> f
10 points of scalar-valued Float64 data
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
```
"""
struct ScalarData{N,T,DT<:AbstractVector} <: PointData{N,T}
    data::DT
    ScalarData{N,T,DT}(a::S) where {N,T,DT<:AbstractVector,S<:PointData} = new{length(a),T,typeof(a.data)}(a.data)
    ScalarData{N,T,DT}(a::S) where {N,T,DT<:AbstractVector,S<:AbstractVector} = new{length(a),T,typeof(a)}(a)
end

#@wraparray ScalarData data 1

function ScalarData(data::AbstractVector{T}) where {T <: Number}
  ScalarData{length(data),T,typeof(data)}(data)
end

ScalarData(n::Int;dtype=Float64) = ScalarData(zeros(dtype,n))

"""
    VectorData <: PointData

A wrapper for a one-dimensional array of two-component vector-valued data. The
resulting wrapper can be indexed as though the first component and second
component are stacked on top of each other.

# Constructors
- `VectorData(d::AbstractVector[,dtype=Float64])` constructs a wrapper for the one-dimensional array of data `d`, splitting `d` into the `u` and `v` components evenly.
- `VectorData(u::AbstractVector,v::AbstractVector)` constructs a wrapper for the vector components data `u` and `v`.
- `VectorData(n::Int)` constructs a wrapper with zeros of length `n` for both components.
- `VectorData(x::PointData)` constructs a wrapper for zero components of the
   same length as that wrapped by `x`.
- `VectorData(n::Int,dtype=ComplexF64)` constructs a wrapper with complex-valued zeros
   of length `n` for both components.

# Example

```jldoctest
julia> f = VectorData(10,dtype=ComplexF64);

julia> f.v[1:5] = 1:5;

julia> f
10 points of vector-valued Complex{Float64} data
10×2 Array{Complex{Float64},2}:
 0.0+0.0im  1.0+0.0im
 0.0+0.0im  2.0+0.0im
 0.0+0.0im  3.0+0.0im
 0.0+0.0im  4.0+0.0im
 0.0+0.0im  5.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im

julia> f[7] = 1im; f[18] = 0.2;

julia> f
10 points of vector-valued Complex{Float64} data
10×2 Array{Complex{Float64},2}:
 0.0+0.0im  1.0+0.0im
 0.0+0.0im  2.0+0.0im
 0.0+0.0im  3.0+0.0im
 0.0+0.0im  4.0+0.0im
 0.0+0.0im  5.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+1.0im  0.0+0.0im
 0.0+0.0im  0.2+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
```
"""
struct VectorData{N,T,DT<:AbstractVector} <: PointData{N,T}
    data::DT
    u::ScalarData{N,T}
    v::ScalarData{N,T}
end

function VectorData(data::AbstractVector{T}) where {T <: Number}
  nc = NDIM
  @assert mod(length(data),nc) == 0
  numpts = div(length(data),nc)
  u = ScalarData(view(data,1:numpts))
  v = ScalarData(view(data,numpts+1:length(data)))
  VectorData{numpts,T,typeof(data)}(data,u,v)
end

function VectorData(u::AbstractVector{T},v::AbstractVector{T}) where {T <: Number}
  @assert length(u) == length(v)
  VectorData([u;v])
end

VectorData(x::Tuple{AbstractVector{T},AbstractVector{T}}) where {T <: Number} = VectorData(x...)
VectorData(n::Int;dtype=Float64) = VectorData(zeros(dtype,NDIM*n))

#=
A note on the TensorData data type: if we were to write the components as a matrix,
then it would be:

du/dx  dv/dx
du/dy  dv/dy

or, (grad v)_ij = partial v_j/partial i. This ensures that divergence of
the gradient correctly leads to the Laplacian.
=#

"""
    TensorData <: PointData

A wrapper for a one-dimensional array of 2x2 tensor-valued data, with fields
`dudx`, `dudy`, `dvdx`, `dvdy`. The resulting wrapper can be indexed as though these four components are stacked
on top of each other.

# Constructors
- `TensorData(d::AbstractVector[,dtype=Float64])` constructs a wrapper for the one-dimensional array of data `d`, splitting `d` into the four components evenly.
- `TensorData(dudx,dudy,dvdx,dvdy)` constructs a wrapper for the tensor components data, each of type `AbstractVector`
- `TensorData(n::Int)` constructs a wrapper with zeros of length `n` for all components.
- `TensorData(x::PointData[,dtype=Float64])` constructs a wrapper for zero components of the
   same length as that wrapped by `x`.

# Example

"""
struct TensorData{N,T,DT<:AbstractVector} <: PointData{N,T}
    data::DT
    dudx::ScalarData{N,T}
    dudy::ScalarData{N,T}
    dvdx::ScalarData{N,T}
    dvdy::ScalarData{N,T}
end

function TensorData(data::AbstractVector{T}) where {T <: Number}
  nc = NDIM*NDIM
  @assert mod(length(data),nc) == 0
  numpts = div(length(data),nc)
  dudx = ScalarData(view(data,         1:numpts))
  dudy = ScalarData(view(data,  numpts+1:2*numpts))
  dvdx = ScalarData(view(data,2*numpts+1:3*numpts))
  dvdy = ScalarData(view(data,3*numpts+1:4*numpts))
  TensorData{numpts,T,typeof(data)}(data,dudx,dudy,dvdx,dvdy)
end

function TensorData(dudx::AbstractVector{T},dudy::AbstractVector{T},
                dvdx::AbstractVector{T},dvdy::AbstractVector{T}) where {T <: Number}
  @assert length(dudx) == length(dudy) == length(dvdx) == length(dvdy)
  TensorData([dudx;dudy;dvdx;dvdy])
end

TensorData(x::NTuple{NDIM*NDIM,AbstractVector{T}}) where {T <: Number} = TensorData(x...)
TensorData(n::Int;dtype=Float64) = TensorData(zeros(dtype,n*NDIM*NDIM))

for f in (:ScalarData,:VectorData,:TensorData)
  @eval (::Type{$f{N,T,DT}})() where {N,T<:Number,DT<:AbstractVector} = $f(N,dtype=T)
  @eval $f(::PointData{N,T};dtype=T) where {N,T<:Number} = $f(N,dtype=dtype)
  @eval similar(::$f{N,T,DT};element_type=T) where {N,T<:Number,DT<:AbstractVector} = $f(N,dtype=element_type)
end

for f in (:VectorData,:TensorData)
  # if argument u is an AbstractVector, then make this new instance a point to it
  @eval (::Type{$f{N,T,DT}})(u::AbstractVector) where {N,T<:Number,DT<:AbstractVector} = $f(u)
  # if argument u is the same PointData type, then make this new instance
  # a pointer to the `data` field of u
  @eval (::Type{$f{N,T,DT}})(u::$f{N,T,DT}) where {N,T<:Number,DT<:AbstractVector} = $f(u.data)
end

parent(A::PointData) = A.data
parentindices(A::PointData) = parentindices(A.data)

@propagate_inbounds Base.getindex(A::PointData,i::Int) = getindex(A.data,i)
@propagate_inbounds Base.setindex!(A::PointData, v, i::Int) = setindex!(A.data,v,i)

size(A::PointData{N,T},d::Int) where {N,T} = size(A.data,d)
size(A::PointData{N,T}) where {N,T} = size(A.data)

#=
"""
    size(A::PointData,d::Int) -> Int

Return twice the number of vector data points if `d` is 1 (the sum of the length of the `u` and
`v` vectors) and 1 if `d` is 2. This is consistent with the interpretation of VectorData as
a stacked pair of columns, corresponding to the `u` and `v` components, respectively.
"""
Base.size(::VectorData{N,T},d::Int) where {N,T} = d == 1 ? NDIM*N : 1

"""
    size(A::VectorData) -> Tuple

Return a tuple of the number of vector data points by the number of dimensions.
"""
Base.size(A::VectorData) = (size(A,1),)

#@propagate_inbounds Base.getindex(A::VectorData{N,T},i::Int) where {N,T} =
#   i > N ? A.data[i-N] : A.u[i]
#@propagate_inbounds Base.setindex!(A::VectorData{N,T}, v, i::Int) where {N,T} =
#   i > N ? A.v[i-N] = convert(T, v) : A.u[i] = convert(T, v)

"""
    size(A::TensorData,d::Int) -> Int

Return four times the number of tensor data points if `d` is 1 (the sum of the length of the four components)
and 1 if `d` is 2. This is consistent with the interpretation of TensorData as
a stacked set of columns.
"""
Base.size(::TensorData{N,T},d::Int) where {N,T} = d == 1 ? NDIM*NDIM*N : 1

"""
    size(A::TensorData) -> Tuple

Return a tuple of the number of tensor data points by the number of dimensions.
"""
Base.size(A::TensorData) = (size(A,1),)
=#

function show(io::IO, m::MIME"text/plain", pts::ScalarData{N,T}) where {N,T}
  println(io,"$N points of scalar-valued $T data")
  show(io,m,pts.data)
end

function show(io::IO, m::MIME"text/plain", pts::VectorData{N,T}) where {N,T}
  println(io,"$N points of vector-valued $T data")
  show(io,m,pts.data)
end

function show(io::IO, m::MIME"text/plain", pts::TensorData{N,T}) where {N,T}
  println(io,"$N points of tensor-valued $T data")
  show(io,m,pts.data)
end

include("basicpointoperations.jl")

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids