text:
__precompile__()

module MultipleScattering

## Symmetries
export Symmetry, AbstractSymmetry, AbstractPlanarSymmetry, AbstractAzimuthalSymmetry
export WithoutSymmetry, PlanarSymmetry, PlanarAzimuthalSymmetry, AzimuthalSymmetry, RadialSymmetry

## Shapes
export Shape, Circle, Rectangle, Box, EmptyShape, Halfspace, Plate, TimeOfFlight, TimeOfFlightFromPoint

export outer_radius, volume, name, iscongruent, (≅), congruent, in, issubset, origin, shape, Sphere, (==), isequal, show
export boundary_functions, boundary_points, boundary_data, bounding_box, corners
export points_in_shape, bottomleft, topright

## Physical mediums
export PhysicalMedium, outgoing_basis_function, regular_basis_function, outgoing_radial_basis, regular_radial_basis, outgoing_translation_matrix, regular_translation_matrix, estimate_regular_basisorder, estimate_outgoing_basisorder, basisorder_to_basislength, basis_coefficients, basislength_to_basisorder, internal_field, check_material
export spatial_dimension, field_dimension

## Electromagnetic
export Electromagnetic

## Acoustics
export Acoustic, AcousticCapsule, sound_hard, hard, rigid, zero_neumann, sound_soft, soft, zero_dirichlet, pressure_release, impedance

## Particles
export AbstractParticle, Particle, CapsuleParticle, AbstractParticles

## RegularSources
export AbstractSource, RegularSource, source_expand, regular_spherical_coefficients, self_test, constant_source, point_source, plane_source, regular_spherical_source, (*), (+)
export PlaneSource

## Main simulation and results
export Simulation, run, TimeSimulation, forcing, field
export SimulationResult, FrequencySimulation, FrequencySimulationResult, size

## Impulses for time response
export ContinuousImpulse, TimeDiracImpulse, FreqDiracImpulse, GaussianImpulse
export DiscreteImpulse, continuous_to_discrete_impulse, DiscreteTimeDiracImpulse, DiscreteGaussianImpulse
export TimeSimulationResult, frequency_to_time, time_to_frequency
export ω_to_t, t_to_ω, firstnonzero

export random_particles, statistical_moments

export t_matrix, get_t_matrices
export scattering_matrix

import Printf: @printf

using StaticArrays: SVector
using OffsetArrays: OffsetArray

using SpecialFunctions: besselj, hankelh1
using WignerSymbols, GSL
using Random, LinearAlgebra, RecipesBase, Statistics
using ProgressMeter

# Generic machinery common to all physical models
include("types.jl")
include("shapes/shapes.jl")

include("physics/special_functions.jl") # Special functions missing from Base library
include("physics/physical_medium.jl")

include("particle.jl")
include("source.jl")
include("result.jl")
include("simulation.jl")
include("impulse.jl")
include("time_simulation.jl")
include("random/random.jl")
include("t_matrix.jl")
include("scattering_matrix.jl")

# Specific physical models
include("physics/diffbessel.jl")
include("physics/acoustics/export.jl")
include("physics/electromagnetism.jl")

#Plotting and graphics
include("plot/plot.jl")

# Precompile hints
include("precompile.jl")

end # module


text:
# Character Sets
#
# Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

@api public  CharSet, UniPlusCharSet
@api develop charset_types

struct CharSet{CS}   end

"""List of installed character sets"""
const charset_types = []

CharSet(s)  = CharSet{Symbol(s)}()

# List of basic character sets

show(io::IO, ::Type{CharSet{S}}) where {S} = print(io, "CharSet{:$S}")

const UniPlusCharSet = CharSet{:UniPlus}
push!(charset_types, UniPlusCharSet)

for lst in cse_info
    length(lst) > 2 && continue
    nam = lst[1]
    csnam = symstr(nam, "CharSet")
    @eval const $csnam = CharSet{$(quotesym(nam))}
    @eval show(io::IO, ::Type{$csnam}) = print(io, $(string(csnam)))
    @eval push!(charset_types, $csnam)
    @eval @api $(String(nam)[1] == '_' ? :develop : :public) $csnam
end

# Handle a few quirks
charset(::Type{<:AbstractChar}) = UTF32CharSet
charset(::Type{UInt8})          = BinaryCharSet  # UInt8 instead of "BinaryChr"
charset(::Type{Char})           = UniPlusCharSet # Char instead of "UniPlusChr"

ext: jl


text:
function save_mesh(file::AbstractString, m::AirfoilMesh, format::Symbol; kwargs...)
    open(file, "w") do io
        save_mesh(io, m, format; kwargs...)
    end
end

const _valid_formats = (:TECPLOT, :FLUENT, :PLOT3D)

function save_mesh(io::IO, m::AirfoilMesh, format::Symbol; kwargs...)
    @assert format in _valid_formats "invalid format :$format, must be one of $(_valid_formats)"
    save_mesh(Val(format), io, m; kwargs...)
end

include("tecplot.jl")
include("fluent.jl")
include("plot3d.jl")


text:
using Logging
global_logger(ConsoleLogger(stderr, Logging.Debug)) # modify the debugging level
# this file is best used interactively
using Distributed
addprocs(4)
using ScatteringTransform
@everywhere using Revise
@everywhere using ScatteringTransform
using Plots
using Statistics
# simple example; a sum of sine waves that has a discontinuous jump
t = 0.001:6π/100:6π
f = 2 .* max.(0,-(t.-3*π))./(t.-3*π) .* (sin.(2π*t)) .+ 10 .* max.(0,t.-3*π)./(t.-3*π) .+ max.(0,t.-3*π)./(t.-3*π).* sin.(4π*t.+π)
plot(t, f)

# First create the transform configuration
layers = stParallel(2, length(f))

# two major versions:
###########################################################################################################
#####    the first outputs a datatype specifically for working with Scattering transforms, and has both all internal results and the output. It is not particulary appropriate for use on large datasets
###########################################################################################################

@time output = st(f, layers, absType(), thin=false) # run twice to get the time. The first is compile time

###########################################################################################################
####    the second outputs a single array where the last dimension corresponds to the transformed data. Useful for large datasets, and more optimized
###########################################################################################################
@time thinOp = st(f, layers, absType(), outputSubsample=(-1,3), thin=true); # run twice to get the time. The first is compile time

# computing the morlet wavelet coefficients subsampled for comparison
waves = wavelet(WT.Morlet(5.0), averagingType = WT.Father())
cwtf = cwt(f, waves)
firstLayer = zeros(eltype(cwtf), 25, size(cwtf, 2))
for i=1:size(firstLayer,2)
  firstLayer[:,i] = ScatteringTransform.bspline(ScatteringTransform.bspline(abs.(cwt(f, waves)[:,i]), 2), 2)
end

plot(abs.(firstLayer)[:,end,1])
# let's take a look at some of the coefficients
# first, for comparison, the cwt using the same Morlet wavelets is
heatmap(abs.(firstLayer)[:,2:end,1]', xaxis="Space",yaxis="Frequency", title="Morlet coefficients of f")
# where the discontinuity is clearly visible, as is the difference in frequency
# For comparison, the first layer (the difference in magnitude is due to normalization)
heatmap(output[2][:,:,1]', xlabel="space", ylabel="Frequency",title="First Layer Output (smoothed version of the previous)")
# the discontinuity is less visible here, but it remains in the next layer

# let's take a closer look at the 6th frequency, which is where the sinusoid above is represented
p = pathType([6]) # an object that we can use as an index into the output
plot(output[p], title="6th frequency coefficients") # biased towards the left, where the oscillation is present
# now it's children
p = pathType([6,:])
heatmap(output[p][:,:,1]', xlabel="space", ylabel="Frequency", title="6th frequency variation")
# now the highest frequency, which should reflect the discontinuity
p = pathType([15,:])
heatmap(output[p][:,:,1]', xlabel="space", ylabel="Frequency", title="15th frequency decomposition")
p = pathType([14,:])
heatmap(output[p][:,:,1]', xlabel="space", ylabel="Frequency", title="15th frequency decomposition")

# An example from a sonar dataset
using JLD2
@load "examples/resp1999.0_10.2.jld2"
resp = permutedims(resp, (1,3,2))[:,:,1:45]
# this example is the reflection off of a sharkfin whose interior liquid has speed 1999.0 at distance 10.2 meters
# it is now ordered as (range)×(cross-range)×(angle)

layers = stParallel(2, size(resp,1))
size(resp)
# since this is going to be pretty large, we should subsample heavily (see the docs via ?st for options)
@time tmpOutput = st(resp[:,:,1], layers, absType(), outputSubsample=(-1,5))
output = zeros(size(tmpOutput)..., size(resp,3))
for i=1:size(resp,3)
    @time output[:,:,i] = st(resp[:,:,i], layers, absType(), outputSubsample=(-1,5))
end
output
# output is now prepared. At this point, the preprocessing is done, and you could do regression, clustering, etc.


text:
"""
	RigidBody{R<:Rotation} <: LieGroupModel

Abstraction of a dynamical system with free-body dynamics, with a 12 or 13-dimensional state
vector: `[p; q; v; ω]`
where `p` is the 3D position, `q` is the 3 or 4-dimension attitude representation, `v` is the
3D linear velocity, and `ω` is the 3D angular velocity.

# Interface
Any single-body system can leverage the `RigidBody` type by inheriting from it and defining the
following interface:
```julia
forces(::MyRigidBody, x, u)  # return the forces in the world frame
moments(::MyRigidBody, x, u) # return the moments in the body frame
inertia(::MyRigidBody, x, u) # return the 3x3 inertia matrix
mass(::MyRigidBody, x, u)  # return the mass as a real scalar
```

Instead of defining `forces` and `moments` you can define the higher-level `wrenches` function
	wrenches(model::MyRigidbody, z::AbstractKnotPoint)
	wrenches(model::MyRigidbody, x, u)

# Rotation Parameterization
A `RigidBody` model must specify the rotational representation being used. Any `Rotations.Rotation{3}`
can be used, but we suggest one of the following:
* `UnitQuaternion`
* `MRP`
* `RodriguesParam`
"""
abstract type RigidBody{R<:Rotation} <: LieGroupModel end

LieState(::RigidBody{R}) where R = LieState(R, (3,6))

function Base.rand(model::RigidBody{D}) where {D}
    n,m = size(model)
    r = @SVector rand(3)
    q = rand(D)
    v = @SVector rand(3)
    ω = @SVector rand(3)
    x = build_state(model, r, q, v, ω)
    u = @SVector rand(m)  # NOTE: this is type unstable
    return x,u
end

function Base.zeros(model::RigidBody{D}) where D
    n,m = size(model)
    r = @SVector zeros(3)
    q = one(D)
    v = @SVector zeros(3)
    ω = @SVector zeros(3)
    x = build_state(model, r, q, v, ω)
    u = @SVector zeros(m)  # NOTE: this is type unstable
    return x,u
end

@inline rotation_type(::RigidBody{D}) where D = D

"""
    gen_inds(model::RigidBody)

Generate a `NamedTuple` containing the indices of the position (`r`), orientation (`q`),
linear velocity (`v`), and angular velocity (`ω`) from the state vector for `model`.
"""
@generated function gen_inds(model::RigidBody{R}) where R
    iF = SA[1,2,3]
    iM = SA[4,5,6]
    ir, iq, iv, iω = SA[1,2,3], SA[4,5,6], SA[7,8,9], SA[10,11,12]
    if R <: UnitQuaternion
        iq = push(iq, 7)
        iv = iv .+ 1
        iω = iω .+ 1
    end
    quote
        m = control_dim(model)
        iu = $iω[end] .+ SVector{m}(1:m)
        return (r=$ir, q=$iq, v=$iv, ω=$iω, u=iu)
    end
end

# Getters
@inline Base.position(model::RigidBody, x) = x[gen_inds(model).r]
@inline linear_velocity(model::RigidBody, x) = x[gen_inds(model).v]
@inline angular_velocity(model::RigidBody, x) = x[gen_inds(model).ω]

for rot in [RodriguesParam, MRP, RotMatrix, RotationVec, AngleAxis]
    @eval orientation(model::RigidBody{<:$rot}, x::AbstractVector, renorm=false) = ($rot)(x[4],x[5],x[6])
end
function orientation(model::RigidBody{<:UnitQuaternion}, x::AbstractVector,
        renorm=false)
    q = UnitQuaternion(x[4],x[5],x[6],x[7], renorm)
    return q
end

function flipquat(model::RigidBody{<:UnitQuaternion}, x)
    return @SVector [x[1], x[2], x[3], -x[4], -x[5], -x[6], -x[7],
        x[8], x[9], x[10], x[11], x[12], x[13]]
end

"""
    parse_state(model::RigidBody{R}, x, renorm=false)

Return the position, orientation, linear velocity, and angular velocity as separate vectors.
The orientation will be of type `R`. If `renorm=true` and `R <: UnitQuaternion` the quaternion
will be renormalized.
"""
function parse_state(model::RigidBody, x, renorm=false)
    r = position(model, x)
    p = orientation(model, x, renorm)
    v = linear_velocity(model, x)
    ω = angular_velocity(model, x)
    return r, p, v, ω
end

function RBState(model::RigidBody, x, renorm=false)
    RBState(parse_state(model, x, renorm)...)
end

function RBState(model::RigidBody, x::RBState, renorm=false)
    if renorm
        renorm(x)
    else
        x
    end
end

@inline RBState(model::RigidBody, z::AbstractKnotPoint, renorm=false) = RBState(model, state(z), renorm)

"""
    build_state(model::RigidBody{R}, x::RBState) where R
    build_state(model::RigidBody{R}, x::AbstractVector) where R
    build_state(model::RigidBody{R}, r, q, v, ω) where R

Build the state vector for `model` using the `RBState` `x`. If `R <: UnitQuaternion` this
    returns `x` cast as an `SVector`, otherwise it will convert the quaternion in `x` to
    a rotation of type `R`.

Also accepts as arguments any arguments that can be passed to the constructor of `RBState`.
"""
@inline build_state(model::RigidBody{<:UnitQuaternion}, x::RBState) = SVector(x)
function build_state(model::RigidBody{R}, x::RBState) where R <: Rotation
    r = position(x)
    q = Rotations.params(R(orientation(x)))
    v = linear_velocity(x)
    ω = angular_velocity(x)
    SA[
        r[1], r[2], r[3],
        q[1], q[2], q[3],
        v[1], v[2], v[3],
        ω[1], ω[2], ω[3],
    ]
end

function build_state(model::RigidBody{R}, args...) where R <: Rotation
    x_ = RBState(args...)
    build_state(model, x_)
end

@inline function build_state(model::RigidBody,
        r::AbstractVector, q::AbstractVector, v::AbstractVector, ω::AbstractVector)
    @assert length(q) == 3
    SA[
        r[1], r[2], r[3],
        q[1], q[2], q[3],
        v[1], v[2], v[3],
        ω[1], ω[2], ω[3],
    ]
end

@inline function build_state(model::RigidBody{<:UnitQuaternion},
        r::AbstractVector, q::AbstractVector, v::AbstractVector, ω::AbstractVector)
    @assert length(q) == 4
    SA[
        r[1], r[2], r[3],
        q[1], q[2], q[3], q[4],
        v[1], v[2], v[3],
        ω[1], ω[2], ω[3],
    ]
end

function fill_state(model::RigidBody{<:UnitQuaternion}, x::Real, q::Real, v::Real, ω::Real)
    @SVector [x,x,x, q,q,q,q, v,v,v, ω,ω,ω]
end

function fill_state(model::RigidBody, x::Real, q::Real, v::Real, ω::Real)
    @SVector [x,x,x, q,q,q, v,v,v, ω,ω,ω]
end

############################################################################################
#                                DYNAMICS
############################################################################################
function dynamics(model::RigidBody, x, u, t=0)

    r,q,v,ω = parse_state(model, x)

    z = StaticKnotPoint(x, u, 0.0, t)
    ξ = wrenches(model, z)
    F = SA[ξ[1], ξ[2], ξ[3]]  # forces in world frame
    τ = SA[ξ[4], ξ[5], ξ[6]]  # torques in body frame
    m = mass(model)
    J = inertia(model)
    Jinv = inertia_inv(model)

    qdot = Rotations.kinematics(q,ω)
    if velocity_frame(model) == :world
        xdot = v
        vdot = F ./ m
    elseif velocity_frame(model) == :body
        xdot = q * v
        vdot = q\(F ./ m) - ω × v
    end
    ωdot = Jinv*(τ - ω × (J*ω))

    build_state(model, xdot, qdot, vdot, ωdot)
    # [xdot; qdot; vdot; ωdot]
end

# Use the StaticArrays methods since we know the size will always be small enough
@inline function dynamics!(model::RigidBody{D}, xdot, x, u, t=0) where D
    xdot .= dynamics(model, x, u, t)
end

@inline wrenches(model::RigidBody, z::AbstractKnotPoint) = wrenches(model, state(z), control(z), time(z))
function wrenches(model::RigidBody, x, u, t)
    F = forces(model, x, u, t)
    M = moments(model, x, u, t)
    SA[F[1], F[2], F[3], M[1], M[2], M[3]]
end

@inline forces(model::RigidBody, x, u, t) = forces(model, x, u)
@inline moments(model::RigidBody, x, u, t) = moments(model, x, u)

@inline mass(::RigidBody) = throw(ErrorException("Not implemented"))
@inline inertia(::RigidBody)::SMatrix{3,3} = throw(ErrorException("Not implemented"))
@inline inertia_inv(model::RigidBody) = inv(inertia(model))
@inline forces(::RigidBody, x, u)::SVector{3} = throw(ErrorException("Not implemented"))
@inline moments(::RigidBody, x, u)::SVector{3} = throw(ErrorException("Not implemented"))
@inline velocity_frame(::RigidBody) = :world # :body or :world

function jacobian!(model::RigidBody{<:UnitQuaternion}, F, y, x, u, t)
    iF = SA[1,2,3]
    iM = SA[4,5,6]
    ir,iq,iv,iω,iu = gen_inds(model)

    # Extract the info from the state and model
    r,q,v,ω = parse_state(model, x)
    m = mass(model)
    J = inertia(model)
    Jinv = inertia_inv(model)

    ξ = wrenches(model, x, u, t)
    f = SA[ξ[1], ξ[2], ξ[3]]  # forces in world frame
    τ = SA[ξ[4], ξ[5], ξ[6]]  # torques in body frame

    # Calculate the Jacobian wrt the wrench and multiply the blocks according to the sparsity
    F .= 0
    Jw = uview(get_data(F), 8:13, :)
    wrench_jacobian!(Jw, model, z)
    js = wrench_sparsity(model)
    if velocity_frame(model) == :world
        tmp = I * 1/m
    else
        tmp = 1/m * RotMatrix(inv(q))
    end
    js[1,1] && (Jw[iF, ir] .= tmp * Jw[iF, ir])
    js[1,2] && (Jw[iF, iq] .= tmp * Jw[iF, iq])
    js[1,3] && (Jw[iF, iv] .= tmp * Jw[iF, iv])
    js[1,4] && (Jw[iF, iω] .= tmp * Jw[iF, iω])
    js[1,5] && (Jw[iF, iu] .= tmp * Jw[iF, iu])
    js[2,1] && (Jw[iM, ir] .= Jinv * Jw[iM, ir])
    js[2,2] && (Jw[iM, iq] .= Jinv * Jw[iM, iq])
    js[2,3] && (Jw[iM, iv] .= Jinv * Jw[iM, iv])
    js[2,4] && (Jw[iM, iω] .= Jinv * Jw[iM, iω])
    js[2,5] && (Jw[iM, iu] .= Jinv * Jw[iM, iu])

    # Add in the parts of the Jacobian that are not functions of the wrench
    F[iq, iq] .= 0.5 * Rotations.rmult(Rotations.pure_quaternion(ω))
    F[iq, iω] .= 0.5 * Rotations.lmult(q) * Rotations.hmat()
    F[iω, iω] .= Jinv*(skew(J*ω) - skew(ω)*J) + F[iω,iω]
    if velocity_frame(model) == :world
        F[1,8] += 1
        F[2,9] += 1
        F[3,10] += 1
    else
        F[ir,iq] .+= Rotations.∇rotate(q, v)
        F[ir,iv] .+= RotMatrix(q)
        F[iv,iq] .+= 1/m * Rotations.∇rotate(inv(q), f) * Rotations.tmat()
        F[iv,iv] .+= -Rotations.skew(ω)
        F[iv,iω] .+=  Rotations.skew(v)
    end

    return F
end

function wrench_jacobian!(F, model::RigidBody, z)
    function w(x)
        wrenches(model, StaticKnotPoint(z, x))
    end
    ForwardDiff.jacobian!(F, w, z.z)
end

"""
    wrench_sparsity(model::RigidBody)

Specify the sparsity of the wrench Jacobian of `model` as a `js = SMatrix{2,5,Bool,10}`.
The elements of `js` correspond to the block elements of the wrench Jacobian:

```julia
[∂F/∂r ∂F/∂q ∂F/∂v ∂F/∂ω ∂F/∂u;
 ∂M/∂r ∂M/∂q ∂M/∂v ∂M/∂ω ∂M/∂u]
```

where `js[i,j] = false` if the corresponding partial derivative is always zero.

Note that this is only for performance improvement of continuous-time Jacobians of rigid bodies;
specifying the sparsity is completely optional.

# Example
For a fully-actuated satellite where `F = q*u[1:3]` and `M = u[4:6]`, the wrench sparsity
would be

```julia
SA[false true  false false true;
   false false false false true]
```
"""
wrench_sparsity(model::RigidBody) = @SMatrix ones(Bool,2,5)


text:
# TextRow

"""
    TextRow{M, N}(catcols, contcols, categorydict) <: Block

`Block` for table rows with M categorical and N continuous columns. `data`
is valid if it satisfies the `AbstractRow` interface in Tables.jl, values
present in indices for categorical and continuous columns are consistent,
and `data` is indexable by the elements of `catcols` and `contcols`.
"""
struct TextRow{M,N,T} <: Block
    catcols::NTuple{M}
    contcols::NTuple{N}
    categorydict::T
end

function TextRow(catcols, contcols, categorydict)
    TextRow{length(catcols),length(contcols)}(catcols, contcols, categorydict)
end

function checkblock(block::TextRow, x)
    columns = Tables.columnnames(x)
    (
        all(col -> col ∈ columns, (block.catcols..., block.contcols...)) &&
        all(
            col ->
                haskey(block.categorydict, col) &&
                    (ismissing(x[col]) || x[col] ∈ block.categorydict[col]),
            block.catcols,
        ) &&
        all(col -> ismissing(x[col]) || x[col] isa Number, block.contcols)
    )
end

function mockblock(block::TextRow)
    cols = (block.catcols..., block.contcols...)
    vals = map(cols) do col
        col in block.catcols ? rand(block.categorydict[col]) : rand()
    end
    return NamedTuple(zip(cols, vals))
end

"""
    setup(TextRow, data[; catcols, contcols])

Create a `TextRow` block from data container `data::TextDataset`. If the
categorical and continuous columns are not specified manually, try to
guess them from the dataset's column types.
"""
function setup(::Type{TextRow}, data; catcols=nothing, contcols=nothing)
    catcols_, contcols_ = getcoltypes(data)
    catcols = isnothing(catcols) ? catcols_ : catcols
    contcols = isnothing(contcols) ? contcols_ : contcols

    return TextRow(
        catcols,
        contcols,
        gettransformdict(data, DataAugmentation.Categorify, catcols),
    )
end

function Base.show(io::IO, block::TextRow)
    print(io, ShowCase(block, (:catcols, :contcols), show_params=false, new_lines=true))
end

# ## Interpretation

# function showblock!(io, ::ShowText, block::TextRow, obs) end


text:
using Test
using QuantumOptics
using LinearAlgebra, Random

@testset "operators-lazysum" begin

Random.seed!(0)

D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))
D(x1::StateVector, x2::StateVector) = norm(x2-x1)

b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(1)
b3b = GenericBasis(5)

b_l = b1a⊗b2a⊗b3a
b_r = b1b⊗b2b⊗b3b

# Test creation
@test_throws ArgumentError LazySum()
@test_throws AssertionError LazySum([1., 2.], [randoperator(b_l)])
@test_throws bases.IncompatibleBases LazySum(randoperator(b_l, b_r), sparse(randoperator(b_l, b_l)))
@test_throws bases.IncompatibleBases LazySum(randoperator(b_l, b_r), sparse(randoperator(b_r, b_r)))

# Test copy
op1 = 2*LazySum(randoperator(b_l, b_r), sparse(randoperator(b_l, b_r)))
op2 = copy(op1)
@test op1 == op2
@test !(op1 === op2)
op2.operators[1].data[1,1] = complex(10.)
@test op1.operators[1].data[1,1] != op2.operators[1].data[1,1]
op2.factors[1] = 3.
@test op2.factors[1] != op1.factors[1]

# Test dense & sparse
op1 = randoperator(b_l, b_r)
op2 = sparse(randoperator(b_l, b_r))
@test 0.1*op1 + 0.3*dense(op2) == dense(LazySum([0.1, 0.3], [op1, op2]))
@test 0.1*sparse(op1) + 0.3*op2 == sparse(LazySum([0.1, 0.3], [op1, op2]))

# Arithmetic operations
# =====================
op1a = randoperator(b_l, b_r)
op1b = randoperator(b_l, b_r)
op2a = randoperator(b_l, b_r)
op2b = randoperator(b_l, b_r)
op3a = randoperator(b_l, b_r)
op1 = LazySum([0.1, 0.3], [op1a, sparse(op1b)])
op1_ = 0.1*op1a + 0.3*op1b
op2 = LazySum([0.7, 0.9], [sparse(op2a), op2b])
op2_ = 0.7*op2a + 0.9*op2b
op3 = LazySum(op3a)
op3_ = op3a

x1 = Ket(b_r, rand(ComplexF64, length(b_r)))
x2 = Ket(b_r, rand(ComplexF64, length(b_r)))
xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))

# Addition
@test_throws bases.IncompatibleBases op1 + dagger(op2)
@test 1e-14 > D(op1+op2, op1_+op2_)

# Subtraction
@test_throws bases.IncompatibleBases op1 - dagger(op2)
@test 1e-14 > D(op1 - op2, op1_ - op2_)
@test 1e-14 > D(op1 + (-op2), op1_ - op2_)
@test 1e-14 > D(op1 + (-1*op2), op1_ - op2_)

# Test multiplication
@test_throws ArgumentError op1*op2
@test LazySum([0.1, 0.1], [op1a, op2a]) == LazySum(op1a, op2a)*0.1
@test LazySum([0.1, 0.1], [op1a, op2a]) == 0.1*LazySum(op1a, op2a)
@test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2)
@test 1e-11 > D((op1+op2)*(x1+0.3*x2), (op1_+op2_)*(x1+0.3*x2))
@test 1e-12 > D(dagger(x1)*dagger(0.3*op2), dagger(x1)*dagger(0.3*op2_))

# Test division
@test 1e-14 > D(op1/7, op1_/7)

# Test identityoperator
Idense = identityoperator(DenseOperator, b_r)
id = identityoperator(LazySum, b_r)
@test isa(id, LazySum)
@test dense(id) == Idense
@test 1e-11 > D(id*x1, x1)

Idense = identityoperator(DenseOperator, b_l)
id = identityoperator(LazySum, b_l)
@test isa(id, LazySum)
@test dense(id) == Idense
@test 1e-11 > D(xbra1*id, xbra1)

# Test tr and normalize
op1 = randoperator(b_l)
op2 = randoperator(b_l)
op3 = randoperator(b_l)
op = LazySum([0.1, 0.3, 1.2], [op1, op2, op3])
op_ = 0.1*op1 + 0.3*op2 + 1.2*op3

@test tr(op_) ≈ tr(op)
op_normalized = normalize(op)
@test tr(op_) ≈ tr(op)
@test 1 ≈ tr(op_normalized)
op_copy = deepcopy(op)
normalize!(op_copy)
@test tr(op) != tr(op_copy)
@test 1 ≈ tr(op_copy)

# Test partial tr
op1 = randoperator(b_l)
op2 = randoperator(b_l)
op3 = randoperator(b_l)
op123 = LazySum([0.1, 0.3, 1.2], [op1, op2, op3])
op123_ = 0.1*op1 + 0.3*op2 + 1.2*op3

@test 1e-14 > D(ptrace(op123_, 3), ptrace(op123, 3))
@test 1e-14 > D(ptrace(op123_, 2), ptrace(op123, 2))
@test 1e-14 > D(ptrace(op123_, 1), ptrace(op123, 1))

@test 1e-14 > D(ptrace(op123_, [2,3]), ptrace(op123, [2,3]))
@test 1e-14 > D(ptrace(op123_, [1,3]), ptrace(op123, [1,3]))
@test 1e-14 > D(ptrace(op123_, [1,2]), ptrace(op123, [1,2]))

@test_throws ArgumentError ptrace(op123, [1,2,3])

# Test expect
state = Ket(b_l, rand(ComplexF64, length(b_l)))
@test expect(op123, state) ≈ expect(op123_, state)

state = DenseOperator(b_l, b_l, rand(ComplexF64, length(b_l), length(b_l)))
@test expect(op123, state) ≈ expect(op123_, state)

# Permute systems
op1a = randoperator(b1a)
op2a = randoperator(b2a)
op3a = randoperator(b3a)
op1b = randoperator(b1a)
op2b = randoperator(b2a)
op3b = randoperator(b3a)
op1c = randoperator(b1a)
op2c = randoperator(b2a)
op3c = randoperator(b3a)
op123a = op1a⊗op2a⊗op3a
op123b = op1b⊗op2b⊗op3b
op123c = op1c⊗op2c⊗op3c
op = LazySum([0.3, 0.7, 1.2], [op123a, sparse(op123b), op123c])
op_ = 0.3*op123a + 0.7*op123b + 1.2*op123c

@test 1e-14 > D(permutesystems(op, [1, 3, 2]), permutesystems(op_, [1, 3, 2]))
@test 1e-14 > D(permutesystems(op, [2, 1, 3]), permutesystems(op_, [2, 1, 3]))
@test 1e-14 > D(permutesystems(op, [2, 3, 1]), permutesystems(op_, [2, 3, 1]))
@test 1e-14 > D(permutesystems(op, [3, 1, 2]), permutesystems(op_, [3, 1, 2]))
@test 1e-14 > D(permutesystems(op, [3, 2, 1]), permutesystems(op_, [3, 2, 1]))

# Test gemv
op1 = randoperator(b_l, b_r)
op2 = randoperator(b_l, b_r)
op3 = randoperator(b_l, b_r)
op = LazySum([0.1, 0.3, 1.2], [op1, op2, op3])
op_ = 0.1*op1 + 0.3*op2 + 1.2*op3

state = Ket(b_r, rand(ComplexF64, length(b_r)))
result_ = Ket(b_l, rand(ComplexF64, length(b_l)))
result = deepcopy(result_)
operators.gemv!(complex(1.), op, state, complex(0.), result)
@test 1e-13 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemv!(alpha, op, state, beta, result)
@test 1e-13 > D(result, alpha*op_*state + beta*result_)

state = Bra(b_l, rand(ComplexF64, length(b_l)))
result_ = Bra(b_r, rand(ComplexF64, length(b_r)))
result = deepcopy(result_)
operators.gemv!(complex(1.), state, op, complex(0.), result)
@test 1e-13 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemv!(alpha, state, op, beta, result)
@test 1e-13 > D(result, alpha*state*op_ + beta*result_)

# Test gemm
op1 = randoperator(b_l, b_r)
op2 = randoperator(b_l, b_r)
op3 = randoperator(b_l, b_r)
op = LazySum([0.1, 0.3, 1.2], [op1, op2, op3])
op_ = 0.1*op1 + 0.3*op2 + 1.2*op3

state = randoperator(b_r, b_r)
result_ = randoperator(b_l, b_r)
result = deepcopy(result_)
operators.gemm!(complex(1.), op, state, complex(0.), result)
@test 1e-12 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemm!(alpha, op, state, beta, result)
@test 1e-12 > D(result, alpha*op_*state + beta*result_)

state = randoperator(b_l, b_l)
result_ = randoperator(b_l, b_r)
result = deepcopy(result_)
operators.gemm!(complex(1.), state, op, complex(0.), result)
@test 1e-12 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemm!(alpha, state, op, beta, result)
@test 1e-12 > D(result, alpha*state*op_ + beta*result_)

end # testset


text:
function BFS!(G, s)
   v = 1:size(G.AdjMatrix, 1)
   for u in v
     if u != s
        G.V[u].color = 'W'
		G.V[u].d = typemax(Int)
		G.V[u].pred = -1
     end
   end
   G.V[s].color = 'G'
   G.V[s].d = 0
   G.V[s].pred = -1

   Q = Int[]
   push!(Q, s)
   while ~isempty(Q)
        u = Q[1]
	deleteat!(Q, 1)
        for v in findall(G.AdjMatrix[u,:] .==1)
	   if G.V[v].color=='W'
	      G.V[v].color = 'G'
	      G.V[v].d = G.V[u].d + 1
	      G.V[v].pred = u
	      push!(Q, v)
	   end
	end
	G.V[u].color = 'B'
   end
end

function getPath(G, idStart, idEnd)
	path = Int[]
	temp = idEnd;
	while (temp != idStart)
		path = [temp; path]
		temp = G.V[temp].pred;
	end
	return [idStart; path]
end


text:
# Test suite for testing stall logic.
#
# Strategy: Use "TestFifos" to do dependency injection to test the logic of
# the stall suites. This still involves creating a whole "AsapCore", but we
# will have much better control over the state of the various fifos for
# testing.
@testset "Testing Stall Logic" begin
    EmptyFullFifo = AsapSim.EmptyFullFifo

    num_outputs = 8

    # Instantiate a specializd core.
    core = AsapSim.AsapCore{EmptyFullFifo}(
        clock_period = 1,
        fifos = [EmptyFullFifo(), EmptyFullFifo()],
        # Populate all 8 output directions.
        outputs = Dict{Int,EmptyFullFifo}(i => EmptyFullFifo() for i in 1:num_outputs)
    )

    # ------------------------------- #
    # --- Test normal IO stalling --- #
    # ------------------------------- #

    # Set all obuf mask bits to 1. For normal instructions, should not stall
    # right now because all EmptyFullFifo are neither empty nor full.
    core.obuf_mask = trues(8)

    @test AsapSim.stall_check_ibuf(core, 1, false) == false
    @test AsapSim.stall_check_ibuf(core, 2, true) == false

    @test AsapSim.stall_check_obuf(core, 1, false) == false
    @test AsapSim.stall_check_obuf(core, 8, true) == false

    # Make a test instruction reading from both input fifos and writing to
    # an output fifo.
    test_instruction = AsapSim.AsapInstructionKeyword(
        src1 = Loc(AsapSim.IBUF, 0),
        src2 = Loc(AsapSim.IBUF_NEXT, 1),
        dest = Loc(AsapSim.OUTPUT, 0),
    )

    # Since all FIFOs are ready to go, core should not stall.
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.NoStall

    # Start iterating over fifos, first setting inputs to "empty", then output
    # to "full". Make sure that stall logic fires correctly.
    core.fifos[1].empty = true
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.EmptyIbuf

    core.fifos[1].empty = false
    core.fifos[2].empty = true
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.EmptyIbuf

    core.fifos[2].empty = false
    core.outputs[1].full = true
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.FullObuf

    core.outputs[1].full = false

    # New test instruction testing the broadcast mechanism.
    test_instruction = AsapSim.AsapInstructionKeyword(
        dest = AsapSim.OBUF
    )

    # Right now, all outputs are not full, so stall should not be triggered.
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.NoStall

    core.outputs[8].full = true
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.FullObufMask

    core.outputs[8].full = false
    core.outputs[5].full = true
    @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.FullObufMask

    core.outputs[5].full = false

    # Iterate through all output buffers.
    # Set each output to full in turn but deactivate the mask bit for that
    # output.
    #
    # Core should not stall.
    for i in 1:num_outputs
        core.outputs[i].full = true
        core.obuf_mask[i] = false
        @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.NoStall

        # Reset for next test.
        core.outputs[i].full = false
        core.obuf_mask[i] = true
    end

    # Do this again, but deactivate all mask bits first.
    #
    # Then, set each output to "full" and its corresponding mask bit high.
    #
    # Should stall in each case.

    # Set all mask bits to false.
    core.obuf_mask .= false
    for i in 1:num_outputs
        core.outputs[i].full = true
        core.obuf_mask[i] = true
        @test AsapSim.stall_check_io(core, test_instruction) == AsapSim.FullObufMask

        core.outputs[i].full = false
        core.obuf_mask[i] = false
    end

    # ------------------------------ #
    # --- Test STALL instruction --- #
    # ------------------------------ #

    # Reset the state of all fifos.
    for fifo in core.fifos
        fifo.empty = true
    end

    for fifo in values(core.outputs)
        fifo.full = true
    end
    core.obuf_mask .= true

    # Create a mask checking all IO.
    all_mask = 0b10011

    # Stall should return "true" since all inputs are empty and all outputs
    # are full.
    @test AsapSim.stall_check_stall_op(core, all_mask) == AsapSim.StallOp

    # Iterate over inputs, setting each non-empty. Test that the core should
    # not stall. Then do the same for the outputs.
    for fifo in core.fifos
        fifo.empty = false
        @test AsapSim.stall_check_stall_op(core, all_mask) == AsapSim.NoStall
        fifo.empty = true
    end

    for fifo in values(core.outputs)
        fifo.full = false
        @test AsapSim.stall_check_stall_op(core, all_mask) == AsapSim.NoStall
        fifo.full = true
    end

    # Only check outputs, but make inputs non-empty.
    for fifo in core.fifos
        fifo.empty = false
    end

    output_mask = 0b10000
    @test AsapSim.stall_check_stall_op(core, output_mask) == AsapSim.StallOp

    # Do similar checking of the output mask as we did for the IO Stall above.
    core.obuf_mask .= true
    for fifo in values(core.outputs)
        fifo.full = true
    end

    # Ignore non-full fifos. Should always stall.
    for i in 1:num_outputs
        core.outputs[i].full = false
        core.obuf_mask[i] = false
        @test AsapSim.stall_check_stall_op(core, output_mask) == AsapSim.StallOp

        core.outputs[i].full = true
        core.obuf_mask[i] = true
    end

    # Flip the state of everything.
    core.obuf_mask .= false
    for fifo in values(core.outputs)
        fifo.full = true
    end

    for i in 1:num_outputs
        core.outputs[i].full = true
        core.obuf_mask[i] = true
        @test AsapSim.stall_check_stall_op(core, output_mask) == AsapSim.StallOp

        core.outputs[i].full = false
        core.obuf_mask[i] = false
    end
end


text:
__precompile__()

module HyperbolicDiffEq

using Reexport

@reexport using PolynomialBases
import PolynomialBases: integrate, interpolate, interpolate!,
                        compute_coefficients, compute_coefficients!,
                        evaluate_coefficients, evaluate_coefficients!

using Roots
using ArgCheck
using Parameters
using StaticArrays
using MappedArrays

using RecipesBase
using LaTeXStrings

@reexport using DiffEqBase

# interfaces
import DiffEqBase: solve

import Base: *

import StaticArrays: similar_type

# types
"""
An abstract type for a computational mesh.
"""
abstract type AbstractMesh end

# for broadcasting; treat as scalar
Base.broadcastable(mesh::AbstractMesh) = Ref(mesh)

"""
An abstract type representing a balance law in `Dim` space dimensions.
"""
abstract type AbstractBalanceLaw{Dim} end

# for broadcasting; treat as scalar
Base.broadcastable(model::AbstractBalanceLaw) = Ref(model)

"""
An abstract type representing a scalar balance law using `T` as real type in
`Dim` space dimensions.
"""
abstract type ScalarBalanceLaw{T,Dim} <: AbstractBalanceLaw{Dim} end
variables(model::ScalarBalanceLaw{T}) where {T} = T

"""
An abstract type representing a Riemann problem.
"""
abstract type AbstractRiemannProblem end

# for broadcasting; treat as scalar
Base.broadcastable(prob::AbstractRiemannProblem) = Ref(prob)

"""
An abstract type representing the solution of a Riemann problem.
"""
abstract type AbstractRiemannSolution end

# for broadcasting; treat as scalar
Base.broadcastable(sol::AbstractRiemannSolution) = Ref(sol)

"""
An abstract type for a semidiscretisation.
"""
abstract type AbstractSemidiscretisation end

# for broadcasting; treat as scalar
Base.broadcastable(semidisc::AbstractSemidiscretisation) = Ref(semidisc)

"""
An abstract type for a reconstruction.
"""
abstract type AbstractReconstruction end

# for broadcasting; treat as scalar
Base.broadcastable(recons::AbstractReconstruction) = Ref(recons)

"""
An abstract type for a numerical flux.
"""
abstract type NumericalFlux end

# for broadcasting; treat as scalar
Base.broadcastable(flux::NumericalFlux) = Ref(flux)

"""
An abstract type for a dissipation operator.
"""
abstract type DissipationOperator end

# for broadcasting; treat as scalar
Base.broadcastable(op::DissipationOperator) = Ref(op)

include("meshes/meshes1d.jl")
include("meshes/compute_coefficients.jl")

include("finite_volumes/finite_volumes.jl")
include("finite_volumes/numerical_fluxes.jl")
include("finite_volumes/periodic_fv_1d.jl")
include("finite_volumes/reconstructions.jl")
include("finite_volumes/modified_eno.jl")
include("finite_volumes/central_reconstruction.jl")

include("tecno/tecno.jl")

include("flux_difference/uniform_flux_difference.jl")
include("flux_difference/uniform_flux_difference_1d.jl")

include("riemann_problems.jl")
include("balance_laws/constant_linear_advection.jl")
include("balance_laws/burgers.jl")
include("balance_laws/cubic.jl")
include("balance_laws/quartic.jl")
include("balance_laws/quintic.jl")
include("balance_laws/sextic.jl")
include("balance_laws/septic.jl")
include("balance_laws/octic.jl")
include("balance_laws/buckley_leverette.jl")
include("balance_laws/quartic_nonconvex.jl")
include("balance_laws/shallow_water.jl")
include("balance_laws/euler.jl")
include("balance_laws/keyfitz_kranzer.jl")

include("utils.jl")

# models
export ConstantLinearAdvection, Burgers, IntegralQuantitiesBurgers,
       Cubic, Quartic, Quintic, Sextic, Septic, Octic,
       BuckleyLeverette,
       QuarticNonconvex
export ShallowWater, ShallowWaterVar1D
export Euler, EulerVar1D, EulerVar2D, EulerVar3D, IntegralQuantitiesEuler
export KeyfitzKranzer, KeyfitzKranzerVar, IntegralQuantitiesKeyfitzKranzer

export flux, max_abs_speed, naive_max_abs_speed, variables, kinetic_energy,
       entropy, conserved_variables, primitive_variables, entropy_variables,
       satisfies_physical_constraints, flux_potential

export RiemannProblem

# meshes
export UniformMesh1D, UniformPeriodicMesh1D
export compute_coefficients, compute_coefficients!
export evaluate_coefficients, evaluate_coefficients!

# general semidiscretisations
export semidiscretise, max_dt

# finite volume methods
export FirstOrderFV, UniformPeriodicReconstructedFV1D, ScalarUniformPeriodicTecno1D
export CentralReconstruction
export ModifiedENO, ClassicalChoiceENO, BiasedENOChoice,
       MinL2Choice, BiasedMinL2Choice, LexMinLegendreChoice
export GodunovFlux, LocalLaxFriedrichsFlux, HartenLaxVanLeerFlux, HLL, SuliciuFlux,
       KineticFlux
export CentralFlux, MorinishiFlux, DucrosEtAlFlux, KennedyGruberFlux, PirozzoliFlux
export EnergyConservativeFlux, EnergyConservativeFlux1Param, EnergyConservativeFlux2Param,
       IsmailRoeFluxEC, ChandrashekarFluxEC, RanochaFluxECandKEP,
       L2L4ConservativeFlux, L2L2sConservativeFlux
export FluxPlusDissipation, LocalLaxFriedrichsDissipation, ScalarDissipation,
       MatrixDissipation

# flux difference methods
export UniformFluxDiffDisc1D, UniformPeriodicFluxDiffDisc1D,
        UniformPeriodicFluxDiffDisc2D, UniformPeriodicFluxDiffDisc3D

# utilities
export logmean, integrate, order

end # module


text:
using ReTest, AdvancedHMC, ForwardDiff, Zygote

include("common.jl")

@testset "contrib" begin
    @testset "ad" begin
        metric = UnitEuclideanMetric(D)
        h_hand = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
        h_forwarddiff = Hamiltonian(metric, ℓπ, ForwardDiff)
        h_zygote = Hamiltonian(metric, ℓπ, Zygote)
        for x in [rand(D), rand(D, 10)]
            v_hand, g_hand = h_hand.∂ℓπ∂θ(x)
            v_forwarddiff, g_forwarddiff = h_forwarddiff.∂ℓπ∂θ(x)
            v_zygote, g_zygote = h_hand.∂ℓπ∂θ(x)
            @test v_hand ≈ v_forwarddiff
            @test v_hand ≈ v_zygote
            @test g_hand ≈ g_forwarddiff
            @test g_hand ≈ g_forwarddiff
        end
    end
end


text:
include("../../validation/loocv.jl")

"""
Computes logdet hat(ΣθinvX) given logdet(ΣθinvX), where the hat indicates that
a new point has been incorporated.
"""
function extendlogdetΣθ(old_ΣθinvX, new_choleskyΣθ, new_X)
    lin_sys_loocv()
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using CUDA
using GNSSSignals
using Random
using StaticArrays
using Statistics
using StructArrays
using Test
using Tracking
using TrackingLoopFilters

using Acquisition: AcquisitionResults
using Unitful: MHz, kHz, Hz, s, ms, dBHz

include("code_replica.jl")
include("carrier_replica.jl")
include("downconvert.jl")
include("discriminators.jl")
include("gps_l1.jl")
include("gps_l5.jl")
include("galileo_e1b.jl")
include("boc.jl")
include("secondary_code_or_bit_detector.jl")
include("bit_buffer.jl")
include("correlator.jl")
include("tracking_state.jl")
include("tracking_results.jl")
include("tracking_loop.jl")
include("cn0_estimation.jl")

if CUDA.functional()
    include("cuda/bit_buffer.jl")
    include("cuda/boc.jl")
    include("cuda/cn0_estimation.jl")
    include("cuda/discriminators.jl")
    include("cuda/galileo_e1b.jl")
    include("cuda/gps_l1.jl")
    include("cuda/gps_l5.jl")
    include("cuda/secondary_code_or_bit_detector.jl")
    include("cuda/tracking_loop.jl")
    include("cuda/tracking_results.jl")
    include("cuda/tracking_state.jl")
end


text:
function convergence_parameters(p::Param , q0::Int)

    #Base number of nodes
    N0 = 2^q0 + 1

    #Create parameter list
    sim_param = []

    for q in 1:3

        #Number of nodes
        Nx = (2^(q-1))*(N0-1) + 1

        #Name of file
        filename    = string( p.fname , "_Nx_", Nx , "_Ny_", Nx , "_tmax_" , p.tmax)
        #Create simulation parameters
        s   = Param( p ,  xnodes = Nx , ynodes = Nx , fname = filename )

        #Add to list
        push!(sim_param , s)
    end

    return sim_param

end

function test_convergence( p::Param , q0::Int )

    #Create parameter list
    sim_param = convergence_parameters(p, q0)

    for x in sim_param
        solve_wave_equation_2D(x)
    end

end

function get_data_costum(p::Param , itt::Int , q0::Int, q::Int)

    #Get data at a given time
    ttt , xxx , yyy , zzz  , _ =  get_fields(p , itt)

    N0 = 2^q0 + 1

    #Obtain common nodes
    x_out   = [ xxx[ j + ( 2^(q-q0) - 1 )*(j-1)  ] for j in 1:N0  ]
    y_out   = [ yyy[ j + ( 2^(q-q0) - 1 )*(j-1)  ] for j in 1:N0  ]

    ψ = reshape(zzz , (p.xnodes , p.ynodes) )
    ψ_out  = zeros(N0,N0)

    for i in 1:N0 , j in 1:N0
        ψ_out[i,j]  = ψ[ i + ( 2^(q-q0) - 1 )*(i-1)  , j + ( 2^(q-q0) - 1 )*(j-1) ]
    end

    return x_out , y_out , ψ_out
end

function explore_convergence(p::Param , q0::Int)

    sim_param = convergence_parameters(p , q0)

    #Max iteration
    max_iter , _ = get_t_max(sim_param[1])

    #Create Figure
    fig = Figure( resolution = (1000 , 950), font  = "Helvetica")

    ax1 = fig[4,1]      = Axis(fig  , width = 600 , height = 600 )
    ax2 = fig[4,2]      = Axis(fig  , width = 200 , height = 600 , xlabelvisible = :false , yaxisposition = :right)
    ax3 = fig[1:3,1]    = Axis(fig , width = 600 , height = 200 , yticksvisible = true , xaxisposition = :top)

    linkyaxes!(ax1, ax2)
    linkxaxes!(ax1, ax3)

    slider1 = Slider( fig[1,2] , range = 1:1:max_iter , height = 200/3 , width = 200)
    slider2 = Slider( fig[2,2] , range = 1:1:sim_param[1].xnodes  , height = 200/3 ,  width = 200)
    slider3 = Slider( fig[3,2] , range = 1:1:sim_param[1].ynodes  , height = 200/3 , width = 200)

    #Data for plotting

    x , y , _ = get_data_costum(sim_param[1] , 1 , q0 , q0 )

    #Get fields at a certain iteration that is controled by slider 1
    ψ1 =  lift(slider1.value) do v
        get_data_costum(sim_param[1] , v[] , q0 ,q0 )[3]
#            get_fields( sim_param[1] , v[] )[4]
        end

    ψ2 =  lift(slider1.value) do v
        get_data_costum(sim_param[2] , v[] , q0 ,q0+1 )[3]
        end

    ψ3 =  lift(slider1.value) do v
        get_data_costum(sim_param[3] , v[] , q0 ,q0+2 )[3]
        end

    #Get slices and differences for convergence analysis and scale them
    scale_factor = 100
    ψ1_x_slice = lift(slider2.value) do index
        scale_factor*(ψ1[][index,:] -ψ2[][index,:])
    end
    ψ2_x_slice = lift(slider2.value) do index
        scale_factor*2*(ψ2[][index,:] - ψ3[][index,:])
    end

    ψ3_x_slice = lift(slider2.value) do index
        scale_factor*4*(ψ2[][index,:] - ψ3[][index,:])
    end

    ψ4_x_slice = lift(slider2.value) do index
        scale_factor*16*(ψ2[][index,:] - ψ3[][index,:])
    end

    scale_factor = 100
    ψ1_y_slice = lift(slider3.value) do index
        scale_factor*(ψ1[][:,index] -ψ2[][:,index])
    end
    ψ2_y_slice = lift(slider3.value) do index
        scale_factor*2*(ψ2[][:,index] - ψ3[][:,index])
    end

    ψ3_y_slice = lift(slider3.value) do index
        scale_factor*4*(ψ2[][:,index] - ψ3[][:,index])
    end

    ψ4_y_slice = lift(slider3.value) do index
        scale_factor*16*(ψ2[][:,index] - ψ3[][:,index])
    end

    #Get vertical and horizontal lines to see cuts on histogram
    x_coord = lift(slider2.value) do index
        [ x[index] for i in 1:(2^q0+1) ]
    end

    y_coord = lift(slider3.value) do index
        [ y[index] for i in 1:(2^q0+1) ]
    end

    #y3 = lift(slider3.value) do index
    #    y[][:,index]
    #end

    #Plot histogram of lowest resolution and slices
    heatmap!(ax1 , x , y ,ψ1 ,colorrange=(-5,5) , colormap = :vik )
    lines!(ax1 , x_coord  , x )
    lines!(ax1 , x  , y_coord )

    #Plot x-slices with convergence analysis
    l1 = lines!(ax2 , ψ1_x_slice  , x,  linewidth  = 5.0  ,color = RGBA{Float32}(0.0f0, 0.8f0 , 0.9f0 , 1.0f0)) # color = ψ1_x_slice, colormap = :vik)
    l2 = lines!(ax2 , ψ2_x_slice  , x , color = :black , linestyle = :solid )
    l3 = lines!(ax2 , ψ3_x_slice  , x , color = :red   , linestyle = :solid)
    l4 = lines!(ax2 , ψ4_x_slice  , x , color = :green , linestyle = :solid)
    axislegend(ax1, [l1, l2 ,l3 ,l4], ["ψ1-ψ2", "2(ψ2-ψ3)" , "4(ψ2-ψ3)" , "16(ψ2-ψ3)"], position = :rb,
    orientation = :horizontal)

    #Plot y slices
    lines!(ax3 , y , ψ1_y_slice  ,  linewidth  = 5.0  ,color = RGBA{Float32}(0.0f0, 0.8f0 , 0.9f0 , 1.0f0)) # color = ψ1_x_slice, colormap = :vik)
    lines!(ax3 , y , ψ2_y_slice  ,  color = :black , linestyle = :solid )
    lines!(ax3 , y , ψ3_y_slice  ,  color = :red   , linestyle = :solid)
    lines!(ax3 , y , ψ4_y_slice  ,  color = :green , linestyle = :solid)

    #Adjust limits of heatmap
    xlims!(ax1 , p.xmin ,p.xmax)
    ylims!(ax1 , p.ymin ,p.ymax)

    #Adjust limits of vertical
    xlims!(ax2 , -1 ,1)
    ylims!(ax3 , -1 ,1)

    #Do heatmap

    display(fig)

end

function see_energy_covergence(p::Param , q0::Int)

    sim_param = convergence_parameters(p , q0)

    #Max iteration
    max_iter , _ = get_t_max(sim_param[2])

    #Create Figure
    fig = Figure( resolution = (1000 , 950), font  = "Helvetica")
    ax = fig[1,1] = Axis(fig , title = "Test")

    ls = []
cols = [:red , :blue , :black]
    for (ix , x) in enumerate(sim_param)
        t , x = get_energy(x)
        push!(ls , lines!(ax , t , x  , color = cols[ix]) )
    end

    axislegend(ax, ls, ["1" , "2" , "3"], position = :rb)

     #Adjust limits of heatmap
     ylims!(ax , 0,10)
    display(fig)

end

function plot_energy(p::Param , savetrue::Bool , name::String)

    #Max iteration
    max_iter , _ = get_t_max(p)

    #Create Figure
    fig = Figure( resolution = (800 , 450), font  = "Helvetica")
    ax = fig[1,1] = Axis(fig)

    t , x = get_energy(p)

    lines!(ax , t , x )

    ax.ylabel = "Energy Density "
    ax.xlabel = "Time"

    ylims!(ax , 0,10)
    display(fig)

    if(savetrue)
        save( name , fig )
    end

end


text:
module Test_emmsyn

using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Random
using Test

# Testing examples for AFDSYN
println("Test_emmsyn")
@testset "emmsyn" begin
#rand()

##
p = 1; mf = 0; mw = 0;
sysf = fdimodset(rss(1,p,mf+mw),faults = 1:mf,noise = mf.+Vector(1:mw))
sysr = FDFilterIF(sysf.sys,0,0,mf,mw)

# solve an EMMP for a single reference model
@time Q, R, info = emmsyn(sysf,sysr; minimal = false); info
@test iszero(R.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,sysr) == 0

# solve an EMMP for a bank of two reference models
sysr = FDIFilterIF([sysf.sys,sysf.sys],0,0,mf,mw)
@time Q, R, info = emmsyn(sysf,sysr; minimal = false); info
@test iszero(vcat(R.sys...)-vcat(Q.sys...)*sysf.sys, atol = 1.e-7) && fdimmperf(R,sysr) == [0,0]

##
p = 2; mf = 2;
sysf = fdimodset(rss(1,p,mf,stable = true),faults = 1:mf)
sysr = FDFilterIF(sysf.sys,0,0,mf)

# solve an EMMP with sysr = sysf: solution Q = I
@time Q, Rf, info = emmsyn(sysf,sysr,atol = 1.e-7); info
@test iszero(Rf.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(Rf,info.M*sysr) == 0

# solve an EMMP for a bank of two reference models
sysr = FDIFilterIF([sysf.sys,sysf.sys],0,0,mf)
Q, Rf, info = emmsyn(sysf,sysr,atol = 1.e-7); info
@test iszero(vcat(Rf.sys...)-vcat(Q.sys...)*sysf.sys, atol = 1.e-7) &&
      isapprox(fdimmperf(Rf,info.M*sysr), [0,0], atol = 1.e-7)

##
p = 3; mf = 2;
sysf = fdimodset(rss(5,p,mf,stable=true),faults = 1:mf)

sysr = FDFilterIF(sysf.sys,0,0,mf)

# solve an EMMP with sysr = sysf: solution Q = I
@time Q, Rf, info = emmsyn(sysf, sysr, atol = 1.e-7); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rf.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) < 1.e-7 &&
      iszero(R.sys-info.M*sysr.sys, atol = 1.e-7)

# solve an EMMP with sysr = sysf: solution Q = I
@time Q, Rf, info = emmsyn(sysf, sysr, atol = 1.e-7, minimal = true); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rf.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) < 1.e-7

# solve an EMMP for a bank of two reference models
sysr = FDIFilterIF([sysf.sys,sysf.sys],0,0,mf)
@time Q, Rf, info = emmsyn(sysf,sysr,atol = 1.e-7); info
@test iszero(vcat(Rf.sys...)-vcat(Q.sys...)*sysf.sys, atol = 1.e-7) &&
      isapprox(fdimmperf(Rf,info.M*sysr), [0,0], atol = 1.e-7)

##
p = 3; mf = 2;
sysf = fdimodset(rss(1,p,mf,stable=true),faults = 1:mf)
sysr = FDFilterIF(sysf.sys[1:mf,:],0,0,mf)

# solve an EMMP with sysr = sysf: solution Q = [I 0]
@time Q, Rf, info = emmsyn(sysf, sysr, atol = 1.e-7, minimal = false); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rf.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) < 1.e-7 &&
      iszero(R.sys-info.M*sysr.sys, atol = 1.e-7)

# solve an EMMP using the computed design matrix
@time Q1, Rf1, info1 = emmsyn(sysf,sysr, atol = 1.e-7, minimal = false, HDesign = info.HDesign);
@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

# solve an EMMP with sysr = sysf: solution Q = [I 0]
@time Q, Rf, info = emmsyn(sysf, sysr, atol = 1.e-7, minimal = true); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rf.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) < 1.e-7 &&
      iszero(R.sys-info.M*sysr.sys, atol = 1.e-7) && ismissing(info.HDesign)

@time Q1, Rf1, info1 = emmsyn(sysf,sysr, atol = 1.e-7, minimal = true, HDesign = eye(2,3));
@test iszero(Q.sys-Q1.sys,atol=1.e-7) && iszero(Rf.sys-Rf1.sys,atol=1.e-7)

##
p = 3; mf = 2;
sysf = fdimodset(rss(1,p,mf,stable=true),faults = 1:mf)
sysr = FDFilterIF(sysf.sys[1:mf,:],0,0,mf)

# solve an EMMP with sysr = sysf: solution Q = [I 0]
@time Q, Rf, info = emmsyn(sysf, sysr, atol = 1.e-7, minimal = false, simple = true); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rf.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) < 1.e-5 &&
      iszero(R.sys-info.M*sysr.sys, atol = 1.e-7)

##
p = 3; mf = 2; mw = 2;
sysf = fdimodset(rss(1,p,mf+mw,stable=true),faults = 1:mf,noise = mf.+Vector(1:mw))
sysr = FDFilterIF(sysf.sys[1:mf,sysf.faults],0,0,mf)

# solve an EMMP with reference model without noise input
@time Q, Rfw, info = emmsyn(sysf,sysr; atol = 1.e-7, minimal = false); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rfw.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) ≈ fdimmperf(R)

sysr = FDFilterIF(sysf.sys[1:mf,:],0,0,mf,mw)

# solve an EMMP with reference model with noise input
@time Q, Rfw, info = emmsyn(sysf,sysr; atol = 1.e-7, minimal = false); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rfw.sys-Q.sys*sysf.sys, atol = 1.e-7) && fdimmperf(R,info.M*sysr) < 1.e-7

##
p = 3; mu = 2; md = 1; mf = 2; mw = 0;
m = mu+md+mf+mw;
sysf = fdimodset(rss(3,p,m), c = 1:mu,d = mu .+ (1:md),f = (mu+md) .+ (1:mf), n = (mu+md+mf) .+ (1:mw))
sysr = fdimodset(dss(zeros(3,mu+md)),c = 1:mu,d = mu .+ (1:md))

@time Q, Rfw, info = emmsyn(sysf,sysr; atol = 1.e-7); info
R = fdIFeval(Q, sysf; atol = 1.e-7, minimal = true);
@test iszero(Rfw.sys-Q.sys*[sysf.sys[:,[sysf.controls;sysf.disturbances]]; eye(mu,mu+md)], atol = 1.e-7) &&
       fdimmperf(R,info.M*sysr) < 1.e-7

# Example 5.12c - Solution of an EMMP

# enter output and fault vector dimensions
p = 3; mf = 3;
# generate random dimensions for system order and input vectors
nu = Int(floor(1+4*rand())); mu = Int(floor(1+4*rand()));
nd = Int(floor(1+4*rand())); md = Int(floor(1+4*rand()));

# define random Gu(s) and Gd(s) with triplex sensor redundancy
# and Gf(s) for triplex sensor faults
Gu = ones(3,1)*rss(nu,1,mu); # enter Gu(s) in state-space form
Gd = ones(3,1)*rss(nd,1,md); # enter Gd(s) in state-space form
Gf = eye(3);                 # enter Gf(s) for sensor faults
atol = 1.e-7;                # tolerance for rank tests

# build model with faults
#sysf = [Gu Gd Gf]
sysf = fdimodset([Gu Gd Gf], c = 1:mu, d = mu.+(1:md), f = (mu+md).+(1:mf))

# enter reference model
Mr = FDFilterIF(dss([ 0 1 -1; -1 0 1; 1 -1 0]), 0, 0, mf)

# solve an exact model-matching problem using EMMSYN
@time Q, R, info = emmsyn(sysf,Mr; atol); info

# one step solution
# solve Qbar*Ge = Me, where Ge = [Gu Gd Gf; I 0 0] and Me = [0 0 Mr ].
Ge = [sysf.sys; eye(mu,mu+md+mf)]; Me = [zeros(p,mu+md) Mr.sys];
Qbar = glsol(Ge, Me; atol)[1]

# compare solutions
@test iszero(Q.sys-Qbar; atol)

# solve an EMMP without explicit nullspace computation
Mre = FDFilterIF([zeros(mf,mu) Mr.sys],mu,0,mf);

@time Q1, R1, info1 = emmsyn(sysf,Mre; atol); info1
Rt1 = fdIFeval(Q1, sysf; atol = 1.e-7, minimal = true);

@test iszero(Q.sys-Q1.sys; atol) && iszero(R1.sys[:,R1.faults]-R.sys[:,R.faults]; atol)

# Example 5.13c - Solution of an EMMP using EMMSYN

# define s as an improper transfer function
s = rtf('s');
# enter Gu(s), Gf(s) and Mr(s)
Gu = [s/(s^2+3*s+2) 1/(s+2);
     s/(s+1) 0;
      0 1/(s+2)];
Mr = dss(eye(2));                  # enter Mr(s)
p, mu = size(Gu); mf = mu

# build the synthesis model with additive faults
sysf = fdimodset(dss(Gu), c = 1:mu, f = 1:mu);

# enter reference model
Mr = fdimodset(dss(eye(mf)), f = 1:mf)

atol = 1.e-7                # tolerance for rank tests
sdeg = -1                   # set stability degree

# solve an exact model-matching problem using EMMSYN
@time Q, R, info = emmsyn(sysf, Mr; atol, sdeg, minimal = false); info

# check solution
G = [sysf.sys; eye(mu,mu+mf)]; F = [zeros(mf,mu) info.M*Mr.sys];
@test iszero(Q.sys*G-F; atol)

# solve an exact model-matching problem using EMMSYN
@time Q1, R1, info1 = emmsyn(sysf, Mr; atol, sdeg, nullspace = false, minimal = false); info

# check solution
G = [sysf.sys; eye(mu,mu+mf)]; F = [zeros(mf,mu) info1.M*Mr.sys];
@test iszero(Q1.sys*G-F; atol)

# Example 5.4c - Solution of an EFDP using EFDSYN
# define s as an improper transfer function
s = rtf('s');
# define Gu(s) and Gd(s)
Gu = [(s+1)/(s-2); (s+2)/(s-3)];     # enter Gu(s)
Gd = [(s-1)/(s+2); 0];               # enter Gd(s)
Gf = [(s+1)/(s-2) 0; (s+2)/(s-3) 1]; # enter Gf(s)
p = 2; mu = 1; md = 1; mf = 2;       # set dimensions

# setup the synthesis model with faults with Gf(s) = [ Gu(s) [0;1]]
sysf = fdimodset(dss([Gu Gd]),c =1,d = 2,f = 1,fs = 2);

# call of EFDSYN with the options for stability degree -3 and the synthesis
# of a scalar output filter
Qt, Rft, info = efdsyn(sysf, sdeg = -3, rdim = 1);

# normalize Q and Rf to match example
scale = evalfr(Rft.sys[1,1],Inf)[1,1];
Q = FDFilter(Qt.sys/scale,Qt.outputs,Qt.controls);
Rf = FDFilterIF(Rft.sys/scale,faults = Rft.faults);

#  solve an EMMP
@time QM, RM, info = emmsyn(sysf,Rf)
@test iszero(Q.sys-QM.sys,atol=1.e-7) && iszero(Rf.sys-RM.sys,atol=1.e-7)

## Model Niemann 1998, Optim. Appl. Meth.
n = 5; mu = 2; md = 1; mf = 2; m = mu+md+mf; p = 5;
a = [
-0.0782 0.2939 0.0220 0.0208 -0.0291
0.0077 -0.0278 0.0015 -0.0015 0.0026
-1.3251 5.32608 -0.5263 0.2214 -0.4777
1.0809 -4.4452 0.3770 -0.4631 0.4032
2.1532 -8.6386 0.7811 -0.5745 0.7816]*1000;
bd = [
    0.0443
   -0.0042
    0.7910
   -0.6598
   -1.2881]*1000;
bu = [
    -0.007197 0.003005
0.003473 0.000284
1.218502 -0.032866
1.322502 0.020147
-0.082265 0.024388 ];
bf = bu;
c = eye(n); dd = zeros(n,md); du = zeros(n,mu); df = zeros(n,mf);

sys = dss(a,[bu bd bf],c,[du dd df]);
sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf) # determine achievable structure matrix
nb = size(SFDI,1)
@time Q, Rf = efdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1);
@time Q, Rf = efdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1, sdeg = -5, smarg = -0.05,
                                  FDfreq = 0, FDGainTol = 0.0001);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == ones(Int,nb)

#  solve an EMMP
@time QM, RM, info = emmsyn(sysf,Rf; atol = 1.e-7, minimal = true);
# @test iszero(vcat(Q.sys...)-vcat(QM.sys...),atol=1.e-7) && iszero(Rf.sys-RM.sys,atol=1.e-7)

R1 = fdIFeval(QM,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R1.sys...)[:,[R1.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R1.sys...)[:,R1.faults]-vcat(RM.sys...),atol=1.e-7) &&
      order.(QM.sys) == ones(Int,nb) &&
      iszero(vcat((info.M*Rf).sys...)[:,Rf.faults]-vcat(RM.sys...),atol=1.e-7)

## Example 2, Xiong & Saif 2000, Int. J. Robust Nonlinear Control
n = 5; mu = 3; md = 2; mf = 2; m = mu+md+mf; p = 3;
a = [
0 0 -0.0034 0 0
0 -0.041 0.0013 0 0
0 0 -1.1471 0 0
0 0 -0.0036 0 0
0 0.094 0.0057 0 -0.051];
bu = [
    -1 0 0
     0 0 0
     0 0 0.948
     0.916 -1 0
    -0.598 0 0];
bd = [
     0 1
     0.062 -0.132
     0 -7.189
     0 0
     0 0 ];
bf = bu[:,1:2];
c = [
1 0 0 0 0
0 0 1 0 0
0 0 0 1 0];
dd = zeros(p,md); du = zeros(p,mu); df = zeros(p,mu-1);

sys = dss(a,[bu bd bf],c,[du dd df]);

sysf = fdimodset(sys, c = 1:mu, d = mu .+ (1:md) , f = (mu+md) .+ (1:mf));

SFDI = fdigenspec(sysf) # determine achievable structure matrix
nb = size(SFDI,1)
Q, Rf = efdisyn(sysf, SFDI; atol = 1.e-7, rdim = 1, sdeg =-2, smarg = -1,
                                  FDfreq = 0, FDGainTol = 0.0001);

# check synthesis conditions: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rf
R = fdIFeval(Q,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R.sys...)[:,[R.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R.sys...)[:,R.faults]-vcat(Rf.sys...),atol=1.e-7) &&
      order.(Q.sys) == [1, 1, 1]

# check weak and strong fault detectability
@test fditspec(Rf) == fdisspec(Rf) == SFDI

#  solve an EMMP
@time QM, RM, info = emmsyn(sysf,Rf; atol = 1.e-7, minimal = true);
# @test iszero(vcat(Q.sys...)-vcat(QM.sys...),atol=1.e-7) && iszero(Rf.sys-RM.sys,atol=1.e-7)

R1 = fdIFeval(QM,sysf); # form Q*[Gu Gd Gf;I 0 0];
@test iszero(vcat(R1.sys...)[:,[R1.controls;R.disturbances]],atol=1.e-7) &&
      iszero(vcat(R1.sys...)[:,R1.faults]-vcat(RM.sys...),atol=1.e-7) &&
      order.(QM.sys) == ones(Int,nb) &&
      iszero(vcat((info.M*Rf).sys...)[:,Rf.faults]-vcat(RM.sys...),atol=1.e-7)

end # test

end   # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
#=
    Contents: Example: Display last visit with client-side sessions.
=#
using HttpServer
using AccessControl

# Handler
function home!(req, res)
    session = session_read(req)
    if session == ""                             # "id" cookie does not exist...session hasn't started...start a new session.
        session  = session_create!("")           # username = ""
        res.data = "This is your first visit."
	session["lastvisit"] = string(now())
        session_write!(res, session)
    else
        last_visit = session["lastvisit"]
        res.data   = "Welcome back. Your last visit was at $last_visit."
	session["lastvisit"] = string(now())
        session_write!(res, session)
    end
end

# App
function app(req::Request)
    res = Response()
    if req.resource == "/home"
        home!(req, res)
    else
        notfound!(res)
    end
    res
end

### Run the app under HTTPS
# Generate keys/server.key and keys/server.crt if they don't already exist
rel(filename::AbstractString, p::AbstractString) = joinpath(dirname(filename), p)
if !isfile("keys/server.crt")
    @unix_only begin
	run(`mkdir -p $(rel(@__FILE__, "keys"))`)
	run(`openssl req -x509 -nodes -days 365 -newkey rsa:2048 -keyout $(rel(@__FILE__, "keys/server.key")) -out $(rel(@__FILE__, "keys/server.crt"))`)
    end
end

# Define and run the server
server = Server((req, res) -> app(req))
cert   = MbedTLS.crt_parse_file(rel(@__FILE__, "keys/server.crt"))
key    = MbedTLS.parse_keyfile(rel(@__FILE__, "keys/server.key"))
run(server, port = 8000, ssl = (cert, key))


text:
export Characteristic
"""
A domain that is described by its characteristic function.
"""
struct Characteristic{N,T} <: Domain{SVector{N,T}}
    char    ::  Function
    box     ::  Domain{SVector{N,T}}
end

export characteristic
characteristic(char::Function, dom)  = Characteristic(char,boundingbox(dom))

DomainSets.indomain(x, c::Characteristic) = c.char(x)

boundingbox(c::Characteristic) = c.box

show(io::IO, c::Characteristic) = print(io, "a domain described by a characteristic function")


text:
module GravitationalAcceleration

# This implements the gravitational acceleration

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam
import Base.show, GeoParams.param_info
using ..MaterialParameters: MaterialParamsInfo

abstract type AbstractGravity{_T} <: AbstractMaterialParam end

export  compute_gravity,        # calculation routines
        ConstantGravity,        # constant
        param_info

# Constant Gravity -------------------------------------------------------
"""
    GravityConstant(g=9.81m/s^2)

Set a constant value for the gravitational acceleration:
```math
    g  = 9.81 m s^{-2}
```
"""
@with_kw_noshow struct ConstantGravity{_T,U}   <: AbstractGravity{_T}
    g::GeoUnit{_T,U}              =   9.81m/s^2               # gravitational acceleration
end
ConstantGravity(args...) = ConstantGravity(convert.(GeoUnit,args)...)

function param_info(s::ConstantGravity) # info about the struct
    return MaterialParamsInfo(Equation = L"g = 9.81 m s^{-2}" )
end

# Calculation routine
function compute_gravity(s::ConstantGravity{_T}) where _T
    @unpack_val g   = s

    return g
end

# Print info
function show(io::IO, d::ConstantGravity{_T})  where _T
    print(io, "Gravitational acceleration: g=$(UnitValue(d.g))")
end
#-------------------------------------------------------------------------

# Help info for the calculation routines
"""
compute_gravity(s:<AbstractGravity)

Returns the gravitational acceleration

"""
compute_gravity

end


text:
using BinaryBuilder

name = "Luna"
version = v"0.23.0"

sources = [
    GitSource("https://bitbucket.org/remnrem/luna-base.git",
              "6d333a4034f7022ba4d1aa99ce2f8afddfd6832e"),
    DirectorySource("./bundled"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/luna-base
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/fix_it.patch
mkdir -p ${libdir}
mkdir -p ${bindir}
def_windows=""
if [[ ${target} == *-apple-* ]]; then
    suspicious_arch="MAC"
elif [[ ${target} == *-mingw* ]]; then
    suspicious_arch="WINDOWS"
    def_windows="WINDOWS=1"
else
    suspicious_arch="LINUX"
fi
make -j${nproc} ARCH=${suspicious_arch} FFTW=${prefix} PREFIX=${prefix} LIBDIR=${libdir} ${def_windows}
cp "luna${exeext}" "${bindir}/"
cp "libluna.${dlext}" "${libdir}/"
"""

platforms = expand_cxxstring_abis(supported_platforms())

products = [
    ExecutableProduct("luna", :luna),
    LibraryProduct("libluna", :libluna),
]

dependencies = [
    Dependency("FFTW_jll"),
    Dependency("Zlib_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Autogenerated wrapper script for XyceWrapper_jll for x86_64-linux-gnu-cxx03-julia_version+1.6.0
export xycelib

using Xyce_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("XyceWrapper")
JLLWrappers.@declare_library_product(xycelib, "libxycelib.so")
function __init__()
    JLLWrappers.@generate_init_header(Xyce_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        xycelib,
        "lib/libxycelib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DensestCommonSubgraph
using FactCheck
using Gadfly
using MAT
using DataFrames
using ProgressMeter
FactCheck.onlystats(true)
FactCheck.setstyle(:default)
using HDF5, JLD

function test_instances2(m::Integer, n::Integer, p::Real, k::Integer)
  G = Graph[]
  for i=1:m
    cG = planted_clique(n, p, k)
    push!(G , cG)
  end
  return G
end

function test_dcs_lp(nt::Integer=10 )
  m = 1
  Z = zeros(nt, m)
  D = zeros(nt, 4)  # Density of extracted subgraph
  V = zeros(nt, 4)  # LP optimum
  S = cell(nt, 4)   # subgraphs extracted
  I = cell(nt, 1)   # graph instances
  X = DataFrame(GID=Int64[] , RelDensity=Float64[], RelOptVal=Float64[], Jaccard=Float64[], Method=String[])
  p_bar = Progress(nt, 1, "Computing initial pass...", 50)
  for i=1:nt
    n = 20
    p = 0.5
    k = rand(1:ifloor(n/2))
    G = test_instances2(m, n, p, k);
    I[i] = G
    (s_true, d_val ) = dcs_parallel(G, 2);
    d_true = common_subgraph_density(G, s_true)
    for jj=1:m
      Z[i, jj] = subgraph_density(G[jj], s_true)
      println("i: $i j: $jj $(Z[i, jj])")
    end

    # @fact d_true => roughly(d_val/2.0)
    # println("OPT: density: $(lpad(d_true, 20, ' ')) subgraph: $s_true")

    D[i, 1] = d_true
    V[i, 1] = d_val
    S[i, 1] = s_true
    (p_dcs_lp, X_dcs, Y_dcs, t_dcs) = dcs_lp(G);
    s_dcs = extract_primal_solution(Y_dcs)
    d_dcs = common_subgraph_density(G, s_dcs)
    # println("DCS: density: $(lpad(d_dcs, 20, ' ')) subgraph: $s_dcs")

    D[i, 2] = d_dcs
    V[i, 2] = t_dcs
    S[i, 2] = s_dcs
    if t_dcs > d_val
      println("i: $i integrality gap: dcs: $t_dcs opt: $d_val  ratio: $(t_dcs/d_val)")
    end
    X = vcat(X, DataFrame(GID=i, RelDensity = d_dcs/d_true , RelOptVal=t_dcs/d_val, Jaccard=length(intersect(s_dcs , s_true))/length(union(s_dcs, s_true)), Method="DCS"))

    # @fact d_dcs => less_than_or_equal(d_true)

    (p_charikar_lp, Xc, Yc, t_c) = charikar_lp(dcs_intersection_graph(G))
    s_c = extract_primal_solution(Yc)
    d_c = common_subgraph_density(G, s_c)

    D[i, 3] = d_c
    V[i, 3] = t_c
    S[i, 3] = s_c

    X = vcat(X, DataFrame(GID=i, RelDensity = d_c/d_true , RelOptVal=t_c/d_val, Jaccard=length(intersect(s_c , s_true))/length(union(s_c, s_true)), Method="MIN"))

    # println("CHARIKAR: density: $(lpad(d_c, 20, ' ')) subgraph: $s_c")

    @fact d_c => less_than_or_equal(d_true)

    (p_avg_lp, Xavg, Yavg, t_avg) = charikar_lp(dcs_avg_graph(G))
    s_avg = extract_primal_solution(Yavg)
    d_avg = common_subgraph_density(G, s_avg)
    # println("AVG: density: $(lpad(d_avg, 20, ' ')) subgraph: $s_avg")

    D[i, 4] = d_avg
    V[i, 4] = t_avg
    S[i, 4] = s_avg

    X = vcat(X, DataFrame(GID=i, RelDensity = d_avg/d_true , RelOptVal=t_avg/d_val, Jaccard=length(intersect(s_avg , s_true))/length(union(s_avg, s_true)), Method="AVG"))

    # @fact d_avg => less_than_or_equal(d_true)
    println("$i/$nt opt: $(rpad(d_true, 20, ' ')) DCS: $(rpad(d_dcs, 20, ' ')) CHARIKAR: $(rpad(d_c, 20, ' ')) AVG: $(p_avg_lp.optval)\n")
    # println("\n")
    next!(p_bar)
  end
  return (I, X, S, Z)
end

# function evaluate_primal(G::GraphVec, methods::String[]= LP_METHODS)
#   m = length(G)
#   @assert m >= 1
#   n = size(G, 1)
#   S = Array{Int64, 1}[]
#   X = DataFrame(GID=Int64[] , RelDensity=Float64[], RelOptVal=Float64[], Jaccard=Float64[], Method=String[])
#   (s_true, d_val) = dcs_parallel(G, 2);
#   push!(S, s_true)
#   d_true = common_subgraph_density(G, s_true)
#   for m=methods
#     (cS, d, t) = dcs_lp_extract(G; method = m)
#     push!(S, cS)
#     X = vcat(X, DataFrame(GID=gid, RelDensity = d/d_true , RelOptVal=t/d_val, Jaccard=length(intersect(S , s_true))/length(union(cS, s_true)), Method=m))
#   end
#   return
# end

nt = 100;
(I, X, S, Z) = test_dcs_lp(nt)
legend = "I: graph instances, S: subgraphs returned by OPT,DCS,LP_min,LP_avg, X: stats, Z: densities by S_opt"
using HDF5, JLD
@save "output_test_dcs_lp_2.jld" I X S Z legend
p1 = plot(X, xgroup="Method", x="RelOptVal", Geom.subplot_grid(Geom.histogram), Guide.xlabel(nothing), Guide.title("Ratio of LP optimum to density of S_{OPT}"), color="Method")
p2 = plot(X, xgroup="Method", x="RelDensity", Geom.subplot_grid(Geom.histogram), Guide.xlabel(nothing), Guide.title("Ratio of densities of subgraph obtained by LP rounding to S_{OPT}"), color="Method")
draw(PDF("p1a.pdf", 9inch, 6inch), p1)
draw(PDF("p2a.pdf", 9inch, 6inch), p2)


text:
"""
integrate a TaylorNModel with respect to the variable number `which`.
Optionally adds `x0` to the result.
"""
function integrate(f::TaylorNModel, which=1, x0=0)

    p = integrate(f.p, which)  # not necessary if an already complete Taylor series, in which case p2 == f.p

    Δ = integral_bound(f, which)

    g = TaylorNModel(f.n, f.x0, f.I, p, Δ)

    g.p[0] += x0  # constant term

    # for k in 0:(x0.order)
    #     TaylorSeries.add!(g.p, g.p, x0, k)
    # end

    return g

end

doc"""
Bound the integral of a `TaylorNModel` `f` with respect to the variable `which`.
"""
function integral_bound(f::TaylorNModel, which)

    high_order_term = f.p[end]  # a HomogeneousPolynomial

    Δ = ( bound(high_order_term, f.x0, f.I) + f.Δ ) * diam(f.I[which])

    return Δ
end

bound(f::HomogeneousPolynomial, x0, I) = f( [(I - x0)...] )
# applies the hom poly to the bounds


text:
# Translating
# https://proplot.readthedocs.io/en/latest/basics.html#Formatting-subplots
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

f, axs = proplot.subplots(ncols=2, nrows=2, share=0, tight=true, axwidth=1.7)
axs[1].plot(range(1, 10, length=80), cumsum(rand(80, 5) .- 0.5, dims=1))
[ax.format(
    suptitle="Format command demo",
    abc=true, abcloc="ul", abcstyle="a.",
    title="Main", ltitle="Left", rtitle="Right",  # different titles
    urtitle="Title A", lltitle="Title B", lrtitle="Title C",  # extra titles
    collabels=["Column label 1", "Column label 2"], rowlabels=["Row label 1", "Row label 2"],
    xlabel="x-axis", ylabel="y-axis",
    xlim=(1, 10), xticks=1, xscale="log",
    ylim=(-2, 2), yticks=proplot.arange(-2, 2), yticklabels=("a", "bb", "c", "dd", "e"),
    xtickdir="inout", xtickminor=false,
    ytickloc="both", yticklabelloc="both", ygridminor=true,
    linewidth=0.8, gridlinewidth=0.8, gridminorlinewidth=0.5,
) for ax in axs]
f.savefig("figs/Formatting-subplots.svg", transparent=false)

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
## How do you do ?
a=10
b=20
c=a+b
println(c)
function sum(x,y)
    return x+y
end


text:
module Test_SignalTorque

using  Modia3D
import Modia3D.ModiaMath

@signal Signal(;A=10.0) begin
   y1 = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Signal, sim::ModiaMath.SimulationState)
    signal.y1.value  = signal.A
end

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly PendulumDrivenKinematically(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx, m=1.0) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))

   # Pendulum
   body   = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(body; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))

   rev    = Modia3D.Revolute(world, frame1)

   sig    = Signal()
   signal = Modia3D.SignalToFlangeTorque(sig.y1)
   Modia3D.connect(signal, rev)
end

gravField = Modia3D.UniformGravityField(g=0.0)
pendulum = PendulumDrivenKinematically(Lx=1.6, m=0.5, sceneOptions = Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.3))
model    = Modia3D.SimulationModel(pendulum)
# ModiaMath.print_ModelVariables(model)
result   = ModiaMath.simulate!(model, stopTime=0.5)
ModiaMath.plot(result, ["sig.y1", "rev.phi", "rev.tau"] )

println("... success of Test_SignalTorque.jl!")
end


text:
export clipnorm!, clipvalue!

function maxnorm(x::Var, threshold::Float64)
    norm = norm(x.data, p=2)
    if norm <= threshold
        data = 0
    else

    end
    Var
end

function clipnorm!(x::Array{T}, threshold) where {T}
    (ndims(x) == 1 || ndims(x) == 2) || throw("Not supported yet.")
    z = mapreducedim(v -> v*v, +, x, 1)
    for j = 1:length(z)
        norm = sqrt(z[j])
        z[j] = norm <= T(threshold) ? 1 : T(threshold)/norm
    end
    x .*= z
end

function clipvalue!{T}(x::Array{T}, value)
    value = T(value)
    @inbounds for i = 1:length(x)
        x[i] = min(max(x[i],-value), value)
    end
end


text:
using Bokeh

plotfile("display.html")

x = linspace(0, 2pi)

# plot returns a Plot object
myplot = plot(x, sin.(x))
# which can then be displayed with
showplot(myplot)

# if you wish to generate a plot but not display it, you can use genplot
genplot(myplot)

# you can also specify the filename directly to genplot:
genplot(myplot, "display1.html")

# underneither genplot there's renderplot which returns the actual html
html = renderplot(myplot)
println(html[1:100], "...")

# however bokeh.jl also keeps a reference to the current plot, which you can access
plot(x, cos.(x))
showplot(curplot())

# for convienience there's a shorthand for this:
showplot()

# Note: in this example, we've generated display.html and opened it in the browser
# 3 times, however probably all three plots will show the final plot,
# that's because the browser is too slow for julia which has generated the last version of
# display.html before the browser has time to opening it the first two times
# (or so I assume????)


text:
using Base.Test
using ComparativeAutograder: calculate_error, TestCaseResult

t1 = TestCaseResult(
    3,
    nothing, nothing,
    0.0, "", ""
)

t2 = TestCaseResult(
    4,
    nothing, nothing,
    0.0, "", ""
)

@test calculate_error(t1, t2) == 1


text:
# simple Julia APIs over Kubernetes Swagger interface

function sel(label::String, op::Symbol)
    @assert op === :exists
    label
end
sel(label::String, op::Symbol, items::String...) = label * " " * string(op) * " (" * join(items, ",") * ")"
sel(cnd::String...) = join(cnd, ", ")

_delopts(; kwargs...) = Typedefs.MetaV1.DeleteOptions(; preconditions=Typedefs.MetaV1.Preconditions(; kwargs...), kwargs...)

function _get_apictx(ctx::KuberContext, O::Symbol, apiversion::Union{String,Nothing})
    if apiversion !== nothing
        k = Kuber.api_group_type(apiversion)
        apictx = k(ctx.client)
    else
        kapi = ctx.modelapi[O]
        apictx = kapi.api(ctx.client)
    end
    apictx
end

function list(ctx::KuberContext, O::Symbol, name::String, apiversion::Union{String,Nothing}=nothing; namespace::Union{String,Nothing}=ctx.namespace, kwargs...)
    isempty(ctx.apis) && set_api_versions!(ctx)

    apictx = _get_apictx(ctx, O, apiversion)
    namespaced = (namespace !== nothing) && !isempty(namespace)
    allnamespaces = namespaced && (namespace == "*")

    if allnamespaces
        apicall = eval(Symbol("list$(O)ForAllNamespaces"))
        return apicall(apictx, name; kwargs...)
    elseif namespaced
        apicall = eval(Symbol("listNamespaced$O"))
        return apicall(apictx, name, namespace; kwargs...)
    else
        apicall = eval(Symbol("list$O"))
        return apicall(apictx, name; kwargs...)
    end
end

function list(ctx::KuberContext, O::Symbol, apiversion::Union{String,Nothing}=nothing; namespace::Union{String,Nothing}=ctx.namespace, kwargs...)
    isempty(ctx.apis) && set_api_versions!(ctx)

    apictx = _get_apictx(ctx, O, apiversion)
    namespaced = (namespace !== nothing) && !isempty(namespace)
    allnamespaces = namespaced && (namespace == "*")

    if allnamespaces
        apicall = eval(Symbol("list$(O)ForAllNamespaces"))
        return apicall(apictx; kwargs...)
    elseif namespaced
        apicall = eval(Symbol("listNamespaced$O"))
        return apicall(apictx, namespace; kwargs...)
    else
        apicall = eval(Symbol("list$O"))
        return apicall(apictx; kwargs...)
    end
end

function get(ctx::KuberContext, O::Symbol, name::String, apiversion::Union{String,Nothing}=nothing; max_tries::Integer=1, kwargs...)
    isempty(ctx.apis) && set_api_versions!(ctx; max_tries=max_tries)

    apictx = _get_apictx(ctx, O, apiversion)
    try
        apicall = eval(Symbol("read$O"))
        @repeat max_tries try
            return apicall(apictx, name; kwargs...)
        catch e
            @retry if isa(e, IOError)
                @debug("Retrying ", "read$O")
                sleep(2)
            end
        end
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("readNamespaced$O"))
        @repeat max_tries try
            return apicall(apictx, name, ctx.namespace; kwargs...)
        catch e
            @retry if isa(e, IOError)
                @debug("Retrying ", "readNamespaced$O")
                sleep(2)
            end
        end
    end
end

function get(ctx::KuberContext, O::Symbol, apiversion::Union{String,Nothing}=nothing; label_selector=nothing, namespace::Union{String,Nothing}=ctx.namespace, max_tries::Integer=1)
    isempty(ctx.apis) && set_api_versions!(ctx; max_tries=max_tries)

    apictx = _get_apictx(ctx, O, apiversion)
    try
        apiname = "list$O"
        (namespace === nothing) && (apiname *= "ForAllNamespaces")
        apicall = eval(Symbol(apiname))
        @repeat max_tries try
            return apicall(apictx; labelSelector=label_selector)
        catch e
            @retry if isa(e, IOError)
                @debug("Retrying ", apiname)
                sleep(2)
            end
        end
    catch ex
        isa(ex, UndefVarError) || rethrow()
        (namespace === nothing) && rethrow()
        apicall = eval(Symbol("listNamespaced$O"))
        @repeat max_tries try
            return apicall(apictx, namespace; labelSelector=label_selector)
        catch e
            @retry if isa(e, IOError)
                @debug("Retrying ", "listNamespaced$O")
                sleep(2)
            end
        end
    end
end

function put!(ctx::KuberContext, v::T) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    put!(ctx, Symbol(vjson["kind"]), vjson)
end

function put!(ctx::KuberContext, O::Symbol, d::Dict{String,Any})
    isempty(ctx.apis) && set_api_versions!(ctx)

    apictx = _get_apictx(ctx, O, get(d, "apiVersion", nothing))
    try
        apicall = eval(Symbol("create$O"))
        return apicall(apictx, d)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("createNamespaced$O"))
        return apicall(apictx, ctx.namespace, d)
    end
end

function delete!(ctx::KuberContext, v::T; kwargs...) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    kind = vjson["kind"]
    name = vjson["metadata"]["name"]
    delete!(ctx, Symbol(kind), name, get(vjson, "apiVersion", nothing); kwargs...)
end

function delete!(ctx::KuberContext, O::Symbol, name::String, apiversion::Union{String,Nothing}=nothing; kwargs...)
    isempty(ctx.apis) && set_api_versions!(ctx)
    apictx = _get_apictx(ctx, O, apiversion)

    params = [apictx, name]

    try
        apicall = eval(Symbol("delete$O"))
        return apicall(params...; kwargs...)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("deleteNamespaced$O"))
        push!(params, ctx.namespace)
        return apicall(params...; kwargs...)
    end
end

function update!(ctx::KuberContext, v::T, patch, patch_type) where {T<:SwaggerModel}
    vjson = convert(Dict{String,Any}, v)
    kind = vjson["kind"]
    name = vjson["metadata"]["name"]
    update!(ctx, Symbol(kind), name, patch, patch_type, get(vjson, "apiVersion", nothing))
end

function update!(ctx::KuberContext, O::Symbol, name::String, patch, patch_type, apiversion::Union{String,Nothing}=nothing)
    isempty(ctx.apis) && set_api_versions!(ctx)

    apictx = _get_apictx(ctx, O, apiversion)

    try
        apicall = eval(Symbol("patch$O"))
        return apicall(apictx, name, patch; _mediaType=patch_type)
    catch ex
        isa(ex, UndefVarError) || rethrow()
        apicall = eval(Symbol("patchNamespaced$O"))
        return apicall(apictx, name, ctx.namespace, patch; _mediaType=patch_type)
    end
end

"""
Fetch logs of the specified Pod

Parameters:
- name::String (required)

Keyword Args:
- container::String : The container for which to stream logs. Defaults to only container if there is one container in the pod.
- follow::Bool : Follow the log stream of the pod. Defaults to false.
- limitBytes::Int32 : If set, the number of bytes to read from the server before terminating the log output. This may not display a complete final line of logging, and may return slightly more or slightly less than the specified limit.
- pretty::String : If 'true', then the output is pretty printed.
- previous::Bool : Return previous terminated container logs. Defaults to false.
- sinceSeconds::Int32 : A relative time in seconds before the current time from which to show logs. If this value precedes the time a pod was started, only logs since the pod start will be returned. If this value is in the future, no logs will be returned. Only one of sinceSeconds or sinceTime may be specified.
- sinceTime::String : An RFC3339 timestamp from which to show logs. If this value precedes the time a pod was started, only logs since the pod start will be returned. If this value is in the future, no logs will be returned. Only one of sinceSeconds or sinceTime may be specified.
- tailLines::Int32 : If set, the number of lines from the end of the logs to show. If not specified, logs are shown from the creation of the container or sinceSeconds or sinceTime
- timestamps::Bool : If true, add an RFC3339 or RFC3339Nano timestamp at the beginning of every line of log output. Defaults to false.

Returns: String of all log entries, one per line
"""
get_logs(ctx::KuberContext, pod_name::String; kwargs...) = get(ctx, :PodLog, pod_name; kwargs...)

list_namespaced_custom_metrics(ctx::KuberContext, metricname::String; kwargs...) = list(ctx, :MetricValue, "metrics/"*metricname; kwargs...)
list_namespaced_custom_metrics(ctx::KuberContext, objecttype::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/*/" * metricname; kwargs...)
list_namespaced_custom_metrics(ctx::KuberContext, objecttype::String, objectname::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/" * objectname * "/" * metricname; kwargs...)

list_custom_metrics(ctx::KuberContext, objecttype::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/*/" * metricname; namespace=nothing, kwargs...)
list_custom_metrics(ctx::KuberContext, objecttype::String, objectname::String, metricname::String; kwargs...) = list(ctx, :MetricValue, objecttype * "/" * objectname * "/" * metricname; namespace=nothing, kwargs...)


text:
module Tinker

using Gtk.ShortNames, GtkReactive, Graphics, Colors, Images, IntervalSets
using Luxor.isinside, Luxor.Point

img_ctxs = Signal([])

abstract type Shape end

# Rectangle structure
struct Rectangle <: Shape
    x::Number
    y::Number
    w::Number
    h::Number
end

Rectangle() = Rectangle(0,0,-1,-1)
Base.isempty(R::Rectangle) = R.w <= 0 || R.h <= 0

mutable struct ImageContext{T}
    id::Int
    image
    canvas::GtkReactive.Canvas
    zr::Signal{ZoomRegion{T}}
    zl::Int # for tracking zoom level
    mouseactions::Dict{String,Any}
    shape::Signal{<:Shape} # Tracks type of selection in the environment
    points::Signal{<:AbstractArray} # Holds points that define shape outline
    rectview::Signal{<:AbstractArray} # Holds bounding box of selection
end

ImageContext() = ImageContext(-1, nothing, canvas(),
Signal(ZoomRegion((1:10, 1:10))), -1, Dict("dummy"=>Signal(false),"dummy2"=>Signal(4)),
Signal(Rectangle()), Signal([]), Signal([]))

# Returns a view of an image with the given bounding region
function get_view(image,x_min,y_min,x_max,y_max)
    xleft,yleft = Int(floor(Float64(x_min))),Int(ceil(Float64(y_min)))
    xright,yright = Int(floor(Float64(x_max))),Int(ceil(Float64(y_max)))
    min_x, max_x = map(x->x, extrema(indices(image,2)))
    min_y, max_y = map(y->y, extrema(indices(image,1)))
    (xleft < min_x) && (xleft = min_x)
    (yleft < min_y) && (yleft = min_y)
    (xright > max_x) && (xright = max_x)
    (yright > max_y) && (yright = max_y)
    return view(image, yleft:yright, xleft:xright)
end

# Calculates tolerance based on zoom level
function get_tolerance(zr::ZoomRegion)
  tol = 5*(IntervalSets.width(zr.currentview.x)/IntervalSets.width(zr.fullview.x))
  return tol
end

function get_tolerance(ctx::ImageContext)
  zr = value(ctx.zr)
  get_tolerance(zr)
end

# Creates a Rectangle out of any two points
function Rectangle(p1::XY,p2::XY)
    x, w = min(p1.x, p2.x), abs(p2.x - p1.x)
    y, h = min(p1.y, p2.y), abs(p2.y - p1.y)
    return Rectangle(x, y, w, h)
    (p1.x == p2.x) || (p1.y == p2.y) && return Rectangle()
end

# rectangle draw function
function drawrect(ctx, rect, color, width)
    set_source(ctx, color)
    set_line_width(ctx, width)
    rectangle(ctx, rect.x, rect.y, rect.w, rect.h)
    stroke(ctx)
end;

# Polygon shape
struct Polygon <: Shape
    pts::AbstractArray
end

Polygon() = Polygon([])
Base.isempty(P::Polygon) = isempty(P.pts)

# Handles modify rectangles
struct Handle <: Shape
    r::Shape
    pos::String # refers to which side or corner of rectangle handle is on
    x::Float64
    y::Float64
end

Handle() = Handle(Rectangle(),"",0,0)
Base.isempty(H::Handle) = isempty(H.r)

# Creates handle given a Rectangle and a position
function Handle(r::Rectangle, pos::String)
    # Position of handle refers to center coordinate of handle based on rect
    position_coord = Dict("tlc"=>(r.x,r.y),"ts"=>(r.x+(r.w/2),r.y),
                          "trc"=>(r.x+r.w,r.y),"rs"=>(r.x+r.w,r.y+(r.h/2)),
                          "brc"=>(r.x+r.w,r.y+r.h),"bs"=>(r.x+(r.w/2),r.y+r.h),
                          "blc"=>(r.x,r.y+r.h),"ls"=>(r.x,r.y+(r.h/2)))
    xy = get(position_coord, pos, (-Inf,-Inf))
    if xy == (-Inf,-Inf)
        println("Not a valid Handle position.")
        return Handle()
    else
        x = position_coord[pos][1]
        y = position_coord[pos][2]
        return Handle(r,pos,xy[1],xy[2])
    end
end

function Handle(p::Polygon, i::Int)
    pt = p.pts[i]
    return Handle(p, "", pt.x,pt.y)
end

# Draws a handle
function drawhandle(ctx, handle::Handle, d)
  if !isempty(handle)
    rectangle(ctx, handle.x-(d), handle.y-(d),
              d*2, d*2)
    set_source(ctx,colorant"white")
    fill_preserve(ctx)
    set_source(ctx,colorant"black")
    set_line_width(ctx,1.0)
    stroke_preserve(ctx)
  end
end; # like drawrect, but makes x,y refer to center of handle

# A rectangle with handles at all 8 positions
struct RectHandle <: Shape
    r::Rectangle
    h::NTuple{8,Handle}
end

RectHandle() = RectHandle(Rectangle())
Base.isempty(RH::RectHandle) = isempty(RH.r)

# Creates a RectHandle given just a Rectangle
function RectHandle(r::Rectangle)
    # derive all 8 handles from r
    # numbered 1-8, 1=tlc, moving clockwise around rectangle
    h = (Handle(r, "tlc"), Handle(r, "ts"), Handle(r, "trc"), Handle(r, "rs"),
         Handle(r, "brc"), Handle(r, "bs"), Handle(r, "blc"), Handle(r, "ls"))
    return RectHandle(r,h)
end

# A polygon with handles at every vertex
struct PolyHandle <: Shape
    p::Polygon
    h::AbstractArray
end

PolyHandle() = PolyHandle(Polygon([]),[])

function PolyHandle(p::Polygon)
    # makes polyhandle
    h = [Handle(p,1)]
    for i in 2:length(p.pts)
        push!(h,Handle(p,i))
    end
    return PolyHandle(p,h)
end

function PolyHandle(p::AbstractArray)
  return PolyHandle(Polygon(p))
end

# Draws RectHandle
function drawrecthandle(ctx, rh::RectHandle, d, color1, width)
    drawrect(ctx, rh.r, color1, width)
    for n in 1:length(rh.h)
        drawhandle(ctx, rh.h[n], d)
    end
end

# draws PolyHandle
function drawpolyhandle(ctx, ph::PolyHandle, d, color, width)
  drawline(ctx, ph.p.pts, color, width)
  for i in 1:length(ph.h)
    drawhandle(ctx, ph.h[i], d)
  end
end

# Connects an array of points
function drawline(ctx, l, color, width)
    isempty(l) && return
    p = first(l)
    move_to(ctx, p.x, p.y)
    set_source(ctx, color)
    set_line_width(ctx, width)
    for i = 2:length(l)
        p = l[i]
        line_to(ctx, p.x, p.y)
    end
    stroke(ctx)
end

# Set of versatile methods to draw many shapes with the same function
function drawshape(ctx, sh::RectHandle, d, color, width)
  # draw RectHandle
  drawrecthandle(ctx, sh, d, color, width)
end

function drawshape(ctx, sh::Polygon, d, color, width)
  # draw Polygon
  drawline(ctx, sh.pts, color, width)
end

function drawshape(ctx, sh::PolyHandle, d, color, width)
  # draw PolyHandle
  drawpolyhandle(ctx, sh, d, color, width)
end

# Checks if an array of points qualifies as a polygon
function ispolygon(p::AbstractVector)
    length(p) < 4 && return false
    p[1] != p[end] && return false
    for i in 1:(length(p)-1)
        p[i]!=p[i+1] && return true
    end
    return false
end

# Moves polygon to a given location
function move_polygon_to(p::AbstractArray, pt::XY)
    # find diff b/t start & pt; add diff to all in p
    if ispolygon(p)
        diff = XY(pt.x-p[1].x,pt.y-p[1].y)
        map(n -> XY(n.x+diff.x,n.y+diff.y),p)
    else
        #println("Not a polygon")
        return p
    end
end

# Converts an XY to a Point
function Point(p::XY)
    Point(p.x,p.y)
end

include("zoom_interaction.jl")
include("selection_actions.jl")

## Sets up an image in a separate window with the ability to adjust view
function init_image(image::AbstractArray; name="Tinker")
    # set up window
    win = Window(name, size(image,2), size(image,1))
    c = canvas(UserUnit)
    push!(win, c)

    # set up a zoom region
    zr = Signal(ZoomRegion(image))

    # create view
    imagesig = map(zr) do r
        cv = r.currentview
        view(image, UnitRange{Int}(cv.y), UnitRange{Int}(cv.x))
    end;

    # create a view diagram
    viewdim = map(zr) do r
        fvx, fvy = r.fullview.x, r.fullview.y # x, y range of full view
        cvx, cvy = r.currentview.x, r.currentview.y # x, y range of currentview
        xfull, yfull =
            (fvx.right-fvx.left),(fvy.right-fvy.left) # width of full view
        xcurrent, ycurrent =
            (cvx.right-cvx.left),(cvy.right-cvy.left) # width of current view
        # scale
        xsc,ysc = 0.1*(xcurrent/xfull), 0.1*(ycurrent/yfull)
        # offset
        x_off,y_off = cvx.left+(0.01*xcurrent),cvy.left+(0.01*ycurrent)
        # represents full view
        rect1 = Rectangle(x_off, y_off, xsc*xfull, ysc*yfull)
        # represents current view
        rect2 = Rectangle(x_off+(cvx.left*xsc), y_off+(cvy.left*ysc),
                          xsc*xcurrent, ysc*ycurrent)
        return [rect1,rect2]
    end

    # Placeholder dictionary for context
    dummydict = Dict("pandrag"=>Signal(false),"zoomclick"=>Signal(false),
    "select"=>Signal(false),"selmode"=>Signal(rectangle_mode))

    # Context
    imagectx=ImageContext(length(value(img_ctxs))+1,image,c,zr,1,dummydict,
    Signal(Shape, Rectangle()),Signal([]),Signal(view(image,1:size(image,2),
    1:size(image,1))))

    # Mouse actions
    pandrag = init_pan_drag(c, zr)
    zoomclick = init_zoom_click(imagectx)
    select = init_selection_actions(imagectx)

    push!(pandrag["enabled"],false)
    push!(zoomclick["enabled"],false)

    imagectx.mouseactions = Dict("pandrag"=>pandrag["enabled"],"zoomclick"=>
    zoomclick["enabled"],"select"=>select["enabled"],"selmode"=>select["mode"])

    imagectx.rectview = map(imagectx.points) do pts
        if !isempty(pts)
            x_min,x_max = minimum(map(n->n.x,pts)),maximum(map(n->n.x,pts))
            y_min,y_max =minimum(map(n->n.y,pts)),maximum(map(n->n.y,pts))
            get_view(image,x_min,y_min,x_max,y_max)
        else
            view(image,1:size(image,2),1:size(image,1))
        end
    end

    append!(c.preserved, [pandrag,zoomclick,select])

    # draw
    redraw = draw(c, imagesig, zr, viewdim, imagectx.points, imagectx.shape) do cnvs,img,r,vd,pts,sh
        copy!(cnvs, img) # show image on canvas at current zoom level
        set_coordinates(cnvs, r) # set canvas coordinates to zr
        ctx = getgc(cnvs)
        # draw view diagram if zoomed in
        if r.fullview != r.currentview
            drawrect(ctx, vd[1], colorant"blue", 2.0)
            drawrect(ctx, vd[2], colorant"blue", 2.0)
        end
        d = get_tolerance(imagectx)
        if ispolygon(pts)
          # draw shape
          if typeof(sh) == RectHandle
              recthandle = RectHandle(Rectangle(XY(pts[1].x,pts[1].y),XY(pts[3].x,pts[3].y)))
              drawrecthandle(ctx, recthandle, d, colorant"yellow", 1.0)
          elseif typeof(sh) == PolyHandle
              drawpolyhandle(ctx, PolyHandle(pts), d, colorant"yellow", 1.0)
          else
              drawline(ctx, pts, colorant"yellow", 1.0)
          end
        else
          # draw working line
          drawline(ctx, pts, colorant"yellow", 1.0)
          if typeof(sh) == PolyHandle && length(pts) > 0
              drawrect(ctx, Rectangle(pts[1].x-d,pts[1].y-d,2*d,2*d), colorant"yellow", 1.0)
          end
        end
    end

    showall(win);

    # shows measurements
    display_measure(imagectx)

    push!(img_ctxs, push!(value(img_ctxs), imagectx))
    return imagectx
end;

init_image(file::AbstractString) = init_image(load(file); name=file)

@enum Mode zoom_mode rectangle_mode freehand_mode polygon_mode

function set_mode(ctx::ImageContext, mode::Mode)
    push!(ctx.mouseactions["pandrag"], false)
    push!(ctx.mouseactions["zoomclick"], false)
    push!(ctx.mouseactions["select"],false)
    if mode == zoom_mode # turn on zoom controls
        println("Zoom mode")
        push!(ctx.mouseactions["pandrag"], true)
        push!(ctx.mouseactions["zoomclick"], true)
        push!(ctx.mouseactions["select"],false)
    elseif mode == rectangle_mode # turn on rectangular region selection controls
        println("Rectangle mode")
        push!(ctx.mouseactions["select"],true)
        push!(ctx.mouseactions["selmode"],rectangle_mode)
    elseif mode == freehand_mode # freehand select
        println("Freehand mode")
        push!(ctx.mouseactions["select"],true)
        push!(ctx.mouseactions["selmode"],freehand_mode)
    elseif mode == polygon_mode # polygon select
        println("Polygon mode")
        push!(ctx.mouseactions["select"],true)
        push!(ctx.mouseactions["selmode"],polygon_mode)
    end
end

set_mode(sig::Signal, mode) = set_mode(value(sig), mode)

# Sets all contexts
function set_mode_all(mode::Mode)
  for c in value(img_ctxs)
    set_mode(c, mode)
  end
  nothing
end

include("measure.jl")
include("guisetup.jl")

end # module


text:
push!(LOAD_PATH,"../src/")
using Documenter, SeisNoise

makedocs(
    modules = [SeisNoise],
    sitename = "SeisNoise.jl",
    # Uncomment below for local build
    # format = Documenter.HTML(prettyurls = false),
    authors = "Tim Clements",
    pages = Any[
        "Home" => "index.md",
        "Workflow" => "types.md",
        "Using the GPU" => "using_gpus.md",
        "Pre-Processing" => "preprocessing.md",
        "Computing FFTs" => "fft.md",
        "Correlation" => "correlation.md",
        "Extending SeisNoise" => "extend.md",
        "Examples" => "examples.md",
        "Contributer's Guide" => "contributing.md",
        "Function Index" => "func_index.md",
        ],
)

deploydocs(
    repo = "github.com/tclements/SeisNoise.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
)