text:
"""
    update_wake_shedding_locations!(wakes, wake_shedding_locations,
        surfaces, ref, fs, dt, additional_velocity, Vte, nwake, eta)

Update the wake shedding locations.  Also update the first chordwise wake panels
to account for the new wake shedding location

# Arguments
 - `wakes`: Vector of wakes corresponding to each surface, represented by matrices
    of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns) where `nw` is the
    number of chordwise wake panels and `ns` is the number of spanwise panels.
 - `wake_shedding_locations`: Shedding location coordinates for each surface for
    each trailing edge vertex.
 - `surfaces`: Vector of surfaces, represented by matrices of surface panels
    (see [`SurfacePanel`](@ref) of shape (nc, ns) where `nc` is the number of
    chordwise panels and `ns` is the number of spanwise panels
 - `reference`: Reference parameters (see [`Reference`](@ref))
 - `freestream`: Freestream parameters (see [`Freestream`](@ref))
 - `dt`: Time step (seconds)
 - `additional_velocity`: Function defining additional velocity field
 - `Vte`: Velocity experienced at the trailing edge due to surface motion.
 - `nwake`: Number of chordwise wake panels to use from each wake in `wakes`
 - `eta`: Time step fraction used to define separation between trailing
    edge and wake shedding location.  Typical values range from 0.2-0.3.
"""
@inline function update_wake_shedding_locations!(wakes, wake_shedding_locations,
    surfaces, ref, fs, dt, additional_velocity, Vte, nwake, eta)

    # get number of surfaces
    nsurf = length(surfaces)

    # loop through all surfaces
    for isurf = 1:nsurf

        # number of spanwise panels
        ns = size(wakes[isurf], 2)

        # update wake shedding location
        for j = 1:ns+1

            # extract trailing edge coordinate
            if j < ns + 1
                rte = bottom_left(surfaces[isurf][end, j])
            else
                rte = bottom_right(surfaces[isurf][end, j-1])
            end

            # freestream velocity
            V = freestream_velocity(fs)

            # rotational velocity
            V += rotational_velocity(rte, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                V += additional_velocity(rte)
            end

            # velocity due to surface motion
            if !isnothing(Vte)
                V += Vte[isurf][j]
            end

            # update wake shedding location coordinates
            wake_shedding_locations[isurf][j] = rte + eta*V*dt

        end

        if nwake[isurf] > 0
            # loop through first row of wake panels
            for j = 1:ns
                # update wake panel with wake shedding location coordinates
                rtl = wake_shedding_locations[isurf][j]
                rtr = wake_shedding_locations[isurf][j+1]

                # preserve other wake panel coordinates
                rbl = bottom_left(wakes[isurf][1,j])
                rbr = bottom_right(wakes[isurf][1,j])

                # preserve core size
                core_size = get_core_size(wakes[isurf][1,j])

                # preserve circulation strength
                gamma = circulation_strength(wakes[isurf][1,j])

                # replace the old wake panel
                wakes[isurf][1,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)
            end
        end

    end

    return wakes, wake_shedding_locations
end

"""
    get_wake_velocities!(wake_velocities, surfaces, wakes, ref, fs, Γ,
        additional_velocity, Vte, symmetric, repeated_points, nwake,
        surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat)

# Arguments
 - `wake_velocities`: Velocities at the corners of the wake panels in `wakes`
 - `surfaces`: Vector of surfaces, represented by matrices of surface panels
    (see [`SurfacePanel`](@ref) of shape (nc, ns) where `nc` is the number of
    chordwise panels and `ns` is the number of spanwise panels
 - `wakes`: Vector of wakes corresponding to each surface, represented by matrices
    of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns) where `nw` is the
    number of chordwise wake panels and `ns` is the number of spanwise panels.
 - `reference`: Reference parameters (see [`Reference`](@ref))
 - `freestream`: Freestream parameters (see [`Freestream`](@ref))
 - `Γ`: Circulation of all surface panels stored in a single vector
 - `additional_velocity`: Function defining additional velocity field
 - `Vte`: Velocity at the trailing edge vertices on each surface due to surface motion
 - `symmetric`: (required) Flag for each surface indicating whether a mirror
    image across the X-Z plane should be used when calculating induced velocities
 - `repeated_points`: Dictionary of the form `Dict((isurf, i) => [(jsurf1, j1),
    (jsurf2, j2)...]` which defines repeated trailing edge points.  Trailing edge
    point `i` on surface `isurf` is repeated on surface `jsurf1` at point `j1`,
    `jsurf2` at point `j2`, and so forth. See [`repeated_trailing_edge_points`](@ref)
 - `nwake`: Number of chordwise wake panels to use from each wake in `wakes`,
    defaults to all provided wake panels
 - `surface_id`: Surface ID for each surface.  The finite core model is disabled
    when calculating the influence of surfaces/wakes that share the same ID.
 - `wake_finite_core`: Flag for each wake indicating whether the finite core
    model should be enabled when calculating the wake's influence on itself and
    surfaces/wakes with the same surface ID.  Defaults to `true` for each surface.
 - `wake_shedding_locations`: Shedding location coordinates for each surface for
    each trailing edge vertex.
 - `trailing_vortices`: Flags to enable/disable trailing vortices, defaults to
    `true` for each surface
 - `xhat`: Direction in which to shed trailing vortices, defaults to [1, 0, 0]
"""
@inline function get_wake_velocities!(wake_velocities, surfaces, wakes, ref, fs, Γ,
    additional_velocity, Vte, symmetric, repeated_points, nwake,
    surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat)

    # number of surfaces
    nsurf = length(surfaces)

    # loop through all surfaces
    for isurf = 1:nsurf

        # number of chordwise and spanwise panels
        nc, ns = size(surfaces[isurf])
        nw = nwake[isurf]

        # velocity at the wake shedding locations
        for is = 1:ns+1

            # check if this point is a duplicate, skip if it is
            if (isurf, is) in keys(repeated_points)
                for (jsurf, js) in repeated_points[(isurf, is)]
                    # NOTE: we assume that a point is not repeated on the same surface
                    if jsurf < isurf
                        wake_velocities[isurf][1, is] = wake_velocities[jsurf][1, js]
                        continue
                    end
                end
            end

            # get vertex location
            rc = wake_shedding_locations[isurf][is]

            # freestream velocity
            V = freestream_velocity(fs)

            # rotational velocity
            V += rotational_velocity(rc, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                V += additional_velocity(rc)
            end

            # velocity at the trailing edge
            wake_velocities[isurf][1,is] = V
        end

        # velocity at all other wake vertices
        cr = CartesianIndices((2:nw+1, 1:ns+1))

        # loop through all vertices
        for I in cr

            # check if this point is a duplicate, skip if it is
            if (isurf, I[2]) in keys(repeated_points)
                for (jsurf, js) in repeated_points[(isurf, I[2])]
                    # NOTE: we assume that a point is not repeated on the same surface
                    if jsurf < isurf
                        wake_velocities[isurf][I] = wake_velocities[jsurf][I[1], js]
                        continue
                    end
                end
            end

            # get vertex location
            if I[1] <= nw && I[2] <= ns
                rc = top_left(wakes[isurf][I[1], I[2]])
            elseif I[1] == nw + 1 && I[2] <= ns
                rc = bottom_left(wakes[isurf][I[1]-1, I[2]])
            elseif I[1] <= nw && I[2] == ns + 1
                rc = top_right(wakes[isurf][I[1], I[2]-1])
            else # I[1] == nw + 1 && I[2] == ns + 1
                rc = bottom_right(wakes[isurf][I[1]-1, I[2]-1])
            end

            # freestream velocity
            wake_velocities[isurf][I] = freestream_velocity(fs)

            # rotational velocity
            wake_velocities[isurf][I] += rotational_velocity(rc, fs, ref)

            # additional velocity field
            if !isnothing(additional_velocity)
                wake_velocities[isurf][I] += additional_velocity(rc)
            end

            # induced velocity from each surface and wake
            jΓ = 0 # index for accessing Γ
            for jsurf = 1:nsurf

                # number of panels on sending surface
                Ns = length(surfaces[jsurf])

                # check if receiving point is repeated on the sending surface
                if isurf == jsurf
                    # the surfaces are the same
                    same_surface = true
                    # vertex spanwise coordinate
                    js = I[2]
                else
                    # the surfaces are different, but the point could still be a duplicate
                    if (isurf, I[2]) in keys(repeated_points)
                        # the vertex is duplicated

                        # check if the vertex is on the sending surface
                        idx = findfirst(x -> x[1] == jsurf, repeated_points[(isurf, I[2])])

                        if isnothing(idx)
                            # the vertex is not duplicated on the sending surface
                            same_surface = false
                        else
                            # the vertex is duplicated on the sending surface
                            same_surface = true
                            # vertex spanwise coordinate on sending surface
                            js = repeated_points[(isurf, I[2])][idx][2]
                        end
                    else
                        # the vertex is not duplicated

                        # the vertex is not on the sending surface
                        same_surface = false
                    end
                end

                # extract circulation values corresponding to the sending surface
                vΓ = view(Γ, jΓ+1:jΓ+Ns)

                # induced velocity from this surface
                wake_velocities[isurf][I] += induced_velocity(rc, surfaces[jsurf], vΓ;
                    finite_core = surface_id[isurf] != surface_id[jsurf],
                    wake_shedding_locations = wake_shedding_locations[jsurf],
                    symmetric = symmetric[jsurf],
                    trailing_vortices = false,
                    xhat = xhat)

                # add induced velocity from the wake
                if same_surface
                    # vertex location on wake
                    J = CartesianIndex(I[1], js)

                    # induced velocity from wake on its own vertex
                    wake_velocities[isurf][I] += induced_velocity(J, wakes[jsurf];
                        finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf],
                        symmetric = symmetric[jsurf],
                        nc = nwake[jsurf],
                        trailing_vortices = trailing_vortices[jsurf],
                        xhat = xhat)
                else
                    # induced velocity from wake on another wake's vertex
                    wake_velocities[isurf][I] += induced_velocity(rc, wakes[jsurf];
                        finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf],
                        symmetric = symmetric[jsurf],
                        nc = nwake[jsurf],
                        trailing_vortices = trailing_vortices[jsurf],
                        xhat = xhat)
                end

                jΓ += Ns # increment Γ index
            end
        end
    end

    return wake_velocities
end

"""
    translate_wake(panel, wake_velocities, dt)

Return a translated copy of the wake panel `panel` given the wake corner velocities
`wake_velocities` and the time step `dt`

# Arguments
 - `panel`: Wake panel (of type [`WakePanel`](@ref))
 - `wake_velocities`: Matrix containing the velocities at each of the four corners
    of `panel`
 - `dt`: Time step (seconds)
"""
@inline function translate_wake(panel::WakePanel, wake_velocities, dt)

    # extract corners
    rtl = top_left(panel)
    rtr = top_right(panel)
    rbl = bottom_left(panel)
    rbr = bottom_right(panel)

    # get vortex filament length
    lt = norm(rtr - rtl)
    lb = norm(rbl - rbr)
    ll = norm(rtl - rbl)
    lr = norm(rbr - rtr)
    l1 = lt + lb + ll + lr

    # translate corners
    rtl += wake_velocities[1,1]*dt
    rtr += wake_velocities[1,2]*dt
    rbl += wake_velocities[2,1]*dt
    rbr += wake_velocities[2,2]*dt

    # get new vortex filament length
    lt = norm(rtr - rtl)
    lb = norm(rbl - rbr)
    ll = norm(rtl - rbl)
    lr = norm(rbr - rtr)
    l2 = lt + lb + ll + lr

    # use previous core size
    core_size = get_core_size(panel)

    # correct vorticity for vortex stretching
    gamma = circulation_strength(panel)*l1/l2

    return WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)
end

"""
    translate_wake!(wake, wake_velocities, dt; nwake = size(wake, 1))

Translate the wake panels in `wake` given the corner velocities `wake_velocities`
and the time step `dt`.

# Arguments
 - `wake`: Matrix of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns)
    where `nw` is the number of chordwise wake panels and `ns` is the number of
    spanwise panels, defaults to no wake panels
 - `wake_velocities`: Velocities at each of the vertices corresponding to the
    wake panels in `wake`
 - `dt`: Time step

# Keyword Arguments
 - `nwake`: Number of chordwise wake panels to use from `wake`, defaults to all
    provided wake panels
"""
@inline function translate_wake!(wake, wake_velocities, dt; nwake = size(wake, 1))

    nw = nwake
    ns = size(wake, 2)
    cw = CartesianIndices((nw, ns))

    for I in cw

        panel = wake[I]

        vV = view(wake_velocities, I[1]:I[1]+1, I[2]:I[2]+1)

        wake[I] = translate_wake(panel, vV, dt)
    end

    return wake
end

"""
    shed_wake!(wake, wake_shedding_locations, wake_velocities, dt, surface, Γ, nwake)

Shed a new wake panel from the wake shedding locations and translate existing
wake panels.

# Arguments
 - `wake`: Matrix of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns)
    where `nw` is the number of chordwise wake panels and `ns` is the number of
    spanwise panels
 - `wake_shedding_locations`: Vector of length `ns` which stores the coordinates
    where wake panels are shed from the trailing edge of `surface`.
 - `wake_velocities`: Velocities at each of the vertices corresponding to the
    wake panels in `wake`
 - `dt`: Time step (seconds)
 - `surface`: Matrix of surface panels (see [`SurfacePanel`](@ref)) of shape
    (nc, ns) where `nc` is the number of chordwise panels and `ns` is the number
    of spanwise panels
 - `Γ`: Circulation strength of each surface panel in `surface`
 - `nwake`: Number of chordwise wake panels to use from `wake`, defaults to all
    provided wake panels
"""
@inline function shed_wake!(wake::AbstractMatrix, wake_shedding_locations,
    wake_velocities, dt, surface, Γ, nwake)

    nc, ns = size(surface)
    nw = size(wake, 1)
    ls = LinearIndices((nc, ns))

    # replace the last chordwise panels with the newly shed wake panels
    for j = 1:ns

        # shedding location coordinates
        rtl = wake_shedding_locations[j]
        rtr = wake_shedding_locations[j+1]

        # shed coordinates
        rbl = rtl + wake_velocities[1, j]*dt
        rbr = rtr + wake_velocities[1, j+1]*dt

        # use core size from the shedding panel
        core_size = get_core_size(surface[end, j])

        # use circulation strength from the shedding panel
        gamma = Γ[ls[end,j]]

        # replace the oldest wake panel
        wake[end,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)
    end

    # translate all existing wake panels except the most recently shed panels
    translate_wake!(wake, wake_velocities, dt, nwake = min(nwake, nw-1))

    # shift wake panels to make newly shed wake panel first
    rowshift!(wake)

    return wake
end

"""
    shed_wake!(wakes, wake_shedding_locations, wake_velocities, dt, surfaces, Γ, nwake)

Shed a new wake panel from the wake shedding locations and translate existing
wake panels.

# Arguments
 - `wakes`: Vector of wakes corresponding to each surface, represented by matrices
    of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns) where `nw` is the
    number of chordwise wake panels and `ns` is the number of spanwise panels.
 - `wake_shedding_locations`: Shedding location coordinates for each surface for
    each trailing edge vertex.
 - `wake_velocities`: Velocities at each of the vertices corresponding to the
    wake panels in `wake`
 - `dt`: Time step (seconds)
 - `surfaces`: Vector of surfaces, represented by matrices of surface panels
    (see [`SurfacePanel`](@ref) of shape (nc, ns) where `nc` is the number of
    chordwise panels and `ns` is the number of spanwise panels
 - `Γ`: Circulation strength of each surface panel in `surfaces`
 - `nwake`: Number of chordwise wake panels to use from each wake in `wakes`,
    defaults to all provided wake panels
"""
@inline function shed_wake!(wakes::AbstractVector{<:AbstractMatrix}, wake_shedding_locations,
    wake_velocities, dt, surfaces, Γ, nwake)

    iΓ = 0
    for i = 1:length(surfaces)

        N = length(surfaces[i])

        vΓ = view(Γ, iΓ+1:iΓ+N)

        shed_wake!(wakes[i], wake_shedding_locations[i], wake_velocities[i], dt,
            surfaces[i], vΓ, nwake[i])

        iΓ += N
    end

    return wakes
end

"""
    rowshift!(A)

Circularly shifts the rows of a matrix down one row.
"""
function rowshift!(A)

    ni, nj = size(A)

    for j = 1:nj
        tmp = A[ni,j]
        for i = ni:-1:2
            A[i,j] = A[i-1,j]
        end
        A[1,j] = tmp
    end

    return A
end


text:
export dot_similarity, normalize!, cosine_similarity, cosine_distance, angle_distance

"""
Computes the dot product between vectors (it may disapear in a near future in favor of LinearAlgebra.dot)
"""
function dot_similarity(a::AbstractVector{T}, b::AbstractVector{T})::Float64 where {T <: Real}
    s = 0.0
    @fastmath @inbounds @simd for i in eachindex(a)
        s += a[i] * b[i]
    end

    s
end

"""
normalize!

Divides the vector by its norm (in-place operation)
"""
function normalize!(a::AbstractVector{T}) where {T <: Real}
    xnorm = sqrt(dot_similarity(a, a))
    (xnorm <= eps(T)) && error("A valid vector for cosine's normalize! cannot have a zero norm $xnorm -- vec: $a")
    invnorm = 1.0 / xnorm
    @fastmath @inbounds @simd for i in eachindex(a)
        a[i] = a[i] * invnorm
    end

    a
end

"""
normalize!

Divides all column vectors of the matrix by its norm (in-place operation)
"""
function normalize!(a::AbstractMatrix{T}) where {T <: Real}
    for i in size(a, 2)
        normalize!(@view a[:, i])
    end

    a
end

"""
angle_distance

Computes the angle between two vectors, it expects normalized vectors (see normalize! method)
"""
function angle_distance(a::AbstractVector{T}, b::AbstractVector{T})::Float64 where {T <: Real}
    m = max(-1.0, dot_similarity(a, b))
    acos(min(1.0, m))
end

"""
cosine_distance

Computes the cosine distance between two vectors, it expects normalized vectors (see normalize! method).
Please use angle_distance if you are expecting a metric function (instead cosine_distance is a faster
alternative whenever the triangle inequality is not needed)
"""
function cosine_distance(a::AbstractVector{T}, b::AbstractVector{T})::Float64 where {T <: Real}
    1 - dot_similarity(a, b)
end

ext: jl


text:
using Core.Compiler: AbstractInterpreter, CodeInstance, MethodInstance, WorldView, NativeInterpreter
using InteractiveUtils

function infer_function(interp, tt)
    # Find all methods that are applicable to these types
    mthds = _methods_by_ftype(tt, -1, typemax(UInt))
    if mthds === false || length(mthds) != 1
        error("Unable to find single applicable method for $tt")
    end

    mtypes, msp, m = mthds[1]

    # Grab the appropriate method instance for these types
    mi = Core.Compiler.specialize_method(m, mtypes, msp)

    # Construct InferenceResult to hold the result,
    result = Core.Compiler.InferenceResult(mi)

    # Create an InferenceState to begin inference, give it a world that is always newest
    world = Core.Compiler.get_world_counter()
    frame = Core.Compiler.InferenceState(result, #=cached=# true, interp)

    # Run type inference on this frame.  Because the interpreter is embedded
    # within this InferenceResult, we don't need to pass the interpreter in.
    Core.Compiler.typeinf(interp, frame)

    # Give the result back
    return (mi, result)
end

struct ExtractingInterpreter <: Core.Compiler.AbstractInterpreter
    code::Dict{MethodInstance, CodeInstance}
    native_interpreter::NativeInterpreter
    msgs::Vector{Tuple{MethodInstance, Int, String}}
    optimize::Bool
end

ExtractingInterpreter(;optimize=false) = ExtractingInterpreter(
    Dict{MethodInstance, Any}(),
    NativeInterpreter(),
    Vector{Tuple{MethodInstance, Int, String}}(),
    optimize
)

import Core.Compiler: InferenceParams, OptimizationParams, get_world_counter,
    get_inference_cache, code_cache,
    WorldView, lock_mi_inference, unlock_mi_inference, InferenceState
InferenceParams(ei::ExtractingInterpreter) = InferenceParams(ei.native_interpreter)
OptimizationParams(ei::ExtractingInterpreter) = OptimizationParams(ei.native_interpreter)
get_world_counter(ei::ExtractingInterpreter) = get_world_counter(ei.native_interpreter)
get_inference_cache(ei::ExtractingInterpreter) = get_inference_cache(ei.native_interpreter)

# No need to do any locking since we're not putting our results into the runtime cache
lock_mi_inference(ei::ExtractingInterpreter, mi::MethodInstance) = nothing
unlock_mi_inference(ei::ExtractingInterpreter, mi::MethodInstance) = nothing

code_cache(ei::ExtractingInterpreter) = ei.code
Core.Compiler.get(a::Dict, b, c) = Base.get(a,b,c)
Core.Compiler.get(a::WorldView{<:Dict}, b, c) = Base.get(a.cache,b,c)
Core.Compiler.haskey(a::Dict, b) = Base.haskey(a, b)
Core.Compiler.haskey(a::WorldView{<:Dict}, b) =
    Core.Compiler.haskey(a.cache, b)
Core.Compiler.setindex!(a::Dict, b, c) = setindex!(a, b, c)
Core.Compiler.may_optimize(ei::ExtractingInterpreter) = ei.optimize
Core.Compiler.may_compress(ei::ExtractingInterpreter) = false
Core.Compiler.may_discard_trees(ei::ExtractingInterpreter) = false

function Core.Compiler.add_remark!(ei::ExtractingInterpreter, sv::InferenceState, msg)
    @show msg
    push!(ei.msgs, (sv.linfo, sv.currpc, msg))
end

macro code_typed_nocache(ex0...)
    esc(:(@InteractiveUtils.code_typed interp=$(ExtractingInterpreter)() $(ex0...)))
end

export @code_typed_nocache


text:
#=

Binary Hierarchical SoftMax

A hierarchical softmax with a binary tree

B. Piwowarski
April 2015

=#

type BinaryHierarchicalSoftmaxUpdates
  size::Int
  offsets::Vector{Int}
  nodes::Vector{Int}
  values
end

# Holds the different structures
type BinaryHierarchicalSoftmax{D<:Device, F<:Float} <: Layer
  # List of parents, using the sign for denoting which branch
  parents::Vector{Int}

  # The number of leaves
  nleaves::Int

  # The depth of each leave
  depth::Vector{Int}

  weight::matrixParameters(D,F)
  bias::vectorParameters(D,F)

  output::denseRealMatrix(D,F)
  grad_input::denseRealVector(D,F)

  updates::BinaryHierarchicalSoftmaxUpdates

  function BinaryHierarchicalSoftmax(inputSize::Int, parents::Vector{Int})
    self = new()

    self.parents = parents
    self.nleaves = div(size(parents, 1) + 1, 2)

    nbInnerNodes = size(parents, 1) - self.nleaves

    @assert(sum(parents) == 0, "Tree is not balanced")
    @assert(parents[size(parents,1)] == 0, "Root is not the last node")

    # Computes the depth at each inner node
    # (this is useful to store the updates)
    function compute_depth(ix::Int)
        if self.depth[ix] == -1 then
           local pix = abs(self.parents[ix + self.nleaves])
           if pix == 0 then
              self.depth[ix] = 0
           else
              self.depth[ix] = compute_depth(pix - self.nleaves) + 1
           end
        end
        return self.depth[ix]
     end

     self.depth = fill(-1, nbInnerNodes)
     for ix = 1:nbInnerNodes
        compute_depth(ix)
     end

     # Initialize the different fields
     self.weight = MatrixParameters{D,F}(nbInnerNodes, inputSize)
     self.bias = VectorParameters{D,F}(nbInnerNodes)

     #  No gradient for the targets
     self.grad_input = Any[ array(D, F), 0 ]

     self.output = array(D, F, 0, 0)

     #  Stores updates
     self.updates = BinaryHierarchicalSoftmaxUpdates(0, Array(Int, 0), Array(Int, 0), array(D, F, 0))

     reset(self)

     return self
  end

end

export BinaryHierarchicalSoftmax

# Resets
function reset{D<:Device, F<:Float}(self::BinaryHierarchicalSoftmax{D, F}, stdv=0)
  if stdv != 0 then
    stdv = stdv * sqrt(3)
  else
    stdv = 1. / sqrt(size(self.weight.values, 2))
  end

  function r!(x)
    rand!(x)
    broadcast!(-, x, x, 0.5)
    broadcast!(*, x, x, stdv)
  end

  r!(self.weight.values)
  r!(self.bias.values)
end

function forward!{F<:Float}(self::BinaryHierarchicalSoftmax{CPUDevice, F}, inputs::Array)
  const D = CPUDevice
  input::matrixOf(D, F), targets::Array{Int, 2} = inputs

  @assert(ndims(input) == 1 || ndims(input) == 2, "input must be vector or matrix")
  @assert(size(input, 1) == size(targets, 1), "Batch size of inputs [$(size(inputs, 1))] and targets [$(size(targets, 1))] differ")

   nframe::Int = size(input, 1)
   if size(self.output) != (nframe, 1)
     self.output = array(D, F, nframe, 1)
   end

   output::matrixOf(D,F) = self.output

    #  Stores the heavy computation before computing grad_input or updating gradient
    tdepth::Int = 0
    for frame = 1:nframe
       local current::Int = self.parents[targets[frame]]
       local ix::Int = abs(current) - self.nleaves
       tdepth = tdepth + self.depth[ix] + 1
    end

    self.updates.size = nframe
    if size(self.updates.offsets, 1) < nframe + 1 then
      resize!(self.updates.offsets, nframe + 1)
    end
    if size(self.updates.nodes, 1) < tdepth
      resize!(self.updates.nodes, tdepth)
      resize!(self.updates.values, tdepth)
    end

    # Create type stable values
    offset::Int = 1
    @stabilize values::vectorOf(D,F) = self.updates.values
    @stabilize weight::matrixOf(D,F) = self.weight.values
    @stabilize bias::vectorOf(D,F) = self.bias.values
    current_p::F = zero(F)

    # ---
    for frame = 1:nframe
       self.updates.offsets[frame] = offset

       local current::Int = self.parents[targets[frame, 1]]
       local logp::Float64 = 0
       while current != 0
          sign::Float64 = 1.
          if current <= 0 then
             sign = -1.
             current = -current
          end

          local ix::Int = current - self.nleaves
          current_p = 1. + exp(sign * (Base.dot(weight[ix], input[frame]) + bias[ix]))

          self.updates.nodes[offset] = current
          values[offset] = sign * (1. / current_p - 1.)

          logp = logp - log(current_p)
          current = self.parents[current]
          offset = offset + 1

       end

       #  Set the output
       output[frame, 1] = logp
    end

    #  Marks the end
    self.updates.offsets[nframe+1] = offset
    return output
 end

 function compute_inputgradient!{F<:Float}(self::BinaryHierarchicalSoftmax{CPUDevice, F}, inputs, gradOutput)
    input, targets = inputs

    if self.grad_input != 0 then
       local _grad_input0 = self.grad_input[1]
       if size(_grad_input0) != size(input) then
         _grad_input0 = self.grad_input[1] = zeros(D, F, size(input)...)
       end
       grad_input0::DenseMatrix(D,F) = _grad_input0

       values::vectorOf(D,F) = self.updates.values
       weight::matrixOf(D,F) = self.weight.values

       for frame = 1:self.updates.size
          for j = self.updates.offsets[frame]:(self.updates.offsets[frame+1] - 1)
             local current = self.updates.nodes[j]
             local ix = current - self.nleaves
             grad_input0[frame] += gradOutput[frame] * values[j] * weight[ix]
          end
       end

       return self.grad_input
    end
 end

function update_gradient!{D<:Device, F<:Float}(self::BinaryHierarchicalSoftmax{D, F}, inputs, gradOutput, scale::F=1.)
  input::matrixOf(D,F), targets = inputs

  @assert(ndims(input) == 2, "Input must be a matrix (ndims=$(ndims(input)))")

  values::vectorOf(D,F) = self.updates.values
  gradWeight::matrixOf(D,F) = self.weight.gradient
  gradBias::vectorOf(D,F) = self.bias.gradient

  for frame = 1:self.updates.size
    for j = self.updates.offsets[frame]:(self.updates.offsets[frame+1] - 1)
      local current = self.updates.nodes[j]
      local ix = current - self.nleaves
      local c = scale * values[j] * gradOutput[frame]

      gradWeight[ix] += c * input[frame]
      gradBias[ix] += c
    end
  end
 end


text:
cb = pyimport("palettable.colorbrewer.diverging")
cmgood = cb.RdYlBu_11_r.mpl_colormap
rdylbl = cb.RdYlBu_11_r.mpl_colormap

matlab_colors() = [[0, 0.4470, 0.7410, 1],
                   [0.8500, 0.3250, 0.0980, 1],
                   [0.9290, 0.6940, 0.1250, 1],
                   [0.4940, 0.1840, 0.5560, 1],
                   [0.4660, 0.6740, 0.1880, 1],
                   [0.3010, 0.7450, 0.9330, 1],
                   [0.6350, 0.0780, 0.1840, 1]
                   ]

SEISMIC_WIDE = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "seismic_wide",
    [
        (0, 0, 0.3, 1),
        (0, 0, 1, 1),
        (0.6, 0.6, 1, 1),
        (0.9, 0.9, 1, 1),
        (1, 1, 1, 1),
        (1, 0.9, 0.9, 1),
        (1, 0.6, 0.6, 1),
        (1, 0, 0, 1),
        (0.5, 0, 0, 1),
    ],  # Extra white in middle from seismic
    N = 250,
)
plt.register_cmap("seismic_wide", SEISMIC_WIDE)
SEISMIC_WIDE2 = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "seismic_wider",
    [
        (0, 0, 0.3, 1),
        (0, 0, 0.7, 1),
        (0.1, 0.1, 0.9, 1),
        (0.3, 0.3, 0.95, 1),
        (0.6, 0.6, 1, 1),
        (0.85, 0.85, 1, 1),
        (0.92, 0.92, 1, 0.99),
        (0.98, 0.98, 1, 0.98),
        (1, 1, 1, 0.95),
        (1, 0.98, 0.98, 0.98),
        (1, 0.92, 0.92, 0.99),
        (1, 0.85, 0.85, 1),
        (1, 0.6, 0.6, 1),
        (1, 0.3, 0.3, 1),
        (0.9, 0.1, 0.1, 1),
        (0.7, 0, 0, 1),
        (0.3, 0, 0, 1),
    ],
    N = 250,
)
plt.register_cmap("seismic_wider", SEISMIC_WIDE2)

SEISMIC_WIDE_Y = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "seismic_wide_y",
    [
        (0, 0, 0.7, 1),
        (0, 0, 1, 1),
        (0.1, 0.3, 1, 1),
        (0.2, 0.6, 0.95, 1),
        (0.5, 0.8, 0.85, 1),
        (0.8, 0.95, 0.80, 1),
        (0.9, 0.92, 0.78, 1),
        (0.95, 0.95, 0.75, 1),
        (0.95, 0.92, 0.70, 1),
        (0.95, 0.85, 0.65, 1),
        (0.97, 0.65, 0.3, 1),
        (1, 0.3, 0.2, 1),
        (1, 0.1, 0.05, 1),
        (1, 0, 0, 1),
        (0.7, 0, 0, 1),
    ],
    N = 450,
)
plt.register_cmap("seismic_wide_y", SEISMIC_WIDE_Y)

DSY5 = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "discrete_seismic_y5",
    broadcast(
        t -> t ./ 256,
        [
            (5, 113, 176, 256),
            (146, 197, 222, 256),
            # (247, 247, 247, 245),  # To make red-white-blue
            (247, 247, 191, 256),  # To make red-yellow-blue
            (244, 165, 130, 256),
            (202, 0, 32, 256),
        ],
    ),
    N = 5,
)
plt.register_cmap("discrete_seismic_y5", DSY5)

DSY7 = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "discrete_seismic_y7",
    broadcast(
        t -> t ./ 256,
        [
            (69, 117, 199, 256),
            (145, 191, 219, 256),
            (224, 243, 248, 256),
            (255, 255, 191, 256),
            (254, 224, 144, 256),
            (252, 141, 89, 256),
            (215, 48, 39, 256),
        ],
    ),
    N = 7,
)
plt.register_cmap("discrete_seismic_y7", DSY7)

DSY7 = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "discrete_seismic_y7",
    broadcast(
        t -> t ./ 256,
        [
            (69, 117, 199, 256),
            (145, 191, 219, 256),
            (224, 243, 248, 256),
            (255, 255, 191, 256),
            (254, 224, 144, 256),
            (252, 141, 89, 256),
            (215, 48, 39, 256),
        ],
    ),
    N = 7,
)
plt.register_cmap("discrete_seismic_y7", DSY7)

NSY = plt.matplotlib.colors.LinearSegmentedColormap.from_list(
    "seismic_narrow_y",
    broadcast(
        t -> t ./ 256,
        [
            (9, 37, 99, 256),
            (39, 97, 199, 256),
            (105, 121, 219, 256),
            (255, 255, 191, 255),
            (252, 111, 69, 256),
            (215, 38, 19, 256),
            (105, 18, 9, 256),
        ],
    ),
    N = 256,
)
plt.register_cmap("seismic_narrow_y", NSY)


text:
module TestExamples3DEulerGravity

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_3d_dgsem")

@testset "Compressible Euler with self-gravity" begin
  @trixi_testset "elixir_eulergravity_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
      l2   = [0.0004276779201667428, 0.00047204222332596204, 0.00047204222332608705, 0.0004720422233259819, 0.0010987026250960728],
      linf = [0.003496616916238704, 0.003764418290373106, 0.003764418290377103, 0.0037644182903766588, 0.008370424899251105],
      resid_tol = 1.0e-4, tspan = (0.0, 0.2))
  end
end

end # module


text:
module ReadIR
using JSON2
using MLStyle
using CamJulia
using CamJulia.CAM

export load_aeson, aeson_to_ir

function load_aeson(path)
    open(path) do f
        let s = read(f, String)
            JSON2.read(s, NamedTuple)
        end
    end
end

function aeson_to_ir(x::NamedTuple)
    f = @match x.tag begin
        "ComLet"    => com_let
        "ComLetrec" => com_letrec
        "ComIf"     => com_if
        "ComWhile"  => com_while
        "ComMutate" => com_mut
        "ComFun"    => com_fun
        "ComApp"    => com_app
        "ComVar"    => com_var
        "ComBlock"  => com_block
        "ComTuple"  => com_tuple
        "ComProj"   => com_proj
        "ComSymbol" => com_sym
        "ComInt"    => com_int
        "ComBigInt" => com_bigint
        "ComDouble" => com_double
        "ComCh"     => com_ch
        "ComBool"   => com_bool
        "ComStr"    => com_str
        "ComNil"    => nothing
        "ComInternal" => com_internal
        a => throw(a)
    end
    f === nothing ? NilConst() : f(x.contents)
end

@active Load(x) begin
    aeson_to_ir(x)
end

list_conv(xs) = map(aeson_to_ir, xs)

com_let = @λ [name, Load(value), Load(body)] ->
          Let(name, value, body)

com_letrec = @λ [
        [[fst, Load(snd)] for (fst, snd) in bindings],
        Load(body)
    ] -> LetRec(bindings, body)

com_if = @λ [Load(cond), Load(iftrue), Load(iffasle)] ->
            If(cond, iftrue, iffasle)

com_while = @λ [Load(cond), Load(body)] ->
               While(cond, body)

com_mut = @λ [varname, Load(value)] -> Mutate(varname, value)

com_fun = @λ [args, Load(body)] -> Fun(args, body)

com_app = @λ [Load(fn), [Load(arg) for arg in args]] -> App(fn, args)

com_var = Var

com_block  = @λ [Load(expr) for expr in blocks] -> Block(blocks)

com_tuple = @λ [Load(elt) for elt in elts] -> Join(elts)

com_proj = @λ [Load(major), Load(ith)] -> Proj(major, ith)

com_bigint = BigIntConst ∘ IdrisIntegerT

com_sym = SymConst

com_int = IntConst

com_str = StrConst

com_double = DoubleConst

com_ch =  @λ s -> ChConst(s[1])

com_bool = BoolConst

com_internal = Internal

end

max_stars_repo_name: thautwarm/idris-cam

max_issues_repo_name: thautwarm/idris-cam

max_forks_repo_name: thautwarm/idris-cam


text:
function form(f::Function, args...; attrs...) :: HTMLString
  normal_element(f, "form", [args...], attr(attrs...))
end

"""
$TYPEDSIGNATURES
"""
function form(children::Union{String,Vector{String}} = "", args...; attrs...) :: HTMLString
  normal_element(children, "form", [args...], attr(attrs...))
end

"""
$TYPEDSIGNATURES
"""
function attr(attrs...)
  attrs = Pair{Symbol,Any}[attrs...]

  for p in attrs
    p[1] == :enctype && return attrs
  end

  push!(attrs, :enctype => "multipart/form-data")
end


text:
export iauGst00a
"""
Greenwich apparent sidereal time (consistent with IAU 2000
resolutions).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   uta,utb    double    UT1 as a 2-part Julian Date (Notes 1,2)
   tta,ttb    double    TT as a 2-part Julian Date (Notes 1,2)

Returned (function value):
               double    Greenwich apparent sidereal time (radians)

Notes:

   1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

         Part A        Part B

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

   2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession-nutation.  If UT1 is
   used for both purposes, errors of order 100 microarcseconds
   result.

   3. This GAST is compatible with the IAU 2000 resolutions and must be
   used only in conjunction with other IAU 2000 compatible
   components such as precession-nutation.

   4. The result is returned in the range 0 to 2pi.

   5. The algorithm is from Capitaine et al. (2003) and IERS
   Conventions 2003.

Called:
   iauGmst00    Greenwich mean sidereal time, IAU 2000
   iauEe00a     equation of the equinoxes, IAU 2000A
   iauAnp       normalize angle into range 0 to 2pi

References:

   Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
   implement the IAU 2000 definition of UT1", Astronomy &
   Astrophysics, 406, 1135-1149 (2003)

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauGst00a(double uta, double utb, double tta, double ttb)

function iauGst00a(uta::Real, utb::Real, tta::Real, ttb::Real)
   return ccall((:iauGst00a, libsofa_c), Cdouble,
               (Cdouble, Cdouble, Cdouble, Cdouble),
               convert(Float64, uta), convert(Float64, utb),
               convert(Float64, tta), convert(Float64, ttb))
end


text:
export WithMemory
export MemoizingParser

"""
    WithMemory(x) <: AbstractString

String wrapper with memoization of next match states for parsers at indices.
Memoization is sometimes recommended as a way of improving the performance of parser combinators (like state machine optimization and compilation for regular languages).

!!! note
    A snappy performance gain could not be demonstrated so far, probably because the costs of state memory allocation for caching are often greater than recomputing a match.
    If you have a case where your performance benefits with this, let me know!
```
"""
struct WithMemory{S,M} <: StringWrapper
    x::S
    mem::M
    function WithMemory(x::S,mem::M) where {S,M}
        new{S,M}(x,mem)
    end
end
function WithMemory(x)
    WithMemory(x,Dict())
end
Base.show(io::IO, x::WithMemory) =
    print(io,x.x)

"""
    MemoizingParser{P,S,T}

[`WrappedParser`](@ref) memoizing all match states.
For slow parsers with a lot of backtracking this parser can help improve speed.

(Sharing a good example where memoization makes a difference is appreciated.)
"""
@auto_hash_equals struct MemoizingParser{P,S,T} <: WrappedParser{P,S,T}
    parser::P
    function MemoizingParser(p)
        new{typeof(p),state_type(p),result_type(p)}(p)
    end
end

_deepmap_parser(f::Function,mem::AbstractDict,x::MemoizingParser,a...;kw...) =
    MemoizingParser(deepmap_parser(f,mem,x.parser,a...;kw...))

@inline function iterate_state(parser::MemoizingParser, sequence, till, posi,after,state)
    error("for memoizing, wrap sequence in WithMemory. Todo: automize wrapping in root parser with optimize")
    iterate_state(parser.parser, sequence, till,posi,after,state)
end

@inline Base.@propagate_inbounds function iterate_state(parser::MemoizingParser, sequence::WithMemory, till, posi,after,state)
    get!(sequence.mem,(parser.parser,posi,state)) do
        copy(iterate_state(parser.parser, sequence,till,posi,after,state))
    end
end


text:
using Test, Base.CoreLogging
import Base.CoreLogging: Info

import InteractiveUtils
InteractiveUtils.versioninfo(verbose=true)

using GPUCompiler

using LLVM, LLVM.Interop

include("util.jl")

@testset "GPUCompiler" begin

GPUCompiler.reset_runtime()

GPUCompiler.enable_timings()

include("native.jl")
include("ptx.jl")
include("spirv.jl")
include("gcn.jl")
include("bpf.jl")

include("examples.jl")

haskey(ENV, "CI") && GPUCompiler.timings()

end


text:
ENV["PYTHON_JL_RUNTIME_PYTHON"] = Sys.which("python")
dir(x) = fieldnames(typeof(x))
singlify(x::Float64) = Float32(x)

using Pkg
using PyCall
using OrderedCollections
using LinearAlgebra
using SpikeSynchrony
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units
using Distributed
import DataStructures
using JLD
using Metaheuristics
using Plots
unicodeplots()

global ngt_spikes
global opt_vec
global extremum
global extremum_param
global ngt_spikes
global fitness

PARALLEL=false
THREADS=false
if PARALLEL
	if nprocs()==1
		addprocs(8)
	else
		@everywhere include("../utils.jl")
		@everywhere include("../current_search.jl")
		eval(macroexpand(quote @everywhere using SpikingNeuralNetworks end))
		eval(macroexpand(quote @everywhere using PyCall end))

	end
else
	include("../utils.jl")
end

###
global cell_type="ADEXP"
global vecp=false
###

(vmgtv,vmgtt,ngt_spikes,ground_spikes,julia_version_vm) = get_data()
global simdur = last(vmgtt)*1000

println("Ground Truth")
plot(vmgtt[:],vmgtv[:]) |> display

#lower,upper=get_izhi_ranges()

lower,upper=get_adexp_ranges()
ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

function loss(E,ngt_spikes,ground_spikes)
    spikes = get_spikes(E)
    spikes = [s/1000.0 for s in spikes]

	#opt_vec = [i[1] for i in opt_vec]
	#s_a = signal(opt_vec, length(opt_vec)/last(vmgtt))
	#s_b = signal(vmgtv, length(vmgtt)/last(vmgtt))

	maxt = findmax(sort!(unique(vcat(spikes,ground_spikes))))[1]
    if size(spikes)[1]>1
        t, S = SPIKE_distance_profile(spikes, ground_spikes;t0=0,tf = maxt)
        spkdistance = SpikeSynchrony.trapezoid_integral(t, S)/(t[end]-t[1]) # == SPIKE_distance(y1, y2)
    else
        spkdistance = 10.0
    end
	if length(spikes)>1

		custom_raster2(spikes,ground_spikes)
		custom_raster(spikes,ground_spikes)

		#savefig(crp, "aligned_VMs_adexp.png")
		#display(crp)
	end
	spkdistance*=spkdistance

    delta = abs(size(spikes)[1] - ngt_spikes)
    return spkdistance+delta

end

function loss(param)
    current = current_search(cell_type,param,ngt_spikes)
	if cell_type=="IZHI"
        pp = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(;N = 1, param = pp)
    end
    if cell_type=="ADEXP"
        adparam = SNN.ADEXParameter(;a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10])

        E = SNN.AD(;N = 1, param=adparam)
    end

    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms

    E.I = [current*nA]
	#@show(simdur)
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = simdur)
    vecplot(E, :v) |> display
    error = loss(E,ngt_spikes,ground_spikes)
	@show(error)

    error
end
#=
function Evolutionary.value!(::Val{:serial}, fitness, objfun, population::AbstractVector{IT}) where {IT}
	Threads.@threads for i in 1:length(population)
		fitness[i] = value(objfun, population[i])
	end
end
=#
ɛ = 0.125
options = GA(
    populationSize = 10,
    ɛ = ɛ,
    selection = ranklinear(1.5),#ss,
    crossover = intermediate(1.0),#line(1.0),#xovr,
    mutation = uniform(1.0),#domainrange(fill(1.0,ts)),#ms
)
@time result = Evolutionary.optimize(loss,lower,upper, initd, options,
    Evolutionary.Options(iterations=140, successive_f_tol=25, show_trace=true, store_trace=true)#,parallelisation=:thread)
)
import Plots
gr()

fitness = minimum(result)
println("GA: ɛ=$ɛ) => F: $(minimum(result))")# C: $(Evolutionary.iterations(result))")
extremum_param = Evolutionary.minimizer(result)

using SignalAnalysis
opt_vec,opt_spikes = checkmodel(extremum_param,cell_type,ngt_spikes)
trace = result.trace
filename = string("../JLD/PopulationScatter_adexp.jld")
save(filename,"trace",trace,"opt_vec", opt_vec,"extremum_param", extremum_param, "vmgtt", vmgtt, "vmgtv", vmgtv, "ground_spikes", ground_spikes,"opt_spikes",opt_spikes)
evo_loss = [t.value for t in trace]
loss_evolution = []
for (i,l) in enumerate(evo_loss)
	if i>1
		append!(loss_evolution,l)
	end
end
display(plot(loss_evolution))
println("Fittest model parameters")
println(extremum_param)

opt_vec = [i[1] for i in opt_vec]

s_a = signal(opt_vec, length(opt_vec)/last(vmgtt))
s_b = signal(vmgtv, length(vmgtt)/last(vmgtt))
p = plot(s_a)
p2 = plot!(p, s_b)
display(p2)
savefig(p2, "../imgs/aligned_VMs_adexp.png")

crp=custom_raster(opt_spikes,ground_spikes)
savefig(crp, "../imgs/aligned_VMs_adexp.png")
display(crp)

#display(plot(s_a))
#display(plot(s_b))

@show(result.minimizer)
@show(fitness)

#plot(opt_vec)|> display
#plot(vmgtt[:],vmgtv[:])|> display

#=
D = 10
bounds = [3ones(D) 40ones(D)]'

bounds = [lower upper]'
a = view(bounds, 1, 1)
b = view(bounds, 1, 2)

information = Information(f_optimum = 0.0)
options = Options( seed = 1, iterations=10, f_calls_limit =10)

D = size(bounds, 2)
nobjectives=1
methods = [
        SMS_EMOA(N = 20, n_samples=20, options=options),
        NSGA2(options=options),
        #MOEAD_DE(gen_ref_dirs(1, 1), options=Options( seed = 1, iterations = 5)),
        #NSGA3(options=options),
      ]

#for method in [methods[1]]
    #f_calls = 0
result = ( optimize(loss, bounds, NSGA2(options=options)) )
show(IOBuffer(), "text/html", result)
show(IOBuffer(), "text/plain", result.population)
show(IOBuffer(), "text/html", result.population)
show(IOBuffer(), result.population[1])
#end
#, label = "randData", ylabel = "Y axis",color = :red, legend = :topleft, grid = :off, xlabel = "Numbers Rand")
#p = twiny()
#plot!(p,vec, label = "log(x)", legend = :topright, box = :on, grid = :off, xlabel = "Log values") |> display
#return fitness
#end
#end
#break
#end
=#


text:
using Documenter, DemoCards
using MIRT

# examples_templates, examples_theme = cardtheme("grid")
# examples, examples_cb = makedemos("examples", examples_templates)

format = Documenter.HTML(
    edit_link = "master",
    prettyurls = get(ENV, "CI", nothing) == "true",
#   assets = [examples_theme],
)

makedocs(;
    modules = [MIRT],
    format,
    sitename = "MIRT.jl",
    pages = [
        "Home" => "index.md",
        "Table of Contents" => "toc.md",
        # "Examples" => examples,
        "Function References" => "reference.md",
    ],
    authors = "Jeff Fessler and contributors",
)

# examples_cb()

deploydocs(;
    repo = "github.com/JeffFessler/MIRT.jl.git",
    devbranch = "master",
    versions = ["stable" => "v^", "dev" => "dev"],
)


text:
function graphdiff(g1::Graph, g2::Graph)
    new_edges = Set{Pair{Int,Int}}()
    for u in nodes(g1)
        for v in neighbors(g1,u)
            if u<=v
                push!(new_edges, Pair(u,v))
            end
        end
    end
    for u in nodes(g2)
        for v in neighbors(g2,u)
            if u<=v
                delete!(new_edges, Pair(u,v))
            end
        end
    end
    new_edges
end

function graphdiff(graph1, graph2)
    edges_g1 = readdlm(graph1, Int)
    edges_g2 = readdlm(graph2, Int)
    new_edges = Set{Pair{Int,Int}}()
    for i=1:size(edges_g1, 1)
        u = edges_g1[i,1]
        v = edges_g1[i,2]
        if u>v
            u,v = v,u
        end
        push!(new_edges, Pair(u,v))
    end
    for i=1:size(edges_g2, 1)
        u = edges_g2[i,1]
        v = edges_g2[i,2]
        if u>v
            u,v = v,u
        end
        delete!(new_edges, Pair(u,v))
    end
    new_edges
end

function modularity(g, c)
    labels = unique(c)
    length(labels) == 1 && return 0.0
    m = ne(g)
    e = Dict{Int,Int}()
    a = Dict{Int,Int}()
    for u in nodes(g)
        for v in neighbors(g,u)
            if u<v
                c1 = c[u]
                c2 = c[v]
                if c1 == c2
                    e[c1] = get(e,c1,0) + 2
                end
                a[c1] = get(a,c1,0) + 1
                a[c2] = get(a,c2,0) + 1
            end
        end
    end
    modularity_value = 0.0
    if m > 0
        for i in labels
            if i > 0
                tmp = haskey(a,i) ? a[i]/2/m : 0
                modularity_value += haskey(e,i) ? e[i]/2/m : 0
                modularity_value -= tmp*tmp
            end
        end
    end
    modularity_value
end


text:
import TrajectoryOptimization: states
struct CopyModel{K,N,M,L} <: AbstractModel
    model::L
    ix::SVector{N,Int}
    iu::SVector{M,Int}
    xinds::Vector{SVector{N,Int}}
    uinds::Vector{SVector{M,Int}}
end

function CopyModel(model::L, K::Int) where L <: AbstractModel
    n,m = size(model)
    xind = SVector{n}(1:n)
    uind = SVector{m}(1:m)
    xinds = [xind .+ (i-1)*n for i = 1:K]
    uinds = [uind .+ (i-1)*m for i = 1:K]
    ix = SVector{n}(1:n)
    iu = SVector{m}(n .+ (1:m))
    return CopyModel{K,n,m,L}(model, ix, iu, xinds, uinds)
end

Base.size(::CopyModel{K,N,M}) where {K,N,M} = N*K,M*K

function build_state(model::CopyModel{K}, xs...) where K
    @assert length(xs) == K
    x = xs[1]
    for i = 2:K
        x = [x; xs[i]]
    end
    return x
end

function Base.zeros(model::CopyModel{K,N,M}) where {K,N,M}
    x,u = zeros(model.model)
    return repeat(x,K), repeat(u,K)
end

function Base.rand(model::CopyModel{K}) where K
    x,u = rand(model.model)
    for i = 2:K
        x_,u_ = rand(model.model)
        x = [x; x_]
        u = [u; u_]
    end
    return x,u
end

function states(model::CopyModel{K,N,M}, x) where {K,N,M}
    reshape(x,N,K)
    return [x[i] for i in model.xinds]
end

states(model::CopyModel{K,L}, z::KnotPoint, k::Int) where {L,K} =
    states(model, state(z), k)
function states(model::CopyModel{K,L}, x, k::Int) where {L,K}
    return x[model.xinds[k]]
end

function states(model::CopyModel, Z::Traj, k::Int)
    [state(z)[model.xinds[k]] for z in Z]
end

function controls(model::CopyModel{K,N,M}, u) where {K,N,M}
    reshape(u,M,K)
end

controls(model::CopyModel{K,L}, z::KnotPoint, k::Int) where {L,K} =
    controls(model, control(z), k)
function controls(model::CopyModel{K,L}, u, k::Int) where {L,K}
    return u[model.uinds[k]]
end

function controls(model::CopyModel, Z::Traj, k::Int)
    [control(z)[model.uinds[k]] for z in Z]
end

function TrajectoryOptimization.get_trajectory(model::CopyModel, Z::Traj, k::Int)
    xinds = model.xinds
    uinds = model.uinds
    map(Z) do z
        x = state(z)[xinds[k]]
        u = control(z)[uinds[k]]
        KnotPoint(x,u,z.dt,z.t)
    end
end

function TrajectoryOptimization.dynamics(model::CopyModel{K}, x, u, t=0.0) where K
    xinds = model.xinds
    uinds = model.uinds
    xdot = dynamics(model.model, x[xinds[1]], u[uinds[1]], t)
    for i = 2:K
        xdot_ = dynamics(model.model, x[xinds[i]], u[uinds[i]], t)
        xdot = [xdot; xdot_]
    end
    return xdot
end

function TrajectoryOptimization.discrete_dynamics(::Type{Q}, model::CopyModel{K},
        z::KnotPoint) where {K,Q <: TrajectoryOptimization.Implicit}
    xinds = model.xinds
    uinds = model.uinds
    x,u = state(z), control(z)
    xdot = discrete_dynamics(Q, model.model, x[xinds[1]], u[uinds[1]], z.t, z.dt)
    for i = 2:K
        xdot_ = discrete_dynamics(Q, model.model, x[xinds[i]], u[uinds[i]], z.t, z.dt)
        xdot = [xdot; xdot_]
    end
    return xdot
end

function TrajectoryOptimization.jacobian(model::CopyModel{K,N0,M0},
        z::KnotPoint{T,N,M,NM}) where {K,N0,M0,T,N,M,NM}
    A0 = @SMatrix zeros(N0,N0)
    B0 = @SMatrix zeros(N0,M0)
    ix = model.ix
    iu = model.iu

    xinds = model.xinds
    uinds = model.uinds
    x,u = state(z), control(z)

    # Process fist model
    z_ = [x[xinds[1]]; u[uinds[1]]]
    z_ = TrajectoryOptimization.StaticKnotPoint(z_, ix, iu, z.dt, z.t)
    ∇f = jacobian(model.model, z_)
    A,B = ∇f[ix,ix], ∇f[ix,iu]
    # append zeros after
    for i = 2:K
        A = [A  A0]
        B = [B  B0]
    end

    # loop over the rest of the models, appending below the first
    for i = 2:K
        z_ = [x[xinds[i]]; u[uinds[i]]]
        z_ = TrajectoryOptimization.StaticKnotPoint(z_, ix, iu, z.dt, z.t)
        ∇f = jacobian(model.model, z_)
        A_,B_ = ∇f[ix,ix], ∇f[ix,iu]

        # prepend zeros
        for j = 1:i-1
            A_ = [A0 A_]
            B_ = [B0 B_]
        end
        # append zeros
        for j = i+1:K
            A_ = [A_ A0]
            B_ = [B_ B0]
        end
        # stack with original
        A = [A; A_]
        B = [B; B_]
    end
    return [A B]
end

# function TrajectoryOptimization.discrete_jacobian(::Type{Q}, model::CopyModel{K,N0,M0},
#         z::KnotPoint{T,N,M,NM}) where {K,N0,M0,T,N,M,NM,Q<:TrajectoryOptimization.Implicit}
#     A0 = @SMatrix zeros(N0,N0)
#     B0 = @SMatrix zeros(N0,M0)
#     ix = model.ix
#     iu = model.iu
#     it = N0+M0+1
#
#     xinds = model.xinds
#     uinds = model.uinds
#     x,u = state(z), control(z)
#     dt = @SVector [z.dt]
#
#     z_ = [x[xinds[1]]; u[uinds[1]]; dt]
#     ∇f = discrete_jacobian(Q, model.model, z_, z.t, ix, iu)
#     A,B,C = ∇f[ix,ix], ∇f[ix,iu], ∇f[ix,it]
#     for i = 2:K
#         A = [A  A0]
#         B = [B  B0]
#     end
#     for i = 2:K
#         z_ = [x[xinds[i]]; u[uinds[i]]; dt]
#         ∇f = discrete_jacobian(Q, model.model, z_, z.t, ix, iu)
#         A_,B_,C_ = ∇f[ix,ix], ∇f[ix,iu], ∇f[ix,it]
#
#         for j = 1:i-1
#             A_ = [A0 A_]
#             B_ = [B0 B_]
#         end
#         for j = i+1:K
#             A_ = [A_ A0]
#             B_ = [B_ B0]
#         end
#         A = [A; A_]
#         B = [B; B_]
#         C = [C; C_]
#     end
#     return [A B C]
# end

function discrete_jacobian!(::Type{Q}, ∇f, model::Dynamics.CopyModel{K,N0,M0},
		z::KnotPoint{T,N,M,NM}) where {T,N,M,NM,K,N0,M0,Q<:TrajectoryOptimization.Implicit}
    xinds = model.xinds
    uinds = model.uinds
	ix = model.ix
	iu = model.iu
	iz = [ix; iu; @SVector [N0+M0+1]]
	x,u = state(z), control(z)
	for i = 1:K
		z_ = [x[xinds[i]]; u[uinds[i]]]
		z_ = StaticKnotPoint(z_, ix, iu, z.dt, z.t)
		∇f_ = discrete_jacobian(Q, model.model, z_)
		A,B,C = ∇f_[ix,ix], ∇f_[ix,iu], ∇f_[ix,N0+M0+1]

		ix_ = ix .+ (i-1)*N0
		iu_ = iu .+ (i-1)*M0 .+ N .- N0
		it_ = NM+1
		∇f[ix_,ix_] .= A
		∇f[ix_,iu_] .= B
		∇f[ix_,it_] .= C
	end
end

function state_diff_size(model::CopyModel{K}) where K
    return state_diff_size(model.model)*K
end

state_diff(model::CopyModel, x::SVector, x0::SVector) = x-x0
function state_diff(model::CopyModel{K,N,M,L}, x::SVector, x0::SVector) where {K,N,M,L<:RigidBody}
    xinds = model.xinds
    dx = state_diff(model.model, x[xinds[1]], x0[xinds[1]])
    for i = 2:K
        dx_ = state_diff(model.model, x[xinds[i]], x0[xinds[i]])
        dx = [dx; dx_]
    end
    return dx
end

TrajectoryOptimization.state_diff_jacobian!(G, model::CopyModel, Z::Traj) = nothing
function TrajectoryOptimization.state_diff_jacobian!(G,
        model::CopyModel{<:Any,<:Any,<:Any,L}, Z::Traj) where {L<:RigidBody{R}} where R
    for k in eachindex(Z)
        G[k] .= state_diff_jacobian(model, state(Z[k]))
    end
end

state_diff_jacobian(::CopyModel, x::SVector) = I
@generated function state_diff_jacobian(model::CopyModel{K,N,M,L},
        x::SVector) where {K,N,M,L<:RigidBody{R}} where R
    if R <: UnitQuaternion
        if R <: UnitQuaternion{T,IdentityMap} where T
            return :(I)
        else
            G0 = @SMatrix zeros(N,N-1)
        end
    else
        G0 = @SMatrix zeros(N,N)
    end
    quote
        xinds = model.xinds
        G = state_diff_jacobian(model.model, x[xinds[1]])
        for i = 2:K
            G = [G $G0]
        end
        for i = 2:K
            G_ = state_diff_jacobian(model.model, x[xinds[i]])
            for j = 1:i-1
                G_ = [$G0 G_]
            end
            for j = i+1:K
                G_ = [G_ $G0]
            end
            G = [G; G_]
        end
        G
    end
end

@generated function TrajectoryOptimization.∇²differential(model::CopyModel{K,N,M,L},
        x::SVector, dx::SVector) where {K,N,M,L<:RigidBody{R}} where R

    if R <: UnitQuaternion
        if R <: UnitQuaternion{T,IdentityMap} where T
            return :(I*0)
        else
            G0 = @SMatrix zeros(N-1,N-1)
        end
    else
        G0 = @SMatrix zeros(N,N)
    end
    quote
        dx_ = reshape(dx,:,K)
        ix_ = SVector{12}(1:12)
        xinds = model.xinds
        G = TrajectoryOptimization.∇²differential(model.model, x[xinds[1]], dx[ix_])
        for i = 2:K
            G = [G $G0]
        end
        for i = 2:K
            dx_ = dx[ix_ .+ (i-1)*12]
            G_ = TrajectoryOptimization.∇²differential(model.model, x[xinds[i]], dx_)
            for j = 1:i-1
                G_ = [$G0 G_]
            end
            for j = i+1:K
                G_ = [G_ $G0]
            end
            G = [G; G_]
        end
        G
    end
end

function ∇²differential!(G, model::Dynamics.CopyModel{K,N,M,L},
        x::SVector, dx::Vector) where {K,N,M,L<:RigidBody{R}} where R
    ix = SVector{12}(1:12)
    xinds = model.xinds
    for i = 1:K
        ix_ = ix .+ (i-1)*12
        dx_ = dx[ix_]
        G_ = ∇²differential(model.model, x[xinds[i]], dx_)
        G[ix_, ix_] .= G_
    end
end


text:
type Wssmp                              # WSMP symmetric type
    ia::Vector{Cint}
    ja::Vector{Cint}
    avals::Vector{Cdouble}
    diag::Vector{Cdouble}
    perm::Vector{Cint}
    invp::Vector{Cint}
    mrp::Vector{Cint}
    iparm::Vector{Cint}
    dparm::Vector{Cdouble}
end

function Wssmp(ia::Vector{Cint},ja::Vector{Cint},avals::Vector{Cdouble},
                 dd::Vector{Cdouble}=Cdouble[])
    1 == ia[1] || error("Column pointers must be 1-based.  That is, ia[1] == 1")
    (n = length(ia) - 1) > 0 || error("0 × 0 sparse matrices not allowed")
    (ld = length(dd)) == 0 || ld == n ||
        error("length(diag) = $length(diag), should be 0 or $n")
    length(ja) == length(avals) == ia[end] - 1 || throw(DimensionMismatch(""))
    for j in ja
        0 < j ≤ n || error("elements of ja must be in [1,$n]")
    end
    sorted = true
    for i in 1:n         # check for sorted columns and lower triangle
        ia[i] ≤ ia[i+1] || error("elements in ia must be non-decreasing")
        for k in ia[i]:(ia[i+1]-1)
            (ld == 0 ? i ≤ ja[k] : i < ja[k]) || error("matrix is not lower triangular")
            k == ia[i] || ja[k-1] ≤ ja[k] || (sorted = false)
        end
    end
    if !sorted                          # sort columns if necessary
        info = zeros(Cint,1)
        ccall((:ws_sortindices_d__,libwsmp),Void,
              (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint}),
              &n,&n,ia,ja,avals,info)
        info[1] == 0 || error("error code $(info[1]) from ws_sortindices_d__")
    end
    ans = Wssmp(ia,ja,avals,dd,Array(Cint,(n,)),Array(Cint,(n,)),Cint[],
                  zeros(Cint,64),zeros(64))
    wssmp(ans,0)
    ld > 0 && (ans.iparm[4] = 1)
    ans
end

function Wssmp(S::SparseMatrixCSC{Cdouble,Cint}, MSC::Bool=false)
    issym(S) || error("Matrix S is not symmetric")
    A = MSC ? tril(S,-1) : tril(S)
    Wssmp(A.colptr,A.rowval,A.nzval,MSC ? diag(S) : Cdouble[])
end

function wssmp(a::Wssmp, task::Integer)
    0 ≤ task ≤ 3 || error("task = $task should be in the range [0,3]")
    a.iparm[2] ≤ task || error("wssmp called with task = $task but a.iparm[2] = $(a.iparm[2])")
    a.iparm[3] = task
    task ≡ 0 && (a.iparm[1] = 0)        # basic initialization
    task ≡ 3 && a.iparm[32] ≠ 0 && length(a.diag) ≡ 0 && (a.diag = zeros(length(a.perm)))
    ccall((:wssmp_,libwsmp),Void,
          (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},
           Ptr{Cint},Ptr{Cint},Ptr{Void},Ptr{Cint},Ptr{Cint},Ptr{Cint}, Ptr{Cdouble}),
          &(length(a.ia)-1),a.ia,a.ja,a.avals,a.diag,a.perm,a.invp,C_NULL,&1,&0,C_NULL,&0,
          a.mrp,a.iparm,a.dparm)
    if a.iparm[64] ≠ 0
        error("error code $(a.iparm[64]) from ccall to wssmp_ with a.iparm[3] == $(a.iparm[3])")
    end
    a
end

function wssmp(a::Wssmp, b::StridedVecOrMat{Cdouble}, task::Integer=4)
    4 ≤ task ≤ 5 || error("task = $task should be 4 or 5")
    a.iparm[3] = task
    ccall((:wssmp_,libwsmp),Void,
          (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},
           Ptr{Cint},Ptr{Cint},Ptr{Void},Ptr{Cint},Ptr{Cint},Ptr{Cint}, Ptr{Cdouble}),
          &(length(a.ia)-1),a.ia,a.ja,a.avals,a.diag,a.perm,a.invp,
          b,&stride(b,2),&size(b,2),C_NULL,&0,a.mrp,a.iparm,a.dparm)
    if a.iparm[64] ≠ 0
        error("error code $(a.iparm[64]) from ccall to wssmp_ with a.iparm[3] == $(a.iparm[3])")
    end
    b
end

Base.A_ldiv_B!(A::Wssmp, b::StridedVecOrMat{Cdouble}) = (A.iparm[2] = min(A.iparm[2],4);wssmp(A,b,4))
(\)(A::Wssmp, b::StridedVecOrMat{Cdouble}) = (A.iparm[2] = min(A.iparm[2],4);wssmp(A,copy(b),4))

Base.size(A::Wssmp) = (n = length(A.ia) - 1; (n,n))

function Base.size(A::Wssmp,i::Integer)
    i ≤ 0 && throw(BoundsError())
    n = length(A.ia) - 1
    1 ≤ i ≤ 2 && return n
    return 1
end

function *(A::Wssmp,x::Vector{Cdouble})
    A.iparm[4] == 0 || error("Wssmp object A must be in CSC format, not MSC")
    wssmatvec(A.ia,A.ja,A.avals,x,similar(x))
end

## for (nm,sym,op5) in ((:wmmrb,:wmmrb_,3),
##                      (:wkktord,:wkktord_,0))
##     @eval begin
##         function $nm(A::SparseMatrixCSC{Cdouble,Cint})
##             issym(A) || error("Matrix A must be symmetric")
##             adj = Base.SparseMatrix.fkeep!(copy(A), (i,j,x,other)->(i!=j), None)
##             n = size(adj,1)
##             options = zeros(Cint,5)
##             options[1] = 3
##             options[3] = 1
##             options[5] = $op5
##             perm = Array(Cint,n)
##             invp = Array(Cint,n)
##             ccall(($(string(nm,"_")),libwsmp),Void,(Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},
##                                                  Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
##                   &n,adj.colptr,adj.rowval,options,&1,perm,invp,C_NULL,&0)
##             perm,invp
##         end
##     end
## end

function wsmmrb(A::SparseMatrixCSC{Cdouble,Cint},options::Vector{Cint}=Int32[3,0,1,0,3])
    issym(A) || error("Matrix A must be symmetric")
    5 ≤ length(options) || error("options vector must be of length 5 or more")
    adj = Base.SparseMatrix.fkeep!(copy(A), (i,j,x,other)->(i ≠ j), None)
    n = size(adj,1)
    perm = Array(Cint,n)
    invp = Array(Cint,n)
    ccall((:wmmrb_,libwsmp),Void,(Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},
                                  Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
          &n,adj.colptr,adj.rowval,options,&1,perm,invp,C_NULL,&0)
    perm,invp
end

function wskktord(A::SparseMatrixCSC{Cdouble,Cint},options::Vector{Cint}=Int32[3,0,1,0,0])
    issym(A) || error("Matrix A must be symmetric")
    5 ≤ length(options) || error("options vector must be of length 5 or more")
    adj = Base.SparseMatrix.fkeep!(copy(A), (i,j,x,other)->(i ≠ j), None)
    n = size(adj,1)
    perm = Array(Cint,n)
    invp = Array(Cint,n)
    ccall((:wkktord_,libwsmp),Void,(Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},
                                  Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
          &n,adj.colptr,adj.rowval,options,&1,perm,invp,C_NULL,&0)
    perm,invp
end


text:
struct Inserter <: Rewriter
    fs
    positions
end

function Inserter(pairs::AbstractVector)
    Inserter(getindex.(pairs, 1), getindex.(pairs, 2))
end

Base.show(io::IO, inserter::Inserter) = print(io, "Inserter($(inserter.positions))")

function rewritedoc(inserter::Inserter, p::AbstractPath, doc::XNode)
    for (f, pos) in zip(inserter.fs, inserter.positions)
        if isnothing(selectfirst(doc, pos.selector))
            error("Could not find position $pos to insert into.")
        end
        doc = insertfirst(doc, f(p, doc), pos)
    end
    return doc
end

function toccreator(;
    docsel = SelectTag(:article),
    hierarchysels = SelectTag.((:h2, :h3, :h4)))
    return function createtoc(p, doc)
        content = selectfirst(doc, docsel)
        if isnothing(content)
            error("Could not find contents at $docsel to create a table of contents from.")
        end
        return maketoc(content, hierarchysels)
    end
end


text:
"""
    Antenna

Abstract supertype for all antenna-like types.
"""
abstract type Antenna end

"""
    AbstractDipole <: Antenna

Abstract type for dipole antennas.

Subtypes of `AbstractDipole` have an orientation `azimuth_angle` (``ϕ``) and
`inclination_angle` (``γ``) where

| [rad] |     ϕ     |     γ      |
|:-----:|:---------:|:----------:|
|   0   | end fire  | vertical   |
|  π/2  | broadside | horizontal |

!!! note

    Angles describe the orientation of the antenna, not the radiation pattern.
"""
abstract type AbstractDipole <: Antenna end

"""
    azimuth(d::AbstractDipole)

Return the azimuth angle ``ϕ`` of `d` measured towards positive `y` direction from `x`.

See also: [`AbstractDipole`](@ref)
"""
azimuth

"""
    inclination(d::AbstractDipole)

Return the inclination (dip) angle ``γ`` of `d` measured from vertical.

See also: [`AbstractDipole`](@ref)
"""
inclination

"""
    Dipole

Dipole antenna with arbitrary orientation described by `azimuth_angle` and
`inclination_angle` from vertical.

See also: [`AbstractDipole`](@ref)
"""
struct Dipole <: AbstractDipole
    azimuth_angle::Float64
    inclination_angle::Float64
end
azimuth(d::Dipole) = d.azimuth_angle
inclination(d::Dipole) = d.inclination_angle

"""
    VerticalDipole

Dipole antenna with inclination angle ``γ = 0`` from the vertical.
"""
struct VerticalDipole <: AbstractDipole end
azimuth(d::VerticalDipole) = 0.0
inclination(d::VerticalDipole) = 0.0

"""
    HorizontalDipole

Dipole antenna with inclination angle ``γ = π/2`` from the vertical and `azimuth_angle`
orientation ``ϕ`` where:

|  ϕ [rad] | orientation |
|:--------:|:-----------:|
|     0    |  end fire   |
|    π/2   |  broadside  |

"""
struct HorizontalDipole <: AbstractDipole
    azimuth_angle::Float64  # radians
end
azimuth(d::HorizontalDipole) = d.azimuth_angle
inclination(d::HorizontalDipole) = π/2


text:
module forwardmodel
import InfiniteOpt: value, InfiniteModel, supports
using Interpolations: Interpolations
import Parameters: @with_kw
import DataFrames: DataFrame

import jcontrol: Track, upsample, int, ξ, closest_point_idx
import ..control: KinematicsProblem, DynamicsProblem, MTMproblem
import ..bicycle: State, Bicycle

export Solution, run_forward_model, solution2state

int(x) = (Int64 ∘ round)(x)

@with_kw struct Solution
    δt::Float64
    t::Vector{Float64} = Vector{Float64}[]    # time
    s::Vector{Float64} = Vector{Float64}[]    # track progression
    x::Vector{Float64} = Vector{Float64}[]    # position
    y::Vector{Float64} = Vector{Float64}[]
    θ::Vector{Float64} = Vector{Float64}[]    # orientation
    δ::Vector{Float64} = Vector{Float64}[]    # steering angle
    u::Vector{Float64} = Vector{Float64}[]    # velocity  | fpr DynamicsProblem its the longitudinal velocity component
    ω::Vector{Float64} = Vector{Float64}[]    # angular velocity

    n::Vector{Float64} = Vector{Float64}[]    # track errors
    ψ::Vector{Float64} = Vector{Float64}[]

    # Kinematic problem only
    β::Vector{Float64} = Vector{Float64}[]

    # Dynamics problem only
    v::Vector{Float64} = Vector{Float64}[]   # lateral velocity
    Fu::Vector{Float64} = Vector{Float64}[]

    # controls
    u̇::Vector{Float64} = Vector{Float64}[]
    δ̇::Vector{Float64} = Vector{Float64}[]
end

function Solution(df::DataFrame)::Solution
    _keys = filter!(k -> Symbol(k) != :δt, names(df))
    dt = "δt" ∈ names(df) ? df[1, :δt] : 0.0
    return Solution(; Dict(map(k -> Symbol(k)=>df[:, k], _keys))..., δt=dt)
end

"""
To get the bike's trajectory in time and wrt the allocentric
track's reference frame (instead of curvilinear coordinates),
this function takes the error `n` and `ψ` and computes the bike
position wrt to the track's centerline.

Since `n` and `ψ` are defined wrt to s, first the bike's position
is computed at increments of s, then this is upsampled and converted
to the corresponding time increments.

Finally, the bike's velocity and angular velocities are computed.
"""
function run_forward_model(
    problemtype::MTMproblem, track::Track, model::InfiniteModel; δt=0.005
)
    # get model's data and upsample
    n = ξ(value(model[:n]))
    ψ = ξ(value(model[:ψ]))
    s = ξ(value(model[:s]))

    # compute bike position with variable Δt (wrt s)
    svalues = 1:0.01:length(value(model[:n]))
    II() = zeros(Float64, length(svalues))
    Xs, Ys, θs = II(), II(), II()
    for (ni, i) in enumerate(svalues)
        _n, _ψ = n[i], ψ[i]

        # get closest track waypoint
        closest = argmin(abs.(track.S .- s[i]))

        # get position and orientation of track
        xp = track.X[closest]
        yp = track.Y[closest]
        θp = track.θ[closest]

        # get position of bike
        if _n == 0
            Xs[ni] = xp
            Ys[ni] = yp
        else
            norm_angle = θp + π / 2

            Xs[ni] = xp + _n * cos(norm_angle)
            Ys[ni] = yp + _n * sin(norm_angle)
        end

        # get orientation of bike
        θs[ni] = θp + _ψ
    end

    # get other variables
    δs = map(ξ(value(model[:δ])), svalues)
    δ̇s = map(ξ(value(model[:δ̇])), svalues)
    us = map(ξ(value(model[:u])), svalues)
    ωs = map(ξ(value(model[:ω])), svalues)
    ns = map(ξ(value(model[:n])), svalues)
    ψs = map(ξ(value(model[:ψ])), svalues)

    if problemtype isa KinematicsProblem
        βs = map(ξ(value(model[:β])), svalues)
        u̇s = map(ξ(value(model[:u̇])), svalues)
    else
        vs = map(ξ(value(model[:v])), svalues)
        βs = map(
                ξ(atan.(vs ./ (us .+ eps()))), svalues
            )
        Fus = map(ξ(value(model[:Fu])), svalues)
    end

    # get values at regular Δt
    Ts = map(ξ(value(model[:t])), svalues)
    time = Ts[1]:δt:Ts[end]
    I() = zeros(Float64, length(time))
    T, X, Y, θ, δ, δ̇, u̇ = I(), I(), I(), I(), I(), I(), I()
    u, ω, n, ψ, β, v, Fu = I(), I(), I(), I(), I(), I(), I()
    for (i, t) in enumerate(time)
        idx = findfirst(Ts .>= t)
        idx = isnothing(idx) ? 1 : idx

        T[i] = Ts[idx]
        X[i] = Xs[idx]
        Y[i] = Ys[idx]
        θ[i] = θs[idx]
        δ[i] = δs[idx]
        δ̇[i] = δ̇s[idx]
        u[i] = us[idx]
        ω[i] = ωs[idx]
        n[i] = ns[idx]
        ψ[i] = ψs[idx]
        β[i] = βs[idx]

        if problemtype isa KinematicsProblem
            u̇[i] = u̇s[idx]
        else
        v[i] = vs[idx]
            Fu[i] = Fus[idx]
        end
    end

    # get `s`
    S::Vector{Float64} = []
    for i in 1:length(X)
        # get the closest track point
        idx = closest_point_idx(track.X, X[i], track.Y, Y[i])
        push!(S, track.S[idx])
    end

    return Solution(
        δt=δt,
        t = T,
        s = S,
        x = X,
        y = Y,
        θ = θ,
        u = u,
        ω = ω,
        δ = δ,
        δ̇ = δ̇,
        u̇ = u̇,
        n = n,
        ψ = ψ,
        β = β,
        v = v,
        Fu=Fu,
    )
end

function trimsolution(sol::Solution, s0, s1)
    start = findfirst(sol.s .>= s0)
    isnothing(start) && return nothing

    stop = findlast(sol.s .<= s1)
    isnothing(start) && return nothing

    return Solution(
        δt= sol.δt,
        t = sol.t[start:stop],
        s = sol.s[start:stop],
        x = sol.x[start:stop],
        y = sol.y[start:stop],
        θ = sol.θ[start:stop],
        u = sol.u[start:stop],
        ω = sol.ω[start:stop],
        δ = sol.δ[start:stop],
        δ̇ =  sol.δ̇[start:stop],
        # u̇ =  sol.u̇[start:stop],
        n = sol.n[start:stop],
        ψ = sol.ψ[start:stop],
        β = sol.β[start:stop],
        v = sol.v[start:stop],
        Fu= sol.Fu[start:stop],
    )
end

function Base.lastindex(solution::Solution)::Int
    return length(solution.x)
end

function Base.getindex(solution::Solution, state)
    vars =  (:x, :y, :θ, :δ, :δ̇, :ω, :u, :β, :v, :n, :ψ, :Fu)
    return State(;
        Dict(
            map(v -> v=>getfield(solution, v)[state], vars)
        )...
    )
end

"""
Take the value of a forward model solution at a given svalue or timestep
and return a State with the solution values there/then.
"""
function solution2state(svalue::Number, solution::Solution; at=:place)::State
    # get solution sample IDX
    if at == :place
        svalue = svalue < 1 ? svalue * 258 : svalue
        idx = findfirst(solution.s .>= svalue)
        idx = isnothing(idx) ? length(solution.x) : idx
    elseif at == :time
        idx = findfirst(solution.t .> svalue)
        idx = isnothing(idx) ? 1 : idx
    end

    vars =  (:x, :y, :θ, :δ, :δ̇, :ω, :u, :β, :v, :n, :ψ, :Fu)
    state = State(;
        Dict(
            map(v -> v=>getfield(solution, v)[idx], vars)
        )...
    )
    # state.ψ = mod(2π - state.ψ, 2π)
    return state
end

"""
The value of `s` of a solution at a timestep.
"""
function solution2s(t::Number, solution::Solution)::Float64
    # get solution sample IDX
    idx = findfirst(solution.t .> t)
    idx = isnothing(idx) ? length(solution.t) : idx  # if solution is too short

    return solution.s[idx]
end
end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
const _PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))
include_dependency(joinpath(_PKG_ROOT_DIR, "Project.toml"))

const COBREXA_VERSION =
    VersionNumber(Pkg.TOML.parsefile(joinpath(_PKG_ROOT_DIR, "Project.toml"))["version"])

function _print_banner()
    c = Base.text_colors
    n = c[:normal]
    y = c[:bold] * c[:cyan]

    println(
        "
                   $(y)//   $(n) |
        \\\\\\\\\\  // $(y)//    $(n) | $(c[:bold])COBREXA.jl $(c[:normal])
         \\\\ \\\\// $(y)//     $(n) |
          \\\\ \\/ $(y)//      $(n) | $(c[:bold])CO$(c[:normal])nstraint-$(c[:bold])B$(c[:normal])ased $(c[:bold])R$(c[:normal])econstruction
           \\\\  $(y)//       $(n) | and $(c[:bold])EX$(c[:normal])ascale $(c[:bold])A$(c[:normal])nalysis in Julia
           //  $(y)\\\\       $(n) |
          // $(y)/\\ \\\\      $(n) | Documentation: https://lcsb-biocore.github.io/COBREXA.jl
         // $(y)//\\\\ \\\\      $(n)|
        // $(y)//  \\\\\\\\\\     $(n)| v$(COBREXA_VERSION) - $(Dates.year(Dates.today()))
       //                $(n)|
        ",
    )
end


text:
abstract  RandomProcess

doc"""
  SampleInterval(tmin, tmax, npts)

Time interval for evaluluation of randowm process where tmin is the start time, tmax the end time
and npts the number of points in the sample.
"""
immutable SampleInterval{T<:Integer, U<:Real}
  npts::T
  tmax::U
  tmin::U
  SampleInterval(npts::T, tmax::U, tmin::U) = tmax > tmin && npts > U(1) ? new(npts, tmax, tmin) : error("must have tmax > tmin && npts > 0")
end

# outer constructors
SampleInterval{T<:Integer, U<:Real}(npts::T, tmax::U, tmin::U) = SampleInterval{T, U}(npts, tmax, tmin)
SampleInterval(npts::Integer, tmax::Integer, tmin::Integer) = SampleInterval(UInt64(npts), Float64(tmax), Float64(tmin))
SampleInterval(npts::Real, tmax::Real, tmin::Real) = SampleInterval(UInt64(npts), promote(tmax, tmin)...)
SampleInterval(npts::Integer, tmax::Real, tmin::Real) = SampleInterval(npts, promote(tmax, tmin)...)
SampleInterval(npts::Integer, tmax::Real) = SampleInterval(npts, tmax, zero(tmax))

# conversions
convert{T<:Integer, U<:Real, V<:Integer, W<:Real}(::Type{SampleInterval{T, U}}, npts::V, tmax::W, tmin::W) = SampleInterval(T(npts), U(tmax), U(tmin))
convert{T<:Integer, U<:Real, V<:Integer, W<:Real}(::Type{SampleInterval{T, U}}, s::SampleInterval{V, W}) = SampleInterval(T(s.npts), U(s.tmax), U(s.tmin))

# parameters
params(sampleInterval::SampleInterval) = (sampleInterval.npts, sampleInterval.tmax, sampleInterval.tmin)
Δt(sampleInterval::SampleInterval) =  (sampleInterval.tmax - sampleInterval.tmin) / sampleInterval.npts

function samples(sampleInterval::SampleInterval)
    collect(linspace(sampleInterval.tmin, sampleInterval.tmax, sampleInterval.npts))
end

doc"""
BrownianMotion(μ, σ)

Brownian motion with mean μ and standard deviation σ.

https://en.wikipedia.org/wiki/Brownian_motion
"""
immutable BrownianMotion{T<:Real} <: RandomProcess
  μ::T
  σ::T
  BrownianMotion(μ, σ) = σ > zero(σ) ? new(μ, σ) : error("must have σ > 0")
end

# outer constructors
BrownianMotion{T<:Real}(μ::T, σ::T) = BrownianMotion{T}(μ, σ)
BrownianMotion(μ::Integer, σ::Integer) = BrownianMotion(Float64(μ), Float64(σ))
BrownianMotion(μ::Real, σ::Real) = BrownianMotion(promote(μ, σ)...)
BrownianMotion(μ::Real) = BrownianMotion(μ, 1.0)
BrownianMotion() = BrownianMotion(0.0, 1.0)

# conversions
convert{T<:Real, U<:Real}(::Type{BrownianMotion{T}}, μ::U, σ::U) = BrownianMotion(T(μ), T(σ))
convert{T<:Real, U<:Real}(::Type{BrownianMotion{T}}, bm::BrownianMotion{U}) = BrownianMotion(T(bm.μ), T(bm.σ))

# parameters
params(randomProcess::BrownianMotion) = (randomProcess.μ, randomProcess.σ)

# generation
function increments(randomProcess::BrownianMotion, interval::SampleInterval)
  μ = Δt(interval) * randomProcess.μ
  σ = sqrt(Δt(interval)) * randomProcess.σ
  bm = zeros(Float64, interval.npts-1)
  for n = 1:interval.npts-1
    bm[n] = rand(Normal(μ, σ))
  end
  return bm
end

function rand(randomProcess::BrownianMotion, interval::SampleInterval)
  bm = zeros(Float64, interval.npts)
  Δbm = increments(randomProcess, interval)
  for n = 2:interval.npts
    bm[n] = bm[n-1] + Δbm[n-1]
  end
  return bm
end

doc"""
GeometricBrownianMotion(μ, σ)

geometric Brownian motion with mean μ and standard deviation σ.

https://en.wikipedia.org/wiki/Geometric_Brownian_motion
"""

immutable GeometricBrownianMotion{T<:Real} <: RandomProcess
  brownianMotion::BrownianMotion{T}
  s0::T
  GeometricBrownianMotion(brownianMotion, s0) = s0 != zero(s0) ? new{T}(brownianMotion, s0) : error("must have s0 ! 0")
end

# outer constructors
GeometricBrownianMotion{T<:Real}(μ::T, σ::T, s0::T) = GeometricBrownianMotion{T}(BrownianMotion(μ, σ), s0)
GeometricBrownianMotion(μ::Integer, σ::Integer, s0::Integer) = GeometricBrownianMotion(Float64(μ), Float64(σ), Float64(s0))
GeometricBrownianMotion(μ::Real, σ::Real, s0::Real) = GeometricBrownianMotion(promote(μ, σ, s0)...)
GeometricBrownianMotion(μ::Real, σ::Real) = GeometricBrownianMotion(Float64(μ), Float64(σ), 1.0)
GeometricBrownianMotion(μ::Real) = GeometricBrownianMotion(Float64(μ), 1.0, 1.0)
GeometricBrownianMotion() = GeometricBrownianMotion(0.0, 1.0, 1.0)

# conversions
convert{T<:Real, U<:Real}(::Type{GeometricBrownianMotion{T}}, μ::U, σ::U, s0::T) = GeometricBrownianMotion(T(μ), T(σ, t(s0)))
convert{T<:Real, U<:Real}(::Type{GeometricBrownianMotion{T}}, gbm::GeometricBrownianMotion{U}) = GeometricBrownianMotion(T(gbm.μ), T(gbm.σ), T(gbm.s0))

# parameters
params(randomProcess::GeometricBrownianMotion) = (randomProcess.brownianMotion.μ, randomProcess.brownianMotion.σ, randomProcess.s0)

# generation
function increments(randomProcess::GeometricBrownianMotion, interval::SampleInterval)
  Δbm = increments(randomProcess.brownianMotion, interval)
  map(x -> exp(x), Δbm)
end

function rand(randomProcess::GeometricBrownianMotion, interval::SampleInterval)
  gbm = ones(Float64, interval.npts)
  Δgbm = increments(randomProcess, interval)
  gbm[1] = randomProcess.s0
  for n = 1:interval.npts-1
      gbm[n+1] = gbm[n] * Δgbm[n]
  end
  return gbm
end


text:
abstract NodeFunc

type Relu <: NodeFunc ; end
nfunc(a, ::Relu)             = max(0., a)
dnfunc(ds, a, s, ::Relu)     = a < 0. ? 0. : ds

type SoftRelu <: NodeFunc ; end
nfunc(a, ::SoftRelu)         = log(1. + exp(a))
dnfunc(ds, a, s, ::SoftRelu) = exp(a) / exp(s) * ds

type Sigmoid <: NodeFunc ; end
nfunc(a, ::Sigmoid)          = 1. / (1. + exp(-a))
dnfunc(ds, a, s, ::Sigmoid)  = s * (1 - s) * ds

type Tanh <: NodeFunc ; end
nfunc(a, ::Tanh)             = (1. - exp(-2a)) / (1. + exp(-2a))
dnfunc(ds, a, s, ::Tanh)     = (1 - s*s) * ds

# exp(x) / (1+exp(x))

function nfunc!(s::AbstractArray, a::AbstractArray, f::NodeFunc)
    for i in 1:length(s)
        s[i] = nfunc(a[i], f)
    end
end

function dnfunc!(ds::AbstractArray, a::AbstractArray, s::AbstractArray, f::NodeFunc)
    for i in 1:length(s)
        ds[i] = dnfunc(ds[i], a[i], s[i], f)
    end
end

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
using Test
using DFTK
using IntervalArithmetic

include("testcases.jl")

function discretized_hamiltonian(T, testcase)
    Ecut = 10  # Hartree

    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    atoms = [spec => testcase.positions]
    # disable symmetry for interval
    model = model_DFT(Array{T}(testcase.lattice), atoms, [:lda_x, :lda_c_vwn], symmetry=:off)

    # For interval arithmetic to give useful numbers,
    # the fft_size should be a power of 2
    fft_size = nextpow.(2, determine_grid_size(model.lattice, Ecut))
    basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1), fft_size=fft_size)

    ham = Hamiltonian(basis; ρ=guess_density(basis))
end

@testset "Application of an LDA Hamiltonian with Intervals" begin
    T = Float64
    ham = discretized_hamiltonian(T, silicon)
    hamInt = discretized_hamiltonian(Interval{T}, silicon)

    hamk = ham.blocks[1]
    hamIntk = hamInt.blocks[1]

    x = randn(Complex{T}, length(G_vectors(ham.basis.kpoints[1])))
    ref = hamk * x
    res = hamIntk * Interval.(x)

    # Difference between interval arithmetic and normal application less than 1e-10
    @test maximum(mid, abs.(res .- ref)) < 1e-10

    # Maximal error done by interval arithmetic less than
    @test maximum(radius, abs.(res)) < 1e-10
end


text:
export neo_hookean, compute_absorbing_boundary_condition_matrix

@doc raw"""
    neo_hookean(u::Union{Array{Float64, 1}, PyObject},
        μ::Union{Array{Float64, 1}, PyObject},
        λ::Union{Array{Float64, 1}, PyObject},
        mesh::Mesh)
    neo_hookean(u::Array{Float64, 1}, μ::Array{Float64, 1}, λ::Array{Float64, 1}, mesh::Mesh)

Computes the elastic stored energy density for  a common neo-Hookean stored energy model.

$$\psi =  \frac{\mu}{2} (I_{c} - 2) - \frac{\mu}{2} \ln(J) + \frac{\lambda}{8}\ln(J)^{2}$$

Here the deformation gradient $F$ is defined by

$$F = I + \nabla u$$

the right Cauchy-Green tensor $C$ is given by

$$C = F^{T} F$$

and the scalars $J$ and $I_c$

$$\begin{split}J     &= \det(C), \\
I_{c} &= {\rm trace}(C).\end{split}$$

$\mu$ and $\lambda$ are Lame parameters.
"""
function neo_hookean(u::Union{Array{Float64, 1}, PyObject},
      μ::Union{Array{Float64, 1}, PyObject},
      λ::Union{Array{Float64, 1}, PyObject},
      mesh::Mesh)
    @assert length(μ)==length(λ)==get_ngauss(mesh)
    @assert length(u) == 2mesh.ndof
    neo_hookean_ = load_op_and_grad(AdFem.libmfem,"neo_hookean", multiple=true)
    u,mu,lamb = convert_to_tensor(Any[u,μ,λ], [Float64,Float64,Float64])
    psi, indices, vv = neo_hookean_(u,mu,lamb)
    J = RawSparseTensor(indices, vv, 2mesh.ndof, 2mesh.ndof)
    set_shape(psi, (2mesh.ndof,)), J
end

function neo_hookean(u::Array{Float64, 1}, μ::Array{Float64, 1}, λ::Array{Float64, 1}, mesh::Mesh)
    @assert length(μ)==length(λ)==get_ngauss(mesh)
    @assert length(u) == 2mesh.ndof
    N = get_ngauss(mesh)*(2*mesh.elem_ndof)^2
    psi = zeros(2mesh.ndof)
    indices = zeros(Int64, 2*N)
    vv = zeros(N)
    @eval ccall((:NH_forward_Julia, $LIBMFEM), Cvoid,
         (Ptr{Cdouble}, Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            $psi, $indices, $vv, $u, $μ, $λ)
    J = sparse(indices[1:2:end] .+ 1, indices[2:2:end] .+ 1, vv, 2mesh.ndof, 2mesh.ndof)
    psi, J
end

function compute_absorbing_boundary_condition_matrix(
    ρ::Union{PyObject, Array{Float64, 1}},
    cs::Union{PyObject, Array{Float64, 1}},
    cp::Union{PyObject, Array{Float64, 1}},
    bdedge::Array{Int64, 2}, mmesh::Mesh)
    @assert length(ρ)==length(cs)==length(cp)==4*size(bdedge, 1)
    ρ, cs, cp = convert_to_tensor([ρ, cs, cp], [Float64, Float64, Float64])
    n = get_edge_normal(bdedge, mmesh)
    m = [n[:,2] -n[:,1]]
    n1 = constant(repeat(n[:,1], 1, 4)'[:])
    n2 = constant(repeat(n[:,2], 1, 4)'[:])
    m1 = constant(repeat(m[:,1], 1, 4)'[:])
    m2 = constant(repeat(m[:,2], 1, 4)'[:])
    a11 = ρ*cp*n1^2+ρ*cs*m1^2
    a12 = ρ*cp*n1*n2+ρ*cs*m1*m2
    a21 = ρ*cp*n1*n2+ρ*cs*m1*m2
    a22 = ρ*cp*n2^2+ρ*cs*m2^2
    A11 = compute_fem_boundary_mass_matrix1(a11, bdedge, mmesh)
    A12 = compute_fem_boundary_mass_matrix1(a12, bdedge, mmesh)
    A21 = compute_fem_boundary_mass_matrix1(a21, bdedge, mmesh)
    A22 = compute_fem_boundary_mass_matrix1(a22, bdedge, mmesh)
    K = [A11 A21
        A12 A22]
end


text:
#*************************************************************************#
#
# Linear Type
#
#*************************************************************************#
"""
    LinearSparse(ninputs, noutputs, i, j)

Initialize `ParametricModel` representing sparse linear operator `z = A*x`.
Where `A` is a sparse matrix with dimension `(noutputs, ninputs) `
cointaining the the parameter vector `Θ` in the position specified
in the elements of `i` and `j`.
"""
struct LinearSparse <: ParametricModel
    ninputs::Int
    noutputs::Int
    nparams::Int
    i::Vector{Int}
    j::Vector{Int}
end

# Constructor
function LinearSparse(ninputs::Int, noutputs::Int, i::Vector{Int}, j::Vector{Int})
    if (length(i) != length(j))
        throw(ArgumentError("Index vector `i` and `j` should have the same length."))
    end
    nparams = length(i)
    LinearSparse(ninputs, noutputs, nparams, i, j)
end

#*************************************************************************#
#
# Generate Buffers
#
#*************************************************************************#
function dx_buffer(lin::LinearSparse)
    return sparse(lin.i, lin.j, ones(lin.nparams), lin.noutputs, lin.ninputs)
end

function dΘ_buffer(lin::LinearSparse)
    row = lin.i
    column = 1:lin.nparams
    return sparse(row, column, ones(lin.nparams), lin.noutputs, lin.nparams)
end

#*************************************************************************#
#
# Compute Output and its Derivatives
#
#*************************************************************************#
function _evaluate!(lin::LinearSparse,
                    x::AbstractVector{<:Real},
                    Θ::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64},
                    dΘ::AbstractMatrix{Float64})

    # TODO: Improve memory allocation

    # Get Matrix
    matrix = sparse(lin.i, lin.j, Θ, lin.noutputs, lin.ninputs)

    # Compute z
    z .= matrix*x

    # Compute dz
    if !isempty(dx)
        dx .= matrix
    end

    # Compute dΘ
    if !isempty(dΘ)
        row = lin.i
        column = 1:lin.nparams
        value = [x[j] for j in lin.j]
        dΘ .= sparse(row, column, value, lin.noutputs, lin.nparams)
    end
end