text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libiup_mglplot_/include/libiup_mglplot_.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libiup_mglplot_ = "iup_mglplot") : (const libiup_mglplot_ = "libiup_mglplot")  # Name of shared lib.

export IupMglPlotOpen, IupMglPlot, IupMglPlotBegin, IupMglPlotAdd1D, IupMglPlotAdd2D, IupMglPlotAdd3D,
IupMglPlotEnd, IupMglPlotNewDataSet, IupMglPlotInsert1D, IupMglPlotInsert2D, IupMglPlotInsert3D,
IupMglPlotSet1D, IupMglPlotSet2D, IupMglPlotSet3D,
IupMglPlotSetFormula, IupMglPlotSetData, IupMglPlotLoadData,
IupMglPlotSetFromFormula, IupMglPlotTransform, IupMglPlotTransformXYZ, IupMglPlotDrawMark,
IupMglPlotDrawLine, IupMglPlotDrawText, IupMglPlotPaintTo

function IupMglPlotOpen()
  ccall( (:IupMglPlotOpen, libiup_mglplot_), Cvoid, (), )
end

function IupMglPlot()
  ccall( (:IupMglPlot, libiup_mglplot_), Ptr{Ihandle}, (), )
end

function IupMglPlotBegin(ih::Ptr{Cint}, dim::Int)
  ccall( (:IupMglPlotBegin, libiup_mglplot_), Cvoid, (Ptr{Cint}, Cint), ih, dim)
end

function IupMglPlotAdd1D(ih::Ptr{Cint}, name::String, y::Cfloat)
  ccall( (:IupMglPlotAdd1D, libiup_mglplot_), Cvoid, (Ptr{Cint}, Ptr{UInt8}, Cfloat), ih, name, y)
end

function IupMglPlotAdd2D(ih::Ptr{Cint}, x::Cfloat, y::Cfloat)
  ccall( (:IupMglPlotAdd2D, libiup_mglplot_), Cvoid, (Ptr{Cint}, Cfloat, Cfloat), ih, x, y)
end

function IupMglPlotAdd3D(ih::Ptr{Cint}, x::Cfloat, y::Cfloat, z::Cfloat)
  ccall( (:IupMglPlotAdd3D, libiup_mglplot_), Cvoid, (Ptr{Cint}, Cfloat, Cfloat, Cfloat), ih, x, y, z)
end

function IupMglPlotEnd(ih::Ptr{Cint})
  ccall( (:IupMglPlotEnd, libiup_mglplot_), Cint, (Ptr{Cint},), ih)
end

function IupMglPlotNewDataSet(ih::Ptr{Ihandle}, dim::Int)
  ccall( (:IupMglPlotNewDataSet, libiup_mglplot_), Cint, (Ptr{Ihandle}, Cint), ih, dim)
end

function IupMglPlotInsert1D(ih::Ptr{Ihandle}, ds_index::Int, sample_index::Int, names::Ptr{Ptr{UInt8}}, y::Ptr{Cfloat}, count::Int)
  ccall( (:IupMglPlotInsert1D, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Ptr{UInt8}}, Ptr{Cfloat}, Cint), ih, ds_index, sample_index, names, y, count)
end

function IupMglPlotInsert2D(ih::Ptr{Ihandle}, ds_index::Int, sample_index::Int, x::Ptr{Cfloat}, y::Ptr{Cfloat}, count::Int)
  ccall( (:IupMglPlotInsert2D, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Cint), ih, ds_index, sample_index, x, y, count)
end

function IupMglPlotInsert3D(ih::Ptr{Ihandle}, ds_index::Int, sample_index::Int, x::Ptr{Cfloat}, y::Ptr{Cfloat}, z::Ptr{Cfloat}, count::Int)
  ccall( (:IupMglPlotInsert3D, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cint), ih, ds_index, sample_index, x, y, z, count)
end

function IupMglPlotSet1D(ih::Ptr{Ihandle}, ds_index::Int, names::Ptr{Ptr{UInt8}}, y::Ptr{Cfloat}, count::Int)
  ccall( (:IupMglPlotSet1D, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{Ptr{UInt8}}, Ptr{Cfloat}, Cint), ih, ds_index, names, y, count)
end

function IupMglPlotSet2D(ih::Ptr{Ihandle}, ds_index::Int, x::Ptr{Cfloat}, y::Ptr{Cfloat}, count::Int)
  ccall( (:IupMglPlotSet2D, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Cint), ih, ds_index, x, y, count)
end

function IupMglPlotSet3D(ih::Ptr{Ihandle}, ds_index::Int, x::Ptr{Cfloat}, y::Ptr{Cfloat}, z::Ptr{Cfloat}, count::Int)
  ccall( (:IupMglPlotSet3D, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cint), ih, ds_index, x, y, z, count)
end

function IupMglPlotSetFormula(ih::Ptr{Ihandle}, ds_index::Int, formulaX::String, formulaY::String, formulaZ::String, count::Int)
  ccall( (:IupMglPlotSetFormula, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Cint), ih, ds_index, formulaX, formulaY, formulaZ, count)
end

function IupMglPlotSetData(ih::Ptr{Ihandle}, ds_index::Int, data::Ptr{Cfloat}, count_x::Int, count_y::Int, count_z::Int)
  ccall( (:IupMglPlotSetData, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{Cfloat}, Cint, Cint, Cint), ih, ds_index, data, count_x, count_y, count_z)
end

function IupMglPlotLoadData(ih::Ptr{Ihandle}, ds_index::Int, filename::String, count_x::Int, count_y::Int, count_z::Int)
  ccall( (:IupMglPlotLoadData, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{UInt8}, Cint, Cint, Cint), ih, ds_index, filename, count_x, count_y, count_z)
end

function IupMglPlotSetFromFormula(ih::Ptr{Ihandle}, ds_index::Int, formula::String, count_x::Int, count_y::Int, count_z::Int)
  ccall( (:IupMglPlotSetFromFormula, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Ptr{UInt8}, Cint, Cint, Cint), ih, ds_index, formula, count_x, count_y, count_z)
end

function IupMglPlotTransform(ih::Ptr{Ihandle}, x::Cfloat, y::Cfloat, z::Cfloat, ix::Ptr{Cint}, iy::Ptr{Cint})
  ccall( (:IupMglPlotTransform, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cfloat, Cfloat, Cfloat, Ptr{Cint}, Ptr{Cint}), ih, x, y, z, ix, iy)
end

function IupMglPlotTransformXYZ(ih::Ptr{Ihandle}, ix::Int, iy::Int, x::Ptr{Cfloat}, y::Ptr{Cfloat}, z::Ptr{Cfloat})
  ccall( (:IupMglPlotTransformXYZ, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cint, Cint, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}), ih, ix, iy, x, y, z)
end

function IupMglPlotDrawMark(ih::Ptr{Ihandle}, x, y, z)
  ccall( (:IupMglPlotDrawMark, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cfloat, Cfloat, Cfloat), ih, x, y, z)
end

function IupMglPlotDrawLine(ih::Ptr{Ihandle}, x1, y1, z1, x2, y2, z2)
  ccall( (:IupMglPlotDrawLine, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat, Cfloat), ih, x1, y1, z1, x2, y2, z2)
end

function IupMglPlotDrawText(ih::Ptr{Ihandle}, text::String, x, y, z)
  ccall( (:IupMglPlotDrawText, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Ptr{UInt8}, Cfloat, Cfloat, Cfloat), ih, text, x, y, z)
end
# void IupMglPlotPaintTo(Ihandle* ih, const char* format, int w, int h, double dpi, void* data);
function IupMglPlotPaintTo(ih::Ptr{Ihandle}, format::String, w::Int, h::Int, dpi::Int, data::String)
  ccall( (:IupMglPlotPaintTo, libiup_mglplot_), Cvoid, (Ptr{Ihandle}, Ptr{UInt8}, Cint, Cint, Cfloat, Ptr{UInt8}), ih, format, w, h, Cfloat(dpi), data)
end


text:
# push!(LOAD_PATH,"Y:/_raspberry/") # for local run
# push!(LOAD_PATH,"../")

using Documenter, NumDisplay

makedocs(
    sitename = "NumDisplay Documentation",
    modules = [NumDisplay],
    pages = [
        "Home" => "index.md",
        "Methods" => [
            "Direct" => "direct.md",
            "BCD" => "bcd.md",
            "SPI" => "spi.md",
        ],
        "API refs" => "api.md",
    ],
)

deploydocs(
    repo = "github.com/metelkin/NumDisplay.jl.git",
    target = "build",
)

ext: jl


text:
using Documenter, OhMyQSIM

makedocs(sitename="OhMyQSIM",
         pages=Any[
         "Getting Started"    =>    "index.md",
         "Algorithms"         =>    "algorithms.md"
         ])

deploydocs(
    repo = "github.com/nmoran/OhMyQSIM.jl.git",
)


text:
using Documenter, Hop

makedocs()

deploydocs(
    repo = "github.com/mistguy/Hop.jl.git",
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    julia = "1.0"
)


text:
"""
Abstract type representing convex classes of probability distributions ``\\mathcal{G}``.
"""
abstract type ConvexPriorClass end

struct PriorVariable{C<:ConvexPriorClass,V}
    convexclass::C
    finite_param::V
    model::Any
end

Base.broadcastable(prior::PriorVariable) = Ref(prior)

function (priorvariable::PriorVariable)(p::AbstractVector{<:Real})
    convexclass = priorvariable.convexclass
    convexclass(p)
end

function (priorvariable::PriorVariable)()
    priorvariable(JuMP.value.(priorvariable.finite_param))
end

abstract type AbstractSimplexPriorClass <: ConvexPriorClass end
abstract type AbstractMixturePriorClass <: AbstractSimplexPriorClass end

# TODO: implement correct projection onto simplex and check deviation is not too big
# though this should mostly help with minor numerical difficulties
function fix_πs(πs)
    πs = max.(πs, 0.0)
    πs = πs ./ sum(πs)
end

"""
    DiscretePriorClass(support) <: Empirikos.ConvexPriorClass

Type representing the family of all discrete distributions supported on a subset
of `support`, i.e., it represents all `DiscreteNonParametric` distributions with
`support = support` and `probs` taking values on the probability simplex.

Note that `DiscretePriorClass(support)(probs) == DiscreteNonParametric(support, probs)`.

# Examples
```jldoctest
julia> gcal = DiscretePriorClass([0,0.5,1.0])
DiscretePriorClass | support = [0.0, 0.5, 1.0]

julia> gcal([0.2,0.2,0.6])
DiscreteNonParametric{Float64, Float64, Vector{Float64}, Vector{Float64}}(support=[0.0, 0.5, 1.0], p=[0.2, 0.2, 0.6])
```
"""
struct DiscretePriorClass{S} <: AbstractSimplexPriorClass
    support::S
end

DiscretePriorClass(; support=DataBasedDefault()) = DiscretePriorClass(support)

support(convexclass::DiscretePriorClass) = convexclass.support

function (convexclass::DiscretePriorClass)(p::AbstractVector{<:Real})
    DiscreteNonParametric(support(convexclass), fix_πs(p))
end

nparams(convexclass::DiscretePriorClass) = length(support(convexclass))
Distributions.components(convexclass::DiscretePriorClass) = Dirac.(support(convexclass))

function Base.show(io::IO, gcal::DiscretePriorClass)
    print(io, "DiscretePriorClass | support = ")
    show(IOContext(io, :compact => true), support(gcal))
end

function prior_variable!(model, convexclass::AbstractSimplexPriorClass)
    n = nparams(convexclass)
    tmp_vars = @variable(model, [i = 1:n])
    set_lower_bound.(tmp_vars, 0.0)
    con = @constraint(model, sum(tmp_vars) == 1.0)
    PriorVariable(convexclass, tmp_vars, model)
end

# Dirac
function pdf(prior::PriorVariable{<:DiscretePriorClass}, Z::EBayesSample)
    @unpack convexclass, finite_param, model = prior
    pdf_combination = likelihood.(Z, support(convexclass))
    @expression(model, dot(finite_param, pdf_combination))
end

# Dirac
function cdf(prior::PriorVariable{<:DiscretePriorClass}, Z::EBayesSample)
    @unpack convexclass, finite_param, model = prior
    cdf_combination = _cdf.(likelihood_distribution.(Z, support(convexclass)), response(Z))
    @expression(model, dot(finite_param, cdf_combination))
end

function (target::LinearEBayesTarget)(prior::PriorVariable{<:DiscretePriorClass})
    @unpack convexclass, finite_param, model = prior
    linear_functional_evals = target.(support(convexclass))
    @expression(model, dot(finite_param, linear_functional_evals))
end

"""
    MixturePriorClass(components) <: Empirikos.ConvexPriorClass

Type representing the family of all mixture distributions with mixing components equal to
`components`, i.e., it represents all `MixtureModel` distributions with
`components = components` and `probs` taking values on the probability simplex.

Note that `MixturePriorClass(components)(probs) == MixtureModel(components, probs)`.

# Examples

```jldoctest
julia> gcal = MixturePriorClass([Normal(0,1), Normal(0,2)])
MixturePriorClass (K = 2)
Normal{Float64}(μ=0.0, σ=1.0)
Normal{Float64}(μ=0.0, σ=2.0)

julia> gcal([0.2,0.8])
MixtureModel{Normal{Float64}}(K = 2)
components[1] (prior = 0.2000): Normal{Float64}(μ=0.0, σ=1.0)
components[2] (prior = 0.8000): Normal{Float64}(μ=0.0, σ=2.0)
```
"""
struct MixturePriorClass{S} <: AbstractMixturePriorClass
    components::S
end

MixturePriorClass() = MixturePriorClass(nothing)
Distributions.components(mixclass::MixturePriorClass) = mixclass.components
function Distributions.component(mixclass::AbstractSimplexPriorClass, i::Integer)
    Distributions.components(mixclass)[i]
end
nparams(mixclass::AbstractMixturePriorClass) = length(components(mixclass))
Distributions.ncomponents(mixclass::AbstractSimplexPriorClass) = nparams(mixclass)
function (mixclass::AbstractMixturePriorClass)(p::AbstractVector{<:Real})
    MixtureModel(components(mixclass), fix_πs(p))
end

# Similar to show method for MixtureModel
function Base.show(io::IO, gcal::MixturePriorClass)
    K = ncomponents(gcal)
    println(io, "MixturePriorClass (K = $K)")
    Ks = min(K, 8)
    for i = 1:Ks
        println(io, Distributions.component(gcal, i))
    end
    if Ks < K
        println(io, "The rest are omitted ...")
    end
end

function pdf(prior::PriorVariable{<:AbstractMixturePriorClass}, Z::EBayesSample)
    @unpack convexclass, finite_param, model = prior
    pdf_combination = pdf.(components(convexclass), Z)
    @expression(model, dot(finite_param, pdf_combination))
end

function cdf(prior::PriorVariable{<:AbstractMixturePriorClass}, Z::EBayesSample)
    @unpack convexclass, finite_param, model = prior
    cdf_combination = cdf.(components(convexclass), Z)
    @expression(model, dot(finite_param, cdf_combination))
end

function (target::LinearEBayesTarget)(prior::PriorVariable{<:AbstractMixturePriorClass})
    @unpack convexclass, finite_param, model = prior
    linear_functional_evals = target.(components(convexclass))
    @expression(model, dot(finite_param, linear_functional_evals))
end

"""
    GaussianScaleMixtureClass(σs) <: Empirikos.ConvexPriorClass

Type representing the family of mixtures of Gaussians with mean `0` and standard deviations
equal to `σs`. `GaussianScaleMixtureClass(σs)` represents the same class of distributions
as `MixturePriorClass.(Normal.(0, σs))`

```jldoctest
julia> gcal = GaussianScaleMixtureClass([1.0,2.0])
GaussianScaleMixtureClass | σs = [1.0, 2.0]

julia> gcal([0.2,0.8])
MixtureModel{Normal{Float64}}(K = 2)
components[1] (prior = 0.2000): Normal{Float64}(μ=0.0, σ=1.0)
components[2] (prior = 0.8000): Normal{Float64}(μ=0.0, σ=2.0)
```
"""
struct GaussianScaleMixtureClass{S} <: AbstractMixturePriorClass
    σs::S
end

GaussianScaleMixtureClass() = GaussianScaleMixtureClass(DataBasedDefault())

function Base.show(io::IO, gcal::GaussianScaleMixtureClass)
    print(io, "GaussianScaleMixtureClass | σs = ")
    show(IOContext(io, :compact => true), gcal.σs)
end

components(convexclass::GaussianScaleMixtureClass) = Normal.(0, convexclass.σs)


text:
export Bellhop

"""
$(TYPEDEF)
A propagation model based on an external FORTRAN Bellhop executable.
"""
struct Bellhop{T} <: PropagationModel{T}
  env::T
  nbeams::Int
  minangle::Float32
  maxangle::Float32
  gaussian::Bool
  debug::Bool
  function Bellhop(env, nbeams, minangle, maxangle, gaussian, debug)
    nbeams < 0 && (nbeams = 0)
    -π/2 ≤ minangle ≤ π/2 || throw(ArgumentError("minangle should be between -π/2 and π/2"))
    -π/2 ≤ maxangle ≤ π/2 || throw(ArgumentError("maxangle should be between -π/2 and π/2"))
    minangle < maxangle || throw(ArgumentError("maxangle should be more than minangle"))
    new{typeof(env)}(check(Bellhop, env), nbeams, Float32(minangle), Float32(maxangle), gaussian, debug)
  end
end

"""
    Bellhop(env; gaussian=false, debug=false)
    Bellhop(env, nbeams, minangle, maxangle, gaussian, debug)

Create a Bellhop propagation model.
"""
Bellhop(env; gaussian=false, debug=false) = Bellhop(env, 0, -80°, 80°, gaussian, debug)

### interface functions

function check(::Type{Bellhop}, env::Union{<:UnderwaterEnvironment,Missing})
  if env === missing
    mktempdir(prefix="bellhop_") do dirname
      try
        bellhop(dirname, false)
      catch e
        e isa BellhopError && e.details == ["Unable to execute bellhop.exe"] && throw(e)
      end
    end
  else
    seabed(env) isa RayleighReflectionCoef || throw(ArgumentError("Seabed type not supported"))
    seasurface(env) === Vacuum || throw(ArgumentError("Only vacuum seasurface supported"))
  end
  env
end

function arrivals(model::Bellhop, tx1::AcousticSource, rx1::AcousticReceiver)
  mktempdir(prefix="bellhop_") do dirname
    writeenv(model, [tx1], [rx1], "A", dirname)
    bellhop(dirname, model.debug)
    readarrivals(joinpath(dirname, "model.arr"))
  end
end

function transfercoef(model::Bellhop, tx1::AcousticSource, rx::AcousticReceiverGrid2D; mode=:coherent)
  if mode === :coherent
    taskcode = "C"
  elseif mode === :incoherent
    taskcode = "I"
  elseif mode === :semicoherent
    taskcode = "S"
  else
    throw(ArgumentError("Unknown mode :" * string(mode)))
  end
  mktempdir(prefix="bellhop_") do dirname
    writeenv(model, [tx1], rx, taskcode, dirname)
    bellhop(dirname, model.debug)
    readshd(joinpath(dirname, "model.shd"))
  end
end

function transfercoef(model::Bellhop, tx1::AcousticSource, rx1::AcousticReceiver; mode=:coherent)
  if mode === :coherent
    taskcode = "C"
  elseif mode === :incoherent
    taskcode = "I"
  elseif mode === :semicoherent
    taskcode = "S"
  else
    throw(ArgumentError("Unknown mode :" * string(mode)))
  end
  mktempdir(prefix="bellhop_") do dirname
    writeenv(model, [tx1], [rx1], taskcode, dirname)
    bellhop(dirname, model.debug)
    readshd(joinpath(dirname, "model.shd"))[1]
  end
end

function eigenrays(model::Bellhop, tx1::AcousticSource, rx1::AcousticReceiver)
  mktempdir(prefix="bellhop_") do dirname
    writeenv(model, [tx1], [rx1], "E", dirname)
    bellhop(dirname, model.debug)
    readrays(joinpath(dirname, "model.ray"))
  end
end

function rays(model::Bellhop, tx1::AcousticSource, θ::AbstractArray, rmax)
  θ isa AbstractRange || length(θ) == 1 || throw(ArgumentError("Bellhop only supports uniformly spaced angles"))
  all(-π/2 .< θ .< π/2) || throw(ArgumentError("θ must be between -π/2 and π/2"))
  mktempdir(prefix="bellhop_") do dirname
    writeenv(model, [tx1], [AcousticReceiver(rmax, 0.0)], "R", dirname; minangle=minimum(θ), maxangle=maximum(θ), nbeams=length(θ))
    bellhop(dirname, model.debug)
    readrays(joinpath(dirname, "model.ray"))
  end
end

rays(model::Bellhop, tx1::AcousticSource, θ, rmax) = rays(model, tx1, [θ], rmax)[1]

### helper functions

struct BellhopError <: Exception
  details::Vector{String}
end

function Base.show(io::IO, e::BellhopError)
  if length(e.details) == 1
    println(io, e.details[1])
  else
    println(io, "Bellhop said:")
    for s ∈ e.details
      println(io, "  ", s)
    end
  end
end

function bellhop(dirname, debug)
  infilebase = joinpath(dirname, "model")
  outfilename = joinpath(dirname, "output.txt")
  try
    run(pipeline(ignorestatus(`bellhop.exe $infilebase`); stdout=outfilename, stderr=outfilename))
    if debug
      @info "Bellhop run completed in $dirname, press ENTER to delete intermediate files..."
      readline()
    end
  catch
    throw(BellhopError(["Unable to execute bellhop.exe"]))
  end
  err = String[]
  checkerr!(err, outfilename)
  checkerr!(err, joinpath(dirname, "model.prt"))
  if length(err) > 0
    throw(BellhopError(err))
  end
end

function checkerr!(err, filename)
  output = false
  open(filename) do f
    for s in eachline(f)
      if output || occursin("ERROR", uppercase(s))
        push!(err, s)
        output = true
      end
    end
  end
end

function writeenv(model::Bellhop, tx::Vector{<:AcousticSource}, rx::AbstractArray{<:AcousticReceiver}, taskcode, dirname; minangle=model.minangle, maxangle=model.maxangle, nbeams=model.nbeams)
  all(location(tx1)[1] == 0.0 for tx1 ∈ tx) || throw(ArgumentError("Bellhop requires transmitters at (0, 0, z)"))
  all(location(tx1)[2] == 0.0 for tx1 ∈ tx) || throw(ArgumentError("Bellhop 2D requires transmitters in the x-z plane"))
  all(location(rx1)[1] >= 0.0 for rx1 ∈ rx) || throw(ArgumentError("Bellhop requires receivers to be in the +x halfspace"))
  all(location(rx1)[2] == 0.0 for rx1 ∈ rx) || throw(ArgumentError("Bellhop 2D requires receivers in the x-z plane"))
  env = model.env
  name = split(basename(dirname), "_")[end]
  filename = joinpath(dirname, "model.env")
  open(filename, "w") do io
    println(io, "'", name, "'")
    flist = [nominalfrequency(tx1) for tx1 ∈ tx]
    f = sum(flist) / length(flist)
    maximum(abs.(flist .- f))/f > 0.2 && @warn("Source frequency varies by more than 20% from nominal frequency")
    @printf(io, "%0.6f\n", f)
    println(io, "1")
    if length(rx) == 1
      maxr = location(rx[1])[1]
    elseif rx isa AcousticReceiverGrid2D
      maxr = maximum(rx.xrange)
    else
      throw(ArgumentError("Receivers must be on a 2D grid"))
    end
    ss = ssp(env)
    sspi = "S"
    ss isa SampledSSP1D && ss.interp === :linear && (sspi = "C")
    print(io, "'", sspi, "VWT")
    alt = altimetry(env)
    if !(alt isa FlatSurface)
      print(io, "*")
      createadfile(joinpath(dirname, "model.ati"), alt, (p...) -> -altitude(p...), maxr, f)
    end
    println(io, "'")
    bathy = bathymetry(env)
    waterdepth = maxdepth(bathy)
    @printf(io, "1 0.0 %0.6f\n", waterdepth)
    if ss isa IsoSSP
      @printf(io, "0.0 %0.6f /\n", soundspeed(ss, 0.0, 0.0, 0.0), )
      @printf(io, "%0.6f %0.6f /\n", waterdepth, soundspeed(ss, 0.0, 0.0, 0.0))
    elseif ss isa SampledSSP1D
      for i ∈ 1:length(ss.z)
        @printf(io, "%0.6f %0.6f /\n", -ss.z[i], ss.c[i])
      end
    else
      for d ∈ range(0.0, waterdepth; length=recommendlength(waterdepth, f))
        @printf(io, "%0.6f %0.6f /\n", d, soundspeed(ss, 0.0, 0.0, -d))
      end
      floor(waterdepth) != waterdepth && @printf(io, "%0.6f %0.6f /\n", waterdepth, soundspeed(ss, 0.0, 0.0, -waterdepth))
    end
    print(io, "'A")
    if !(bathy isa ConstantDepth)
      print(io, "*")
      createadfile(joinpath(dirname, "model.bty"), bathy, depth, maxr, f)
    end
    println(io, "' 0.0") # bottom roughness = 0
    bed = seabed(env)
    c2 = soundspeed(ss, 0.0, 0.0, -waterdepth) * bed.cᵣ
    α = bed.δ * 40π / log(10)      # based on APL-UW TR 9407 (1994), IV-9 equation (4)
    @printf(io, "%0.6f %0.6f 0.0 %0.6f %0.6f /\n", waterdepth, c2, bed.ρᵣ, α)
    printarray(io, [-location(tx1)[3] for tx1 ∈ tx])
    if length(rx) == 1
      printarray(io, [-location(rx[1])[3]])
      printarray(io, [maxr / 1000.0])
    elseif rx isa AcousticReceiverGrid2D
      printarray(io, reverse(-rx.zrange))
      printarray(io, rx.xrange ./ 1000.0)
    end
    println(io, "'", taskcode, model.gaussian ? "B'" : "'")
    @printf(io, "%d\n", nbeams)
    @printf(io, "%0.6f %0.6f /\n", rad2deg(minangle), rad2deg(maxangle))
    @printf(io, "0.0 %0.6f %0.6f\n", 1.01*waterdepth, 1.01 * maxr / 1000.0)
  end
end

function printarray(io, a::AbstractVector)
  println(io, length(a))
  for a1 ∈ a
    @printf(io, "%0.6f ", a1)
  end
  println(io, "/")
end

function recommendlength(x, f)
  # recommendation based on nominal half-wavelength spacing
  λ = 1500.0 / f
  clamp(round(Int, 2x / λ) + 1, 25, 1000)
end

function createadfile(filename, data, func, maxr, f)
  open(filename, "w") do io
    interp = "L"
    if data isa SampledDepth || data isa SampledAltitude
      x = data.x
      data.interp !== :linear && (interp = "C")
    else
      x = range(0.0, maxr; length=recommendlength(maxr, f))
    end
    println(io, "'", interp, "'")
    println(io, length(x))
    for i ∈ 1:length(x)
      @printf(io, "%0.6f %0.6f\n", x[i]/1000.0, func(data, x[i], 0.0))
    end
  end
end

function readrays(filename)
  rays = RayArrival{Float64,Float64}[]
  open(filename, "r") do io
    [readline(io) for i ∈ 1:7]
    while !eof(io)
      s = strip(readline(io))
      length(s) == 0 && break
      aod = parse(Float64, s)
      pts, sb, bb = parse.(Int, split(strip(readline(io)) ,r" +"))
      raypath = Array{NTuple{3,Float64}}(undef, pts)
      for k ∈ 1:pts
        x, d = parse.(Float64, split(strip(readline(io)) ,r" +"))
        raypath[k] = (x, 0.0, -d)
      end
      push!(rays, RayArrival(NaN64, NaN64, sb, bb, -deg2rad(aod), NaN64, raypath))
    end
  end
  rays
end

function readarrivals(filename)
  arrivals = RayArrival{Float64,Missing}[]
  open(filename, "r") do io
    s = strip(readline(io))
    if occursin("2D", s)
      f = parse(Float64, strip(readline(io)))
      v = split(strip(readline(io)) ,r" +")
      n = parse(Int, v[1])
      txdepth = parse.(Float64, v[2:end])
      n == length(txdepth) || error("Wrong number of txdepth entries in arrivals")
      v = split(strip(readline(io)) ,r" +")
      n = parse(Int, v[1])
      rxdepth = parse.(Float64, v[2:end])
      n == length(rxdepth) || error("Wrong number of rxdepth entries in arrivals")
      v = split(strip(readline(io)) ,r" +")
      n = parse(Int, v[1])
      rxrange = parse.(Float64, v[2:end])
      n == length(rxrange) || error("Wrong number of rxrange entries in arrivals")
    else
      v = split(s ,r" +")
      f = parse(Float64, v[1])
      n1, n2, n3 = parse.(Int, v[2:4])
      txdepth = parse.(Float64, split(strip(readline(io)) ,r" +"))
      rxdepth = parse.(Float64, split(strip(readline(io)) ,r" +"))
      rxrange = parse.(Float64, split(strip(readline(io)) ,r" +"))
      n1 == length(txdepth) || error("Wrong number of txdepth entries in arrivals")
      n2 == length(rxdepth) || error("Wrong number of rxdepth entries in arrivals")
      n3 == length(rxrange) || error("Wrong number of rxrange entries in arrivals")
    end
    for j ∈ 1:length(txdepth)
      readline(io)
      for k ∈ 1:length(rxdepth)
        for m ∈ 1:length(rxrange)
          count = parse(Int, strip(readline(io)))
          for n ∈ 1:count
            v = split(strip(readline(io)) ,r" +")
            length(v) == 8 || error("Wrong number of data entries in arrivals")
            A, ph, t, _, aod, aoa = parse.(Float64, v[1:6])
            sb, bb = parse.(Int, v[7:8])
            push!(arrivals, RayArrival(t, A * cis(deg2rad(ph)), sb, bb, -deg2rad(aod), deg2rad(aoa)))
          end
        end
      end
    end
  end
  sort(arrivals; by = a -> a.time)
end

function readshd(filename)
  open(filename, "r") do io
    r = read(io, UInt32)
    seek(io, 4r)
    b = Array{UInt8}(undef, 10)
    read!(io, b)
    strip(String(b)) == "rectilin" || error("Bad shd file format: incorrect ptype")
    seek(io, 8r)
    nfreq = read(io, UInt32)
    nfreq == 1 || error("Bad shd file format: incorrect nfreq")
    nθ = read(io, UInt32)
    nθ == 1 || error("Bad shd file format: incorrect nθ")
    nsx = read(io, UInt32)
    nsy = read(io, UInt32)
    nsd = read(io, UInt32)
    nsd == 1 || error("Bad shd file format: incorrect nsd")
    nrd = read(io, UInt32)
    nrr = read(io, UInt32)
    pressure = Array{ComplexF32}(undef, nrr, nrd)
    for ird ∈ 0:nrd-1
      recnum = 10 + ird
      seek(io, recnum * 4r)
      temp = Array{ComplexF32}(undef, nrr)
      read!(io, temp)
      pressure[:,ird+1] .= -temp    # negative because Bellhop seems to have a 180° phase inversion
    end
    reverse(pressure; dims=2)
  end
end


text:
using FourierGPE

##Set simulation parameters
    L=(16.,16.,16.);
    N=(64,64,64);
    sim = Sim(L,N);
    @unpack_Sim sim;

## Initialse sim
    # parameters
    μ = 25.0;
    γ = 0.05;
    tf = 4/γ;
    Nt = 200;
    t = LinRange(0.,tf,Nt);

## Run sim
    x,y,z = X;
    ψi = randn(N)+im*randn(N);
    ϕi = kspace(ψi,sim);

    @pack_Sim! sim;

## Evolve in k-space
    sol = runsim(sim); # will take a few minutes to run.

# size(xspace(sol[40], sim))

plot_iso(xspace(sol[150], sim))
psi_chaos1 = xspace(sol[10], sim)
psi_chaos2 = xspace(sol[16], sim)
psi_knots1 = xspace(sol[20], sim)
psi_knots2 = xspace(sol[25], sim)
psi_knots3 = xspace(sol[30], sim)
psi_knots4 = xspace(sol[40], sim)
psi_tubes1 = xspace(sol[50], sim)
psi_tubes2 = xspace(sol[90], sim)
psi_ringtube = xspace(sol[100], sim)

@save "box_vorts.jld2" psi_chaos1 psi_chaos2 psi_knots1 psi_knots2 psi_knots3 psi_knots4 psi_tubes1 psi_tubes2 psi_ringtube X

# function dense(i)
#     ψm = xspace(sol[i],sim)
#     density = abs2.(ψm)
#     pmax = maximum(density)
#     return density/pmax
# end

# function densityfilm(Nt,saveto="media/3dquenchiso.gif")
#     scene = Scene()
#     tindex = Node(1)

#     scene = volume(lift(i -> dense(i), tindex), algorithm = :iso,show_axis=false)

#     record(scene, saveto, 1:Nt-10) do i
#         tindex[] = i
#     end
# end

# densityfilm(Nt)


text:
#=
dtft.jl
discrete-time Fourier transform
The main purpose of this code is for testing the NUFFT code.
2019-06-05, Jeff Fessler, University of Michigan
Note: the multi-dimensional version uses collect(Tuple(idx))
that might be inefficient, but for efficiency one should just use NUFFT.
=#

export dtft_init, dtft, dtft_adj

using FFTW: fft
using Random: seed!
#using LinearAlgebra: norm
#using BenchmarkTools: @btime
using Distributed: @sync, @distributed, pmap
using SharedArrays: SharedArray, sdata
using Test: @test
using LinearMapsAA: LinearMapAA

"""
`X = dtft(w, x ; n_shift=?)` 1D DTFT

``X[m] = sum_{n=0}^{N-1} x[n] exp(-i w[m] (n - n_shift)), m=1,…,M``

in
- `w::AbstractVector{<:Real}`	`[M]` frequency locations ("units" radians/sample)
- `x::AbstractVector{<:Number}`	`[N]` 1D signal

option
- `n_shift::Real` often is N/2; default 0

out
- `X::AbstractVector{ComplexF64}`		`[M]` DTFT
"""
function dtft(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	return dtft_loop_n(w, x ; n_shift=n_shift)
end

"""
`X = dtft(w, x ; n_shift=?)` multi-dimensional DTFT (DSFT)

``X[m] = sum_{n=0}^{N-1} x[n] exp(-i w[m,:] (n - n_shift)), m=1,…,M``
where here `n` is a `CartesianIndex`

in
- `w::AbstractMatrix{<:Real}`	`[M,D]` frequency locations ("units" radians/sample)
- `x::AbstractArray{<:Number}`	`[(Nd)]` D-dimensional signal

option
- `n_shift::AbstractVector{<:Real}` often is N/2; default zeros(D)

out
- `X::AbstractVector{ComplexF64}`		`[M]` DTFT
"""
function dtft(w::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Number}
	; n_shift::AbstractVector{<:Real} = zeros(Int, ndims(x)))
	return dtft_loop_n(w, x ; n_shift=n_shift)
end

"""
`d = dtft_init(w, N ; n_shift=?)` for 1D DTFT

in
- `w::AbstractVector{<:Real}`	`[M]` frequency locations ("units" radians/sample)
- `N::Int` size of signal `x`

option
- `n_shift::Real` often is N/2; default 0
- `T::DataType` default `ComplexF64` for testing NUFFT accuracy

out
- `d::NamedTuple` with fields
	`dtft = x -> dtft(x), adjoint = y -> dtft_adj(y), A=LinearMapAA`
"""
function dtft_init(w::AbstractVector{<:Real}, N::Int
		; n_shift::Real = 0,
		T::DataType = ComplexF64)
	M = length(w)
	forw = x -> dtft(w, x ; n_shift=n_shift)
	back = y -> dtft_adj(w, y, N ; n_shift=n_shift)
	A = LinearMapAA(x -> forw(x), y -> back(y), (M, N),
		(name="dtft1", N=(N,)) ; T = ComplexF64)
	return (dtft=forw, adjoint=back, A=A)
end

"""
`d = dtft_init(w, N ; n_shift=?)` for multi-dimensional DTFT (DSFT)

in
- `w::AbstractMatrix{<:Real}`	`[M,D]` frequency locations ("units" radians/sample)
- `N::Dims` `[D]` dimensions of signal `x`

option
- `n_shift::AbstractVector{<:Real}` often is N/2; default zeros(D)
- `T::DataType` default `ComplexF64` for testing NUFFT accuracy

out
- `d::NamedTuple` with fields
	`dtft = x -> dtft(x), adjoint = y -> dtft_adj(y), A=LinearMapAA`
"""
function dtft_init(w::AbstractMatrix{<:Real}, N::Dims
		; n_shift::AbstractVector{<:Real} = zeros(Int, length(N)))
	D = size(w,2)
	length(N) != D && throw(DimensionMismatch("length(N) vs D=$D"))
	length(n_shift) != D && throw(DimensionMismatch("length(n_shift) vs D=$D"))
	M = size(w,1)
	forw = x -> dtft(w, reshape(x, N) ; n_shift=n_shift)
	back = y -> dtft_adj(w, y, N ; n_shift=n_shift)
	A = LinearMapAA(x -> forw(x), y -> back(y)[:], (M, prod(N)),
		(name="dtft$(length(N))", N=N) ; T = ComplexF64)
	return (dtft=forw, adjoint=back, A=A)
end

"""
`x = dtft_adj(w, X, N ; n_shift=?)` adjoint for 1D DTFT

``x[n] = sum_{m=1}^M X[m] exp(i w[m] (n - n_shift)), n=0,…,N-1``

This is the *adjoint* (transpose) of `dtft`, not an *inverse* DTFT.

in
- `w::AbstractVector{<:Real}`	`[M]` frequency locations ("units" radians/sample)
- `X::AbstractVector{ComplexF64}`	`[M]` spectrum values
- `N::Int` size of signal `x`

option
- `n_shift::Real` often is N/2; default 0

out
- `x::AbstractVector{<:Number}`	signal [N]
"""
function dtft_adj(w::AbstractVector{<:Real}, X::AbstractVector{<:Number},
		N::Int ; n_shift::Real = 0)
	M = length(w)
	out = Array{ComplexF64}(undef, N)
	nshift1 = n_shift + 1
	X = ComplexF64.(X)
	for n=1:N # loop over output signal samples
		out[n] = sum(cis.((n - nshift1) * w) .* X)
	end
	return out
end

"""
`x = dtft_adj(w, X, N ; n_shift=?)` adjoint for multi-dimensional DTFT (DSFT)

``x[n] = sum_{m=1}^M X[m] exp(i w[m,:] (n - n_shift)), n=0,…,N-1``
where here `n` is a `CartesianIndex`

in
- `X::AbstractVector{ComplexF64}`		`[M]` DTFT
- `w::AbstractMatrix{<:Real}`	`[M,D]` frequency locations ("units" radians/sample)
- `N::Dims`	`[D]` dimensions of signal `x`

option
- `n_shift::AbstractVector{<:Real}` often is N/2; default zeros(D)

out
- `x::AbstractArray{<:Number}`	`[(N)]` D-dimensional signal
"""
function dtft_adj(w::AbstractMatrix{<:Real}, X::AbstractVector{<:Number},
		N::Dims ; n_shift::AbstractVector{<:Real} = zeros(Int, ndims(x)))

	M,D = size(w)
	nshift1 = n_shift .+ 1

	T = (eltype(w) == Float64) ? ComplexF64 : ComplexF32
	out = Array{T}(undef, N)

	for n=1:prod(N)
		idx = CartesianIndices(N)[n]
		tmp = cis.(w * (collect(Tuple(idx)) - nshift1))
		tmp = transpose(cis.(w * (collect(Tuple(idx)) - nshift1))) * X
		out[idx] = cis.(-w * (collect(Tuple(idx)) - nshift1))' * X
	end
	return out
end

# 1D loop over signal values
function dtft_loop_n(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	N = length(x)
	out = fill(ComplexF64(0), size(w))
	nshift1 = n_shift + 1
	for n=1:N # tried @simd but not faster
		out .+= x[n] * cis.(-(n-nshift1) * w)
	end
	return out
end

# dD loop over signal values
# * `w [M,D]`
# * `x [(Nd)]`
function dtft_loop_n(w::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Number}
		; n_shift::AbstractVector{<:Real} = zeros(Int, size(w,2)))
	M,D = size(w)
	N = size(x)
	length(N) != D && throw(DimensionMismatch("length(N) vs D=$D"))
	length(n_shift) != D && throw(DimensionMismatch("length(n_shift) vs D=$D"))
	ndims(x) != D && throw(DimensionMismatch("ndims(x) vs D=$D"))
	nshift1 = n_shift .+ 1
	T = (eltype(w) == Float64) ? ComplexF64 : ComplexF32
	out = fill(T(0), M)
	for n=1:prod(N)
		idx = CartesianIndices(N)[n]
		out .+= x[idx] * cis.(-w * (collect(Tuple(idx)) - nshift1))
	end
	return out
end

# 1D loop over frequencies
function dtft_loop_m(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	M = length(w)
	N = length(x)
	out = Array{ComplexF64}(undef, M)
	nn = -((0:(N-1)) .- n_shift)
	x = ComplexF64.(x)
	for m=1:M
		out[m] = sum(cis.(nn * w[m]) .* x)
	end
	return out
end

# 1D matrix * vector
function dtft_matvec(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	N = length(x)
	return cis.(w * (-((0:(N-1)) .- n_shift))') * x
end

# 1D single frequency
function dtft_one_w(w::Real, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	N = length(x)
	return sum(cis.(-((0:(N-1)) .- n_shift) * w) .* x)
end

# 1D pmap
function dtft_pmap_m(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	tmp = w -> dtft_one_w(w, x ; n_shift=n_shift)
	return pmap(tmp, w)
end

# 1D distributed loop over frequencies
function dtft_dist_m(w::AbstractVector{<:Real}, x::AbstractVector{<:Number}
		; n_shift::Real = 0)
	M = length(w)
	N = length(x)
	out = SharedArray{ComplexF64}(M)
	nn = -((0:(N-1)) .- n_shift)
	x = ComplexF64.(x)
	@sync @distributed for m=1:M
		out[m] = sum(cis.(nn * w[m]) .* x)
	end
	return sdata(out)
end

# 1D small test data to verify correctness
function dtft_test1( ; N::Int=2^10)
	seed!(0)
	x = randn(ComplexF64, N)
	w = (2*pi) * (0:N-1) / N

	o0 = fft(x)

	o1 = dtft_loop_n(w, x)
	o2 = dtft_loop_m(w, x)
	o3 = dtft_dist_m(w, x)
	o4 = dtft_pmap_m(w, x)
	o5 = dtft_matvec(w, x)
#	@show norm(o1 - o0) / norm(o0)

	@test isapprox(o1, o0)
	@test isapprox(o2, o0)
	@test isapprox(o3, o0)
	@test isapprox(o4, o0)
	@test isapprox(o5, o0)

	true
end

# 2D small test
function dtft_test2( ; N::Dims = (2^3,2^2))
	seed!(0)
	x = randn(ComplexF64, N)
	w1 = (2*pi) * (0:N[1]-1) / N[1]
	w2 = (2*pi) * (0:N[2]-1) / N[2]
	w1 = repeat(w1, 1, N[2])
	w2 = repeat(w2', N[1], 1)
	w = [w1[:] w2[:]]
	d = dtft_init(w, N)

	o0 = fft(x)
#	o1 = dtft_loop_n(w, x)
	o1 = d.dtft(x)
	o1 = reshape(o1, N)
	@test isapprox(o1, o0)
	true
end

# 1D verify consistency
function dtft_test1c( ; N::Int=2^10, M::Int=2^11, n_shift::Real=7)
	seed!(0)
	x = randn(ComplexF64, N)
	w = randn(M)

	o1 = dtft_loop_n(w, x ; n_shift=n_shift)
	o2 = dtft_loop_m(w, x ; n_shift=n_shift)
	o3 = dtft_dist_m(w, x ; n_shift=n_shift)
	o4 = dtft_pmap_m(w, x ; n_shift=n_shift)
	o5 = dtft_matvec(w, x ; n_shift=n_shift)

	@test isapprox(o2, o1)
	@test isapprox(o3, o1)
	@test isapprox(o4, o1)
	@test isapprox(o5, o1)

	d = dtft_init(w, N ; n_shift=n_shift)
	o6 = d.dtft(x)
	@test isapprox(o6, o1)

	o7 = d.A * x
	@test isapprox(o7, o1)

	b1 = dtft_adj(w, o1, N ; n_shift=n_shift)
	b2 = d.adjoint(o1)
	@test isapprox(b2, b1)
	b3 = d.A' * o1
	@test isapprox(b3, b1)

	A = d.A
	@test A.name == "dtft1"
	@test A.N == (N,)

#=
	# time DTFT
	@btime fft($x)
	@btime dtft_loop_n($w, $x)
	@btime dtft_loop_m($w, $x)
	@btime dtft_dist_m($w, $x)
	@btime dtft_pmap_m($w, $x)
	@btime dtft_matvec($w, $x)

	# jf28 results:
	# 12.351 μs (51 allocations: 19.08 KiB)
	# 38.608 ms (6146 allocations: 80.34 MiB)
	# 43.119 ms (2051 allocations: 32.30 MiB)
	# 44.179 ms (11965 allocations: 32.64 MiB)
	# 51.865 ms (20518 allocations: 32.67 MiB)
	# 53.091 ms (7 allocations: 48.03 MiB)

	# conclusion: dtft_loop_n is the fastest
	# (also amenable to parfor=pmap but unhelpful!?)
=#

	true
end

# 2D verify consistency
function dtft_test2c( ;
		M::Int=2^9,
		N::Dims=(2^6,2^4),
		n_shift::AbstractVector{<:Real} = zeros(Int, length(N)))

	seed!(0)
	x = randn(ComplexF64, N)
	w = (rand(M,2) .- 0.5) * 2 * pi

	o1 = dtft_loop_n(w, x ; n_shift=n_shift)
	sd = dtft_init(w, N ; n_shift=n_shift)
	o2 = sd.dtft(x)
	@test isequal(o2, o1)
	o3 = sd.A * x[:]
	@test isequal(o3, o1)

	b1 = dtft_adj(w, o1, N ; n_shift=n_shift)
	b2 = sd.adjoint(o1)
	@test isequal(b2, b1)

	b3 = sd.A' * o1
	b3 = reshape(b3, N)
	@test isequal(b3, b1)

	A = sd.A
	@test A.name == "dtft2"
	@test A.N == N
	true
end

"""
`dtft_test1_adj()`
test adjoint
"""
function dtft_test1_adj( ; N::Int=20, M::Int=30, n_shift::Int=5)
	seed!(0)
	w = randn(M)
	_, _, A = dtft_init(w, N ; n_shift = n_shift)
	isapprox(Matrix(A)', Matrix(A'))
end

"""
`dtft_test2_adj()`
test adjoint
"""
function dtft_test2_adj( ;
		M::Int=2^4,
		N::Dims=(2^3,2^2),
		n_shift::AbstractVector{<:Real} = [6,1.7])
	seed!(0)
	w = (rand(M,2) .- 0.5) * 2 * pi
	sd = dtft_init(w, N ; n_shift=n_shift)
	A = sd.A
	@test A isa LinearMapAA
	isapprox(Matrix(A)', Matrix(A'))
end

"""
`dtft(:test)`
self test
"""
function dtft(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))

	@test dtft_test1()
	@test dtft_test1c()
	@test dtft_test1_adj()

	@test dtft_test2()
	@test dtft_test2c()
	@test dtft_test2_adj()
	true
end


text:
module HearThere

# currently using my hacked version of OSC.jl to enable 1-based indexing and
# the path method. PR has been filed
using OSC
using DataFrames
using PyPlot
using Quaternions
using LightXML
using DataFramesMeta
using Compat

export getrawdata, getcookeddata

# full-scale range in g's (based on 16-bit signed ADC word)
const accelFSR = 16
const accelScale = accelFSR / (2^15)
# full-scale range in degrees / s (based on 16-bit signed ADC word)
const gyroFSR = 2000
const gyroScale = gyroFSR / (2^15)
# full-scale range in μT (based on 14-bit signed ADC word)
const magFSR = 4800
const magScale = magFSR / (2^13)

# magnetic field strength in Cambridge in μT
const magField = 52

"""
Listen on the given port for OSC messages of the form:

    /rawimu iiiiiiiii <accel_x> <accel_y> <accel_z> <gyro_x> <gyro_y> <gyro_z> <mag_x> <mag_y> <mag_z>

and return a DataFrame of the timestamped values.
"""
function getrawdata(duration::Real, oscport::Integer=10001)
    timestamps = Float64[]
    data = [Int32[] for _ in 1:9]
    oscpath = "/rawimu"
    sock = UdpSocket()
    if !bind(sock, ip"0.0.0.0", oscport)
        println("Couldn't bind to port $oscport")
        return makerawdf(timestamps, data)
    end
    try
        # wait for the first message
        println("Waiting for first message...")
        while true
            msg = OscMsg(recv(sock))
            if path(msg) == oscpath
                break
            end
        end
        println("Capturing Data...")
        starttime = time()
        now = time()
        lastreported = 0
        while(now < starttime + duration)
            if(now - lastreported >= 5)
                println("$(round(starttime + duration - now)) seconds left")
                lastreported = now
            end
            msg = OscMsg(recv(sock))
            now = time()
            if path(msg) == oscpath
                push!(timestamps, now)
                for i in 1:9
                    push!(data[i], msg[i])
                end
            end
        end
    finally
        close(sock)
    end
    makerawdf(timestamps, data)
end

"""
Listen on the given port for OSC messages of the form:

    /body_orientation ihffff id timestamp q1 q2 q3 q0
    /body_location  ihfff id timestamp x y z
and
    /ranges ffff anchor0 anchor1 anchor2 anchor3

and return a tuple of DataFrames of the timestamped values.
"""
function getoptidata(duration::Real, oscport::Integer=10001)
    orientationpath = "/body_orientation"
    locationpath = "/body_location"
    orientationdf = DataFrame(
        id=Int32[],
        timestamp=Float64[],
        q0=Float32[],
        q1=Float32[],
        q2=Float32[],
        q3=Float32[])
    locationdf = DataFrame(
        id=Int32[],
        timestamp=Float64[],
        x=Float32[],
        y=Float32[],
        z=Float32[])
    sock = UdpSocket()
    if !bind(sock, ip"0.0.0.0", oscport)
        println("Couldn't bind to port $oscport")
        return (locationdf, orientationdf)
    end
    try
        # wait for the first message
        println("Waiting for first message...")
        while true
            msg = OscMsg(recv(sock))
            if path(msg) in [orientationpath, locationpath]
                break
            else
                println("got unrecognized path: $(path(msg))")
            end
        end
        println("Capturing Data...")
        starttime = time()
        now = time()
        lastreported = 0
        while(now < starttime + duration)
            if(now - lastreported >= 5)
                println("$(round(starttime + duration - now)) seconds left")
                lastreported = now
            end
            msg = OscMsg(recv(sock))
            now = time()
            # opti-track data has an inverted z, so we need to convert
            if path(msg) == orientationpath
                push!(orientationdf, [msg[1], now, msg[6], -msg[3], -msg[4], msg[5]])
            elseif path(msg) == locationpath
                push!(locationdf, [msg[1], now, msg[3], msg[4], -msg[5]])
            else
                println("Got unrecognized path $(path(msg))")
            end
        end
    catch InterruptException
        println("Interrupted.")
    end
    close(sock)
    (locationdf, orientationdf)
end

"""
Listen on the given port for OSC messages of the form:

    /orientation    ffff    q0 q1 q2 q3
    /ranges         ffff     anchor0 anchor1 anchor2 anchor3
    /gpsgeo         ddddd   latitude longitude acceleration horiz_accuracy, vert_accuracy
    /gpsxyz         fff     x y z
    /uwblocalxyz    ffffff  x y z std_x std_y std_z
    /uwbglobalxyz   ffff    x y z confidence
    /fusedxyz       fff     x y z

and return a dictionary of DataFrames of the timestamped values.
"""
function getcookeddata(duration::Real, fileprefix::AbstractString, oscport::Integer=10001)
    dfs = @compat Dict{UTF8String, Any}(
        "/orientation" => DataFrame(
            timestamp=Float64[],
            q0=Float32[], q1=Float32[], q2=Float32[], q3=Float32[]),
        "/ranges" => DataFrame(
            timestamp=Float64[],
            anchor0=Float32[], anchor1=Float32[], anchor2=Float32[], anchor3=Float32[]),
        "/gpsgeo" => DataFrame(
            timestamp=Float64[],
            latitude=Float64[], longitude=Float64[], elevation=Float64[],
            horizAccuracy=Float32[],
            vertAccuracy=Float32[]),
        "/gpsxyz" => DataFrame(
            timestamp=Float64[],
            x=Float32[], y=Float32[], z=Float32[]),
        "/uwblocalxyz" => DataFrame(
            timestamp=Float64[],
            x=Float32[], y=Float32[], z=Float32[],
            stdX=Float32[], stdY=Float32[], stdZ=Float32[]),
        "/uwbglobalxyz" => DataFrame(
            timestamp=Float64[],
            x=Float32[], y=Float32[], z=Float32[],
            confidence=Float32[]),
        "/fusedxyz" => DataFrame(
            timestamp=Float64[],
            x=Float32[], y=Float32[], z=Float32[])
        )

    # returns the filled-in dataframes
    collect_data(dfs, duration, oscport)

    for metric in keys(dfs)
        writetable("$(fileprefix)_$(metric[2:end]).tsv", dfs[metric])
    end

    dfs
end

function collect_data{T<:AbstractString}(dfs::Dict{T, DataFrame}, duration, oscport)
    sock = UdpSocket()
    if !bind(sock, ip"0.0.0.0", oscport)
        println("Couldn't bind to port $oscport")
        return dfs
    end
    try
        # wait for the first message
        println("Waiting for first message...")
        while true
            msg = OscMsg(recv(sock))
            if path(msg) in keys(dfs)
                break
            else
                println("Got unrecognized message path \"$(path(msg))\"")
            end
        end
        println("Capturing Data...")
        starttime = time()
        now = time()
        lastreported = 0
        while(now < starttime + duration)
            if(now - lastreported >= 5)
                println("$(round(starttime + duration - now)) seconds left")
                lastreported = now
            end
            msg = OscMsg(recv(sock))
            now = time()
            msgPath = path(msg)
            df = DataFrame()
            try
                df = dfs[msgPath]
            catch e
                isa(e, KeyError) || rethrow(e)
                println("Got unrecognized message path \"$(path(msg))\"")
                continue
            end

            args = [msg[i] for i in 1:(size(df, 2)-1)]
            row = [now, args...]
            push!(df, [now, args...])
        end
    catch e
        isa(e, InterruptException) || rethrow(e)
        println("Interrupted.")
    finally
        close(sock)
    end

    dfs
end

function rotate(q::Quaternion, v::Vector)
    q = normalize(q)

    imag(q*Quaternion(v)*conj(q))
end

function plotpath(data)
    PyPlot.plot3D(convert(Array, data[:x], NaN),
        convert(Array, data[:y], NaN),
        convert(Array, data[:z], NaN),
        alpha=0.5)
    xlabel("x")
    ylabel("y")
    zlabel("z")
end

function plotorientation(data)
    orientations = Quaternion{Float32}[Quaternion(d...) for d in zip(data[:q0], data[:q1], data[:q2], data[:q3])]
    # this is probably really slow
    path = hcat([rotate(o, Float32[0, 0, 1]) for o in orientations]...)
    PyPlot.plot3D(vec(path[1, :]), vec(path[2, :]), vec(path[3, :]))
    xlim([-1, 1])
    ylim([-1, 1])
    zlim([-1, 1])
    xlabel("x")
    ylabel("y")
    zlabel("z")
end

"""
Takes a set of turntable data collections and returns the x, y, and z
scaling coefficients. Each collection should be with the relevant axis
pointing down, i.e. rotation should be positive. bias_data should be
a measurement when the device is not moving.
"""
function calibrate_gyro_scale(x33, x45, y33, y45, z33, z45, bias_data)
    # convert rpm into gyro ADC units
    units33 = 33 * 360 / 60 / gyroScale
    units45 = 45 * 360 / 60 / gyroScale
    xs = [0, units33, units45]

    # as in y = mx + b
    as = Float64[]
    bs= Float64[]
    x0mean = mean(bias_data[:gyro_x])
    x33mean = mean(x33[:gyro_x])
    x45mean = mean(x45[:gyro_x])
    a, b = linreg(xs, [x0mean, x33mean, x45mean])

    # linreg gives us a, b s.t. x ~= a + b*x', (minimizing the squared error)
    # where x is the measured value and x' is the actual value. Reversing that
    # gives us x' = (x - a) / b. We define our conversion values in terms of addition
    # and multiplication respectively, so we invert them here.
    push!(as, -a)
    push!(bs, 1/b)

    y0mean = mean(bias_data[:gyro_y])
    y33mean = mean(y33[:gyro_y])
    y45mean = mean(y45[:gyro_y])
    a, b = linreg(xs, [y0mean, y33mean, y45mean])
    push!(as, -a)
    push!(bs, 1/b)

    z0mean = mean(bias_data[:gyro_z])
    z33mean = mean(z33[:gyro_z])
    z45mean = mean(z45[:gyro_z])
    a, b = linreg(xs, [z0mean, z33mean, z45mean])
    push!(as, -a)
    push!(bs, 1/b)

    (as, bs)
end

"""
A very simplistic ellipsoid fit of the given data. It just uses the min/max values
of each axis to generate a containing ellipse. It also ignores any rotation.

It returns a tuple of 3-element vectors, the first is the offset α and the 2nd is the
scaling β.

They are meant to be applied as x′ = β(x+α), where x′ is the corrected value
and x is the raw value. We apply the offset first so it can be done in the
integer domain more efficiently.

"""
function calibrate_mag(data)
    xmax = maximum(data[:mag_x])
    xmin = minimum(data[:mag_x])
    ymax = maximum(data[:mag_y])
    ymin = minimum(data[:mag_y])
    zmax = maximum(data[:mag_z])
    zmin = minimum(data[:mag_z])

    α = [
        (xmin + xmax) / -2,
        (ymin + ymax) / -2,
        (zmin + zmax) / -2]
    β = [
         2magField / ((xmax - xmin) * magScale),
         2magField / ((ymax - ymin) * magScale),
         2magField / ((zmax - zmin) * magScale)]

    (α, β)
end

function correct_mag!(data, α, β)
    data[:mag_x] = (data[:mag_x] + α[1]) * β[1]
    data[:mag_y] = (data[:mag_y] + α[2]) * β[2]
    data[:mag_z] = (data[:mag_z] + α[3]) * β[3]
end

function plot_mag(data)
    scatter3D(data[:mag_x], data[:mag_y], data[:mag_z])
    xlabel("x")
    ylabel("y")
    zlabel("z")

    range = maximum(abs([
        minimum(data[:mag_x]),
        minimum(data[:mag_y]),
        minimum(data[:mag_z]),
        maximum(data[:mag_x]),
        maximum(data[:mag_y]),
        maximum(data[:mag_z])]))
    xlim([-range, range])
    ylim([-range, range])
    zlim([-range, range])
end

function plot_gyro(data, bias=[0.0, 0.0, 0.0])
    plot(data[:gyro_x] + bias[1], label="x")
    plot(data[:gyro_y] + bias[2], label="y")
    plot(data[:gyro_z] + bias[3], label="z")
end

"""
Make a dataframe from the timestamps and collected data
"""
function makerawdf(timestamps, data)
    DataFrame(
        timestamp=timestamps,
        accel_x=data[1],
        accel_y=data[2],
        accel_z=data[3],
        gyro_x=data[4],
        gyro_y=data[5],
        gyro_z=data[6],
        mag_x=data[7],
        mag_y=data[8],
        mag_z=data[9])
end

function parseoptitrack(filename)
    bodies = Dict{Int, AbstractString}()
    rootElement = root(parse_file(filename))
    for bodyDesc in get_elements_by_tagname(rootElement, "RigidBodyDesc")
        name = content(find_element(bodyDesc, "Name"))
        id = int(content(find_element(bodyDesc, "ID")))
        bodies[id] = name
    end

    df = DataFrame(
        timestamp=Float64[],
        id=Int[],
        x=Float64[],
        y=Float64[],
        z=Float64[],
        q0=Float64[],
        q1=Float64[],
        q2=Float64[],
        q3=Float64[]
    )
    for dfElement in get_elements_by_tagname(rootElement, "DataFrame")
        # the Latency field apparently means the capture computer's timestamp. See
        # https://forums.naturalpoint.com/viewtopic.php?p=57338
        timestamp = float(content(find_element(dfElement, "Latency")))
        bodydataElement = find_element(dfElement, "RigidBodyData")
        # the rigid body properties aren't contained inside an XML tag, they're all
        # just repeated, so we iterate through the zipped list
        # we're also switching qz and qy as described here:
        # https://forums.naturalpoint.com/viewtopic.php?p=42839
        for (id, x, y, z, q0, q1, q2, q3) in zip(
                get_elements_by_tagname(bodydataElement, "RigidBodyID"),
                get_elements_by_tagname(bodydataElement, "x"),
                get_elements_by_tagname(bodydataElement, "y"),
                get_elements_by_tagname(bodydataElement, "z"),
                get_elements_by_tagname(bodydataElement, "qw"),
                get_elements_by_tagname(bodydataElement, "qx"),
                get_elements_by_tagname(bodydataElement, "qz"),
                get_elements_by_tagname(bodydataElement, "qy"))
            # the optitrack data is 1/2 what it should be because of a calibration
            # error. Also the Z-axis is reversed
            push!(df, [
                timestamp,
                int(content(id)),
                float(content(x)) * 2,
                float(content(y)) * 2,
                -float(content(z)) * 2,
                float(content(q0)),
                float(content(q1)),
                float(content(q2)),
                -float(content(q3))])
        end
    end
    # The OptiTrack reports a location of 0, 0, 0 if the trackable isn't present
    @byrow! df begin
        if :x == 0.0 && :y == 0.0 && :z == 0.0
            :x = NA
            :y = NA
            :z = NA
            :q0 = NA
            :q1 = NA
            :q2 = NA
            :q3 = NA
        end
    end
    df
end

"""
Calculate the ranges for the given path (with x, y, z columns) and add them
to the dataframe
"""
function addranges!(optipath, anchors)
    for (anc, i) in [(:anchor0, 1), (:anchor1, 2), (:anchor2, 3), (:anchor3, 4)]
        optipath[anc] = sqrt(
            (anchors[:x][i] - optipath[:x]).^2 +
            (anchors[:y][i] - optipath[:y]).^2 +
            (anchors[:z][i] - optipath[:z]).^2)
    end
    optipath
end

"""
Takes 2 points (a1 and a2) and two ranges to those points (r1 and r2) and finds
the 0, 1 or 2 points that satisfy the ranges. The are returned as a 2xN matrix,
where N is the number of solutions.
"""
function trilaterate(a1, a2, r1, r2)
    # we first solve in a coordinate system where a1 is at the origin
    x2 = a2[1] - a1[1]
    y2 = a2[2] - a1[2]

    α = r1^2 - r2^2 + x2^2
    a = 1 + x2^2/y2^2
    b = -(α*x2+x2*y2^2)/y2^2
    c = (α^2 + 2α*y2^2+y2^4)/4y2^2 - r1^2

    det = b^2-4a*c
    if det < 0
        sol = Array(Float64, 2, 0)
    elseif det == 0.0
        x0 = -b/2a
        y0 = (r1^2-r2^2+x2^2-2x2*x0+y2^2)/2y2
        sol = [x0, y0]
    else
        # possibly 2 solutions to quadratic
        x0 = [(-b + √det)/2a, (-b - √det)/2a]
        y0 = (r1^2-r2^2+x2^2-2x2*x0+y2^2)/2y2
        sol = vcat(x0', y0')
    end
    # now re-apply the offset. Note this will broadcast in the 2x2 case
    sol .+ a1
end

function cleanranges(ranges)
    for anc in [:anchor0, :anchor1, :anchor2, :anchor3]
        ranges[ranges[anc] .== -1.0, anc] = NA
    end
    ranges
end

function getanchorlocations(rangefile, locationfile)
    anchorDF = DataFrame(
        anchor=0:3,
        x=zeros(4),
        y=zeros(4),
        z=zeros(4),
    )

    locdata = readtable(locationfile)
    counts = zeros(Int, 4)
    for (id, x, y, z) in zip(locdata[:anchor], locdata[:x], locdata[:y], locdata[:z])
        # the values are doubled because the optitrack was calibrated incorrectly
        # so all the reported values are half what they should be. We also
        # divide by 1000 to convert mm to m. Also the Z axis is reversed
        anchorDF[id+1, :x] += 2x / 1000
        anchorDF[id+1, :y] += 2y / 1000
        anchorDF[id+1, :z] += -2z / 1000
        counts[id+1] += 1
    end

    anchorDF[:x][3:4] = anchorDF[:x][3:4] ./ counts[3:4]
    anchorDF[:y][3:4] = anchorDF[:y][3:4] ./ counts[3:4]
    anchorDF[:z][3:4] = anchorDF[:z][3:4] ./ counts[3:4]

    ranges = readtable(rangefile)
    r20 = mean(dropna(ranges[:r20]))
    r21 = mean(dropna(ranges[:r21]))
    r30 = mean(dropna(ranges[:r30]))
    r31 = mean(dropna(ranges[:r31]))

    # we start out considering anchor 3 (index 4) to be the origin
    a0 = trilaterate(
        [anchorDF[4, :x], anchorDF[4, :z]],
        [anchorDF[3, :x], anchorDF[3, :z]],
        r30, r20)
    @assert size(a0) == (2, 2)
    a1 = trilaterate(
        [anchorDF[4, :x], anchorDF[4, :z]],
        [anchorDF[3, :x], anchorDF[3, :z]],
        r31, r21)
    @assert size(a1) == (2, 2)

    # we should get 2 solutions and we happend to know that the correct one is
    # in front of (more positive Z) a3 and a2
    for i in 1:2
        if a0[2, i] > anchorDF[4, :z]
            anchorDF[1, :x] = a0[1, i]
            anchorDF[1, :y] = anchorDF[4, :y]
            anchorDF[1, :z] = a0[2, i]
        end
        if a1[2, i] > anchorDF[3, :z]
            anchorDF[2, :x] = a1[1, i]
            anchorDF[2, :y] = anchorDF[3, :y]
            anchorDF[2, :z] = a1[2, i]
        end
    end

    anchorDF
end

"""
Takes a dataframe of (possibly irregularly) timestamped data and resamples to the given resolution as
if the data had been sampled at a regular interval. Returns an array
of the resampled data. This assumes the data is sorted by timestamp and timestamp.
is in seconds.
"""
function samplereg(data::DataFrame, resolution::AbstractFloat, timefield=:timestamp)
    colnames = filter(n -> n != timefield, names(data))
    endidx = length(data[timefield])
    out = DataFrame()
    out[timefield] = data[timefield][1]:resolution:data[timefield][endidx]
    n = length(out[timefield])
    for name in colnames
        out[name] = Array(Float64, n)
    end

    quat = false
    quatcols = [:q0, :q1, :q2, :q3]
    # we treat quaternions columns specially if present
    if all([col in colnames for col in quatcols])
        quat = true
        # remove the quaternions columns from the regular list
        colnames = filter(n -> !(n in quatcols), colnames)
        for col in quatcols
            out[col] = Array(Float64, n)
        end
    end

    # set up the current interval start/end times
    intstart = data[timefield][1]
    intend = data[timefield][2]

    idx = 1 # input index
    odx = 1 # output index
    while odx <= n
        sampletime = out[odx, timefield]
        # find the data timestamps that bracket our current sample time
        while sampletime > intend
            idx += 1
            intstart = intend
            intend = data[timefield][idx+1]
        end
        if intend == intstart
            α = 0.0
        else
            α = (sampletime - intstart) / (intend - intstart)
        end
        for name in colnames
            if isna(data[idx, name]) || isna(data[idx+1, name])
                out[odx, name] = NA
            else
                out[odx, name] = lerp(data[idx, name], data[idx+1, name], α)
            end
        end
        if quat
            if any([isna(data[idx+off, c]) for off in [0, 1], c in quatcols])
                for col in quatcols
                    out[odx, col] = NA
                end
            else
                p = Quaternion([data[idx, col] for col in quatcols]...)
                q = Quaternion([data[idx+1, col] for col in quatcols]...)
                lerped = nlerp(p, q, α)
                out[odx, :q0] = lerped.s
                out[odx, :q1] = lerped.v1
                out[odx, :q2] = lerped.v2
                out[odx, :q3] = lerped.v3
            end
        end
        odx += 1
    end
    out
end

Quaternion([1, 0, 0, 0]...)
# HearThere.lerp(10, 15, 1)
# p = Quaternion(1, 0, 0, 0.0)
# q = Quaternion(0, 1, 0, 0.0)
# HearThere.nlerp(p, q, 0.99)

"""
returns the lag (in samples) that vector a needs to be delayed relative to
vector b to maximize their cross-correlation, or alternatively the number of
sample of vector b that need to be removed to align the vectors. Assumes the
vectors are at the same sample rate.
"""
function getlag(a, b)
    xc = xcorr(a - mean(a), b - mean(b))
    _, i = findmax(xc)

    max(length(a), length(b)) - i
end

"""
Takes two DataFrames at the same sample rate and returns new ones that are
time-aligned and the same length. The intersection of the column names
excluding the time field are used for the cross-correlation if no offset
is given.
"""
function align(a, b; offset=nothing, timefield=:timestamp)
    if offset == nothing
        acolnames = filter(n -> n != timefield, names(a))
        bcolnames = filter(n -> n != timefield, names(b))
        colnames = intersect(acolnames, bcolnames)

        lagsum = 0
        for colname in colnames
            lagsum += getlag(convert(Array, a[colname], 0.0), convert(Array, b[colname], 0.0))
        end
        avglag = int(round(lagsum / length(colnames)))
    else
        avglag = offset
    end
    if avglag >= 0
        # b needs to be trimmed at the beginning
        bstart = avglag + 1
        astart = 1
    else
        println("WARNING: negative lag ($avglag) untested")
        # a needs to be trimmed at the beginning
        astart = -avglag + 1
        bstart = 1
    end

    n = min(size(a[astart:end, :], 1), size(b[bstart:end, :], 1))

    a = a[astart:astart+n-1, :]
    b = b[bstart:bstart+n-1, :]

    # we don't know what the correct times are anyways, so align them starting
    # at zero
    a[timefield] -= a[1, timefield]
    b[timefield] = a[timefield]
    (a, b)
end

function combine_range_calib{T <: Real, S <: AbstractString}(
    distances::Array{T},
    measurementfiles::Array{S},
    undobias=false)
    # these are the horizontal offsets of the anchors during calibration in meters
    offsets = @compat Dict(
        :anchor0 => -0.165,
        :anchor1 => -0.055,
        :anchor2 => 0.055,
        :anchor3 => 0.165
    )
    df = DataFrame(
        timestamp=Float64[],
        anchor=UTF8String[],
        actual=Float64[],
        measured=Float64[]
    )
    for (r, f) in zip(distances, measurementfiles)
        measuredf = readtable(f)
        for anc in [:anchor0, :anchor1, :anchor2, :anchor3]
            actual = sqrt(offsets[anc]^2 + r^2)
            anchordf = DataFrame(
                timestamp=measuredf[:timestamp],
                anchor=string(anc),
                actual=ones(size(measuredf, 1)) * actual,
                measured=measuredf[anc]
            )
            @byrow! anchordf begin
                if :measured == -1.0 || :measured > 1000
                    :measured = NA
                elseif undobias
                    # reverse the bias that the headtracker is doing internally
                    :measured = reversebias(:measured)
                end
            end
            df = vcat(df, anchordf)
        end
    end

    df
end

"""
Takes 2 dataframes and returns a new dataframe with the errors between
matching columns
"""
function geterr(a, b, timefield=:timestamp)
    acolnames = filter(n -> n != timefield, names(a))
    bcolnames = filter(n -> n != timefield, names(b))
    colnames = intersect(acolnames, bcolnames)

    errdf = DataFrame()
    for colname in colnames
        errdf[colname] = a[colname] - b[colname]
    end
    errdf
end

# this is ported over from the decawave code and is used to correct the range
# measurements
function getrangebias(range)
    # bias table from the decawave driver, used in the tag code to correct the distance
    biastable = [ 1, 2, 4, 5, 6, 8, 9, 10, 12, 13, 15, 18, 20, 22, 24, 27,
                  29, 32, 35, 38, 41, 44, 47, 51, 55, 58, 62, 66, 71, 78, 85,
                  96, 111, 135, 194, 240, 255 ]
    biasoff = -23

    range25 = range * 4 # in units of 25cm
    if range25 > 255
        range25 = 255
    end

    i = 1
    while range25 > biastable[i]
        i += 1
    end
    bias = i - 1 + biasoff # bias is in cm
    bias * 0.01 # convert to m
end

function adjustrange(range)
    range - getrangebias(range)
end

"""
Try to reverse the bias applied in the tag. This can have an error of 1cm
"""
function reversebias(range)
    range + getrangebias(range)
end

"""
Convert a quaternion to the Aerospace Euler Angle representation (Z-Y'-X'').
This representation considers X to be forward (towards the nose), Y is to the
right, and Z is down. This sequence is also known as Yaw-Pitch-Roll. Note that
the primes indicate that the rotations are about the object's local coordinate
system, which is transformed at each step. In the global coordinate system
the same rotation is achieved by rotating in the reverse order, X-Y-Z. Also
note that rotations follow the right-hand rule.

See Kuipers "Quaternions and Rotation Sequences" (1999) for more info.
"""
function aerospace(q::Quaternion)
    # TODO: handle gimbal lock situations
    q = normalize(q)
    # test = 2(q.v1*q.v3+q.s*q.v2)
    test = 2(-q.v1*q.v3+q.s*q.v2)
    if test > 0.998
        ϕ = 0;
        θ = pi/2
        ψ = 2 * atan2(q.v1, q.s);
    elseif test < -0.998
        ϕ = 0;
        θ = -pi/2
        ψ = -2 * atan2(q.v1, q.s);
    else
        ϕ = -atan2(2(q.s*q.v1 + q.v2*q.v3), 1 - 2(q.v1^2+q.v2^2))
        θ = asin(2(q.s*q.v2 - q.v3*q.v1))
        ψ = -atan2(2(q.s*q.v3 + q.v1*q.v2), 1 - 2(q.v2^2+q.v3^2))
    end

    # yaw, pitch, roll
    (ψ, θ, ϕ)
end

function addaerospace!(df::DataFrame)
    yaw = DataArray(Float64, 0)
    pitch = DataArray(Float64, 0)
    roll = DataArray(Float64, 0)
    for i in 1:size(df, 1)
        # convert to aerospace axis convention
        q0 = df[i, :q0]
        q1 = -df[i, :q3]
        q2 = -df[i, :q1]
        q3 = df[i, :q2]

        if isna(q0) || isna(q1) || isna(q2) || isna(q3)
            y = NA
            p = NA
            r = NA
        else
            y, p, r = aerospace(Quaternion(q0, q1, q2, q3, true))
        end
        push!(yaw, y)
        push!(pitch, p)
        push!(roll, r)
    end
    df[:yaw] = yaw
    df[:pitch] = pitch
    df[:roll] = roll

    nothing
end

"""
Linearly interpolate between two values.
"""
lerp(lhs, rhs, t::Real) = (1-t)*lhs + t*rhs

"""
Dot product of two quaternions, as if they're vectors in 4D
"""
dot(p::Quaternion, q::Quaternion) = p.s*q.s + p.v1*q.v1 + p.v2*q.v2 + p.v3*q.v3

"""
Interpolate between two quaternions by linearly interpolating their
components and ensuring it's still a valid rotation by normalizing. This is
still guaranteed to be on the great-circle path between the two quaternions,
(like slerp), but will not travel at a constant rotational velocity with respect
to t.
"""
function nlerp(p::Quaternion, q::Quaternion, t::Real)
    if dot(p, q) < 0
        # negating the quaternion doesn't change the rotation it represents
        # but this ensures we take the short way around instead of the long way
        q = -q
    end
    res = Quaternion(
        lerp(p.s, q.s, t),
        lerp(p.v1, q.v1, t),
        lerp(p.v2, q.v2, t),
        lerp(p.v3, q.v3, t))

    normalize(res)
end

"""
Calculates the error between two unit quaternions representing rotations.
Uses metric Φ6 from Huynh's "Metrics for 3D Rotations: Comparison and Analysis",
which has units of radians.
"""
err(p::Quaternion, q::Quaternion) = 2acos(abs(dot(p, q)))

"""
Calculates the quaternion error for the given dataframes. Assumes they're the
same length
"""
function err(df1::DataFrame, df2::DataFrame)
    @assert(size(df1, 1) == size(df2, 1))
    len = size(df1, 1)
    er = DataFrame(timestamp=Array(Float64, len), err=Array(Float64, len))
    quatcols = [:q0, :q1, :q2, :q3]
    for i in 1:len
        er[i, :timestamp] = df1[i, :timestamp]
        if any(Bool[isna(df[i, col]) for df in DataFrame[df1, df2], col in quatcols])
            er[i, :err] = NA
        else
            p = Quaternion([df1[i, col] for col in quatcols]...)
            q = Quaternion([df2[i, col] for col in quatcols]...)
            er[i, :err] = err(p, q)
        end
    end
    er
end

end # module


text:
"""
from Arb docs:

An arb_t represents a ball over the real numbers, that is, an interval [m±r]≡[m−r,m+r]
where the midpoint m and the radius r are (extended) real numbers and r is nonnegative (possibly infinite).
The result of an (approximate) operation done on arb_t variables is a ball which contains
the result of the (mathematically exact) operation applied to any choice of points in the input balls.
In general, the output ball is not the smallest possible.

The precision parameter passed to each function roughly indicates the precision to which calculations
on the midpoint are carried out (operations on the radius are always done using a fixed, small precision.)
""" ArbReal

# ArbReal structs hold balls given as midpoint, radius

mutable struct ArbReal{P} <: Real     # P is the precision in bits
                    #      midpoint
    mid_exp::Int    # fmpz         exponent of 2 (2^exp)
    mid_size::UInt  # mp_size_t    nwords and sign (lsb holds sign of significand)
    mid_d1::UInt    # significand  unsigned, immediate value or the initial span
    mid_d2::UInt    #   (d1, d2)   the final part indicating the significand, or 0
                    #      radius
    rad_exp::Int    # fmpz       exponent of 2 (2^exp)
    rad_man::UInt   # mp_limb_t  radius, unsigned by definition

    function ArbReal{P}() where {P}
        z = new{P}()
        ccall(@libarb(arb_init), Cvoid, (Ref{ArbReal},), z)
        finalizer(arb_clear, z)
        return z
    end
end

arb_clear(x::ArbReal{P}) where {P} = ccall(@libarb(arb_clear), Cvoid, (Ref{ArbReal},), x)

ArbReal{P}(x::ArbReal{P}) where {P} = x
ArbReal(x::ArbReal{P}) where {P} = x

ArbReal{P}(x::Missing) where {P} = missing
ArbReal(x::Missing) = missing

@inline sign_bit(x::ArbReal{P}) where {P} = isodd(x.mid_size)

ArbReal(x, prec::Int) = prec>=MINIMUM_PRECISION ? ArbReal{workingbits(prec)}(x) : throw(DomainError("bit precision $prec < $MINIMUM_PRECISION"))

function ArbReal(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number}
    if base === 10
        digits = digits > 0 ? bits4digits(digits) : (bits > 0 ? bits : DEFAULT_PRECISION.x)
    elseif base === 2
        digits = bits > 0 ? bits : (digits > 0 ? digits : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbReal(x, digits)
end

swap(x::ArbReal{P}, y::ArbReal{P}) where {P} = ccall(@libarb(arb_swap), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, y)

function copy(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function ArbReal{P}(x::Int64) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_si), Cvoid, (Ref{ArbReal}, Clong), z, x)
    return z
end
ArbReal{P}(x::T) where {P, T<:Union{Int8, Int16, Int32}} = ArbReal{P}(Int64(x))

function ArbReal{P}(x::UInt64) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_ui), Cvoid, (Ref{ArbReal}, Culong), z, x)
    return z
end
ArbReal{P}(x::T) where {P, T<:Union{UInt8, UInt16, UInt32}} = ArbReal{P}(UInt64(x))

function ArbReal{P}(x::Float64) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_d), Cvoid, (Ref{ArbReal}, Cdouble), z, x)
    return z
end
ArbReal{P}(x::T) where {P, T<:Union{Float16, Float32}} = ArbReal{P}(Float64(x))

function ArbReal{P}(x::BigFloat) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_interval_mpfr), Cvoid, (Ref{ArbReal}, Ref{BigFloat}, Ref{BigFloat}, Clong), z, x, x, P)
    return z
end
ArbReal{P}(x::BigInt) where {P} = ArbReal{P}(BigFloat(x))
ArbReal{P}(x::Rational{T}) where {P, T<:Signed} = ArbReal{P}(BigFloat(x))

BigFloat(x::ArbReal{P}) where {P} = BigInt(trunc(BigFloat(ArbFloat{P}(x))))
BigInt(x::ArbReal{P}) where {P} = BigInt(trunc(BigFloat(ArbFloat{P}(x))))

function Base.Integer(x::ArbReal{P}) where {P}
    if isinteger(x)
       abs(x) <= typemax(Int64) ? Int64(x) : BigInt(x)
    else
       throw(InexactError("$x"))
    end
end

function ArbReal{P}(x::Irrational{S}) where {P,S}
    y = ArbFloat{P}(x)
    z = ArbReal{P}(y)
    return z
end

function ArbReal{P}(x::ArbFloat{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, x)
    return z
end

for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil))
    @eval begin
        function $F(x::ArbReal{P}) where {P}
            z = ArbFloat{P}()
            ccall(@libarb($A), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
            return ArbReal{P}(z)
        end
        $F(::Type{T}, x::ArbReal{P}) where {P,T} = T($F(x))
    end
end

trunc(x::ArbReal{P}) where {P} = signbit(x) ? lowerbound(x) : upperbound(x)
trunc(::Type{T}, x::ArbReal{P}) where {P, T} = T(trunc(x))

function midpoint(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function radius(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), z, x)
    return z
end

function midpoint_byref(x::ArbReal{P}) where {P}
    w = ArbReal{P}()
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), w, x)
    z = ArbFloat{P}(w)
    return z
end

function radius_byref(x::ArbReal{P}) where {P}
    w = ArbReal{P}()
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), w, x)
    z = ArbFloat{P}(w)
    return z
end

midpoint(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P} = midpoint_byref(x)
radius(x::ArbReal{P}, ::Type{Mag}) where {P} = radius_byref(x)

function radius(x::ArbReal{P}, ::Type{ArbFloat{P}}) where {P}
    z = ArbFloat{P}()
    mag = radius_byref(x)
    ccall(@libarb(arf_set_mag), Cvoid, (Ref{ArbFloat}, Ref{Mag}), z, mag)
    return z
end


text:
#### EDMF model
using DocStringExtensions

#### Entrainment-Detrainment model

"""
    EntrainmentDetrainment

An Entrainment-Detrainment model for EDMF, containing
all related model and free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct EntrainmentDetrainment{FT <: AbstractFloat}
    "Entrainment TKE scale"
    c_λ::FT = 0.3
    "Entrainment factor"
    c_ε::FT = 0.13
    "Detrainment factor"
    c_δ::FT = 0.52
    "Turbulent Entrainment factor"
    c_t::FT = 0.1
    "Detrainment RH power"
    β::FT = 2
    "Logistic function scale ‵[1/s]‵"
    μ_0::FT = 0.0004
    "Updraft mixing fraction"
    χ::FT = 0.25
    "Minimum updraft velocity"
    w_min::FT = 0.1
    "Exponential area limiter scale"
    lim_ϵ::FT = 0.0001
    "Exponential area limiter amplitude"
    lim_amp::FT = 10
    "Minimum value for turb entr"
    εt_min::FT = 1e-4
end

"""
    SubdomainModel

A subdomain model for EDMF, containing
all related model and free parameters.

TODO: `a_max` is valid for all subdomains,
    but it is insufficient to ensure `a_en`
    is not negative. Limits can be imposed
    for updrafts, but this is a limit
    is dictated by 1 - Σᵢ aᵢ, which must be
    somehow satisfied by regularizing prognostic
    source terms.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct SubdomainModel{FT <: AbstractFloat}
    "Minimum area fraction for any subdomain"
    a_min::FT
    "Maximum area fraction for any subdomain"
    a_max::FT
end

function SubdomainModel(
    ::Type{FT},
    N_up;
    a_min::FT = 0.001,
    a_max::FT = 1 - N_up * a_min,
) where {FT}
    return SubdomainModel(; a_min = a_min, a_max = a_max)
end

"""
    SurfaceModel

A surface model for EDMF, containing all boundary
values and parameters needed by the model.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct SurfaceModel{FT <: AbstractFloat, SV}
    "Temperature ‵[k]‵"
    T::FT = 300.4
    "Liquid water potential temperature ‵[k]‵"
    θ_liq::FT = 299.1
    "Specific humidity ‵[kg/kg]‵"
    q_tot::FT = 22.45e-3
    "Sensible heat flux ‵[w/m^2]‵"
    shf::FT = 9.5
    "Latent heat flux ‵[w/m^2]‵"
    lhf::FT = 147.2
    "Area"
    a::FT
    "Scalar coefficient"
    scalar_coeff::SV = 0
    "Friction velocity"
    ustar::FT = 0.28
    "Monin - Obukhov length"
    obukhov_length::FT = -100
    "Surface covariance stability coefficient"
    ψϕ_stab::FT = 8.3
    "Square ratio of rms turbulent velocity to friction velocity"
    κ_star²::FT = 3.75
    "Height of the lowest level"
    zLL::FT = 60
end

"""
    SurfaceModel{FT}(N_up) where {FT}

Constructor for `SurfaceModel` for EDMF, given:
 - `N_up`, the number of updrafts
"""
function SurfaceModel{FT}(N_up;) where {FT}
    a_surf::FT = 0.1

    surface_scalar_coeff = SVector(
        ntuple(N_up) do i
            percentile_bounds_mean_norm(
                1 - a_surf + (i - 1) * FT(a_surf / N_up),
                1 - a_surf + i * FT(a_surf / N_up),
                1000,
            )
        end,
    )
    SV = typeof(surface_scalar_coeff)
    return SurfaceModel{FT, SV}(;
        scalar_coeff = surface_scalar_coeff,
        a = a_surf,
    )
end

"""
    PressureModel

A pressure model for EDMF, containing
all related model and free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct PressureModel{FT <: AbstractFloat}
    "Pressure drag"
    α_d::FT = 10.0
    "Pressure advection"
    α_a::FT = 0.1
    "Pressure buoyancy"
    α_b::FT = 0.12
    "Default updraft height"
    H_up::FT = 500
end

"""
    MixingLengthModel

A mixing length model for EDMF, containing
all related model and free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct MixingLengthModel{FT <: AbstractFloat}
    "dissipation coefficient"
    c_d::FT = 0.22
    "Eddy Viscosity"
    c_m::FT = 0.14
    "Static Stability coefficient"
    c_b::FT = 0.63
    "Empirical stability function coefficient"
    a1::FT = 0.2
    "Empirical stability function coefficient"
    a2::FT = 100
    "Von Karmen constant"
    κ::FT = 0.4
    "Maximum mixing length"
    max_length::FT = 1e6
    "Prandtl number empirical coefficient"
    ω_pr::FT = 53.0 / 13.0
    "Prandtl number scale"
    Pr_n::FT = 1 #0.74
    "Critical Richardson number"
    Ri_c::FT = 0.25
    "Random small number variable that should be addressed"
    random_minval::FT = 1e-9
    "smooth minimum's fractional upper bound"
    smin_ub::FT = 0.1
    "smooth minimum's regularization minimum"
    smin_rm::FT = 1.5
end

abstract type AbstractStatisticalModel end
struct SubdomainMean <: AbstractStatisticalModel end
struct GaussQuad <: AbstractStatisticalModel end
struct LogNormalQuad <: AbstractStatisticalModel end

"""
    MicrophysicsModel

A microphysics model for EDMF, containing
all related model and free parameters and
assumed subdomain distributions.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct MicrophysicsModel{FT <: AbstractFloat, SM}
    "Subdomain statistical model"
    statistical_model::SM
end

"""
    MicrophysicsModel(
        FT;
        statistical_model = SubdomainMean()
    )

Constructor for `MicrophysicsModel` for EDMF, given:
 - `FT`, the float type used
 - `statistical_model`, the assumed environmental distribution
            of thermodynamic variables.
"""
function MicrophysicsModel(FT; statistical_model = SubdomainMean())
    args = (statistical_model,)
    return MicrophysicsModel{FT, typeof(statistical_model)}(args...)
end

"""
    Environment <: BalanceLaw
A `BalanceLaw` for the environment subdomain arising in EDMF.
"""
Base.@kwdef struct Environment{FT <: AbstractFloat, N_quad} <: BalanceLaw end

"""
    Updraft <: BalanceLaw
A `BalanceLaw` for the updraft subdomains arising in EDMF.
"""
Base.@kwdef struct Updraft{FT <: AbstractFloat} <: BalanceLaw end

"""
    EDMF <: TurbulenceConvectionModel

A turbulence convection model for the EDMF
scheme, containing all closure models and
free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct EDMF{FT <: AbstractFloat, N, UP, EN, ED, P, S, MP, ML, SD} <:
                   TurbulenceConvectionModel
    "Updrafts"
    updraft::UP
    "Environment"
    environment::EN
    "Entrainment-Detrainment model"
    entr_detr::ED
    "Pressure model"
    pressure::P
    "Surface model"
    surface::S
    "Microphysics model"
    micro_phys::MP
    "Mixing length model"
    mix_len::ML
    "Subdomain model"
    subdomains::SD
end

"""
    EDMF(
        FT, N_up, N_quad;
        updraft = ntuple(i -> Updraft{FT}(), N_up),
        environment = Environment{FT, N_quad}(),
        entr_detr = EntrainmentDetrainment{FT}(),
        pressure = PressureModel{FT}(),
        surface = SurfaceModel{FT}(N_up),
        micro_phys = MicrophysicsModel(FT),
        mix_len = MixingLengthModel{FT}(),
        subdomain = SubdomainModel(FT, N_up),
    )
Constructor for `EDMF` subgrid-scale scheme, given:
 - `FT`, the AbstractFloat type used
 - `N_up`, the number of updrafts
 - `N_quad`, the quadrature order. `N_quad^2` is
        the total number of quadrature points
        used for environmental distributions.
 - `updraft`, a tuple containing N_up updraft BalanceLaws
 - `environment`, the environment BalanceLaw
 - `entr_detr`, an `EntrainmentDetrainment` model
 - `pressure`, a `PressureModel`
 - `surface`, a `SurfaceModel`
 - `micro_phys`, a `MicrophysicsModel`
 - `mix_len`, a `MixingLengthModel`
 - `subdomain`, a `SubdomainModel`
"""
function EDMF(
    FT,
    N_up,
    N_quad;
    updraft = ntuple(i -> Updraft{FT}(), N_up),
    environment = Environment{FT, N_quad}(),
    entr_detr = EntrainmentDetrainment{FT}(),
    pressure = PressureModel{FT}(),
    surface = SurfaceModel{FT}(N_up),
    micro_phys = MicrophysicsModel(FT),
    mix_len = MixingLengthModel{FT}(),
    subdomain = SubdomainModel(FT, N_up),
)
    args = (
        updraft,
        environment,
        entr_detr,
        pressure,
        surface,
        micro_phys,
        mix_len,
        subdomain,
    )
    return EDMF{FT, N_up, typeof.(args)...}(args...)
end

import ClimateMachine.TurbulenceConvection: turbconv_sources, turbconv_bcs

"""
    EDMFBCs <: TurbConvBC

Boundary conditions for EDMF.
"""
struct EDMFBCs <: TurbConvBC end
n_updrafts(m::EDMF{FT, N_up}) where {FT, N_up} = N_up
n_updrafts(m::TurbulenceConvectionModel) = 0
turbconv_filters(m::TurbulenceConvectionModel) = ()
turbconv_filters(m::EDMF) = (
    "turbconv.environment.ρatke",
    "turbconv.environment.ρaθ_liq_cv",
    "turbconv.environment.ρaq_tot_cv",
    "turbconv.updraft",
)
n_quad_points(m::Environment{FT, N_quad}) where {FT, N_quad} = N_quad
turbconv_sources(m::EDMF) = (turbconv_source!,)
turbconv_bcs(::EDMF) = EDMFBCs()

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using Documenter, NDReducibles

makedocs(;
    modules=[NDReducibles],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        hide("internals.md"),
    ],
    repo="https://github.com/tkf/NDReducibles.jl/blob/{commit}{path}#L{line}",
    sitename="NDReducibles.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    strict = v"1.0" <= VERSION < v"1.2-",
)

deploydocs(;
    repo="github.com/tkf/NDReducibles.jl",
)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at

#   http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

## Blur estimator pseudocoede

# Algorithm for image quality measure
# Input: Image I of size M×N.
# Output: Image Quality measure (FM) where FM stands for Frequency Domain Image Blur Measure
# Step 1: Compute F which is the Fourier Transform representation of image I
# Step 2: Find Fc which is obtained by shifting the origin of F to centre.
# Step 3: Calculate AF = abs (Fc) where AF is the absolute value of the centered Fourier transform of image I.
# Step 4: Calculate M = max (AF) where M is the maximum value of the frequency component in F.
# Step 5: Calculate TH = the total number of pixels in F whose pixel value > thres, where thres = M/1000.
# Step 6: Calculate Image Quality measure (FM) from equation (1).

#   equation 1: measure = TH / (M x N)

Function blurriness(image)
    # Find M and N
    local   F = fft(I)
    local  Fc = fft(shiftOrigin(I))
    local  AF = abs(Fc) # Absolute value of centered fourier tranformof image I
    local  Mf = max(AF) # Maximum value of frequency component in F
    local  TH = getNoOfPixelsAboveTreshold(Mf/1000)
    local  measure = TH / (M*N)
end

# The idea is then to use a window around every pixel (e.g. 16x16) and then calculate
# a value for every pixel.  This value is then used to select from which image plane
# we should get a pixel (we should always get the sharpest pixel .-)


text:
using SoftRobots
using PyPlot
using MultiPoly

function copy_dict_values(dict::Dict)
    Dict(zip(keys(dict), deepcopy(values(dict))))
end

function test_two_snakes()

    world = SoftRobots.World2D()
    world.friction = 1.5
    world_state = SoftRobots.WorldState()

    snake1, snake1_state = SoftRobots.snake()
    push!(world.objects, snake1)
    world_state[snake1] = snake1_state

    snake2, snake2_state = SoftRobots.snake()
    for i in 1:length(snake2_state.positions)
        snake2_state.positions[i] -= [0, 0.15]
    #     snake2_state.velocities[i] += [-0.2, 4]
    end
    push!(world.objects, snake2)
    world_state[snake2] = snake2_state

    x, y = generators(MPoly{Float64}, :x, :y)
    terrain = SoftRobots.FixedObject()
    # terrain_state = SoftRobots.FixedObjectState(-(0.1 + 0.5x^2 - 1y))
    terrain_state = SoftRobots.FixedObjectState(-0.1 + y)
    push!(world.objects, terrain)
    world_state[terrain] = terrain_state

    vis = SoftRobots.PyPlotVisualizer(world)

    history = typeof(world_state)[]

    function playback(history::Array{SoftRobots.WorldState}, dt=1)
        for i = 1:dt:length(history)
            SoftRobots.draw(vis, history[i])
            xlim([0,1])
            ylim([0,1])
            PyPlot.draw()
        end
    end

    pygui(true)
    PyPlot.ion()

    for j = 0:2000
        SoftRobots.update!(world, world_state, 0.001)
        push!(history, copy_dict_values(world_state))
        if mod(j, 100) == 0
            SoftRobots.draw(vis, world_state)
            xlim([0,1])
            ylim([0,1])
        end
    end

    playback(history, 10)
end

test_two_snakes()


text:
using DSPopt

include("farmer_models.jl")

m = farmer_blocks()
DSPopt.classifyConstrs(m)

"""
NOTE:
  JuMP v0.21 changes the order of indexing variables.
  This affects some of the test lines below, which were based on the older variable indexing.
"""

@testset "Parent model" begin
    @test length(m.variables) == 27
    @test length(m.constraints) == 6
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(m)
    @test length(start) == length(m.constraints) + 1
    @test start[end] == length(index)
    @test start == [0, 2, 4, 6, 8, 10, 12]
    # @test index == [0, 1, 1, 2, 3, 4, 4, 5, 6, 7, 7, 8] # will fail with JuMP v0.21
    @test value == [1., -1., 1., -1., 1., -1., 1., -1., 1., -1., 1., -1.]
    @test rlbd == zeros(6)
    @test rubd == zeros(6)
    # @test obj == [
    #     probability[1]*Cost[1], probability[2]*Cost[1], probability[3]*Cost[1],
    #     probability[1]*Cost[2], probability[2]*Cost[2], probability[3]*Cost[2],
    #     probability[1]*Cost[3], probability[2]*Cost[3], probability[3]*Cost[3],
    #     probability[1]*Purchase[1], probability[2]*Purchase[1], probability[3]*Purchase[1],
    #     probability[1]*Purchase[2], probability[2]*Purchase[2], probability[3]*Purchase[2],
    #     -probability[1]*Sell[1], -probability[2]*Sell[1], -probability[3]*Sell[1],
    #     -probability[1]*Sell[2], -probability[2]*Sell[2], -probability[3]*Sell[2],
    #     -probability[1]*Sell[3], -probability[2]*Sell[3], -probability[3]*Sell[3],
    #     -probability[1]*Sell[4], -probability[2]*Sell[4], -probability[3]*Sell[4]] # will fail with JuMP v0.21
    @test clbd == zeros(27)
    @test cubd == zeros(27) .+ Inf
    @test ctype == "IIIIIIIIICCCCCCCCCCCCCCCCCC"
end

@testset "Child model $i" for (i,subm) in m.children
    @test i > 0
    @test length(subm.variables) == 0
    @test length(subm.constraints) == 5
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(subm, i)
    @test length(start) == length(subm.constraints) + 1
    @test start[end] == length(index)
    @test start == [0, 3, 6, 9, 12, 13]
    # @test index == [0, 3, 6,  0, 9, 15,  3, 12, 18,  6, 21, 24,  21] .+ (i-1) # will fail with JuMP v0.21
    @test value == [1., 1., 1.,  Yield[i,1], 1., -1.,  Yield[i,2], 1., -1.,  Yield[i,3], -1., -1., 1.]
    @test rlbd == [-Inf; Minreq; -Inf]
    @test rubd == [Budget, Inf, Inf, Inf, 6000]
    @test obj == []
    @test clbd == []
    @test cubd == []
    @test ctype == ""
end

@testset "optimize!: $j" for j in [DSPopt.ExtensiveForm, DSPopt.DW]
    status = DSPopt.optimize!(m, solve_type = j, is_stochastic = false)
    @test status == MOI.OPTIMAL
    @test isapprox(objective_value(m), -108390., rtol=0.01)
    @test isapprox(dual_objective_value(m), -108390., rtol=0.01)

    primsol = value()
    dualsol = dual()
    for s = 1:3
        @test primsol[s] == []
    end
    @test dualsol == []

    x = m[:x]
    y = m[:y]
    w = m[:w]
    @test isapprox(value(x[1,1]), 170.0)
    @test isapprox(value(x[1,2]), 170.0)
    @test isapprox(value(x[1,3]), 170.0)
    @test isapprox(value(x[2,1]), 80.0)
    @test isapprox(value(x[2,2]), 80.0)
    @test isapprox(value(x[2,3]), 80.0)
    @test isapprox(value(x[3,1]), 250.0)
    @test isapprox(value(x[3,2]), 250.0)
    @test isapprox(value(x[3,3]), 250.0)
    @test isapprox(value(y[1,1]), 0.0)
    @test isapprox(value(y[1,2]), 0.0)
    @test isapprox(value(y[1,3]), 0.0)
    @test isapprox(value(y[2,1]), 0.0)
    @test isapprox(value(y[2,2]), 0.0)
    @test isapprox(value(y[2,3]), 48.0)
    @test isapprox(value(w[1,1]), 310.0)
    @test isapprox(value(w[1,2]), 225.0)
    @test isapprox(value(w[1,3]), 140.0)
    @test isapprox(value(w[2,1]), 48.0, rtol=0.01)
    @test isapprox(value(w[2,2]), 0.0)
    @test isapprox(value(w[2,3]), 0.0)
    @test isapprox(value(w[3,1]), 6000.0)
    @test isapprox(value(w[3,2]), 5000.0)
    @test isapprox(value(w[3,3]), 4000.0)
    @test isapprox(value(w[4,1]), 0.0)
    @test isapprox(value(w[4,2]), 0.0)
    @test isapprox(value(w[4,3]), 0.0)
    DSPopt.freeSolver(dsp)
    @testset "freeModel" begin
        DSPopt.freeModel(dsp)
        @test dsp.p != C_NULL
        @test length(dsp.numRows) == 0
        @test length(dsp.numCols) == 0
        @test isnan(dsp.primVal)
        @test isnan(dsp.dualVal)
        @test length(dsp.colVal) == 0
        @test length(dsp.rowVal) == 0
        @test dsp.nblocks == -1
        @test dsp.block_ids == []
        @test dsp.is_stochastic == false
        @test dsp.solve_type == DSPopt.Dual
        @test length(dsp.quadConstrs) == 0
        @test length(dsp.linConstrs) == 0
    end
end


text:
using SafeTestsets
const LONGER_TESTS = false

@time @safetestset "Stan_String" begin include("test_stan_string.jl") end
@time @safetestset "Stan" begin include("test_stan.jl") end


text:
######### Stan program example  ###########

using StanSample, MCMCChains

ProjDir = @__DIR__

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
"

datatheta = joinpath(ProjDir, "bernoulli.data.R")
inittheta = joinpath(ProjDir , "bernoulli.init.R")

sm = SampleModel("bernoulli", bernoullimodel);

rc = stan_sample(sm, data=datatheta, init=inittheta)

if success(rc)
  chn = read_samples(sm; output_format=:mcmcchains)
  show(chn)
end


text:
println("Di) Test of the RayleighCompton module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
#
setDefaults("print summary: open", "zzz-RayleighCompton.sum")

if  false
    # Green function for Rayleigh scattering on the ground-configuration levels of B-like neon
    grid             = Radial.Grid(Radial.Grid(true), rnt = 4.0e-6, h = 5.0e-2, rbox = 10.0)
    name             = "Rayleigh scattering on Be-like 1s^2 2s^2 2p neon"
    refConfigs       = [Configuration("[He] 2s^2 2p")]
    levelSymmetries  = [LevelSymmetry(1//2, Basics.plus), LevelSymmetry(3//2, Basics.minus)]
    greenSettings    = GreenSettings(3, [0, 1], 0.01, true, LevelSelection())
    greenRep         = Representation(name, Nuclear.Model(10.), Radial.Grid(true), refConfigs,
                                      GreenExpansion( AtomicState.DampedSpaceCI(), Basics.DeExciteSingleElectron(), levelSymmetries, 3, greenSettings) )
    greenOut         = generate(greenRep, output=true)
    green            = greenOut["Green channels"]

elseif true
    # Green function for Rayleigh scattering on the ground-configuration levels of B-like neon
    grid             = Radial.Grid(Radial.Grid(true), rnt = 4.0e-6, h = 5.0e-2, rbox = 10.0)
    asfSettings      = AsfSettings(AsfSettings(), scField=Basics.HSField())  # not used
    rayleighSettings = RayleighCompton.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], [10., 100.], green, true, true, true,
                                                ExpStokes(), [SolidAngle(1.0, 0.0)], LineSelection())
    setDefaults("unit: rate", "1/s")

    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(10.),
                            initialConfigs=[Configuration("1s^2 2s^2 2p")],
                            finalConfigs  =[Configuration("1s^2 2s^2 2p")],
                            processSettings = rayleighSettings )
    @show wa
    wb = perform(wa)
end

setDefaults("print summary: close", "")


text:
#=
    Contents: Non app-specific handlers.
=#

function notfound!(res)
    res.status = 404
    res.data   = "Requested resource not found."
end

function badrequest!(res)
    res.data   = "Bad request"
    res.status = 400
end

"Redirect user to destination_path."
function redirect!(res, destination_path::AbstractString)
    res.status = 303
    res.headers["Location"] = destination_path
end

# EOF


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using ValueShapes, Distributions, ForwardDiff

@testset "test_distribution_transform" begin
    #=
    using Cuba
    function integrate_over_unit(density::AbstractDensity)
        vs = varshape(density)
        f_cuba(source_x, y) = y[1] = exp(logvalof(density)(vs(source_x)))
        Cuba.vegas(f_cuba, 1, 1).integral[1]
    end
    =#
end

@testset "bat_transform_defaults" begin
    mvn = @inferred(product_distribution([Normal(-1), Normal(), Normal(1)]))
    uniform_prior = @inferred(product_distribution([Uniform(-3, 1), Uniform(-2, 2), Uniform(-1, 3)]))

    posterior_uniform_prior = @inferred(PosteriorDensity(mvn, uniform_prior))
    posterior_gaussian_prior = @inferred(PosteriorDensity(mvn, mvn))

    @test @inferred(bat_transform(PriorToGaussian(), posterior_uniform_prior)).result.prior.dist == @inferred(BAT.StandardMvNormal(3))
    @test @inferred(bat_transform(PriorToUniform(), posterior_gaussian_prior)).result.prior.dist == @inferred(BAT.StandardMvUniform(3))
    @test @inferred(bat_transform(NoDensityTransform(), posterior_uniform_prior)).result.prior.dist == uniform_prior
    pd = @inferred(product_distribution([Uniform() for i in 1:3]))
    density = @inferred(BAT.DistributionDensity(pd))
    @test @inferred(bat_transform(NoDensityTransform(), density)).result.dist == density.dist

    # ToDo: Improve comparison for bounds so `.dist` is not required here:
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardUvUniform()))).result.dist == convert(AbstractDensity, BAT.StandardUvUniform()).dist
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardMvUniform(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvUniform(4)).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardUvNormal()))).result.dist == convert(AbstractDensity, BAT.StandardUvNormal()).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardMvNormal(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvNormal(4)).dist
end


text:
using FinancialToolbox
if !(VERSION.major==0&&VERSION.minor<=6)
	using Test,Dates
else
	using Base.Test
end
function yearFractionTester(StartDate,EndDate,MatlabResults,TestToll=1e-14)
for i=0:FinancialToolbox.currMaxImplemented
	@test(abs(MatlabResults[i+1]-yearfrac(StartDate,EndDate,i))<TestToll)
end
return ;
end

StartDate=Date(1992,12,14);
EndDate=Date(1996,2,28);
@test(-yearfrac(StartDate,StartDate,0)==0)
@test(-yearfrac(StartDate,EndDate,0)==yearfrac(EndDate,StartDate,0))
MatlabResults=[3.20821917808219;3.20555555555556;3.25277777777778;3.20821917808219;3.20555555555556;3.20555555555556;3.20555555555556;3.20821917808219;3.20821917808219;3.25277777777778;3.20821917808219;3.20555555555556;3.20765027322409];
yearFractionTester(StartDate,EndDate,MatlabResults);

EndDate2=Date(1996,2,29);
MatlabResults2=[3.21095890410959;3.20833333333333;3.25555555555556;3.21095890410959;3.20833333333333;3.20833333333333;3.20833333333333;3.21095890410959;3.21095890410959;3.25555555555556;3.21095890410959;3.20833333333333;3.21038251366122];
yearFractionTester(StartDate,EndDate2,MatlabResults2);

EndDate3=Date(2026,2,28);
MatlabResults3=[33.2301369863014;33.2055555555556;33.6916666666667;33.2301369863014;33.2055555555556;33.2055555555556;33.2055555555556;33.2082191780822;33.2301369863014;33.6916666666667;33.2301369863014;33.2055555555556;33.2080844374580];
TestToll=1e-13;
yearFractionTester(StartDate,EndDate3,MatlabResults3,TestToll);

StartDate2=Date(1992,2,29);
EndDate4=Date(1996,2,29);
MatlabResults4=[3.99180327868852;4;4.05833333333333;4.00273972602740;3.99722222222222;4;4;4;3.99180327868852;4.05833333333333;4.00273972602740;4;4]
yearFractionTester(StartDate2,EndDate4,MatlabResults4);

StartDate3=Date(1992,3,30);
EndDate5=Date(1996,12,31);
MatlabResults5=[4.75890410958904;4.75000000000000;4.82500000000000;4.75890410958904;4.75000000000000;4.75000000000000;4.75000000000000;4.75616438356164;4.75890410958904;4.82500000000000;4.75890410958904;4.75000000000000;4.75409836065569]
yearFractionTester(StartDate3,EndDate5,MatlabResults5);

StartDate4=Date(1992,3,31);
EndDate6=Date(1996,12,30);
MatlabResults6=[4.75342465753425;4.75000000000000;4.81944444444445;4.75342465753425;4.75000000000000;4.75000000000000;4.75000000000000;4.75068493150685;4.75342465753425;4.81944444444445;4.75342465753425;4.75000000000000;4.74863387978144]
yearFractionTester(StartDate4,EndDate6,MatlabResults6);

@test_throws(ErrorException,yearfrac(StartDate,EndDate3,FinancialToolbox.currMaxImplemented+1));
@test_throws(ErrorException,yearfrac(StartDate,EndDate3,-1));

@test(daysact(Date(1992,12,14),Date(1992,12,13))==-1)

##Excel test
@test(fromExcelNumberToDate(33952)==Date(1992,12,14))

println("Test Dates Passed")