text:
"""
Backend independent enums which
represent keyboard buttons.
"""
module Keyboard

    using ..AbstractPlotting: INSTANCES # import the docstring extensions
    """
        Keyboard.Button

    Enumerates all keyboard buttons.
    See the implementation for details.
    """
    @enum(Button,
            unknown            = -1,
            # printable keys,
            space              = 32,
            apostrophe         = 39,  # ',
            comma              = 44,  # ,,
            minus              = 45,  # -,
            period             = 46,  # .,
            slash              = 47,  # /,
            _0                  = 48,
            _1                  = 49,
            _2                  = 50,
            _3                  = 51,
            _4                  = 52,
            _5                  = 53,
            _6                  = 54,
            _7                  = 55,
            _8                  = 56,
            _9                  = 57,
            semicolon          = 59,  # ;,
            equal              = 61,  # =,
            a                  = 65,
            b                  = 66,
            c                  = 67,
            d                  = 68,
            e                  = 69,
            f                  = 70,
            g                  = 71,
            h                  = 72,
            i                  = 73,
            j                  = 74,
            k                  = 75,
            l                  = 76,
            m                  = 77,
            n                  = 78,
            o                  = 79,
            p                  = 80,
            q                  = 81,
            r                  = 82,
            s                  = 83,
            t                  = 84,
            u                  = 85,
            v                  = 86,
            w                  = 87,
            x                  = 88,
            y                  = 89,
            z                  = 90,
            left_bracket       = 91,  # [,
            backslash          = 92,  # ,
            right_bracket      = 93,  # ],
            grave_accent       = 96,  # `,
            world_1            = 161, # non-us #1,
            world_2            = 162, # non-us #2,
            # function keys,
            escape             = 256,
            enter              = 257,
            tab                = 258,
            backspace          = 259,
            insert             = 260,
            delete             = 261,
            right              = 262,
            left               = 263,
            down               = 264,
            up                 = 265,
            page_up            = 266,
            page_down          = 267,
            home               = 268,
            _end               = 269,
            caps_lock          = 280,
            scroll_lock        = 281,
            num_lock           = 282,
            print_screen       = 283,
            pause              = 284,
            f1                 = 290,
            f2                 = 291,
            f3                 = 292,
            f4                 = 293,
            f5                 = 294,
            f6                 = 295,
            f7                 = 296,
            f8                 = 297,
            f9                 = 298,
            f10                = 299,
            f11                = 300,
            f12                = 301,
            f13                = 302,
            f14                = 303,
            f15                = 304,
            f16                = 305,
            f17                = 306,
            f18                = 307,
            f19                = 308,
            f20                = 309,
            f21                = 310,
            f22                = 311,
            f23                = 312,
            f24                = 313,
            f25                = 314,
            kp_0               = 320,
            kp_1               = 321,
            kp_2               = 322,
            kp_3               = 323,
            kp_4               = 324,
            kp_5               = 325,
            kp_6               = 326,
            kp_7               = 327,
            kp_8               = 328,
            kp_9               = 329,
            kp_decimal         = 330,
            kp_divide          = 331,
            kp_multiply        = 332,
            kp_subtract        = 333,
            kp_add             = 334,
            kp_enter           = 335,
            kp_equal           = 336,
            left_shift         = 340,
            left_control       = 341,
            left_alt           = 342,
            left_super         = 343,
            right_shift        = 344,
            right_control      = 345,
            right_alt          = 346,
            right_super        = 347,
            menu               = 348,
    )

    """
        Keyboard.Action

    Enumerates all key states/actions in accordance with the GLFW spec.

    $(INSTANCES)
    """
    @enum Action begin
        release = 0
        press   = 1
        repeat  = 2
    end
end

"""
Backend independent enums and fields which
represent mouse actions.
"""
module Mouse
    using ..AbstractPlotting: INSTANCES # import the docstring extensions

    """
        Mouse.Button

    Enumerates all mouse buttons in accordance with the GLFW spec.

    $(INSTANCES)
    """
    @enum Button begin
        left = 0
        middle = 2
        right = 1 # Conform to GLFW
        none = -1 # for convenience
    end

    """
        Mouse.Action

    Enumerates all mouse states/actions in accordance with the GLFW spec.

    $(INSTANCES)
    """
    @enum Action begin
        press   = 1
        release = 0
    end

    # deprecated, remove eventually
    """
        Mouse.DragEnum
    Enumerates the drag states of the mouse.
    $(INSTANCES)
    """
    @enum DragEnum begin
        down
        up
        pressed
        notpressed
    end
end

# Void for no button needs to be pressed,
const ButtonTypes = Union{Nothing, Mouse.Button, Keyboard.Button}

struct KeyEvent
    key::Keyboard.Button
    action::Keyboard.Action
    # mod::ButtonModifier
end

struct MouseButtonEvent
    button::Mouse.Button
    action::Mouse.Action
    # mod::ButtonModifier
end


text:
# functions for integrating over faces

"""
  Does an integral of the form:

  \sum \int f dGamma

  over the interior faces.  This is different than the regular face integrals
  used in the weak form, which are conservative (resL -= integral,
  resR += integral, so the sum is zero).  This is useful for computing
  functionals.  Up to `numDofPerNode` functionals can be computed at the
  same time.

  This function requires parallel communication to already have been completed
  before it is called.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * flux_func: computes the integrand `f` at each face node, `FluxType`

  **Inputs/Outputs

   * val: vector of length `numDofPerNode`, one entry for each functional,
          to be overwritten
"""
function integrateFaceQuantity(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                               flux_func::FluxType, val::Vector) where {Tsol, Tres}

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)

  flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  interfaces = mesh.interfaces
  fill!(val, 0)

  for i=1:length(interfaces)
    iface_i = interfaces[i]

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j)
      qR_j = ro_sview(q_faceR, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(eqn.params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = ro_sview(mesh.nrm_face, :, j, i)
      flux_j = sview(flux_face, :, j)

      flux_func(eqn.params, qL_j, qR_j, aux_vars, nrm_xy, flux_j)

      for k=1:mesh.numDofPerNode
        val[k] += mesh.sbpface.wface[j]*flux_j[k]
      end
    end  # end loop j
  end

  if opts["parallel_data"] == PARALLEL_DATA_FACE
    integrateSharedFaceQuantity(mesh, sbp, eqn, opts, flux_func, val)
  elseif opts["parallel_data"] == PARALLEL_DATA_ELEMENT
    integrateSharedFaceQuantity_element(mesh, sbp, eqn, opts, flux_func, val)
  else
    error("unrecognized parallel data: $(opts["parallel_data"])")
  end

  return nothing
end

"""
  Integrate a quantity over the shared faces and sum.  Assumes parallel
  communication has already been done
"""
function integrateSharedFaceQuantity(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                          flux_func::FluxType, val::Vector) where {Tsol, Tres}

  params = eqn.params

  for idx=1:length(eqn.shared_data)
    data = eqn.shared_data[idx]

    if getParallelData(data) != PARALLEL_DATA_FACE
      throw(ErrorException("cannot use integrateSharedFaceQuantity without parallel face data"))
    end

    # only one process should do the integral for a given face.
    # Arbitrarily choose the lower ranked process
    if data.peernum < mesh.myrank
      continue
    end

    # calculate the flux
    idx = data.peeridx
    interfaces = data.interfaces
    qL_arr = data.q_send
    qR_arr = data.q_recv
    aux_vars_arr = eqn.aux_vars_sharedface[idx]
    nrm_arr = mesh.nrm_sharedface[idx]
    flux_j = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

    for j=1:length(interfaces)
      interface_i = interfaces[j]
      eL = interface_i.elementL
      fL = interface_i.faceL

      # compute the flux
      for k=1:mesh.numNodesPerFace
        qL = ro_sview(qL_arr, :, k, j)
        qR = ro_sview(qR_arr, :, k, j)
        aux_vars = ro_sview(aux_vars_arr, :, k, j)
        parent(aux_vars)[1] = calcPressure(params, qL)
        nrm_xy = ro_sview(nrm_arr, :, k, j)
        flux_k = sview(flux_j, :, k)
        flux_func(params, qL, qR, aux_vars, nrm_xy, flux_k)

        for p=1:mesh.numDofPerNode
          val[p] += mesh.sbpface.wface[k]*flux_k[p]
        end
      end
    end
  end

  return nothing
end

function integrateSharedFaceQuantity_element(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                          flux_func::FluxType, val::Vector) where {Tsol, Tres}
  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)

  for idx=1:length(eqn.shared_data)

    data = eqn.shared_data[idx]

    if getParallelData(data) != PARALLEL_DATA_ELEMENT
      throw(ErrorException("cannot use integrateSharedFaceQuantity without parallel element data"))
    end

    # only one process should do the integral for a given face.
    # Arbitrarily choose the lower ranked process
    if data.peernum < mesh.myrank
      continue
    end

    interfaces = data.interfaces
    qR_arr = data.q_recv
    nrm_arr = mesh.nrm_sharedface[idx]
    aux_vars_arr = eqn.aux_vars_sharedface[idx]
    flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

    start_elnum = mesh.shared_element_offsets[idx]

    for j=1:length(interfaces)
      iface_j = interfaces[j]
      fL = iface_j.faceL

      # interpolate to face
      qL = ro_sview(eqn.q, :, :, iface_j.elementL)
      el_r = iface_j.elementR - start_elnum + 1
      qR = ro_sview(qR_arr, :, :, el_r)
      interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

      # calculate flux
      for k=1:mesh.numNodesPerFace
        qL_k = ro_sview(q_faceL, :, k)
        qR_k = ro_sview(q_faceR, :, k)
        aux_vars = ro_sview(aux_vars_arr, :, k, j)
        nrm_xy = ro_sview(nrm_arr, :, k, j)
        flux_k = sview(flux_face, :, k)

        parent(aux_vars)[1] = calcPressure(eqn.params, qL_k)

        flux_func(eqn.params, qL_k, qR_k, aux_vars, nrm_xy, flux_k)

        for p=1:mesh.numDofPerNode
          val[p] += mesh.sbpface.wface[k]*flux_k[p]
        end
      end  # end k
    end  # end j
  end # end idx

  return nothing
end

#------------------------------------------------------------------------------
# revq

function integrateFaceQuantity_revq(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                               flux_func::FluxType_revq, val_bar::Vector) where {Tsol, Tres}

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)

  q_faceL_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  interfaces = mesh.interfaces

  for i=1:length(interfaces)
    iface_i = interfaces[i]
    fill!(q_faceL_bar, 0); fill!(q_faceR_bar, 0)

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j); qR_j = ro_sview(q_faceR, :, j)
      qL_bar_j = sview(q_faceL_bar, :, j); qR_bar_j = sview(q_faceR_bar, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(eqn.params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = ro_sview(mesh.nrm_face, :, j, i)
      flux_j_bar = sview(flux_face_bar, :, j)

      for k=1:mesh.numDofPerNode
        flux_j_bar[k] = mesh.sbpface.wface[j]*val_bar[k]
      end

      flux_func(eqn.params, qL_j, qL_bar_j, qR_j, qR_bar_j,  aux_vars, nrm_xy, flux_j_bar)

    end  # end loop j

    qL_bar = sview(eqn.q_bar, :, :, iface_i.elementL)
    qR_bar = sview(eqn.q_bar, :, :, iface_i.elementR)
    interiorFaceInterpolate_rev!(mesh.sbpface, iface_i, qL_bar, qR_bar,
                                 q_faceL_bar, q_faceR_bar)

  end

  if opts["parallel_data"] == PARALLEL_DATA_FACE
    error("parallel face data not supported")
  elseif opts["parallel_data"] == PARALLEL_DATA_ELEMENT
    func = (mesh, sbp, eqn, opts, data, data_bar) -> integrateSharedFaceQuantity_element_revq(mesh, sbp, eqn, opts, data, data_bar, flux_func, val_bar)

  else
    error("unrecognized parallel data: $(opts["parallel_data"])")
  end

  # do each set of shared faces and send q_bar back to the owning process
  exchangeData_rev(mesh, sbp, eqn, opts, eqn.shared_data, eqn.shared_data_bar,
                    func)

  # wait for communication to finish
  finishSolutionBarExchange(mesh, sbp, eqn, opts)

  return nothing
end

function integrateSharedFaceQuantity_element_revq(mesh, sbp,
                  eqn::EulerData{Tsol, Tres}, opts, data::SharedFaceData,
                  data_bar::SharedFaceData,
                  flux_func::FluxType_revq,
                  val_bar::Vector) where {Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use integrateSharedFaceQuantity without parallel element data"))
  end

  # only one process should do the integral for a given face.
  # Arbitrarily choose the lower ranked process
  if data.peernum < mesh.myrank
    fill!(data_bar.q_recv, 0)
    return nothing
  end

  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceL_bar = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR_bar = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)

  idx = data_bar.peeridx

  interfaces = data.interfaces
  qR_arr = data.q_recv
  qR_arr_bar = data_bar.q_recv
  nrm_arr = mesh.nrm_sharedface[idx]
  aux_vars_arr = eqn.aux_vars_sharedface[idx]

  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  start_elnum = mesh.shared_element_offsets[idx]

  for j=1:length(interfaces)
    iface_j = interfaces[j]
    fill!(q_faceL_bar, 0); fill!(q_faceR_bar, 0)
    fL = iface_j.faceL

    # interpolate to face
    qL = ro_sview(eqn.q, :, :, iface_j.elementL)
    el_r = iface_j.elementR - start_elnum + 1
    qR = ro_sview(qR_arr, :, :, el_r)
    interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

    # calculate flux
    for k=1:mesh.numNodesPerFace
      qL_k = ro_sview(q_faceL, :, k); qR_k = ro_sview(q_faceR, :, k)
      qL_k_bar = sview(q_faceL_bar, :, k); qR_k_bar = sview(q_faceR_bar, :, k)
      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
      flux_k_bar = sview(flux_face_bar, :, k)

      parent(aux_vars)[1] = calcPressure(eqn.params, qL_k)

      for p=1:mesh.numDofPerNode
        flux_k_bar[p] = mesh.sbpface.wface[k]*val_bar[p]
      end

      flux_func(eqn.params, qL_k, qL_k_bar, qR_k, qR_k_bar, aux_vars, nrm_xy, flux_k_bar)
    end  # end k

    qL_bar = sview(eqn.q_bar, :, :, iface_j.elementL)
    qR_bar = sview(qR_arr_bar, :, :, el_r)
    interiorFaceInterpolate_rev!(mesh.sbpface, iface_j, qL_bar, qR_bar,
                                 q_faceL_bar, q_faceR_bar)
  end  # end j

  return nothing
end

#------------------------------------------------------------------------------
# revm

function integrateFaceQuantity_revm(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                               flux_func::FluxType_revm, val_bar::Vector) where {Tsol, Tres}

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)

  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  interfaces = mesh.interfaces

  for i=1:length(interfaces)
    iface_i = interfaces[i]

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j); qR_j = ro_sview(q_faceR, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(eqn.params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = ro_sview(mesh.nrm_face, :, j, i)
      nrm_bar = sview(mesh.nrm_face_bar, :, j, i)
      flux_j_bar = sview(flux_face_bar, :, j)

      for k=1:mesh.numDofPerNode
        flux_j_bar[k] = mesh.sbpface.wface[j]*val_bar[k]
      end

      flux_func(eqn.params, qL_j, qR_j,  aux_vars, nrm_xy, nrm_bar, flux_j_bar)

    end  # end loop j
  end

  # no need to do parallel communication here because we update
  # mesh.nrm_sharedface_bar, and the mesh reverse mode is responsible for its
  # reverse mode
  if opts["parallel_data"] == PARALLEL_DATA_FACE
    error("parallel face data not supported")
  elseif opts["parallel_data"] == PARALLEL_DATA_ELEMENT
    integrateSharedFaceQuantity_element_revm(mesh, sbp, eqn, opts, flux_func,
                                             val_bar)
  else
    error("unrecognized parallel data: $(opts["parallel_data"])")
  end

  return nothing
end

function integrateSharedFaceQuantity_element_revm(mesh, sbp,
                  eqn::EulerData{Tsol, Tres}, opts,
                  flux_func::FluxType_revm,
                  val_bar::Vector) where {Tsol, Tres}

  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)

  for idx=1:length(eqn.shared_data)
    data = eqn.shared_data[idx]

    if getParallelData(data) != PARALLEL_DATA_ELEMENT
      throw(ErrorException("cannot use integrateSharedFaceQuantity without parallel element data"))
    end

    # only one process should do the integral for a given face.
    # Arbitrarily choose the lower ranked process
    if data.peernum < mesh.myrank
      continue
    end

    interfaces = data.interfaces
    qR_arr = data.q_recv
    nrm_arr = mesh.nrm_sharedface[idx]
    nrm_arr_bar = mesh.nrm_sharedface_bar[idx]
    aux_vars_arr = eqn.aux_vars_sharedface[idx]

    flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

    start_elnum = mesh.shared_element_offsets[idx]

    for j=1:length(interfaces)
      iface_j = interfaces[j]
      fL = iface_j.faceL

      # interpolate to face
      qL = ro_sview(eqn.q, :, :, iface_j.elementL)
      el_r = iface_j.elementR - start_elnum + 1
      qR = ro_sview(qR_arr, :, :, el_r)
      interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

      # calculate flux
      for k=1:mesh.numNodesPerFace
        qL_k = ro_sview(q_faceL, :, k); qR_k = ro_sview(q_faceR, :, k)
        aux_vars = ro_sview(aux_vars_arr, :, k, j)
        nrm_xy = ro_sview(nrm_arr, :, k, j)
        nrm_bar = sview(nrm_arr_bar, :, k, j)
        flux_k_bar = sview(flux_face_bar, :, k)

        parent(aux_vars)[1] = calcPressure(eqn.params, qL_k)

        for p=1:mesh.numDofPerNode
          flux_k_bar[p] = mesh.sbpface.wface[k]*val_bar[p]
        end

        flux_func(eqn.params, qL_k, qR_k, aux_vars, nrm_xy, nrm_bar, flux_k_bar)
      end  # end k
    end  # end j
  end  # end idx

  return nothing
end

ext: jl


text:
using Base.Iterators: product
using Dates
using Documenter: doctest
using Infinity: Infinite, InfExtendedReal, InfExtendedTime, InfMinusInfError, ∞
using Intervals
using Intervals: isfinite
using Serialization: deserialize
using Test
using TimeZones

const BOUND_PERMUTATIONS = product((Closed, Open), (Closed, Open))

include("test_utils.jl")

@testset "Intervals" begin
    include("inclusivity.jl")
    include("endpoint.jl")
    include("interval.jl")
    include("anchoredinterval.jl")
    include("comparisons.jl")
    include("plotting.jl")

    # Note: The output of the doctests currently requires a newer version of Julia
    # https://github.com/JuliaLang/julia/pull/34387
    if VERSION >= v"1.5.0-DEV.163"
        doctest(Intervals)
    else
        @warn "Skipping doctests"
    end
end


text:
function linspace(start::T, stop::T, np::I) where {T<:AbstractFloat, I<:Integer}
    d = (stop-start)/(np-1)
    p = Vector{T}(undef, np)
    p[1] = start
    p[end] = stop
    for i = 2:np-1
        p[i] = start + (i-1)*d
    end
    return p
end

function cosspace(start::T, stop::T, np::I) where {T<:AbstractFloat, I<:Integer}
    p = linspace(zero(T), convert(T,π), np)
    r = (stop-start)/2
    for i = 1:np
        p[i] = start + r*(one(T)-cos(p[i]))
    end
    return p
end

function gen_airfoil(a::AbstractAirfoil{T}, np::Int=65; cos_space::Bool = true) where {T<:AbstractFloat}
    if cos_space
        xc = cosspace(zero(T), one(T), np)
    else
        xc = linspace(zero(T), one(T), np)
    end
    xu = Vector{T}(undef, np)
    yu = Vector{T}(undef, np)
    xl = Vector{T}(undef, np)
    yl = Vector{T}(undef, np)
    @inbounds for i = 1:np
        xu[i] = x_upper(a, xc[i])
        yu[i] = y_upper(a, xc[i])
        xl[i] = x_lower(a, xc[i])
        yl[i] = y_lower(a, xc[i])
    end
    return xu,yu,xl,yl
end

function n_upper(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T<:AbstractFloat}
    θ = atan(dy_upper(a, xc))
    return -sin(θ), cos(θ)
end

function n_lower(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T<:AbstractFloat}
    θ = atan(dy_lower(a, xc))
    return sin(θ), -cos(θ)
end

@inline n_upper(a::AbstractAirfoil{T}, xc::R) where {T<:AbstractFloat, R<:Real} = n_upper(a, convert(T, xc))
@inline n_lower(a::AbstractAirfoil{T}, xc::R) where {T<:AbstractFloat, R<:Real} = n_lower(a, convert(T, xc))

function t_upper(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T<:AbstractFloat}
    dy = dy_upper(a,xc)
    if isinf(dy)
        return zero(T), sign(dy)
    else
        n = sqrt(dy^2+oneunit(T))
        return 1/n, dy/n
    end
end

function t_lower(a::AbstractAirfoil{T}, xc::T)::Tuple{T,T} where {T<:AbstractFloat}
    dy = dy_lower(a,xc)
    if isinf(dy)
        return zero(T), sign(dy)
    else
        n = sqrt(dy^2+oneunit(T))
        return 1/n, dy/n
    end
end

@inline t_upper(a::AbstractAirfoil{T}, xc::R) where {T<:AbstractFloat, R<:Real} = t_upper(a, convert(T, xc))
@inline t_lower(a::AbstractAirfoil{T}, xc::R) where {T<:AbstractFloat, R<:Real} = t_lower(a, convert(T, xc))

for fname in [:x_lower, :x_upper, :y_upper, :y_lower, :dy_upper, :dy_lower]
    @eval begin
        @inline function ($fname)(af::AbstractAirfoil, xc::AbstractVector)
            y = similar(xc)
            @inbounds for i in eachindex(xc)
                y[i] = ($fname)(af, xc[i])
            end
            y
        end
    end
end

for fname in [:fy_upper, :fy_lower, :n_upper, :n_lower, :t_upper, :t_lower]
    @eval begin
        @inline function ($fname)(af::AbstractAirfoil, xc::AbstractVector)
            y = similar(xc, 2, length(xc))
            @inbounds for i in eachindex(xc)
                y[:,i] .= ($fname)(af, xc[i])
            end
            y
        end
    end
end


text:
function insertionsort!(x::AbstractArray)
    n = length(x)
    for i ∈ 1:n
        @inbounds key = x[i]
        j = i-1
        @inbounds while j > 0 && x[j] > key
            @inbounds x[j+1] = x[j]
            j-=1
        end
        @inbounds x[j+1] = key
    end
    return nothing
end


text:
@testset "Quadratic and hermitian forms" begin
  include("QuadForm/Basic.jl")
end


text:
## ASCIIChar ##

abstract ASCIIChar <: Integer

## Conversions ##

convert{T<:ASCIIChar}(::Type{T}, x::Uint8)  = box(T,unbox(Uint8,x))
convert{T<:ASCIIChar}(::Type{T}, x::Int8)   = box(T,unbox(Uint8,x))
convert{T<:ASCIIChar}(::Type{T}, x::Int16)  = box(T,trunc8(unbox(Int16,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Uint16) = box(T,trunc8(unbox(Uint16,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Int32)  = box(T,trunc8(unbox(Int32,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Char)  = box(T,trunc8(unbox(Char,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Uint32) = box(T,trunc8(unbox(Uint32,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Int64)  = box(T,trunc8(unbox(Int64,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Uint64) = box(T,trunc8(unbox(Uint64,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Int128) = box(T,trunc8(unbox(Int128,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Uint128)= box(T,trunc8(unbox(Uint128,x)))
convert{T<:ASCIIChar}(::Type{T}, x::Float32)= box(T,trunc8(checked_fptoui32(unbox(Float32,x))))
convert{T<:ASCIIChar}(::Type{T}, x::Float64)= box(T,trunc8(checked_fptoui64(unbox(Float64,x))))

convert{T<:ASCIIChar}(::Type{Uint8}, x::T)  = box(Uint8,unbox(T,x))
convert{T<:ASCIIChar}(::Type{Int8}, x::T)   = box(Int8,unbox(T,x))
convert{T<:ASCIIChar}(::Type{Int16}, x::T)  = box(Int16,zext16(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Uint16}, x::T) = box(Uint16,zext16(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Int32}, x::T)  = box(Int32,zext32(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Char}, x::T)  = box(Char,zext32(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Uint32}, x::T) = box(Uint32,zext32(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Int64}, x::T)  = box(Int64,zext64(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Uint64}, x::T) = box(Uint64,zext64(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Int128}, x::T) = box(Int128,zext_int(Uint128,unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Uint128}, x::T)= box(Uint128,zext_int(Uint128,unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Float32}, x::T)= box(Float32,uitofp32(unbox(T,x)))
convert{T<:ASCIIChar}(::Type{Float64}, x::T)= box(Float64,uitofp64(unbox(T,x)))

integer(x::ASCIIChar) = int(x)
unsigned(x::ASCIIChar) = uint(x)

## Operations ##

-{T<:ASCIIChar}(x::T, y::Int) = convert(T,( int(x) - y ))
+{T<:ASCIIChar}(x::T, y::Int) = convert(T,( int(x) + y ))
-{T<:ASCIIChar}(x::Int, y::T) = y - x
+{T<:ASCIIChar}(x::Int, y::T) = y + x

## Comparisons ##

=={T<:ASCIIChar}(x::T, y::T) = eq_int(unbox(T,x),unbox(T,y))
!={T<:ASCIIChar}(x::T, y::T) = ne_int(unbox(T,x),unbox(T,y))
<{T<:ASCIIChar}(x::T, y::T) = ult_int(unbox(T,x),unbox(T,y))
<={T<:ASCIIChar}(x::T, y::T) = ule_int(unbox(T,x),unbox(T,y))

## Promotion rules ##

promote_rule{Tb<:ASCIIChar,T<:Number}(::Type{Tb}, ::Type{T}) = T
promote_rule{T<:ASCIIChar}(::Type{T}, ::Type{None}) = T

## Bitwise Operations ##

~{T<:ASCIIChar}(x::T) = box(T,not_int(unbox(T,x)))
(&){T<:ASCIIChar}(x::T, y::T) = box(T,and_int(unbox(T,x),unbox(T,y)))
|{T<:ASCIIChar}(x::T, y::T) = box(T,or_int(unbox(T,x),unbox(T,y)))
($){T<:ASCIIChar}(x::T, y::T) = box(T,xor_int(unbox(T,x),unbox(T,y)))
<<{T<:ASCIIChar}(x::T, y::Int32) = box(T,shl_int(unbox(T,x),unbox(Int32,y)))
>>{T<:ASCIIChar}(x::T, y::Int32) = box(T,lshr_int(unbox(T,x),unbox(Int32,y)))
>>>{T<:ASCIIChar}(x::T, y::Int32) = box(T,lshr_int(unbox(T,x),unbox(Int32,y)))

## Characteristics ##

typemin{T<:ASCIIChar}(::Type{T}) = convert(T,0x00)
typemax{T<:ASCIIChar}(::Type{T}) = convert(T,0xff)

## Show as Character ##

show(io::IO,x::ASCIIChar) = print(io,char(x))

## Bits representation ##

bits(x::ASCIIChar) = bin(reinterpret(Uint8,x),8)

## Array & Strings ##

convert{T<:Union(Uint8,ASCIIChar)}(::Type{Vector{T}}, x::Vector{Union(Uint8,ASCIIChar)}) = reinterpret(T,x)
convert{T<:Union(Uint8,ASCIIChar)}(::Type{Matrix{T}}, x::Matrix{Union(Uint8,ASCIIChar)}) = reinterpret(T,x)

convert{T<:ASCIIChar}(::Type{ASCIIString}, s::Vector{T}) = ASCIIString(reinterpret(Uint8,copy(s)))
convert{T<:ASCIIChar}(::Type{Vector{T}}, s::ASCIIString) = reinterpret(T,copy(s.data))

convert{T<:ASCIIChar}(::Type{Vector{T}}, x::Matrix{T}) = reshape(x,length(x))
convert{T<:ASCIIChar}(::Type{Matrix{T}}, x::Vector{T}) = reshape(x,1,length(x))

bytestring{T<:ASCIIChar}(s::Vector{T}) = bytestring(reinterpret(Uint8,s))

## libc character class testing functions ##

iswascii(c::ASCIIChar) = c < 0x80

for f = (:iswalnum, :iswalpha, :iswblank, :iswcntrl, :iswdigit,
         :iswgraph, :iswlower, :iswprint, :iswpunct, :iswspace,
         :iswupper, :iswxdigit )
    @eval ($f{T<:ASCIIChar})(c::T) = bool(ccall($(expr(:quote,f)), Int32, (T,), c))
end


text:
"""
    PointState

Holds the state variables for a point

# Fields:
 - `u`: Displacement variables for the point (in the global coordinate frame)
 - `theta`: Wiener-Milenkovic rotational displacement variables for the point
 - `F`: Externally applied forces on the point (in the global coordinate frame)
 - `M`: Externally applied moments on the point (in the global coordinate frame)
"""
struct PointState{TF}
    u::SVector{3, TF}
    theta::SVector{3, TF}
    F::SVector{3, TF}
    M::SVector{3, TF}
end

"""
    ElementState

Holds the state variables for an element

# Fields:
 - `u`: Displacement variables for the element (in the global coordinate frame)
 - theta: Wiener-Milenkovic rotational displacement variables for the element
 - `F`: Resultant forces for the element (in the deformed beam coordinate frame)
 - `M`: Resultant moments for the element (in the deformed beam coordinate frame)
 - `P`: Linear momenta of the element (in the deformed beam coordinate frame)
 - `H`: Angular momenta of the element (in the deformed beam coordinate frame)
"""
struct ElementState{TF}
    u::SVector{3, TF}
    theta::SVector{3, TF}
    F::SVector{3, TF}
    M::SVector{3, TF}
    P::SVector{3, TF}
    H::SVector{3, TF}
end

"""
    AssemblyState{TF, TP<:AbstractVector{PointState{TF}}, TE<:AbstractVector{ElementState{TF}}}

Struct for storing state variables for the points and elements in an assembly.

# Fields:
 - `points::TP`: Array of `PointState`s for each point in the assembly
 - `elements::TE`: Array of `ElementState`s for each element in the assembly
"""
struct AssemblyState{TF, TP<:AbstractVector{PointState{TF}}, TE<:AbstractVector{ElementState{TF}}}
    points::TP
    elements::TE
end
Base.eltype(::AssemblyState{TF, TP, TE}) where {TF, TP, TE} = TF

"""
    AssemblyState(system, assembly, x = system.x;
        prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

Post-process the system state given the solution vector `x`.  Return an object
of type `AssemblyState` that defines the state of the assembly for the time step.

If `prescribed_conditions` is not provided, all point state variables are assumed
to be displacements/rotations, rather than their actual identities as used in the
analysis.
"""
function AssemblyState(system, assembly, x = system.x;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}())

    irow_pt = system.irow_pt
    irow_beam = system.irow_beam
    irow_beam1 = system.irow_beam1
    irow_beam2 = system.irow_beam2
    icol_pt = system.icol_pt
    icol_beam = system.icol_beam

    TF = eltype(x)

    istep = system.current_step[]

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    points = Vector{PointState{TF}}(undef, npoint)
    for ipoint = 1:npoint

        if icol_pt[ipoint] <= 0
            # point variables are not state variables, solve for point variables

            # find the first beam connected to the point, checking both sides
            ibeam = findfirst(x -> x == ipoint, assembly.start)
            if !isnothing(ibeam)
                side = -1
            else
                ibeam = findfirst(x -> x == ipoint, assembly.stop)
                side = 1
            end

            beam = assembly.elements[ibeam]

            # get beam state variables
            icol_b = icol_beam[ibeam]
            u_b = SVector(x[icol_b   ], x[icol_b+1 ], x[icol_b+2 ])
            θ_b = SVector(x[icol_b+3 ], x[icol_b+4 ], x[icol_b+5 ])
            F_b = SVector(x[icol_b+6 ], x[icol_b+7 ], x[icol_b+8 ]) .* FORCE_SCALING
            M_b = SVector(x[icol_b+9 ], x[icol_b+10], x[icol_b+11]) .* FORCE_SCALING

            # get beam element properties
            ΔL_b = beam.L
            Ct_b = get_C(θ_b)'
            Cab_b = beam.Cab
            γ_b = element_strain(beam, F_b, M_b)
            κ_b = element_curvature(beam, F_b, M_b)

            # solve for the displacements/rotations of the point
            u = u_b + side*ΔL_b/2*(Ct_b*Cab_b*(e1 + γ_b) - Cab_b*e1)
            theta = θ_b + side*ΔL_b/2*get_Qinv(θ_b)*Cab_b*κ_b
            F = @SVector zeros(3)
            M = @SVector zeros(3)
        else
            # point variables are state variables, extract point state variables
            prescribed = haskey(prescribed_conditions, ipoint)
            if prescribed
                u, theta, F, M = point_variables(x, icol_pt[ipoint],
                    prescribed_conditions[ipoint], istep)
            else
                u, theta, F, M = point_variables(x, icol_pt[ipoint])
            end
        end

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = rotation_parameter_scaling(theta)
        theta *= scaling

        points[ipoint] = PointState{TF}(u, scaling*theta, F, M)
    end

    # all element variables are state variables
    elements = Vector{ElementState{TF}}(undef, nbeam)
    for ibeam = 1:nbeam
        icol = icol_beam[ibeam]
        static = irow_beam[ibeam] <= 0
        u = SVector{3, TF}(x[icol], x[icol+1], x[icol+2])
        theta = SVector{3, TF}(x[icol+3], x[icol+4], x[icol+5])
        F = SVector{3, TF}(x[icol+6], x[icol+7], x[icol+8]) .* FORCE_SCALING
        M = SVector{3, TF}(x[icol+9], x[icol+10], x[icol+11]) .* FORCE_SCALING
        P = ifelse(static, zero(u), SVector{3, TF}(x[icol+12], x[icol+13], x[icol+14]))
        H = ifelse(static, zero(u), SVector{3, TF}(x[icol+15], x[icol+16], x[icol+17]))

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = rotation_parameter_scaling(theta)
        theta *= scaling

        elements[ibeam] = ElementState{TF}(u, theta, F, M, P, H)
    end

    return AssemblyState(points, elements)
end

"""
    write_vtk(name, assembly::Assembly; kwargs...)
    write_vtk(name, assembly::Assembly, state::AssemblyState; kwargs...)
    write_vtk(name, assembly::Assembly, history::Vector{<:AssemblyState}], dt; kwargs...)

Write the deformed geometry (and associated data) to a VTK file for visualization
using ParaView.

The `state` parameter may be omitted to write the original geometry to a VTK file
without any associated data.

If the solution time `history` is provided, the time step must also be provided

# Keyword Arguments
 - `scaling=1.0`: Parameter to scale the deflections (only valid if state is provided)
 - `metadata=Dict()`: Dictionary of metadata for the file(s)
"""
function write_vtk(name, assembly; metadata=Dict())

    # get problem dimensions
    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # extract point locations
    points = Matrix{eltype(assembly)}(undef, 3, npoint)
    for ipoint = 1:npoint
        for i = 1:3
            points[i,ipoint] = assembly.points[ipoint][i]
        end
    end

    # create cells
    cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nbeam]

    # write vtk file
    vtk_grid(name, points, cells) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # add local axis data
        axis_name = ["x-axis", "y-axis", "z-axis"]
        axis_vector = [e1, e2, e3]
        for i = 1:3
            data = Matrix{eltype(assembly)}(undef, 3, nbeam)
            for ibeam = 1:nbeam
                data[:, ibeam] .= assembly.elements[ibeam].Cab*axis_vector[i]
            end
            vtkfile[axis_name[i], VTKCellData()] = data
        end

    end

    return nothing
end

function write_vtk(name, assembly, state; scaling=1.0,
    metadata=Dict())

    # get problem dimensions
    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    # extract point locations
    points = Matrix{eltype(assembly)}(undef, 3, npoint)
    for ipoint = 1:npoint
        for i = 1:3
            points[i,ipoint] = assembly.points[ipoint][i] + scaling*state.points[ipoint].u[i]
        end
    end

    # create cells
    cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nbeam]

    # write vtk file
    vtk_grid(name, points, cells) do vtkfile

        vtkfile["scaling"] = scaling

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # add local axis data
        axis_name = ["x-axis", "y-axis", "z-axis"]
        axis_vector = [e1, e2, e3]
        for i = 1:3
            data = Matrix{eltype(assembly)}(undef, 3, nbeam)
            for ibeam = 1:nbeam
                CtCab = wiener_milenkovic(state.elements[ibeam].theta)'*assembly.elements[ibeam].Cab
                data[:, ibeam] .= CtCab * axis_vector[i]
            end
            vtkfile[axis_name[i], VTKCellData()] = data
        end

        # add point data
        for field in fieldnames(PointState)
            data = Matrix{eltype(assembly)}(undef, 3, npoint)
            for ipoint = 1:npoint
                data[:, ipoint] .= getproperty(state.points[ipoint], field)
            end
            vtkfile[string(field), VTKPointData()] = data
        end

        # add cell data
        for field in fieldnames(ElementState)
            data = Matrix{eltype(assembly)}(undef, 3, nbeam)
            for ibeam = 1:nbeam
                data[:, ibeam] .= getproperty(state.elements[ibeam], field)
            end
            vtkfile[string(field), VTKCellData()] = data
        end

    end

    return nothing
end

function write_vtk(name, assembly, history, dt; scaling=1.0,
    metadata=Dict())

    # get problem dimensions
    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    paraview_collection(name) do pvd

        t = 0.0

        for (current_step, state) in enumerate(history)

            # extract point locations
            points = Matrix{eltype(assembly)}(undef, 3, npoint)
            for ipoint = 1:npoint
                for i = 1:3
                    points[i,ipoint] = assembly.points[ipoint][i] + scaling*state.points[ipoint].u[i]
                end
            end

            # create cells
            cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nbeam]

            # write vtk file
            vtkfile = vtk_grid(name*"-step$current_step", points, cells)

            # add metadata
            vtkfile["scaling"] = scaling
            vtkfile["time"] = t
            for (key, value) in pairs(metadata)
                vtkfile[string(key)] = value
            end

            # add local axis data
            axis_name = ["x-axis", "y-axis", "z-axis"]
            axis_vector = [e1, e2, e3]
            for i = 1:3
                data = Matrix{eltype(assembly)}(undef, 3, nbeam)
                for ibeam = 1:nbeam
                    CtCab = wiener_milenkovic(state.elements[ibeam].theta)' * assembly.elements[ibeam].Cab
                    data[:, ibeam] .= CtCab * axis_vector[i]
                end
                vtkfile[axis_name[i], VTKCellData()] = data
            end

            # add point data
            for field in fieldnames(PointState)
                data = Matrix{eltype(assembly)}(undef, 3, npoint)
                for ipoint = 1:npoint
                    data[:, ipoint] .= getproperty(state.points[ipoint], field)
                end
                vtkfile[string(field), VTKPointData()] = data
            end

            # add cell data
            for field in fieldnames(ElementState)
                data = Matrix{eltype(assembly)}(undef, 3, nbeam)
                for ibeam = 1:nbeam
                    data[:, ibeam] .= getproperty(state.elements[ibeam], field)
                end
                vtkfile[string(field), VTKCellData()] = data
            end

            pvd[t] = vtkfile

            # increment time
            t += dt
        end
    end

    return nothing
end

"""
    write_vtk(name, assembly::Assembly, [state::AssemblyState, ]λ::Number, eigenstate::AssemblyState;
    scaling=1.0, mode_scaling=1.0, cycles=1, steps=100)

Write a series of files corresponding to the elastic motion of the `assembly`
about the deformed state encoded in `state` defined by the eigenvalue `λ` and
the eigenvector encoded in `eigenstate` over the time period specified by `time`.

The steady-state deflections can be scaled with `scaling` and the eigenmode
deflections can be scaled using `mode_scaling`.

The current time is encoded in the metadata tag "time"
"""
function write_vtk(name, assembly, state, λ, eigenstate;
    scaling=1.0,
    mode_scaling=1.0,
    cycles = 1,
    steps = 100)

    npoint = length(assembly.points)
    nbeam = length(assembly.elements)

    damping = real(λ)
    period = 2*pi/imag(λ)

    if damping <= 0.0
        start = -period*cycles
        stop = 0.0
    else
        start = 0.0
        stop = period*cycles
    end

    time=range(start, stop, length=steps)

    paraview_collection(name) do pvd

        for (current_step, t) in enumerate(time)

            # extract point locations
            points = Matrix{eltype(assembly)}(undef, 3, npoint)
            for ipoint = 1:npoint
                for i = 1:3
                    points[i,ipoint] = assembly.points[ipoint][i] +
                        scaling*state.points[ipoint].u[i] +
                        mode_scaling*real(eigenstate.points[ipoint].u[i]*exp(λ*t))
                end
            end

            # create cells
            cells = [MeshCell(PolyData.Lines(), [assembly.start[i], assembly.stop[i]]) for i = 1:nbeam]

            # write vtk file
            vtkfile = vtk_grid(name*"-step$(current_step)", points, cells)

            # add metadata
            vtkfile["scaling"] = scaling
            vtkfile["time"] = t
            vtkfile["phase"] = 2*pi*t/period

            # add local axis data
            axis_name = ["x-axis", "y-axis", "z-axis"]
            axis_vector = [e1, e2, e3]
            for i = 1:3
                data = Matrix{eltype(assembly)}(undef, 3, nbeam)
                for ibeam = 1:nbeam
                    CtCab = wiener_milenkovic(real(state.elements[ibeam].theta .* exp(λ*t)))' *
                        assembly.elements[ibeam].Cab
                    data[:, ibeam] .= CtCab * axis_vector[i]
                end
                vtkfile[axis_name[i], VTKCellData()] = data
            end

            # add point data
            for field in fieldnames(PointState)
                data = Matrix{eltype(assembly)}(undef, 3, npoint)
                for ipoint = 1:npoint
                    data[:, ipoint] .= getproperty(state.points[ipoint], field) +
                         mode_scaling*real(getproperty(eigenstate.points[ipoint], field) .* exp(λ*t))
                end
                vtkfile[string(field), VTKPointData()] = data
            end

            # add cell data
            for field in fieldnames(ElementState)
                data = Matrix{eltype(assembly)}(undef, 3, nbeam)
                for ibeam = 1:nbeam
                    data[:, ibeam] .= getproperty(state.elements[ibeam], field) +
                         mode_scaling*real(getproperty(eigenstate.elements[ibeam], field).* exp(λ*t))
                end
                vtkfile[string(field), VTKCellData()] = data
            end

            pvd[t] = vtkfile
        end
    end

    return nothing
end

"""
    left_eigenvectors(system, λ, V)
    left_eigenvectors(K, M, λ, V)

Compute the left eigenvector matrix `U` for the `system` using inverse power
iteration given the eigenvalues `λ` and the corresponding right eigenvector
matrix `V`.

The complex conjugate of each left eigenvector is stored in each row of the
matrix `U`

Left and right eigenvectors satisfy the following M-orthogonality condition:
 - u'*M*v = 1 if u and v correspond to the same eigenvalue
 - u'*M*v = 0 if u and v correspond to different eigenvalues
This means that U*M*V = I

This function assumes that `system` has not been modified since the eigenvalues
and right eigenvectors were computed.
"""
left_eigenvectors(system, λ, V) = left_eigenvectors(system.K, system.M, λ, V)

function left_eigenvectors(K, M, λ, V)

    # problem type and dimensions
    TC = eltype(V)
    nx = size(V,1)
    nev = size(V,2)

    # allocate storage
    U = rand(TC, nev, nx)
    u = Vector{TC}(undef, nx)
    tmp = Vector{TC}(undef, nx)

    # get entries in M
    iM, jM, valM = findnz(M)

    # compute eigenvectors for each eigenvalue
    for iλ = 1:nev

        # factorize (K + λ*M)'
        KmλMfact = factorize(K' - λ[iλ]'*M')

        # initialize left eigenvector
        for i = 1:nx
            u[i] = U[iλ,i]
        end

        # perform a few iterations to converge the left eigenvector
        for ipass = 1:3
            # get updated u
            mul!(tmp, M, u)
            ldiv!(u, KmλMfact, tmp)
            # normalize u
            unorm = zero(TC)
            for k = 1:length(valM)
                unorm += conj(u[iM[k]])*valM[k]*V[jM[k],iλ]
            end
            rdiv!(u, conj(unorm))
        end

        # store conjugate of final eigenvector
        for i = 1:nx
            U[iλ,i] = conj(u[i])
        end
    end

    return U
end

"""
    correlate_eigenmodes(C)

Return the permutation and the associated corruption index vector which associates
eigenmodes from the current iteration with those of the previous iteration given
the correlation matrix `C`.

The correlation matrix can take one of the following forms (in order of preference):
 - `C = U_p*M*V`
 - `C = U*M_p*V_p`
 - `C = V_p'*V`
 - `C = V'*V_p`
where `U` is a matrix of conjugated left eigenvectors, `M` is the system mass
matrix, `V` is a matrix of right eigenvectors, and `()_p` indicates a variable
from the previous iteration.

Note that the following two forms of the correlation matrix seem to be significantly
inferior to their counterparts listed above: `C = U*M*V_p` and `C = U_p*M_p*V`.
This is likely due to the way in which the left eigenvector matrix is calculated.

The corruption index is the largest magnitude in a given row of `C`
that was not chosen divided by the magnitude of the chosen eigenmode.  It is most
meaningful when using one of the forms of the correlation matrix that uses left
eigenvectors since correct eigenmodes will have magnitudes close to 1 and
incorrect eigenmodes will have magnitudes close to 0.

If the new mode number is already assigned, the next highest unassigned mode
number is used.  In this case a corruption index higher than 1 will be returned,
otherwise the values of the corruption index will always be bounded by 0 and 1.

See "New Mode Tracking Methods in Aeroelastic Analysis" by Eldred, Vankayya, and
Anderson.
"""
function correlate_eigenmodes(C)

    # get row permutation that puts maximum values on the diagonals
    nev = size(C, 1)
    tmp = Vector{real(eltype(C))}(undef, nev)
    perm = zeros(Int, nev)
    corruption = Vector{real(eltype(C))}(undef, nev)
    for i = 1:nev

        # rank each value in this row
        for j = 1:nev
            tmp[j] = abs(C[i,j])
        end
        ranked_modes = sortperm(tmp, rev=true)

        # choose the best fit that is not yet assigned
        i1 = ranked_modes[findfirst((x) -> !(x in perm), ranked_modes)]
        i2 = i1 == ranked_modes[1] ? ranked_modes[2] : ranked_modes[1]

        # assign best eigenmode fit, create corruption index
        perm[i] = i1
        corruption[i] = tmp[i2]/tmp[i1]
    end

    return perm, corruption
end


text:
#=##############################################################################
# DESCRIPTION
    Simulation interface connecting QVLMVehicle vehicles and maneuvers.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Apr 2020
  * License   : MIT
=###############################################################################

function solve(self::Simulation{V, M, R}, Vinf::Function,
                pfield::vpm.ParticleField, wake_coupled::Bool,
                dt::Real, rlx::Real, sigma::Real, rho::Real,
                speedofsound::Real, staticpfield::vpm.ParticleField;
                init_sol::Bool=false, vlm_fsgm=1
                ) where {V<:QVLMVehicle, M<:AbstractManeuver, R}

    # Bring the rotors back to their initial positions at beginning of this
    # time step before `precalculations`. This is needed to avoid double
    # accounting for the RPM when the solver is run
    rotate_rotors(self, -dt)

    vhcl = self.vehicle
    mnvr = self.maneuver
    t = self.t

    # TIME-STEPPING PROCEDURE:
    # -1) Solve one particle field time step
    # 0) Translate and rotate systems
    # 1) Recalculate horseshoes with kinematic velocity
    # 2) Paste previous Gamma solution to new system position after translation
    # 3) Shed semi-infinite wake after translation
    # 4) Calculate wake-induced velocity on VLM and Rotor system
    # 5) Solve VLM and Rotor system
    # 6) Shed unsteady-loading wake after new solution
    # 7) Save new solution as prev solution
    # Iterate
    #
    # On the first time step (pfield.nt==0), it only does steps (5) and (7),
    # meaning that the unsteady wake of the first time step is never shed.
    if self.nt==0
        # NOTE: VLMs and Rotors solutions are losely coupled
        # NOTE: VLMs surfaces and their wakes are also losely coupled

        # Set Vinf on system to get horseshoes
        vlm.setVinf(vhcl.system, Vinf)

        # Solve VLMs
        # TODO: Add Rotor-on-VLM induced velocity
        # NOTE: Since there isn't a kinematic velocity in the first time step
        #   (t==0), the semi-infinite wake would start from a back initial
        #   solution if calculated in this step, so this step is rather going
        #   to waste and we wait for the next step to initiate the wake properly
        if init_sol
            # Here I use the freestream-wake for the first step
            vlm.solve(vhcl.vlm_system, Vinf; t=t, keep_sol=true)
        else
            # Here I completely ignore the wake
            vlm.solve(vhcl.vlm_system, Vinf; t=t, keep_sol=true,
                                                    vortexsheet=(X,t)->zeros(3))
        end

        # Calculate induced velocities to use in rotor solver
        ## Points where to calculate induced velocities
        Xs = _get_midXs(vhcl.rotor_systems)

        ## Evaluate VLM-on-Rotor induced velocity
        Vinds = [vlm.Vind(vhcl.vlm_system, X; t=t, ign_infvortex=init_sol)
                                                                    for X in Xs]

        # Solve Rotors
        for (si, rotors) in enumerate(vhcl.rotor_systems)

            # RPM of this rotor system
            RPM = self.RPMref*get_RPM(mnvr, si, t/self.ttot)

            for (ri, rotor) in enumerate(rotors)

                # Get velocities induced at every blade of this rotor
                this_Vinds = _parse_midXs(vhcl.rotor_systems,  Vinds, si, ri)

                VindVkin = _format_blades(this_Vinds, vhcl.rotor_systems, si, ri)

                vlm.solvefromCCBlade(rotor, Vinf, RPM, rho; _Vinds=VindVkin,
                                              t=t, sound_spd=speedofsound)
            end
        end

    else

        # Initialize the wake system as a VLM of only the wake system on the
        # kinematic velocity
        if self.nt==1
            vlm.solve(vhcl.wake_system, Vinf; t=t, extraVinf=_extraVinf1,
                                keep_sol=true, vortexsheet=(X,t)->zeros(3))
        end

        # Previous rotor system for kinematic and induced velocity calculation
        prev_rotor_systems = _get_prev_rotor_systems(vhcl)

        # ---------- 4) Calculate VPM velocity on VLM and Rotor system -----

        # Control points of VLM
        Xs = _get_Xs(vhcl.vlm_system, "CP"; t=t)

        # Wake-induced velocity at every VLM control point
        Vinds = [vlm.Vind(vhcl.wake_system, X; t=t, only_infvortex=true)
                                                                    for X in Xs]

        # Rotor-induced velocity at every VLM control point (only blade-induced
        # velocities---this ignores the rotor slipstream)
        for (si, rotors) in enumerate(vhcl.rotor_systems)
            for (ri, rotor) in enumerate(rotors)

                vlm.getHorseshoe(rotor, 1)          # Force HS calculation
                vlm._addsolution(rotor, "Gamma",    # Give previous solution
                                prev_rotor_systems[si][ri]._wingsystem.sol["Gamma"])

                Vinds .+= [vlm.Vind(rotor, X; t=t, ign_infvortex=true)
                                                                    for X in Xs]
            end
        end

        vlm._addsolution(vhcl.vlm_system, "Vind", Vinds; t=t)

        # Add dummy VPM velocity for `calc_aerodynamicforce()`
        Vdummy = zeros(3)
        Vdummies = [Vdummy for i in 1:vlm.get_m(vhcl.vlm_system)]
        vlm._addsolution(vhcl.vlm_system, "Vvpm_ApA", Vdummies; t=t)
        vlm._addsolution(vhcl.vlm_system, "Vvpm_AB", Vdummies; t=t)
        vlm._addsolution(vhcl.vlm_system, "Vvpm_BBp", Vdummies; t=t)

        # Calculate induced velocities to use in rotor solver
        ## Points where to calculate induced velocities
        Xs = _get_midXs(vhcl.rotor_systems)

        ## Evaluate wake-on-Rotor induced velocity
        Vinds = [vlm.Vind(vhcl.wake_system, X; t=t, only_infvortex=true)
                                                                    for X in Xs]
        ## Evaluate VLM-on-Rotor induced velocity
        Vinds .+= [vlm.Vind(vhcl.vlm_system, X; t=t, ign_infvortex=true)
                                                                    for X in Xs]
        ## Evaluate Rotor-on-Rotor (and self) induced velocity
        # NOTE: This include the blade-induced velocity of the rotor on itself,
        #   which is already accounted for in the momentum-balance part of BEMT.
        #   Remember to substract this later or the solver will get on a
        #   positive-feedback instability.
        for rotors in vhcl.rotor_systems
            for rotor in rotors
                Vinds .+= [vlm.Vind(rotor, X; t=t, ign_infvortex=true)
                                                                    for X in Xs]
            end
        end

        # ---------- 5) Solve VLM system -----------------------------------
        vlm.solve(vhcl.vlm_system, Vinf; t=t, extraVinf=_extraVinf3,
                                keep_sol=true, vortexsheet=(X,t)->zeros(3))

        # Relaxes (rlx->1) or stiffens (rlx->0) the VLM solution
        if rlx >= 0
            rlxd_Gamma = rlx*vhcl.vlm_system.sol["Gamma"] +
                                (1-rlx)*_get_prev_vlm_system(vhcl).sol["Gamma"]
            vlm._addsolution(vhcl.vlm_system, "Gamma", rlxd_Gamma)
        end

        # ---------- 5) Solve Rotor system ---------------------------------
        for (si, rotors) in enumerate(vhcl.rotor_systems)

            # RPM of this rotor system
            RPM = self.RPMref*get_RPM(mnvr, si, t/self.ttot)

            for (ri, rotor) in enumerate(rotors)

                # Calculate kinematic velocities
                Vkin = _Vkinematic_rotor(vhcl.rotor_systems,
                                         prev_rotor_systems, si, ri, dt)

                # Get velocities induced at every blade of this rotor
                this_Vinds = _parse_midXs(vhcl.rotor_systems, Vinds, si, ri)

                # Calculate and substract the blade-induced velocity of the
                # rotor on itself
                this_Xs = _parse_midXs(vhcl.rotor_systems, Xs, si, ri)
                this_Vinds .-= [vlm.Vind(rotor, X; t=t, ign_infvortex=true)
                                                               for X in this_Xs]

                VindVkin = _format_blades(this_Vinds + Vkin, vhcl.rotor_systems,
                                                                        si, ri)

                vlm.solvefromCCBlade(rotor, Vinf, RPM, rho; _Vinds=VindVkin,
                                            t=t, sound_spd=speedofsound)
            end
        end

        # Rotate rotors
        # NOTE: this means that the rotational velocity hadn't been included in
        # the kinematic velocitie up to this point
        rotate_rotors(self, dt)
    end
end


text:
"""
    ticr!(net, D::DataFrame, optimizeBL::Bool)
    ticr!(net, D::DataCF,    optimizeBL::Bool)
    ticr(D::DataCF)

Goodness-of-fit test for the adequacy of the multispecies network coalescent,
to see if a given population or species network explains the
quartet concordance factor data adequately (see Stenz et al 2015
and [addendum](http://www.stat.wisc.edu/~ane/publis/2015Stenz_TICR_addendum.pdf)
for the method on trees)
Obviously the acronym TICR (Tree Incongruence Checking with R) becomes outdated:
with a method extended to networks and implemented in Julia. Oh well :smiley:.

The tree / network needs to have branch lengths in coalescent units,
must be fully resolved, and must be of level 1.

The model assumes a Dirichlet distribution for the observed quartet concordance
factor, with concentration parameter estimated from the data. An outlier p-value
is calculated for each four-taxon set. Four-taxon sets are then binned into
categories according to their p-values: `0-0.01`, `0.01-0.05`, `0.05-0.10`, and `0.10-1`.
Finally, a chi-square goodness-of-fit test is performed on these binned
frequency with 3 degrees of freedom, to determine if they departs from the
expected proportions `(0.01, 0.04, 0.05, 0.90)`.

- The first version takes a DataFrame object where each row corresponds to a given
  four-taxon set. The DataFrame is modified by having an additional another column
  containing the p-values corresponding to each four-taxon set.
- The second version takes a DataCF object and modifies it by updating
  the expected concordance factors stored in that object.
- The last version (which all others call) assumes that the expected concordance
  factors in the DataCF object are correctly calculated from the test network.

`optimizeBL`: when false, the loglik field of `net` is updated;
when `true`, a copy of `net` with updated branch lengths (in coalescent units)
and update loglik is returned.

output:

- p-value of the χ^2 test
- χ^2 statistic
- value of the pseudo likelihood
- value of the concentration parameter α
- a vector of outlier p-values, one for each four-taxon set
- network (first and second versions):
  `net` with loglik field updated if `optimizeBL` is false;
   copy of `net` with optimized branch lengths and loglik if `optimizeBL` is true

# References

NWM Stenz, B Larget, DA Baum and C Ané (2015).
Exploring tree-like and non-tree-like patterns using genome sequences:
An example using the inbreeding plant species *Arabidopsis thaliana* (L.) Heynh.
Systematic Biology, 64(5):809-823. doi: 10.1093/sysbio/syv039
"""
function ticr!(net, D::DataFrame, optimizeBL::Bool)
    d = readTableCF(D);
    res = ticr!(net, d, optimizeBL); # = (chisq_pval, chisq, pseudolik, alpha, pval)
    D[:p_value] = res[5] # order of value in results: (chisq_pval, chisq, pseudolik, alpha, pval)
    return res
end

function ticr!(net, D::DataCF, optimizeBL::Bool)
    if optimizeBL
        net = topologyMaxQPseudolik!(net,D);
    else
        topologyQPseudolik!(net,D);
    end
    res = ticr(D);
    return (res..., net) # (chisq_pval, chisq, pseudolik, alpha, pval, net)
end

@doc (@doc ticr!) ticr
function ticr(D::DataCF)
    res_alpha = ticr_optimalpha(D)
    alpha = res_alpha[2][1]
    pseudolik = res_alpha[1]
    N = D.numQuartets
    pval = Float64[]
    for i in 1:N
        phat = D.quartet[i].obsCF
        p = D.quartet[i].qnet.expCF
        max = findmax(p)
        p_max = max[1]
        max_idx = max[2]
        p_max_hat = getindex(phat,max_idx)
        p_sort = sort(D.quartet[i].qnet.expCF)
        abs(p_max-p_sort[end-1]) > 1e-6 || warn("Check the network for major quartet")
        d = abs(p_max_hat - p_max)
        temp = [1-(1-p_max)*alpha/2, 0.0]
        shapeAdd = maximum(temp)
        ipval = StatsFuns.betacdf(alpha*p_max+shapeAdd, alpha*(1-p_max)+2*shapeAdd, p_max-d)+
        StatsFuns.betaccdf(alpha*p_max+shapeAdd, alpha*(1-p_max)+2*shapeAdd, p_max+d)
        push!(pval,ipval)
    end
    pcat = CategoricalArrays.cut(pval,[0, 0.01, 0.05, 0.1, 1]) # CategoricalArrays is required by DataFrames
    count = StatsBase.countmap(pcat)
    c = Float64[]
    e = [0.01,0.04,0.05,0.90]*N
    interval = ["[0.0, 0.01)","[0.01, 0.05)","[0.05, 0.1)","[0.1, 1.0)"]
    for i in 1:length(interval)
        if haskey(count,interval[i])
            push!(c,count[interval[i]])
        else
            push!(c,0)
        end
    end
    chisq = 0.0
    for i in 1:length(c)
        chisq += (c[i]-e[i])^2 / e[i]
    end
    chisq_pval = StatsFuns.chisqccdf(3,chisq)
    return (chisq_pval, chisq, pseudolik, alpha, pval)
end

"""
    ticr_optimalpha(D::DataCF)

Find the concentration parameter α by maximizing the pseudo-log-likelihood
of observed quartet concordance factors. The model assumes a Dirichlet distribution
with mean equal to the expected concordance factors calculated from a
phylogenetic network (under ILS and reticulation). These expected CFs
are assumed to be already calculated, and stored in D.

When calculating the pseudo-log-likelihood, this function will check the
observed concordance factors for any values equal to zero: they cause a problem
because the Dirichlet density is 0 at 0 (for concentrations > 1).
Those 0.0 observed CF values are re-set to the minimum of:
- the minimum of all expected concordance factors, and
- the minimum of all nonzero observed concordance factors.

output:

- maximized pseudo-loglikelihood
- value of α where the pseudo-loglikelihood is maximized
- return code of the optimization

The optimization uses NLOpt, with the `:LN_BOBYQA` method.
Optional arguments can tune the optimization differently:
`NLoptMethod`, `xtol_rel` (1e-6 by default),
starting α value `x_start` (1.0 by default).
"""
function ticr_optimalpha(D::DataCF; x_start=1.0::Float64,
        NLoptMethod=:LN_BOBYQA::Symbol, xtol_rel=1e-6::Float64)
    M = D.numQuartets
    logCFtilde = 0.0
    minExpCF = minimum([minimum(q.qnet.expCF) for q in D.quartet])
    minNonZeroObsCF = minimum([minimum(filter(!iszero,q.obsCF)) for q in D.quartet])
    minObsCF = min(minExpCF,minNonZeroObsCF)
    for q in D.quartet
        q.obsCF[find(q.obsCF .== 0)] = minObsCF
    end
    for i in 1:M
        lobsCF = log.(D.quartet[i].obsCF)
        p = D.quartet[i].qnet.expCF
        for i in 1:length(lobsCF)
            logCFtilde += lobsCF[i]*p[i]
        end
    end
    function obj(x::Vector, grad::Vector)
        a = x[1]
        sum = 0
        logCFbar = 0.0
        for i in 1:M
            lobsCF = log.(D.quartet[i].obsCF)
            p = D.quartet[i].qnet.expCF
            p_min = minimum(p)
            if a >= 1/p_min
                b = 0
            else
                b = 1 - a * p_min
            end
            logCFbar += mean(lobsCF)*(1-b)
            sum += lgamma(a+3*b) - lgamma(a*p[1]+b) - lgamma(a*p[2]+b) - lgamma(a*p[3]+b)
        end
        PseudologLik = sum + a*logCFtilde-3*logCFbar
        return PseudologLik
    end
    opt = NLopt.Opt(NLoptMethod,1)
    NLopt.lower_bounds!(opt, 0)
    NLopt.upper_bounds!(opt, 10^5)
    NLopt.maxeval!(opt,1000)
    NLopt.xtol_rel!(opt,xtol_rel)
    NLopt.max_objective!(opt,obj)
    fmax, xmax, ret = NLopt.optimize(opt,[x_start])
    return fmax, xmax, ret
end


text:
# # 1D Hawkes process with delayed alpha kernel

#=
In this example, I simulate either a 1D or a 2D Hawkes process, with a
delayed-alpha interaction kernel

```math
g(t) = H(t-t_{\text{delay}}) \,  \frac{(t-t_{\text{delay}})}{\tau^2} \,
 \exp\left(- \frac{(t-t_{\text{delay}})}{\tau} \right)
```

where $H(x)$ is the Heaviside function: $H(x)=0$ for $x<0$, $H(x)=1$ for $t\geq 0$.

Kernels are always normalized so that their integral is 1.
=#

# ## Initialization
using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:default) #; plotlyjs()
using FFTW

using ProgressMeter
using Random
Random.seed!(0)

using HawkesSimulator; const global H = HawkesSimulator

##
# ## Define and visualize the kernel

mytau = 0.5
mydelay = 2.0
myw = fill(0.85,(1,1)) # 1x1 matrix
myinput = [0.5,]       # 1-dim vector

ker = H.KernelAlphaDelay(mytau,mydelay)

plt = let ts = range(-0.5,8;length=150)
  y = [H.interaction_kernel(_t,ker) for _t in ts]
  plt=plot(ts,y ; linewidth=3,leg=false,xlabel="time (s)",
      ylabel="interaction kernel")
  ymax=ylims()[2]
  plot!(plt,[0,0],[0,ymax];linestyle=:dash,linecolor=:black)
end;
plot(plt)

# Note that the kernel starts after zero, according to the delay indicated.

# As a side note: in order to simulate Hawkes processes, one always
# needs to define a *non-increasing upper limit* to the kernel.
# This is what it looks like for this kernel.

plt = let  ts = range(-0.5,8;length=150)
  y = [H.interaction_kernel(_t,ker) for _t in ts]
  yu = [H.interaction_kernel_upper(_t,ker) for _t in ts]
  plt = plot(ts,y ; linewidth=3,xlabel="time (s)",
      ylabel="interaction kernel", label="true kernel")
  plot!(plt, ts,yu ; linewidth=2, label="upper limit", linestyle=:dash)
  ymax=ylims()[2]
  plot!(plt,[0,0],[0,ymax];linestyle=:dash,linecolor=:black,label="")
end;
plot(plt)

# the closer the upper limit is to the true kernel,
# the more efficient the simulation.

##      #src
# ## Build the network and run it
# I compare the final rate with what I expect from the analytic solution
# (see first example file)

pops = H.PopulationState(ker,1)
ntw = H.RecurrentNetwork(pops,myw,myinput)

function run_simulation!(network,num_spikes,
    t_flush_trigger=300.0,t_flush=100.0)
  t_now = 0.0
  H.reset!(network) # clear spike trains etc
  @showprogress "Running Hawkes process..." for _ in 1:num_spikes
    t_now = H.dynamics_step_singlepopulation!(t_now,network)
    H.flush_trains!(network,t_flush_trigger;Tflush=t_flush)
  end
  H.flush_trains!(network) # flush everything into history
  return t_now
end

n_spikes = 80_000
Tmax = run_simulation!(ntw,n_spikes);

ratenum = H.numerical_rates(pops)[1]
rate_analytic = (I-myw)\myinput
rate_analytic = rate_analytic[1] # from 1-dim vector to scalar

@info "Mean rate -  numerical $(round(ratenum;digits=2)), analytic  $(round(rate_analytic;digits=2))"

##
# ## Covariance density

# First, compute covariance density numerically for a reasonable time step
mytrain = pops.trains_history[1]
mydt = 0.1
myτmax = 60.0
mytaus = H.get_times(mydt,myτmax)
ntaus = length(mytaus)
cov_num = H.covariance_self_numerical(mytrain,mydt,myτmax);

# now compute covariance density
#  analytically (as in Hawkes models), at higher resolution,
# and compare analytic and numeric
#
# Note that the high resolution is not just for a better plot, but also to
# ensure the result is more precise when we move from frequency domain
# (Fourier transforms) to time domain.

function four_high_res(dt::Real,Tmax::Real) # higher time resolution, longer time
  k1,k2 = 2 , 0.01
  myτmax = Tmax * k1
  dt *= k2
  mytaus = H.get_times(dt,myτmax)
  nkeep = div(length(mytaus),k1)
  myfreq = H.get_frequencies_centerzero(dt,myτmax)
  gfou = myw[1,1] .* H.interaction_kernel_fourier.(myfreq,Ref(ker)) |> ifftshift
  ffou = let r=rate_analytic
    covf(g) = r/((1-g)*(1-g'))
    map(covf,gfou)
  end
  retf = real.(ifft(ffou))
  retf[2:end] ./= dt # first element is rate
  return mytaus[1:nkeep],retf[1:nkeep]
end

taush,covfou=four_high_res(mydt,myτmax)

plt= let plt = plot(xlabel="time delay (s)",ylabel="Covariance density")
  plot!(plt,mytaus[2:end], cov_num[2:end] ; linewidth=3, label="numerical" )
  plot!(plt,taush[2:end],covfou[2:end]; label="analytic",linewidth=3,linestyle=:dash)
end;

plot(plt)

# Analytics and numerics match quite well.

# **THE END**

# Literate.markdown("examples/alphadelay.jl","docs/src";documenter=true,repo_root_url="https://github.com/dylanfesta/HawkesSimulator.jl/blob/master") #src


text:
#
# Types supporting Parameters and their connections
#

abstract type ModelParameter <: MimiStruct end

# TBD: rename as ScalarParameter, ArrayParameter, and AbstractParameter?

mutable struct ScalarModelParameter{T} <: ModelParameter
    value::T

    function ScalarModelParameter{T}(value::T) where T
        new(value)
    end

    function ScalarModelParameter{T1}(value::T2) where {T1, T2}
        try
            new(T1(value))
        catch err
            error("Failed to convert $value::$T2 to $T1")
        end
    end
end

mutable struct ArrayModelParameter{T} <: ModelParameter
    values::T
    dim_names::Vector{Symbol} # if empty, we don't have the dimensions' name information

    function ArrayModelParameter{T}(values::T, dims::Vector{Symbol}) where T
        new(values, dims)
    end
end

ScalarModelParameter(value) = ScalarModelParameter{typeof(value)}(value)

Base.convert(::Type{ScalarModelParameter{T}}, value::Number) where {T} = ScalarModelParameter{T}(T(value))

Base.convert(::Type{T}, s::ScalarModelParameter{T}) where {T} = T(s.value)

ArrayModelParameter(value, dims::Vector{Symbol}) = ArrayModelParameter{typeof(value)}(value, dims)

# Allow values to be obtained from either parameter type using one method name.
value(param::ArrayModelParameter)  = param.values
value(param::ScalarModelParameter) = param.value

Base.copy(obj::ScalarModelParameter{T}) where T = ScalarModelParameter(obj.value)
Base.copy(obj::ArrayModelParameter{T}) where T = ArrayModelParameter(obj.values, obj.dim_names)

dim_names(obj::ArrayModelParameter) = obj.dim_names
dim_names(obj::ScalarModelParameter) = []

abstract type AbstractConnection <: MimiStruct end

struct InternalParameterConnection <: AbstractConnection
    src_comp_path::ComponentPath
    src_var_name::Symbol
    dst_comp_path::ComponentPath
    dst_par_name::Symbol
    ignoreunits::Bool
    backup::Union{Symbol, Nothing} # a Symbol identifying the external param providing backup data, or nothing
    offset::Int

    function InternalParameterConnection(src_path::ComponentPath, src_var::Symbol,
                                         dst_path::ComponentPath, dst_par::Symbol,
                                         ignoreunits::Bool, backup::Union{Symbol, Nothing}=nothing; offset::Int=0)
        self = new(src_path, src_var, dst_path, dst_par, ignoreunits, backup, offset)
        return self
    end
end

struct ExternalParameterConnection  <: AbstractConnection
    comp_path::ComponentPath
    param_name::Symbol      # name of the parameter in the component
    external_param::Symbol  # name of the parameter stored in external_params
end

# Converts symbol to component path
function ExternalParameterConnection(comp_name::Symbol, param_name::Symbol, external_param::Symbol)
    return ExternalParameterConnection(ComponentPath(comp_name), param_name, external_param)
end

Base.pathof(obj::ExternalParameterConnection) = obj.comp_path
Base.nameof(obj::ExternalParameterConnection) = obj.param_name


text:
# Test Chrome_result.jl

#include("../src/CGP.jl")
using Main.CGP
MyInt = UInt8
include("../src/Chrome_result.jl")

randgoal(numoutputs) = [rand(collect(0x0:0xf)) for _ in 1:numoutputs]  # Assumes MyInt = UInt8
randgoallist(len,numouputs) = [randgoal(numoutputs) for _ = 1:len]

function mut_evolve(p,funcs)
  c = random_chromosome(p,funcs)
  gl = randgoallist(6,3)
  (worse,same,c,output,matched_goals,matched_goals_list) = mut_evolve_subgoal( c, gl, funcs, 500 )
  (c,gl,output,matched_goals,matched_goals_list)
  @assert output == output_values(c)
end
#goal1 = Main.CGP.MyInt[ 0xb, 0xe ]
#goal2 = Main.CGP.MyInt[ 0x6, 0x9 ]
#goal3 = Main.CGP.MyInt[ 0xa, 0xe ]
#goallist = [goal1, goal2,goal3]

#goallist = [[0x0,0x1],[0x07, 0x01],[0x01, 0x0d],[0x07, 0x0d]]
#output =  [0x07, 0x0d]

numcomponents = 6
routput = rand([0x0,0x1,0x2,0x3],numcomponents)
println("routput: ",routput)
rgoal =   rand([0x0,0x1,0x2,0x3],numcomponents)
println("rgoal: ",rgoal)
cm = components_matched( routput, rgoal )
println("number_components_matched: ",length(cm))
println("components_matched: ",cm)
check = map(x->(routput[cm[x][2]],rgoal[cm[x][3]]), collect(1:length(cm)))
println("check: ",check)
println(map(x->(x[1]==x[2]),check),"  all components should be 1")

numgoals = 4
rgoallist = [rand([0x0,0x1,0x2,0x3],numcomponents) for _=1:numgoals]
gm = goals_matched(routput,rgoallist)
println("goals matched: ",gm)

#=
numinputs = 2
funcs = default_funcs(numinputs)
p = Parameters( numinputs=numinputs, numoutputs=2, numinteriors=6, numlevelsback=6 )
context = construct_context(numinputs)
c = random_chromosome( p, funcs )
max_steps = 20
#result = mut_evolve_exact( random_chromosome( p, funcs ), context, goallist, max_steps )
=#