text:
# run tests in parallel with 2 processes

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))

using PDESolver
#using Base.Test
using Base.Test
using ODLCommonTools
using PdePumiInterface  # common mesh interface - pumi
using SummationByParts  # SBP operators
using Utils
using EulerEquationMod
using LinearSolvers
using NonlinearSolvers   # non-linear solvers
using OptimizationInterface
using ArrayViews
import ArrayViews.view
import MPI
using Input

#------------------------------------------------------------------------------
# define test list
using TestSystem
include("../tags.jl")

global const EulerTests = TestList()
# define global const tags here

# NOTE: For adding parallel tests to TICON, Write them in a separate file and
# include them below `global const EulerTests = TestList()`.

include("test_parallel_derivatives.jl")

"""
  Run the parallel tests and compare against serial results run as part of
  the serial tests
"""
function test_parallel2()
  @testset "----- Testing Parallel -----" begin

    start_dir = pwd()

    # test rk4
    cd("./rk4/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")

    @test isapprox( datas[1], datap[1]) atol=1e-13
    @test isapprox( datas[2], datap[2]) atol=1e-13

    # test staggered_parallel
    cd("../staggered_parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../staggered_serial/error_calc.dat")
    datap = readdlm("error_calc.dat")

    @test isapprox( datas[1], datap[1]) atol=1e-13
    @test isapprox( datas[2], datap[2]) atol=1e-13

    cd("../../")

    cd("./lserk/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")

    @test isapprox( datas[1], datap[1]) atol=1e-13
    @test isapprox( datas[2], datap[2]) atol=1e-13

    cd("../../")

    # test newton
    cd("./newton/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("./error_calc.dat")
    @test isapprox( datas[1], datap[1]) atol=1e-13

    fill!(eqn.q, 1.0)
    fill!(eqn.res, 0.0)

    for i=1:mesh.npeers
      fill!(eqn.shared_data[i].q_send, 1.0)
      fill!(eqn.shared_data[i].q_recv, 1.0)
    end
    evalHomotopy(mesh, sbp, eqn, opts, eqn.res)

    @test isapprox( vecnorm(eqn.res), 0.0) atol=1e-13

    cd(start_dir)

  end

  return nothing
end

function test_parallel_nopre()

  @testset "----- Testing parallel nopre -----" begin
    start_dir = pwd()

    # test rk4
    cd("./rk4/parallel")
    fname = "input_vals_parallel.jl"

    mesh, sbp, eqn, opts = solvePDE(fname)

    EulerEquationMod.evalResidual(mesh, sbp, eqn, opts)
    res_orig = copy(eqn.res)

    opts["precompute_face_flux"] = false
    fill!(eqn.res, 0.0)

    EulerEquationMod.evalResidual(mesh, sbp, eqn, opts)

    @test isapprox( norm(vec(eqn.res - res_orig)), 0.0) atol=1e-13

    cd(start_dir)

    # test_newton
    cd("./newton/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    fill!(eqn.res, 0.0)
    EulerEquationMod.evalResidual(mesh, sbp, eqn, opts)
    res_orig = copy(eqn.res)

    opts["precompute_face_flux"] = false
    fill!(eqn.res, 0.0)

    EulerEquationMod.evalResidual(mesh, sbp, eqn, opts)

    @test isapprox( norm(vec(eqn.res - res_orig)), 0.0) atol=1e-13

    cd(start_dir)

  end  # end facts block

  return nothing
end

function test_restart()
  @testset "----- Testing restart -----" begin
    start_dir = pwd()

    # test rk4
    cd("./rk4/parallel")
    fname = "input_vals_restart"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")

    @test isapprox( datas[1], datap[1]) atol=1e-13
    @test isapprox( datas[2], datap[2]) atol=1e-13

    cd("../../lserk/parallel")
    fname = "input_vals_restart"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")

    @test isapprox( datas[1], datap[1]) atol=1e-13
    @test isapprox( datas[2], datap[2]) atol=1e-13

    cd(start_dir)
  end

  return nothing
end

#test_parallel2()
add_func1!(EulerTests, test_parallel2, [TAG_SHORTTEST])
add_func1!(EulerTests, test_parallel_nopre, [TAG_SHORTTEST])
add_func1!(EulerTests, test_restart, [TAG_SHORTTEST])

#------------------------------------------------------------------------------
# run tests
@testset "----- Running Euler 2 process tests -----" begin
  runTestSystem(EulerTests, solvePDE, ARGS)
end


text:
# ---
# title: 1658. Minimum Operations to Reduce X to Zero
# id: problem1658
# author: AquaIndigo
# date: 2020-11-23
# difficulty: Medium
# categories: Two Pointers, Binary Search, Greedy
# link: <https://leetcode.com/problems/minimum-operations-to-reduce-x-to-zero/description/>
# hidden: true
# ---
#
# You are given an integer array `nums` and an integer `x`. In one operation,
# you can either remove the leftmost or the rightmost element from the array
# `nums` and subtract its value from `x`. Note that this **modifies** the array
# for future operations.
#
# Return _the **minimum number** of operations to reduce _`x` _to **exactly**_
# `0` _if it 's possible_ _, otherwise, return_`-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,4,2,3], x = 5
#     Output: 2
#     Explanation: The optimal solution is to remove the last two elements to reduce x to zero.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [5,6,7,8,9], x = 4
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: nums = [3,2,20,1,1,3], x = 10
#     Output: 5
#     Explanation: The optimal solution is to remove the last three elements and the first two elements (5 operations in total) to reduce x to zero.
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 105`
#   * `1 <= nums[i] <= 104`
#   * `1 <= x <= 109`
#
#
## @lc code=start
using LeetCode

function min_operations(nums::Vector{Int}, x::Int)
    target = sum(nums) - x
    target < 0 && return -1
    println(target)
    left, right, maxn = 0, 0, -1
    len = length(nums)
    while right < len
        target -= nums[right += 1]
        while target < 0
            target += nums[left += 1]
        end
        (target == 0) && (maxn = max(maxn, right - left))
    end
    maxn == -1 ? maxn : len - maxn
end
## @lc code=end

ext: jl


text:
struct InnerConstructor end
macro wrap_type(julia_name, c_name, destructor)
    lower_constructor_name = Symbol("_$julia_name")

    # Automatically add the "Lib" prefix if required.
    if isa(c_name, Symbol)
        c_name = :(Lib.$c_name)
    end

    return esc(
        quote
            # Type definition
            mutable struct $julia_name
                handle::$c_name
                $julia_name(::InnerConstructor) = new($c_name())
            end

            # Use a trick of Lower and Higher constructeors.
            # Lower constructors should have the name `_$julia_name` and not attach
            # finalizers.
            #
            # The higher constructor will simply forward to the lower constructor but
            # attach a finalizer before returning.
            function $julia_name(args...)
                x = $lower_constructor_name(args...)
                attach_finalizer!(x)
                return x
            end

            # Finalizer
            destroy(x::$julia_name) = @apicall $destructor(x)
            attach_finalizer!(x::$julia_name) = finalizer(destroy, x)

            # Conversion functions
            Base.unsafe_convert(::Type{$c_name}, x::$julia_name) = x.handle
            function Base.unsafe_convert(::Type{Ptr{$c_name}}, x::$julia_name)
                return Base.unsafe_convert(Ptr{$c_name}, Base.pointer_from_objref(x))
            end
        end,
    )
end

destroy(x, y...) = (destroy(x); destroy(y...))
destroy(x) = error("Define `destroy` for type $(typeof(x))")

#####
##### Engine
#####

@wrap_type Engine dnnl_engine_t dnnl_engine_destroy
# Lower Constructor
function _Engine(kind::Lib.dnnl_engine_kind_t, index = 0)
    engine = Engine(InnerConstructor())
    @apicall dnnl_engine_create(engine, kind, index)
    return engine
end

@wrap_type Stream dnnl_stream_t dnnl_stream_destroy
# Lower Constructor
function _Stream(engine::Engine)
    stream = Stream(InnerConstructor())
    @apicall dnnl_stream_create(stream, engine, Lib.dnnl_stream_default_flags)
    return stream
end

@wrap_type Attributes dnnl_primitive_attr_t dnnl_primitive_attr_destroy
# Lower Constructor
function _Attributes()
    attributes = Attributes(InnerConstructor())
    @apicall dnnl_primitive_attr_create(attributes)
    @apicall dnnl_primitive_attr_set_scratchpad_mode(
        attributes, Lib.dnnl_scratchpad_mode_user
    )
    return attributes
end

@wrap_type PostOps dnnl_post_ops_t dnnl_post_ops_destroy
function _PostOps()
    postops = PostOps(InnerConstructor())
    @apicall dnnl_post_ops_create(postops)
    return postops
end

@wrap_type PrimitiveDescriptor dnnl_primitive_desc_t dnnl_primitive_desc_destroy
# Lower Constructor
# Automatically apply "Lib.dnnl_primitive_desc_create" if the first argument is
# not already a function.
function _PrimitiveDescriptor(args::Vararg{Any,N}) where {N}
    return _PrimitiveDescriptor(Lib.dnnl_primitive_desc_create, args...)
end

function _PrimitiveDescriptor(f::F, args::Vararg{Any,N}) where {F<:Function,N}
    descriptor = PrimitiveDescriptor(InnerConstructor())
    @apicall f(descriptor, args...)
    return descriptor
end

function Base.copy(x::PrimitiveDescriptor)
    y = PrimitiveDescriptor(InnerConstructor())
    @apicall dnnl_primitive_desc_clone(y, x)
    attach_finalizer!(y)
    return y
end

noforward() = Ptr{Lib.dnnl_primitive_desc}()

@wrap_type Primitive dnnl_primitive_t dnnl_primitive_destroy
# Lower Constructor
function _Primitive(descriptor::PrimitiveDescriptor)
    primitive = Primitive(InnerConstructor())
    @apicall dnnl_primitive_create(primitive, descriptor)
    return primitive
end

@wrap_type MemoryPtr dnnl_memory_t dnnl_memory_destroy
# Lower Constructor
function _MemoryPtr(A::AbstractArray, desc = memorydesc(A))
    return _MemoryPtr(convert(Ptr{Nothing}, pointer(A)), desc)
end

function _MemoryPtr(ptr::Ptr{Nothing}, desc)
    memory = MemoryPtr(InnerConstructor())
    @apicall dnnl_memory_create(memory, desc, global_engine(), ptr)
    return memory
end

# Define an extra conversion function.
@inline function Base.convert(
    ::Type{T}, memory::MemoryPtr
) where {T<:MaybePtr{Lib.dnnl_memory_t}}
    return Base.unsafe_convert(T, memory)
end

#####
##### Functions on Types
#####

function eltwise!(postops::PostOps, f::F, scale = 1) where {F}
    @apicall dnnl_post_ops_append_eltwise(postops, scale, forward_expand(f)...)
    return nothing
end

# Specialize identity to do nothing
eltwise!(postops::PostOps, ::typeof(identity), scale = 1) = nothing

# Attach Post Ops to Attributes
function Base.append!(a::Attributes, p::PostOps)
    @apicall Lib.dnnl_primitive_attr_set_post_ops(a, p)
end

appendsum!(p::PostOps, scale = 1) = @apicall dnnl_post_ops_append_sum(p, scale)
function query_md(
    descriptor::Union{PrimitiveDescriptor,Lib.dnnl_primitive_desc_t}, kind, index = 0
)
    # TODO: Do we need to perform a null check here?
    ptr = Lib.dnnl_primitive_desc_query_md(descriptor, kind, index)
    return unsafe_load(ptr)
end

macro query(sym)
    fullsym = Symbol("dnnl_query_$(sym)_md")
    return :(Lib.$fullsym)
end

#####
##### We manage our own global scratchpad here for now.
#####

const SCRATCHPAD = UInt8[]
const SCRATCHPAD_MEMORY = Ref{Lib.dnnl_memory_t}()

function unsafe_primitive_descriptor(primitive::Primitive)
    _pd = Ref{Lib.dnnl_primitive_desc_t}()
    @apicall dnnl_primitive_get_primitive_desc(primitive, _pd)
    return _pd[]
end

function execute!(primitive::Primitive, args::Tuple{<:AbstractArray,<:Arguments}; kw...)
    return execute!(primitive, args[2]; scratchpad_source = args[1], kw...)
end

function execute!(
    primitive::Primitive,
    _args;
    scratchpad_source = nothing,
    descriptor = nothing,
    wait = true,
)
    # Get the primitive descriptor if it isn't passed in alread.
    if descriptor === nothing
        _descriptor = unsafe_primitive_descriptor(primitive)
    else
        _descriptor = descriptor
    end

    # Create a scratchpad
    md = query_md(_descriptor, Lib.dnnl_query_scratchpad_md)
    bytes = getbytes(md)

    if scratchpad_source === nothing && bytes > length(SCRATCHPAD)
        resize!(SCRATCHPAD, bytes)
    end

    if scratchpad_source === nothing
        scratchpad = SCRATCHPAD
    else
        scratchpad = similar(scratchpad_source, UInt8, bytes)
    end

    GC.@preserve scratchpad begin
        @apicall dnnl_memory_create(
            SCRATCHPAD_MEMORY, md, global_engine(), pointer(scratchpad)
        )
        args = append(
            _args, Lib.dnnl_exec_arg_t(Lib.DNNL_ARG_SCRATCHPAD, SCRATCHPAD_MEMORY[])
        )

        # Finally, call the primitive
        @apicall dnnl_primitive_execute(primitive, global_stream(), length(args), args)
        @apicall dnnl_memory_destroy(SCRATCHPAD_MEMORY[])
        wait && @apicall(dnnl_stream_wait(global_stream()))
    end
    return nothing
end

function _execute!(primitive, args, scratchpad_md)
    _memory = _MemoryPtr(SCRATCHPAD, md)
    @time args = append(_args, Lib.dnnl_exec_arg_t(Lib.DNNL_ARG_SCRATCHPAD, _memory))

    # Finally, call the primitive
    @apicall dnnl_primitive_execute(primitive, global_stream(), length(args), args)
    return destroy(_memory)
end

# Automatically apply recursively to tuples.
kernel_exit_hook(x) = x
kernel_exit_hook(x::Union{<:Tuple,<:NamedTuple}) = map(kernel_exit_hook, x)

function temp_primitive(f::F, args::Vararg{Any,N}) where {F,N}
    desc = _PrimitiveDescriptor(args...)
    primitive = _Primitive(desc)
    ret = f(primitive, desc)

    # If `desc` or `primitive` show up in the returned items from `f`, then attach a
    # finalizer and don't eagerly destroy.
    #
    # If they don't show up, then assume it is safe to eagerly destroy them.
    wrapped = tuplewrap(ret)
    cleanup(desc, wrapped)
    cleanup(primitive, wrapped)
    return kernel_exit_hook(ret)
end

cleanup(x, ys::Tuple) = recursive_any(i -> i === x, ys) ? attach_finalizer!(x) : destroy(x)

tuplewrap(x) = (x,)
tuplewrap(x::Tuple) = x
tuplewrap(x::NamedTuple) = Tuple(x)

# Slightly better type inference then the generic `any`.
recursive_any(f::F, x::Tuple) where {F} = f(x[1]) ? true : recursive_any(f, Base.tail(x))
recursive_any(f::F, x::Tuple{}) where {F} = false


text:
# test extending type and that various generic methods are defined
# to fulfill the API of AbstractYieldCurve
@testset "generic methods and type extensions" begin
    struct MyYield <: Yields.AbstractYieldCurve
        rate
    end

    Yields.discount(c::MyYield,to) = exp(-c.rate * to)
    # Base.zero(c::MyYield,to) = Continuous(c.rate)

    my = MyYield(0.05)
    @test zero(my,1) ≈ Continuous(0.05)
    @test Yields.forward(my,1) ≈ Continuous(0.05)
    @test Yields.forward(my,1,2) ≈ Continuous(0.05)
    @test discount(my,1,2) ≈ exp(-0.05*1)
    @test accumulation(my,1,2) ≈ exp(0.05*1)
    @test accumulation(my,1) ≈ exp(0.05*1)
    @test Yields.__ratetype(my) == Yields.Rate{Float64,Continuous}
    @test Yields.CompoundingFrequency(my) == Continuous()

    @test Yields.par(my,1) |> Yields.rate > 0
end


text:
#=
jinc.jl
=#

export jinc

using SpecialFunctions: besselj1 # grab bessel function

"""
    jinc(x)

Return `jinc(x) = J1(pi*x)/(2x)`,
where `J1` is a Bessel function of the first kind.

Units of `x` are typically cycles/m.

Return type is `promote_type(typeof(x), Float32)`.
"""
function jinc(x::Real)
    T = promote_type(typeof(x), Float32)
    if (x == 0) return convert(T, pi/4) end
    y = abs(x)
    convert(T, besselj1(pi*y) / (2*y))
end


text:
module Test_analysis

using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

println("Test_analysis")
@testset "zeros, rank, poles" begin

# test example for SLICOT subroutine AB08ND
A = [
   1.0   0.0   0.0   0.0   0.0   0.0
   0.0   1.0   0.0   0.0   0.0   0.0
   0.0   0.0   3.0   0.0   0.0   0.0
   0.0   0.0   0.0  -4.0   0.0   0.0
   0.0   0.0   0.0   0.0  -1.0   0.0
   0.0   0.0   0.0   0.0   0.0   3.0];
E = I;
B = [
   0.0  -1.0
  -1.0   0.0
   1.0  -1.0
   0.0   0.0
   0.0   1.0
  -1.0  -1.0];
C = [
   1.0   0.0   0.0   1.0   0.0   0.0
   0.0   1.0   0.0   1.0   0.0   1.0
   0.0   0.0   1.0   0.0   0.0   1.0];
D = [
   0.0   0.0
   0.0   0.0
   0.0   0.0];
sys = dss(A,E,B,C,D);

# zeros
@time val  = gzero(sys)
@test sort(real(val)) ≈ [-1, 2, Inf,Inf]

@time val  = gzero(sys[:,1])
@test sort(real(val)) ≈ [-1, 2, Inf]

@time val, info  = gzeroinfo(sys)
@test sort(real(val)) ≈ [-1, 2, Inf,Inf]  &&
      info.nfz == 2 && info.niev == 4 && info.nisev == 2 && info.niz == 2 &&
      info.nfsz == 1 && info.nfsbz == 0 && info.nfuz == 1 &&
      info.nrank == 8 && info.miev == [2, 2] && info.miz == [1, 1] &&
      info.rki == Int64[] && info.lki == [2] && !info.regular && !info.stable

@time val, info  = gzeroinfo(dss(A,E,B,C,D,Ts=-1))
@test sort(real(val)) ≈ [-1, 2, Inf,Inf]  &&
      info.nfz == 2 && info.niev == 4 && info.nisev == 2 && info.niz == 2 &&
      info.nfsz == 0 && info.nfsbz == 1 && info.nfuz == 1 &&
      info.nrank == 8 && info.miev == [2, 2] && info.miz == [1, 1] &&
      info.rki == Int64[] && info.lki == [2] && !info.regular && !info.stable

# poles
@time val  = gpole(sys)
@test sort(real(val)) ≈ sort([1.0, 1.0, 3.0, -4.0, -1.0, 3.0])

@time val  = gzero(dss(A = sys.A, E = sys.E))
@test sort(real(val)) ≈ sort([1.0, 1.0, 3.0, -4.0, -1.0, 3.0])

@test !isstable(sys)
@test isstable(dss(A = sys.A-5I))

@time val, info  = gpoleinfo(sys)
@test sort(real(val)) ≈ sort([1.0, 1.0, 3.0, -4.0, -1.0, 3.0]) &&
      info.nfev == 6 && info.nfsev == 2 && info.nfuev == 4

# rank
@test gnrank(sys, fastrank = true) == 2 && gnrank(sys, fastrank = false) == 2

@test gnrank(sys[:,1], fastrank = true) == 1 && gnrank(sys[:,1], fastrank = false) == 1

# output decoupling zeros
@time val  = gzero(dss(A = sys.A,E = sys.E,C = sys.C))
@test val ≈ [-1]

# input decoupling zeros
@time val  = gzero(dss(A = sys.A,E = sys.E,B = sys.B))
@test val ≈ [-4]

# test example for SLICOT subroutine AB08ND
A = [
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E = [
    0     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B = [
    -1     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C = [
    0     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];
D = [
    1     2    -2
    0    -1    -2
    0     0     0];

sys = dss(A,E,B,C,D);

# zeros
@time val  = gzero(sys)
@test val ≈ [1, Inf,Inf]

@time val, info  = gzeroinfo(sys)
@test val ≈ [1, Inf,Inf]  &&
      info.nfz == 1 && info.niev == 7 && info.nisev == 1 && info.niz == 2 &&
      info.nfsz == 0 && info.nfsbz == 0 && info.nfuz == 1 &&
      info.nrank == 11 && info.miev == [1, 1, 1, 1, 3] && info.miz == [2] &&
      info.rki == [2] && info.lki == [1] && !info.regular && !info.stable

@time val, info  = gzeroinfo(dss(A,E,B,C,D,Ts=-1))
@test val ≈ [1, Inf,Inf]  &&
      info.nfz == 1 && info.niev == 7 && info.nisev == 1 && info.niz == 2 &&
      info.nfsz == 0 && info.nfsbz == 1 && info.nfuz == 0 &&
      info.nrank == 11 && info.miev == [1, 1, 1, 1, 3] && info.miz == [2] &&
      info.rki == [2] && info.lki == [1] && !info.regular && !info.stable

# poles
@time val  = gpole(sys,check_reg = true)
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf]

@time val  = gzero(dss(A = sys.A, E = sys.E))
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf]

@test !isstable(sys)
@test isstable(dss(A = sys.E, E = sys.A, Ts=-1))

@time val, info  = gpoleinfo(sys)
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf] &&
      info.nfev == 0 && info.nfsev == 0 && info.nfuev == 0 &&
      info.niev == 9 && info.nip == 6 && info.miev == [3, 3, 3] && info.mip == [2, 2, 2] &&
      info.regular && !info.proper && !info.stable

@test gnrank(sys, atol = 1.e-7, fastrank = true) == 2 && gnrank(sys, atol = 1.e-7, fastrank = false) ==2

@test gnrank(sys[:,1], fastrank = true) == 1 && gnrank(sys[:,1], fastrank = false) == 1

# output decoupling zeros
@time val  = gzero(dss(A = sys.A,E = sys.E,C = sys.C))
@test val ≈ [Inf, Inf, Inf, Inf]

# input decoupling zeros
@time val  = gzero(dss(A = sys.A,E = sys.E,B = sys.B))
@test val ≈ Float64[]

# stability tests
@test isstable(rdss(5,3,3, stable = true))
@test isstable(rdss(5,3,3, id = ones(Int,3), stable = true), atol = 1.e-7)

# infinite structure test
n = 5; m = 2; p = 3;
sys = rdss(n,p,m,id=[ones(Int,3);2*ones(Int,2)]);
@time val, info  = gpoleinfo(sys,atol=1.e-7)
@test count(isinf.(val)) == 2 &&
      info.nfev == n && info.niev == 7 && info.nisev == 3 && info.nfsev + info.nfuev == n &&
      info.nip == 2 && info.miev == [1, 1, 1, 2, 2] && info.mip == [1, 1] &&
      info.regular && !info.proper && isregular(sys)

A = rand(n,n); E = rand(n,n);
@time val, info  = gpoleinfo(dss(A=A,E=E))
@test count(isinf.(val)) == 0 &&
      info.nfev == n && info.niev == 0 && info.nisev == 0 && info.nfsev + info.nfuev == n &&
      info.nip == 0 && info.miev == [] && info.mip == [] &&
      info.regular && info.proper

As, Es, = schur(A,E)
@time val, info  = gpoleinfo(dss(A=As,E=Es))
@test count(isinf.(val)) == 0 &&
      info.nfev == n && info.niev == 0 && info.nisev == 0 && info.nfsev + info.nfuev == n &&
      info.nip == 0 && info.miev == [] && info.mip == [] &&
      info.regular && info.proper

@time val, info  = gpoleinfo(dss(A=A,E=triu(E)))
@test count(isinf.(val)) == 0 &&
      info.nfev == n && info.niev == 0 && info.nisev == 0 && info.nfsev + info.nfuev == n &&
      info.nip == 0 && info.miev == [] && info.mip == [] &&
      info.regular && info.proper

end

@testset "gl2norm & gh2norm" begin

sys = rdss(0,0,0);
@time l2norm = gl2norm(sys)
@test l2norm == 0

sys = rdss(0,2,2);
@time l2norm = gl2norm(sys)
@test l2norm == Inf

sys = rdss(0,2,2,disc=true);
@time l2norm = gl2norm(sys)
@test l2norm == norm(sys.D)

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time l2norm = gh2norm(sys)
@test l2norm ≈ 1.419727086450068e+01

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(2,2);
sys = dss(e*a,e,e*b,c,d);
@time l2norm = gh2norm(sys)
@test l2norm ≈ 1.419727086450068e+01

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time l2norm = gl2norm(sys)
@test l2norm ≈ 1.419727086450068e+01

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time l2norm = opnorm(sys,2)
@test l2norm ≈ 1.419727086450068e+01

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d);
@time l2norm = gl2norm(sys)
@test l2norm ≈ 1.419727086450068e+01

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time l2norm = gl2norm(sys')
@test l2norm ≈ 1.419727086450068e+01

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d);
@time l2norm = gl2norm(sys')
@test l2norm ≈ 1.419727086450068e+01

# discrete standard & descriptor
a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d,Ts = 1);
@time l2norm = gh2norm(sys)
@test l2norm ≈ 3.438689619923066e+01

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d,Ts = 1);
@time l2norm = gh2norm(sys)
@test l2norm ≈ 3.438689619923066e+01

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d,Ts = 1);
@time l2norm = gl2norm(sys)
@test l2norm ≈ 3.438689619923066e+01

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d,Ts = 1);
@time l2norm = gl2norm(sys)
@test l2norm ≈ 3.438689619923066e+01

a = [-1 2;0 0]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d);
@time l2norm = gl2norm(sys)
@test l2norm ≈ Inf

a = [-1 2;0 0]; e = [1 0; 0 2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,e,b,c,d);
@time l2norm = gl2norm(sys)
@test l2norm ≈ Inf

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d,Ts = 1);
@time l2norm = gl2norm(sys')
@test l2norm ≈ 3.438689619923066e+01

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d,Ts = 1);
@time l2norm = gl2norm(sys')
@test l2norm ≈ 3.438689619923066e+01

z = Polynomial([0, 1],'z')
g = [z^2+z+1 4*z^2+3*z+2 2*z^2-2;
    z 4*z-1 2*z-2;
    z^2 4*z^2-z 2*z^2-2*z];

sys = dss(g,minimal = true, atol = 1.e-7,Ts = 1)
@time l2norm = gl2norm(sys, atol = 1.e-7)
@test l2norm ≈ gl2norm(sys', atol = 1.e-7)

n = 5; m = 3; p = 2;

Ty = Float64; fast = true;
for Ty in (Float64, Complex{Float64})

for fast in (true, false)

sys = rss(n,p,m,T = Ty, stable = true);
sys.D[:,:] = zeros(Ty,p,m);
@test gh2norm(sys) ≈ gl2norm(sys')

sys = rss(n,p,m,T = Ty, stable = true);
@test gh2norm(sys) ≈ gl2norm(sys')

sys = rss(n,p,m,T = Ty, stable = true, disc = true);
@test gh2norm(sys) ≈ gl2norm(sys')

sys = rdss(n,p,m,T = Ty, stable = true);
sys.D[:,:] = zeros(Ty,p,m);
@test gh2norm(sys) ≈ gl2norm(sys')

sys = rdss(n,p,m,T = Ty, stable = true);
@test gh2norm(sys) ≈ gl2norm(sys')

sys = rdss(n,p,m,T = Ty, stable = true, disc = true);
@test gh2norm(sys) ≈ gl2norm(sys')

sys = rdss(n,p,m, T = Ty, stable = true, id=ones(Int,3));
@test gh2norm(sys) ≈ gl2norm(sys')

end # fast
end # Ty
end # gl2norm

@testset "ghanorm" begin

sys = rdss(0,0,0);
@time hanorm, hsv = ghanorm(sys)
@test hanorm == 0 && hsv == Float64[]

n = 5; m = 3; p = 2;

Ty = Float64; fast = true;
for Ty in (Float64, Complex{Float64})

for fast in (true, false)

sys = rss(n,p,m,stable = true);

try
  @time hanorm, hsv = ghanorm(sys,fast = fast)
  @test true
catch
  @test false
end

try
    @time hanorm, hsv = ghanorm(sys',fast = fast)
    @test false
catch
    @test true
end

sys = rss(n,p,m,T = Ty, stable = true);

@time hanorm, hsv = ghanorm(sys-sys,fast = fast)
@test hanorm < 1.e-7

@time hanorm, hsv = ghanorm(sys+sys,fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

@time hanorm, hsv = ghanorm([sys sys],fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

sys = rss(n,p,m,T = Ty, disc = true, stable = true);
try
    @time hanorm, hsv = ghanorm(sys,fast = fast)
    @test true
catch
    @test false
end

try
    @time hanorm, hsv = ghanorm(sys',fast = fast)
    @test false
catch
    @test true
end
@time hanorm, hsv = ghanorm(sys-sys,fast = fast)
@test hanorm < 1.e-7

@time hanorm, hsv = ghanorm(sys+sys,fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

@time hanorm, hsv = ghanorm([sys sys],fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

sys = rdss(n,p,m,T = Ty, stable = true,id=ones(Int,3));
try
  @time hanorm, hsv = ghanorm(sys,atol=1.e-7)
  @test true
catch
  @test false
end

try
    @time hanorm, hsv = ghanorm(sys',fast = fast)
    @test false
catch
    @test true
end

sys = rdss(n,p,m,T = Ty, stable = true);
@time hanorm, hsv = ghanorm(sys-sys,fast = fast)
@test hanorm < 1.e-7

@time hanorm, hsv = ghanorm(sys+sys,fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

@time hanorm, hsv = ghanorm([sys sys],fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

sys = rdss(n,p,m,T = Ty, disc = true,stable = true,id=ones(Int,3));
try
    @time hanorm, hsv = ghanorm(sys,fast = fast)
    @test true
catch
    @test false
end

try
    @time hanorm, hsv = ghanorm(sys',fast = fast)
    @test false
catch
    @test true
end

sys = rdss(n,p,m,T = Ty, disc = true,stable = true);
@time hanorm, hsv = ghanorm(sys-sys,fast = fast)
@test hanorm < 1.e-7

@time hanorm, hsv = ghanorm(sys+sys,fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

@time hanorm, hsv = ghanorm([sys sys],fast = fast)
@test norm(hsv[n+1:end]) < 1.e-7

end # fast
end # Ty
end # ghanorm

@testset "glinfnorm & ghinfnorm" begin

sys = rdss(0,0,0);
@time linfnorm, fpeak = glinfnorm(sys)
@test linfnorm == 0 && fpeak == 0

sys = rdss(0,2,2);
@time linfnorm, fpeak = glinfnorm(sys)
@test linfnorm == opnorm(sys.D) && fpeak == 0

sys = rdss(0,2,2,disc=true);
@time linfnorm, fpeak = glinfnorm(sys)
@test linfnorm == opnorm(sys.D) && fpeak == 0

a = [0 1;-1 0]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time hinfnorm, fpeak = ghinfnorm(sys)
@test isinf(hinfnorm) && isnan(fpeak)

a = [0 1;-1 0]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time linfnorm, fpeak = glinfnorm(sys)
@test isinf(linfnorm) && fpeak ≈ 1

a = [0 1;-1 0]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d, Ts = 2);
@time linfnorm, fpeak = glinfnorm(sys)
@test isinf(linfnorm) && fpeak ≈ pi/4

a = [-1 2;0 -2]; e = [0 1;0 0]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,e,b,c,d);
@time linfnorm, fpeak = glinfnorm(sys)
@test isinf(linfnorm) && isinf(fpeak)

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time hinfnorm, fpeak = ghinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ hinfnorm &&
      round(hinfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

sys = dss(a+im*a,b+im*b,c+im*c,d+im*d);
@time hinfnorm, fpeak = glinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ hinfnorm &&
      abs(hinfnorm - 56.237554184492573) < 1.e-6 &&  # compare with Matlab results
      abs(fpeak+3.952242830041314) < 0.001

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(2,2);
sys = dss(e*a,e,e*b,c,d);
@time hinfnorm, fpeak = ghinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ hinfnorm &&
      round(hinfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(Complex{Float64},2,2);
sys = dss(e*(a+im*a),e,e*(b+im*b),c+im*c,d+im*d);
@time hinfnorm, fpeak = glinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ hinfnorm &&
      abs(hinfnorm - 56.237554184492573) < 1.e-6 &&  # compare with Matlab results
      abs(fpeak+3.952242830041314) < 0.001

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time linfnorm, fpeak = glinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

@time linfnorm, fpeak = opnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(2,2);
sys = dss(e*a,e,e*b,c,d);
@time linfnorm, fpeak = glinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
sys = dss(a,b,c,d);
@time linfnorm, fpeak = glinfnorm(sys',rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

a = [-1 2;-3 -2]; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = zeros(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d);
@time linfnorm, fpeak = glinfnorm(sys',rtolinf = 0.0000001)
@test opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(1.163398400218353e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.695926597795647,digits=3)

# discrete standard & descriptor
a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d,Ts = 1);
@time hinfnorm, fpeak = ghinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,exp(im*fpeak))) ≈ hinfnorm &&
      round(hinfnorm, digits=6) ≈ round(4.133509781281479e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.753904640692288,digits=3)

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d,Ts = 1);
@time hinfnorm, fpeak = ghinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,exp(im*fpeak))) ≈ hinfnorm &&
      round(hinfnorm, digits=6) ≈ round(4.133509781281479e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.753904640692288,digits=3)

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d,Ts = 1);
@time linfnorm, fpeak = glinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(4.133509781281479e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.753904640692288,digits=3)

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d,Ts = 1);
@time linfnorm, fpeak = glinfnorm(sys,rtolinf = 0.0000001)
@test opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(4.133509781281479e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.753904640692288,digits=3)

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
sys = dss(a,b,c,d,Ts = 1);
@time linfnorm, fpeak = glinfnorm(sys',rtolinf = 0.0000001)
@test opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(4.133509781281479e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.753904640692288,digits=3)

a = [-1 2;-3 -2]/10; b = [2 3 4; 1 2 3]; c = [1 4; 2 2; 1 3]; d = ones(3,3);
e = rand(2,2)
sys = dss(e*a,e,e*b,c,d,Ts = 1);
@time linfnorm, fpeak = glinfnorm(sys',rtolinf = 0.0000001)
@test opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(4.133509781281479e+01,digits=6) &&
      round(fpeak,digits = 3) ≈ round(2.753904640692288,digits=3)

z = Polynomial([0, 1],'z')
g = [z^2+z+1 4*z^2+3*z+2 2*z^2-2;
    z 4*z-1 2*z-2;
    z^2 4*z^2-z 2*z^2-2*z];

sys = dss(g,minimal = true, atol = 1.e-7,Ts = 1)
@time linfnorm, fpeak = glinfnorm(sys, atol = 1.e-7)
@test glinfnorm(sys, atol = 1.e-7)[1] ≈ glinfnorm(sys', atol = 1.e-7)[1] &&
      opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm &&
      round(linfnorm, digits=6) ≈ round(1.048808848170152e+01,digits=6) &&
      fpeak == 0

n = 50; m = 30; p = 20;

# continuous-time standard, real
Ty = Float64; fast = true;
sys = rss(n,p,m,T = Ty,disc=false);
@time linfnorm, fpeak = glinfnorm(sys, fast = fast, rtolinf = 0.0000000001)
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(real(ev)))[1] < 1.e-5

# discrete-time standard, real
Ty = Float64; fast = true;
sys = rss(n,p,m,T = Ty,disc=true);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(abs.(ev).-1))[1] < 1.e-5

# continuous-time descriptor, real
Ty = Float64; fast = true;
sys = rdss(n,p,m,T = Ty,disc=false);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(real(ev)))[1] < 1.e-5

# discrete-time descriptor, real
Ty = Float64; fast = true;
sys = rdss(n,p,m,T = Ty,disc=true);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(abs.(ev).-1))[1] < 1.e-5

# continuous-time standard, complex
Ty = Complex{Float64}; fast = true;
sys = rss(n,p,m,T = Ty,disc=false);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(real(ev)))[1] < 1.e-5 && opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm

# discrete-time standard, complex
Ty = Complex{Float64}; fast = true;
sys = rss(n,p,m,T = Ty,disc=true);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(abs.(ev).-1))[1] < 1.e-5 && opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm

# continuous-time descriptor, complex
Ty = Complex{Float64}; fast = true;
sys = rdss(n,p,m,T = Ty,disc=false);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(real(ev)))[1] < 1.e-5 && opnorm(evalfr(sys,im*fpeak)) ≈ linfnorm

# discrete-time descriptor, complex
Ty = Complex{Float64}; fast = true;
sys = rdss(n,p,m,T = Ty,disc=true);
@time linfnorm, fpeak = glinfnorm(sys,fast = fast, rtolinf = 0.0000000001);
syst = (linfnorm)^2*I-sys'*sys; ev=gpole(inv(syst),atol=1.e-7);
@test sort(abs.(abs.(ev).-1))[1] < 1.e-5 && opnorm(evalfr(sys,exp(im*fpeak))) ≈ linfnorm

n = 5; m = 3; p = 2;

Ty = Float64; fast = true;
Ty = Complex{Float64}; fast = true;
rtolinf = 0.000001
for Ty in (Float64, Complex{Float64})

for fast in (true, false)

sys = rss(n,p,m,T = Ty, stable = true);
sys.D[:,:] = zeros(Ty,p,m);
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

sys = rss(n,p,m,T = Ty, stable = true);
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

sys = rss(n,p,m,stable = true, disc = true);
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

sys = rdss(n,p,m,T = Ty, stable = true);
sys.D[:,:] = zeros(Ty,p,m);
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

sys = rdss(n,p,m,T = Ty, stable = true);
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

sys = rdss(n,p,m,T = Ty, stable = true, disc = true);
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

sys = rdss(n,p,m,T = Ty, stable = true, id=ones(Int,3));
@time hinf, hfpeak = ghinfnorm(sys,fast = fast,rtolinf = rtolinf)
@time linf, lfpeak = glinfnorm(sys',fast = fast,rtolinf = rtolinf)
@test abs(hinf-linf)/hinf < 2*rtolinf && (hfpeak ≈ lfpeak || abs(hfpeak-lfpeak) <= 0.01*abs(lfpeak))

end # fast
end # Ty
end # glinfnorm

@testset "gnugap" begin

## Examples from the Vinnicombe's book, p. 115
s = rtf('s'); N = 3;
sys = (dss(100/(2*s+1)), dss(100/(2*s-1)), dss(100/(s+1)^2))
distgap = zeros(N);
for j = 1:N
    distgap[j], = gnugap(sys[1],sys[j],atol=1.e-7);
end
@test norm(distgap-[0,0.02,0.8988]) <0.001

## Examples from the Vinnicombe's book, p. 125
s = rtf('s');
sys1 = dss(1/(s+0.2));                               # P1
sys2 = dss((s+0.1)^2/(s-0.2)/(s^2+0.1*s+0.01));      # P2
nugap, fpeak = gnugap(sys2,sys1,atol=1.e-7)
@test abs(nugap-0.403) < 0.0001 && abs(fpeak-0.32) < 0.01

## Examples from the Vinnicombe's book, p. 160 with wno ~= 0
s = rtf('s'); N = 3;
sys = (dss(0.3*(s-1)/(s^2+1)), dss(0.3*(s-1)/(s^2+2)), dss(0.3*(s-1)*(s-2)/(s^2+2)/(s-sqrt(2))))
distgap = zeros(N);
for j = 1:N
    distgap[j], = gnugap(sys[1],sys[j],atol=1.e-7);
end

@test norm(distgap-[0,1,1]) < 1.e-6

## Examples from the Vinnicombe's book, p. 201
s = rtf('s');
sys1 = dss(1/s^2);    # P1
sys2 = dss(1/s);      # P2
nugap, fpeak = gnugap(sys2,sys1,atol=1.e-7)
@test nugap ≈ 1/sqrt(2) && fpeak ≈ 1

## Examples from the Vinnicombe's book, p. 202  nugap = 1/sqt(2) for epsi>0
s = rtf('s'); epsi = 0.01;
sys1 = dss([1/s; 1/(s+epsi)]);    # P1
sys2 = dss([1/s; 1/s]);           # P2
nugap, = gnugap(sys2,sys1,atol=1.e-7)
@test nugap ≈ 1/sqrt(2)

## Examples from the Vinnicombe's book, p. 202  modified
#  realization of sys2 is not stabilizable
s = rtf('s'); epsi = 0.01;
sys1 = dss([1/s^2; 1/(s+epsi)]);    # P1
sys2 = dss([1/s^2; 1/s]);           # P2
nugap, = gnugap(sys2,sys1,atol=1.e-7)
@test nugap < epsi

## Examples from the Vinnicombe's book, p. 202  modified
#  realization of sys2 is not stabilizable
s = rtf('s'); epsi = 0.01;
sys1 = dss([s^2+epsi; (s+epsi)]);    # P1
sys2 = dss([s^2; s]);           # P2
nugap, = gnugap(sys2,sys1,atol=1.e-7)
@test nugap < sqrt(2)*epsi

##
z = rtf('z');
sys1 = dss(1/z); sys2 = dss(0.5/z);
nugap, = gnugap(sys1,sys2,atol=1.e-6)
@test nugap ≈ 0.316227766016838

##
z = rtf('z');
sys1 = dss(1/z); sys2 = dss(0.5/z^2);
nugap, = gnugap(sys1,sys2,atol=1.e-6)
@test nugap ≈ 0.948683298050514

## ν-gap = 1
s = rtf('s');
sys1 = dss((s+1)/(s+2)); sys2 = dss(-1);
nugap, = gnugap(sys1,sys2,atol=1.e-7)
@test nugap == 1

## gapmetric produces the wrong result nugap = 0
sys1 = dss(1); sys2 = dss(-1);
nugap, = gnugap(sys1,sys2,atol=1.e-5)
@test nugap == 1

## ν-gap = 1
sys1 = rss(3,2,2); sys1.D[:,:] = eye(2); sys2 = -inv(sys1)';
nugap, fpeak = gnugap(sys1,sys2,atol=1.e-6)
@test nugap == 1

n = 5; m = 3; p = 2;

Ty = Float64; fast = true;
Ty = Complex{Float64}; fast = true;
rtolinf = 0.000001
for Ty in (Float64, Complex{Float64})

for fast in (true, false)

## nugap = 0
sys = rss(n,p,m, T = Ty);
@time nugap, fpeak = gnugap(sys,sys,fast = fast, atol = 1.e-7)
@test nugap < 1.e-10

sys = rss(n,p,m,T = Ty, disc=true);
@time nugap, fpeak = gnugap(sys,sys, fast = fast, atol = 1.e-7)
@test nugap < 1.e-10

sys = rdss(n,p,m,T = Ty);
@time nugap, fpeak = gnugap(sys,sys,fast = fast, atol = 1.e-7)
@test nugap < 1.e-10

sys = rdss(n,p,m,T = Ty, disc=true);
@time nugap, fpeak = gnugap(sys,sys, fast = fast, atol = 1.e-7)
@test nugap < 1.e-10

sys = rss(n,p,m, T = Ty);
@time nugap, fpeak = gnugap(sys,sys,fast = fast, atol = 1.e-7, freq = collect(range(0.1, 0.3, length=5)))
@test maximum(nugap) < 1.e-10

sys = rss(n,p,m,T = Ty, disc=true);
@time nugap, fpeak = gnugap(sys, sys, fast = fast, atol = 1.e-7, freq = collect(range(0.1, 0.3, length=5)))
@test maximum(nugap) < 1.e-10

sys = rdss(n,p,m,T = Ty);
@time nugap, fpeak = gnugap(sys, sys,fast = fast, atol = 1.e-7, freq = collect(range(0.1, 0.3, length=5)))
@test maximum(nugap) < 1.e-10

sys = rdss(n,p,m,T = Ty, disc=true);
@time nugap, fpeak = gnugap(sys, sys, fast = fast, atol = 1.e-7, freq = collect(range(0.1, 0.3, length=5)))
@test maximum(nugap) < 1.e-10

sys1 = rss(n,p,m,T = Ty); sys1d = sys1*0.001/glinfnorm(sys1)[1]; sys2 = sys1+sys1d;
@time nugap, fpeak = gnugap(sys2,sys1,fast = fast, atol = 1.e-7)
@time nugap1, fpeak1 = gnugap(sys2,sys1,fast = fast, atol = 1.e-7,freq=fpeak)
@test nugap ≈ nugap1[1] &&  nugap < 0.01

sys1 = rss(n,p,m,T = Ty); sys1d = rss(3,2,3,stable=true,T = Ty); sys2 = sys1+sys1d*0.001/glinfnorm(sys1d)[1];
@time nugap, fpeak = gnugap(sys2,sys1,fast = fast, atol = 1.e-7)
@time nugap1, fpeak1 = gnugap(sys2,sys1,fast = fast, atol = 1.e-7,freq=fpeak)
@test nugap ≈ nugap1[1] &&  nugap < 0.01

@time nugap, fpeek = gnugap(sys1,sys2,fast = fast, atol=1.e-7)
sysd1 = gbilin(sys1,rtfbilin("c2d")[1])[1]; sysd2 = gbilin(sys2,rtfbilin("c2d")[1])[1];
@time nugap1, fpeak1 = gnugap(sysd2,sysd1,fast = fast, atol=1.e-7)
@test abs(nugap-nugap1) < 0.001

sysc1 = gbilin(sysd1,rtfbilin("d2c")[1])[1]; sysc2 = gbilin(sysd2,rtfbilin("d2c")[1])[1];
nugap2, fpeak2 = gnugap(sysc2,sysc1,fast = fast, atol=1.e-7)
@test abs(nugap-nugap2) < 0.0001

end # fast
end # Ty
end # gnugap

end # module


text:
using Parameters
using Unitful: s, ms, ustrip

# question: should we somehow limit the unit response???
# e.g. with:
sig(x) = 1/(1+exp(-10(x-0.5)))

@with_kw struct AdaptMI{S,I}
  c_x::Float64 = 1.0
  τ_x::typeof(1.0s) = 300ms

  shape_y::S = identity

  c_a::Float64 = 5
  τ_a::typeof(1.0s) = 1.5s

  c_m::Float64 = 10
  τ_m::typeof(1.0s) = 50ms
  W_m::I = identity

  τ_σ::typeof(1.0s) = 100ms
  c_σ::Float64 = 0.3
end

################################################################################
# generic adaptation and mutual-inhibition operation
adaptmi(x;progressbar=true,kw...) = adaptmi(x,AdaptMI(;kw...),progressbar)
function adaptmi(x,params::AdaptMI,progressbar=true)
  @assert :time ∈ axisnames(x)
  time = Axis{:time}

  y = similar(x)

  shape_y = params.shape_y
  τ_x = params.τ_x; c_x = params.c_x
  τ_a = params.τ_a; c_a = params.c_a
  τ_m = params.τ_m; c_m = params.c_m; W_m = params.W_m

  # the (a)daptation (m)utual inhibition
  a = similar(y) # a = adaptation
  m = similar(y) # m = mutual inhibition
  x̃ = similar(y) # (this is the smoothed version of x)

  # temporary, single time slice variables
  x̃_t = zero(y[time(1)])
  y_t = copy(x̃_t)
  a_t = copy(x̃_t)
  m_t = copy(x̃_t)

  dt_x = eltype(a_t)(Δt(x) / τ_x)
  dt_a = eltype(a_t)(Δt(x) / τ_a)
  dt_m = eltype(a_t)(Δt(x) / τ_m)

  progress = progressbar ? Progress(desc="Adapt/Inhibit: ",ntimes(y)) : nothing
  for ti in Base.axes(times(y),1)
    t = time(ti)

    @. begin # across all indices...
      # smooth input
      x̃_t += (c_x*x[t] - x̃_t)*dt_x

      # apply adaptation and inhibition
      y_t = shape_y((1 - c_a*a_t)*x̃_t - c_m*m_t)

      # update adaptation and inhibition
      a_t += (y_t - a_t)*dt_a
      m_t += ($W_m(y_t) - m_t)*dt_m
    end
    x̃[t] = x̃_t
    y[t] = y_t
    a[t] = a_t
    m[t] = m_t

    next!(progress)
  end
  y,a,m,x̃
end

################################################################################
# generic drifting noise function

drift(x,along_axes...;progressbar=true,kw...) =
  drift(x,AdaptMI(;kw...),along_axes,progressbar)
function drift(x,params::AdaptMI,along_axes=typeof.(AxisArayys.axes(x)),
    progressbar=true)

  τ_σ, c_σ = params.τ_σ, log(1+params.c_σ)
  time = Axis{:time}

  y = similar(x)
  σ_t = zero(x[time(1),(ax(1) for ax in along_axes)...])
  dims = size(σ_t)
  progress = progressbar ? Progress(desc="Drift: ",ntimes(y)) : nothing
  for t in Base.axes(times(y),1)
    ya_t = x[time(t)]

    σ_t .+= -σ_t.*(Δt(x)/τ_σ) .+ randn(dims).*(c_σ*sqrt(2Δt(x)/τ_σ))
    ya_t .*= exp.(σ_t)

    y[time(t)] = ya_t

    next!(progress)
  end
  y
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
function alilqr_allocs(solver)
    cost_tol = solver.opts.cost_tolerance
    grad_tol = solver.opts.gradient_tolerance
    al_iter = 1

    allocs = @allocated Altro.reset!(solver)
    allocs += @allocated Altro.set_tolerances!(solver, al_iter, cost_tol, grad_tol)

    cons = Altro.get_constraints(solver)
    Altro.settraj!(cons, solver.ilqr.Z̄)
    let solver = solver.ilqr
        allocs += @allocated J_prev = TO.cost(solver, solver.Z̄)
        allocs += @allocated Altro.errstate_jacobians!(solver.model, solver.G, solver.Z̄)
        allocs += @allocated Altro.dynamics_expansion!(solver, solver.Z̄)
        allocs += @allocated Altro.cost_expansion!(solver.obj, solver.Efull, solver.Z̄)
        allocs += @allocated Altro.error_expansion!(solver.model, solver.Eerr, solver.Efull, solver.G, solver.Z̄)
        allocs += @allocated Altro.backwardpass!(solver)
        allocs += @allocated Jnew = Altro.forwardpass!(solver, J_prev)
        allocs += @allocated copyto!(solver.Z, solver.Z̄)
        dJ = J_prev - Jnew
        allocs += @allocated grad = Altro.gradient!(solver)
        Altro.record_iteration!(solver, Jnew, dJ, grad)
        allocs += @allocated Altro.evaluate_convergence(solver)
    end
    Z̄ = solver.ilqr.Z̄
    conset = get_constraints(solver)
    allocs += @allocated J = TO.cost(solver, Z̄)
    allocs += @allocated c_max = Altro.max_violation(conset)
    allocs += @allocated μ_max = Altro.max_penalty(conset)
    Altro.record_iteration!(solver, J, c_max, μ_max)

    allocs += @allocated Altro.dualupdate!(conset)
    allocs += @allocated Altro.penaltyupdate!(conset)
    allocs += @allocated Altro.reset!(solver.ilqr)
    allocs
end

@testset "AliLQR Basic Tests" begin
prob, opts = Problems.Quadrotor()
al = Altro.ALSolver(prob, copy(opts), use_static=Val(true))
conset = get_constraints(al)
ilqr = Altro.get_ilqr(al)
@test conset[1].E === ilqr.Efull
@test conset[1].Z[1] === ilqr.Z
@test conset[1].opts.solveropts === ilqr.opts === al.opts
@test conset[1].cost === get_objective(ilqr).alcost

@testset "ALiLQR Allocation test" begin
    solver = Altro.ALSolver(Problems.Quadrotor()...)
    solver.opts.verbose = 0
    alilqr_allocs(solver)
    !Sys.iswindows() && @test alilqr_allocs(solver) == 0
end
end


text:
export iauS2pv
"""
Convert position/velocity from spherical to Cartesian coordinates.

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

Status:  vector/matrix support function.

Given:
   theta    double          longitude angle (radians)
   phi      double          latitude angle (radians)
   r        double          radial distance
   td       double          rate of change of theta
   pd       double          rate of change of phi
   rd       double          rate of change of r

Returned:
   pv       double[2][3]    pv-vector

This revision:  2013 June 18

SOFA release 2018-01-30

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

# void iauS2pv(double theta, double phi, double r,
#              double td, double pd, double rd,
#              double pv[2][3])

function iauS2pv(theta::Real, phi::Real, r::Real,
                td::Real, pd::Real, rd::Real)
   pv = zeros(Float64, 3, 2)

   ccall((:iauS2pv, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Cdouble,
       Cdouble, Cdouble, Cdouble,
       Ptr{Cdouble}),
       convert(Float64, theta),
       convert(Float64, phi),
       convert(Float64, r),
       convert(Float64, td),
       convert(Float64, pd),
       convert(Float64, rd),
       pv)

   return SMatrix{2,3}(pv')
end


text:
using Test, OpenEXR, FileIO, Colors, FixedPointNumbers

# helpers for testing channelwise bit-depth truncation
function channelwise_max_diff(a::RGBA, b::RGBA)
    max(abs(a.r - b.r), abs(a.g - b.g), abs(a.b - b.b), abs(a.alpha - b.alpha))
end

function channelwise_max_diff(a::RGB, b::RGB)
    max(abs(a.r - b.r), abs(a.g - b.g), abs(a.b - b.b))
end

function channelwise_max_diff(a::GrayA, b::GrayA)
    max(abs(a.val - b.val), abs(a.alpha - b.alpha))
end

function channelwise_max_diff(a::Gray, b::Gray)
    abs(a.val - b.val)
end

@testset "RoundTrip" begin

    @testset "Identical" begin
        for typ in (RGBA{Float16}, RGB{Float16}, GrayA{Float16}, Gray{Float16})
            img = rand(typ, 256, 512)
            fn = File{DataFormat{:EXR}}(tempname())
            OpenEXR.save(fn, img)
            try
                loaded_img = OpenEXR.load(fn)
                @test typeof(loaded_img) === Array{typ,2}
                @test img == loaded_img
            finally
                rm(fn.filename)
            end
        end
    end

    @testset "Lossless with type conversion" begin
        for (save_type, load_type) in (
            (RGBA{N0f8}, RGBA{Float16}),
            (RGB{N0f8}, RGB{Float16}),
            (GrayA{N0f8}, GrayA{Float16}),
            (Gray{N0f8}, Gray{Float16}),
        )
            img = rand(save_type, 256, 512)
            fn = File{DataFormat{:EXR}}(tempname())
            OpenEXR.save(fn, img)
            try
                loaded_img = OpenEXR.load(fn)
                @test typeof(loaded_img) === Array{load_type,2}
                converted_img = (c -> convert(save_type, c)).(loaded_img)
                @test converted_img == img
            finally
                rm(fn.filename)
            end
        end
    end

    @testset "Lossy with type conversion" begin

        @testset "Bit depth truncation" begin
            for (save_type, load_type) in (
                (RGBA{Float32}, RGBA{Float16}),
                (RGB{Float32}, RGB{Float16}),
                (GrayA{Float32}, GrayA{Float16}),
                (Gray{Float32}, Gray{Float16}),
            )
                img = rand(save_type, 256, 512)
                fn = File{DataFormat{:EXR}}(tempname())
                OpenEXR.save(fn, img)
                try
                    loaded_img = OpenEXR.load(fn)

                    # return type is as expected
                    @test typeof(loaded_img) === Array{load_type,2}

                    # differences are bounded channelwise
                    diffs = map(channelwise_max_diff, img, loaded_img)
                    @test maximum(diffs) <= 2.5e-4
                finally
                    rm(fn.filename)
                end
            end
        end

        @testset "Colorspace conversion" begin
            for save_type in
                (HSV{Float16}, HSL{Float16}, Lab{Float16}, LCHab{Float16}, YIQ{Float16})
                img = rand(save_type, 256, 512)
                fn = File{DataFormat{:EXR}}(tempname())
                OpenEXR.save(fn, img)
                try
                    loaded_img = OpenEXR.load(fn)

                    # return type is as expected
                    @test typeof(loaded_img) === Array{RGB{Float16},2}

                    # differences are bounded channelwise
                    converted_img = (c -> convert(RGB{Float16}, c)).(img)
                    diffs = map(channelwise_max_diff, converted_img, loaded_img)
                    @test maximum(diffs) <= 2.5e-4
                finally
                    rm(fn.filename)
                end
            end
        end
    end
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using BenchmarkTools, BlockDiagonals, FillArrays, LinearAlgebra, Random, Stheno,
    TemporalGPs, Zygote

using DataFrames, DrWatson, PGFPlotsX

using TemporalGPs: AV, AM, Separable, RectilinearGrid, LGSSM, GaussMarkovModel

# const n_blas_threads = Sys.CPU_THREADS;
const n_blas_threads = 4;
BLAS.set_num_threads(n_blas_threads);

const exp_dir_name = "lgssm"

#
# Utilities
#

function load_settings()
    return load(joinpath(datadir(), exp_dir_name, "meta", "settings.bson"))[:settings]
end

load_results() = collect_results(joinpath(datadir(), exp_dir_name))

function dense_dynamics_constructor(rng, N_space, N_time, N_blocks)

    # Build kernel with block-diagonal transition dynamics.
    ft = block_diagonal_dynamics_constructor(rng, N_space, N_time, N_blocks)

    # Densify dynamics and construct new LGSSM.
    return LGSSM(
        GaussMarkovModel(
            Fill(collect(first(ft.gmm.A)), length(ft)),
            ft.gmm.a,
            Fill(collect(first(ft.gmm.Q)), length(ft)),
            ft.gmm.H,
            ft.gmm.h,
            ft.gmm.x0,
        ),
        ft.Σ,
    )
end

function block_diagonal_dynamics_constructor(rng, N_space, N_time, N_blocks)

    # Construct kernel.
    k = Separable(SEKernel(), Matern52Kernel())
    for n in 1:(N_blocks - 1)
        k += k
    end

    f = to_sde(GP(k, GPC()))
    t = range(-5.0, 5.0; length=N_time)
    x = randn(rng, N_space)
    return f(RectilinearGrid(x, t), 0.1)
end

#
# Save settings.
#

wsave(
    joinpath(datadir(), exp_dir_name, "meta", "settings.bson"),
    :settings => dict_list(
        Dict(

            # Method to execute `predict` and construct / evaluate its pullback.
            :implementation => [
                (
                    name = "naive",
                    dynamics_constructor = dense_dynamics_constructor,
                ),
                (
                    name = "dense",
                    dynamics_constructor = dense_dynamics_constructor,
                ),
                (
                    name = "block-diagonal",
                    dynamics_constructor = block_diagonal_dynamics_constructor,
                ),
            ],

            # Number of observations in space.
            :N_space => [247],

            # Number of time points.
            # :N_time => [2, 5, 10, 50, 100, 500, 1_000, 50_000],
            :N_time => [25, 50, 75, 100],

            # Number of things to sum in the kernel.
            :N_blocks => [1, 2, 3],
        ),
    ),
)

#
# Run benchmarking experiments.
#

let
    settings = load_settings()
    N_settings = length(settings)
    for (n, setting) in enumerate(settings)

        # Display current iteration.
        impl = setting[:implementation]
        N_space = setting[:N_space]
        N_time = setting[:N_time]
        N_blocks = setting[:N_blocks]
        println(
            "$n / $N_settings: name=$(impl.name), N_space=$N_space, " *
            "N_time=$N_time, N_blocks = $N_blocks",
        )

        # Build dynamics model.
        rng = MersenneTwister(123456)
        ft = impl.dynamics_constructor(rng, N_space, N_time, N_blocks)
        y = rand(rng, ft)

        # Benchmark rand evaluation.
        rand_results = @benchmark rand($rng, $ft)

        # Benchmark logpdf evaluation and gradient evaluation.
        logpdf_results = @benchmark logpdf($ft, $y)
        logpdf_gradient_results = @benchmark Zygote.gradient(logpdf, $ft, $y)

        # Save results to disk.
        wsave(
            joinpath(
                datadir(),
                exp_dir_name,
                savename(
                    Dict(
                        :name => impl.name,
                        :N_space => N_space,
                        :N_time => N_time,
                        :N_blocks => N_blocks,
                    ),
                    "bson",
                ),
            ),
            Dict(
                :setting => setting,
                :rand_results => rand_results,
                :logpdf_results => logpdf_results,
                :logpdf_gradient_results => logpdf_gradient_results,
                :n_blas_threads => n_blas_threads,
                :n_julia_threads => Threads.nthreads(),
            ),
        )
    end
end

#
# Process results.
#

let
    # Define some plotting settings for each of the implementation types.
    colour_map = Dict(
        "naive" => "black",
        "dense" => "blue",
        "block-diagonal" => "red",
    )
    marker_map = Dict(
        "naive" => "square",
        "dense" => "triangle",
        "block-diagonal" => "*",
    )

    try
        mkdir(plotsdir())
    catch
    end

    try
        mkdir(joinpath(plotsdir(), exp_dir_name))
    catch
    end

    # Load the results and compute some additional columns.
    results = DataFrame(
        map(eachrow(load_results())) do row
            row = merge(
                copy(row), # No need to copy in following blocks.
                (
                    implementation_name = row.setting[:implementation].name,
                    N_space = row.setting[:N_space],
                    N_time = row.setting[:N_time],
                    N_blocks = row.setting[:N_blocks],
                ),
            )
            return row
        end
    )

    # Plot timing curves for each block size for each implementation. This is used to assess
    # after what value of N_time the time-per-iteration saturates.
    by(results, [:N_blocks, :implementation_name]) do group

        # Specify plots to produce.
        N_blocks = first(group.N_blocks)
        impl_name = first(group.implementation_name)
        setting_name = "N_blocks=$N_blocks-impl_name=$impl_name"
        plots = [
            (
                plot = (@pgf Axis(
                    {
                        title = "rand",
                        legend_pos = "north west",
                        xmode="log",
                    }
                )),
                result_name = :rand_results,
                save_name = "rand_$setting_name.tex",
            ),
            (
                plot = (@pgf Axis(
                    {
                        title = "logpdf",
                        legend_pos = "north west",
                        xmode="log",
                    },
                )),
                result_name = :logpdf_results,
                save_name = "logpdf_$setting_name.tex",
            ),
            (
                plot = (@pgf Axis(
                    {
                        title = "logpdf gradient",
                        legend_pos = "north west",
                        xmode="log",
                    },
                )),
                result_name = :logpdf_gradient_results,
                save_name = "logpdf_gradient_$setting_name.tex",
            ),
        ]

        by(group, :N_space) do space_group

            space_group = sort(space_group, :N_time)
            for plt in plots
                push!(
                    plt.plot,
                    (@pgf Plot(
                        {
                            color = colour_map[impl_name],
                            mark_options={fill=colour_map[impl_name]},
                            mark=marker_map[impl_name],
                        },
                        Coordinates(
                            space_group.N_time,
                            time.(space_group[!, plt.result_name]) ./ space_group.N_time,
                        ),
                    )),
                )
                push!(plt.plot, LegendEntry(string(first(space_group.N_space))))
            end
        end

        # Save plots.
        for plt in plots
            pgfsave(
                joinpath(plotsdir(), exp_dir_name, plt.save_name),
                plt.plot;
                include_preamble=true,
            )
        end
    end
end

#
# Hacked together benchmarks for playing around.
#

rng = MersenneTwister(123456);
Ts = [1, 10, 100, 1_000];
N_space = 500;
N_blocks = 1;

model = dense_dynamics_constructor(rng, N_space, 100, N_blocks);
y = rand(rng, model);
@benchmark rand($rng, $model)
@benchmark logpdf($model, $y)

using Profile, ProfileView

@profview logpdf(model, y)
@profview logpdf(model, y)

# Test simple things quickly.
rng = MersenneTwister(123456);
T = 1_000_000;
x = range(0.0; step=0.3, length=T);
f = GP(Matern52Kernel() + Matern52Kernel() + Matern52Kernel() + Matern52Kernel(), GPC());
fx_sde_dense = to_sde(f)(x);
fx_sde_static = to_sde(f, SArrayStorage(Float64))(x);

y = rand(fx_sde_static);
@benchmark logpdf($fx_sde_dense, $y)
@benchmark logpdf($fx_sde_static, $y)

using Profile, ProfileView

@profview logpdf(fx_sde_dense, y)
@profview logpdf(fx_sde_dense, y)


text:
"""
    MultiNLI

A corpus of 433k sentence pairs for NLI.

For details, see the [MultiNLI home page](https://www.nyu.edu/projects/bowman/multinli/)
or read the [2018 paper](https://www.nyu.edu/projects/bowman/multinli/paper.pdf) "A Broad-Coverage Challenge Corpus for Sentence Understanding through Inference" by Adina Williams, NIkita Nangia, and Samuel R. Bowman.

Included data:
```julia
MultiNLI.train_tsv()
MultiNLI.train_jsonl()
MultiNLI.dev_matched_tsv()
MultiNLI.dev_matched_jsonl()
MultiNLI.dev_mismatched_tsv()
MultiNLI.dev_mismatched_jsonl()
```
"""
module MultiNLI

import ...register_data
using DataDeps: @datadep_str

multinli_file(filename) = joinpath(datadep"MultiNLI", "multinli_1.0", filename)

train_tsv()            = multinli_file("multinli_1.0_train.txt")
train_jsonl()          = multinli_file("multinli_1.0_train.jsonl")
dev_matched_tsv()      = multinli_file("multinli_1.0_dev_matched.txt")
dev_matched_jsonl()    = multinli_file("multinli_1.0_dev_matched.jsonl")
dev_mismatched_tsv()   = multinli_file("multinli_1.0_dev_mismatched.txt")
dev_mismatched_jsonl() = multinli_file("multinli_1.0_dev_mismatched.jsonl")

function __init__()
    register_data(
        "MultiNLI",
        """
        The Multi-Genre Natural Language Inference (MultiNLI) corpus
        is a crowd-sourced collection of 433k sentence pairs annotated
        with textual entailment information.

        https://www.nyu.edu/projects/bowman/multinli/

        The dataset is distributed as a 227MB zip file.

        For license details, see details the data description paper:
        https://www.nyu.edu/projects/bowman/multinli/paper.pdf
        """,
        "https://www.nyu.edu/projects/bowman/multinli/multinli_1.0.zip",
        "049f507b9e36b1fcb756cfd5aeb3b7a0cfcb84bf023793652987f7e7e0957822")
end

end # module


text:
using Dates
using UUIDs

export SolverStatus

"""
A simple structure for monitoring a running solver.
"""
mutable struct SolverStatus
    id::String
    host::String
    isAttached::Bool
    userId::String
    robotId::String
    sessionId::String
    iteration::Int64
    currentStep::String
    startedUpTimestamp::String
    attachedSessionTimestamp::String
    detachedSessionTimestamp::String
    lastUpdatedTimestamp::String
    lastIterationDurationSeconds::Float64
    result::String
    SolverStatus() = new(string(uuid4()), "", false, "", "", "", 0, "", string(unix2datetime(time())), "", "", "", 0, "")
    SolverStatus(id::String, host::String, isAttached::Bool, userId::String, robotId::String, sessionId::String, iteration::Int64, currentStep::String, startedUpTimestamp::String, attachedSessionTimestamp::String, detachedSessionTimestamp::String, lastUpdatedTimestamp::String, lastIterationDurationSeconds::Float64, result::String) = new(id, host, isAttached, userId, robotId, sessionId, iteration, currentStep, startedUpTimestamp, attachedSessionTimestamp, detachedSessionTimestamp, lastUpdatedTimestamp, lastIterationDurationSeconds, result)
end


text:
using Test
using GridCreation

@testset "fft alignment" begin
    xfft, yfft = mkCartVecs(1., 4, center=CFFT);
    # alignment
    @test xfft[3] ≈ 0
    @test yfft[3] ≈ 0

    # spacing
    @test xfft[1] ≈ -2
    @test yfft[1] ≈ -2
end

@testset "intersample alignment" begin
    xinter, yinter = mkCartVecs(1., 4, center=CInterSample);
    @test xinter[2] ≈ -0.5
    @test xinter[3] ≈  0.5
    @test yinter[2] ≈ -0.5
    @test yinter[3] ≈  0.5
end

@testset "first sample centering" begin
# test sample index centering is properly aligned
    x0, y0 = mkCartVecs(1., 4, center=1);
    @test x0[1] ≈ 0
    @test y0[1] ≈ 0
end

@testset "polar grid" begin
    xfft, yfft = mkCartVecs(1., 4, center=CFFT);
    ρ, θ = cartVecsToPolarGrid(xfft, yfft);
    @test ρ[3,3] ≈ 0
    @test ρ[3,2] ≈ 1
    @test ρ[2,3] ≈ 1
    @test ρ[3,4] ≈ 1
    @test ρ[4,3] ≈ 1
    @test ρ[4,4] ≈ √2
    @test ρ[2,2] ≈ √2
    @test ρ[2,4] ≈ √2
    @test ρ[4,2] ≈ √2
    @test θ[3,3] ≈ 0
    @test θ[3,4] ≈ 0
    @test θ[3,2] ≈ π
    @test θ[3,1] ≈ π
    @test θ[2,3] ≈ -π/2
    @test θ[4,3] ≈ π/2
end

max_stars_repo_name: JuliaOptics/Gridding

max_issues_repo_name: JuliaOptics/Gridding

max_forks_repo_name: JuliaOptics/Gridding


text:
using Plots, Test, Dates

@testset "Limits" begin
    y=[1.0*i*i for i in 1:10]
    x=[Date(2019,11,i) for i in 1:10]

    rx=[x[3],x[5]]

    p = plot(x,y, widen = false)
    vspan!(p, rx, label="", alpha=0.2)

    ref_ylims = (y[1], y[end])
    ref_xlims = (x[1].instant.periods.value, x[end].instant.periods.value)
    @test Plots.ylims(p) == ref_ylims
    @test Plots.xlims(p) == ref_xlims
    #@static if (haskey(ENV, "APPVEYOR") || haskey(ENV, "CI"))
    @static if haskey(ENV, "APPVEYOR")
        @info "Skipping display tests on AppVeyor"
    else
        @test isa(display(p), Nothing) == true
        closeall()
    end
end # testset

@testset "Date xlims" begin
    y=[1.0*i*i for i in 1:10]
    x=[Date(2019,11,i) for i in 1:10]
    span = (Date(2019,10,31), Date(2019,11,11))

    ref_xlims = map(date->date.instant.periods.value, span)

    p = plot(x,y, xlims=span, widen = false)

    @test Plots.xlims(p) == ref_xlims
    #@static if (haskey(ENV, "APPVEYOR") || haskey(ENV, "CI"))
    @static if haskey(ENV, "APPVEYOR")
        @info "Skipping display tests on AppVeyor"
    else
        @test isa(display(p), Nothing) == true
        closeall()
    end
end # testset

@testset "DateTime xlims" begin
    y=[1.0*i*i for i in 1:10]
    x=[DateTime(2019,11,i,11) for i in 1:10]
    span = (DateTime(2019,10,31,11,59,59), DateTime(2019,11,11,12,01,15))

    ref_xlims = map(date->date.instant.periods.value, span)

    p = plot(x,y, xlims=span, widen = false)
    @test Plots.xlims(p) == ref_xlims
    #@static if (haskey(ENV, "APPVEYOR") || haskey(ENV, "CI"))
    @static if haskey(ENV, "APPVEYOR")
        @info "Skipping display tests on AppVeyor"
    else
        @test isa(display(p), Nothing) == true
        closeall()
    end
end # testset


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "OsiBuilder"
version = v"0.107.9"

# Collection of sources required to build OsiBuilder
sources = [
    "https://github.com/coin-or/Osi/archive/releases/0.107.9.tar.gz" =>
    "e2c8a0ee4a2a0abe7475d67f7f98230e8bfbbcb6e74487877e757c996bfd6d30",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd Osi-releases-0.107.9/
update_configure_scripts
# temporary fix
for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done
mkdir build
cd build/
export CXXFLAGS="-std=c++11"

## STATIC BUILD START
if [ $target = "x86_64-apple-darwin14" ]; then
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
--with-lapack-lib="-L${prefix}/lib -lcoinlapack" \
--with-blas-lib="-L${prefix}/lib -lcoinblas"
## STATIC BUILD END

## DYNAMIC BUILD START
#../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
#--with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
#--with-lapack="-L${prefix}/lib -lcoinlapack" \
#--with-blas="-L${prefix}/lib -lcoinblas"
## DYNAMIC BUILD END

make -j${nproc}
make install

"""
#export LDFLAGS="-L${prefix}/lib -lCoinUtils -lcoinglpk"
#--with-glpk-lib="-L${prefix}/lib -lcoinglpk" --with-glpk-incdir="$prefix/include/coin/ThirdParty" \
# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libOsi", :libOsi),
    LibraryProduct(prefix, "libOsiCommonTests", :libOsiCommonTests)
]
platforms = expand_gcc_versions(platforms)
# To fix gcc4 bug in Windows
push!(platforms, Windows(:i686,compiler_abi=CompilerABI(:gcc6)))
push!(platforms, Windows(:x86_64,compiler_abi=CompilerABI(:gcc6)))

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaOpt/CoinUtilsBuilder/releases/download/v2.10.14-1-static/build_CoinUtilsBuilder.v2.10.14.jl",
    "https://github.com/JuliaOpt/COINBLASBuilder/releases/download/v1.4.6-1-static/build_COINBLASBuilder.v1.4.6.jl",
    "https://github.com/JuliaOpt/COINLapackBuilder/releases/download/v1.5.6-1-static/build_COINLapackBuilder.v1.5.6.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
import MXNet
import POMDPs: Solver

type DeepCrossentropyPolicySolver <: Solver
    epochs::Int64
    n_samples::Int64
    percentile::Float64
    max_frame_iterations::Int64
    nnet_spec::Vector{Int64}
    print_every_n::Int64
end

DeepCrossentropyPolicySolver(epochs = 100, n_samples = 100, percentile = 0.5; max_frame_iterations = 100, nnet_spec = [20, 10], print_every_n = 10) =
    DeepCrossentropyPolicySolver(epochs, n_samples, percentile, max_frame_iterations, nnet_spec, print_every_n)

function solve(solver::DeepCrossentropyPolicySolver, pomdp::MDP; verbose = true)

    policy = DeepCrossentropyPolicy(pomdp; nnet_spec = solver.nnet_spec)
    optimizer = mx.ADAM(lr=0.01)

    for i=1:solver.epochs

        # println(n_states(pomdp), solver.n_samples * solver.max_frame_iterations)

        batch_states  = zeros(Float64, n_states(pomdp), solver.n_samples * solver.max_frame_iterations) # a list of lists of states in each session
        batch_actions = zeros(Int64, solver.n_samples * solver.max_frame_iterations) # a list of lists of actions in each session
        batch_rewards_states = zeros(Float64, solver.n_samples * solver.max_frame_iterations)
        batch_rewards = zeros(Float64, solver.n_samples) # a list of floats - total rewards at each session

        batch_records_index = 0

        for j=1:solver.n_samples

            if verbose && (j == 1 || j % solver.print_every_n == 0)
                println("Epoch: $i, Episode: $j")
            end

            session_states, session_actions, session_reward, session_states_new, session_rewards = run_experiment(pomdp, policy; max_frame_iterations = solver.max_frame_iterations, keep_history = true)

            index_start = batch_records_index + 1
            index_end = batch_records_index + size(session_states)[2]

            batch_states[:, index_start:index_end] = session_states
            batch_actions[index_start:index_end] = session_actions
            batch_rewards_states[index_start:index_end] = session_reward
            batch_rewards[j] = session_reward

            batch_records_index = index_end
        end

        # Identify threshold which filters out best performing games
        reward_threshold = quantile(batch_rewards, solver.percentile);
        reward_valid_indices = batch_rewards_states .>= reward_threshold;

        filtered_states  = batch_states[:, reward_valid_indices];
        filtered_actions = batch_actions[reward_valid_indices];

        if length(filtered_actions) > 0
            dataprovider = mx.ArrayDataProvider(:data => filtered_states, :label => filtered_actions)
            mx.fit(policy.net, optimizer, dataprovider, n_epoch = 2, verbosity = 0)
        end

        println("threshold: $(round(reward_threshold, 2)), mean score: $(round(mean(batch_rewards), 2))")

    end

    return policy

end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
module LeadLagAssistant

# Function to call the GUI
export LLAGUI

# Include the main file .fl
include("mainGUI.jl")
end # module


text:
export Diag

# Diag can have either Vector storage, in which case
# it is a general Diag tensor, or scalar storage,
# in which case the diagonal has a uniform value
struct Diag{T} <: TensorStorage
  data::T
  Diag{T}(data) where {T} = new{T}(data)
end

const NonuniformDiag{T} = Diag{T} where {T<:AbstractVector}
const UniformDiag{T} = Diag{T} where {T<:Number}

Base.@propagate_inbounds Base.getindex(D::NonuniformDiag,i::Int)= data(D)[i]
Base.getindex(D::UniformDiag,i::Int) = data(D)

Base.@propagate_inbounds Base.setindex!(D::Diag,val,i::Int)= (data(D)[i] = val)
Base.setindex!(D::UniformDiag,val,i::Int)= error("Cannot set elements of a uniform Diag storage")

Base.fill!(D::Diag,v) = fill!(data(D),v)

# convert to complex
# TODO: this could be a generic TensorStorage function
Base.complex(D::Diag{T}) where {T} = Diag{complex(T)}(complex(data(D)))

Base.copy(D::Diag{T}) where {T} = Diag{T}(copy(data(D)))

Base.eltype(::Diag{T}) where {T} = eltype(T)
Base.eltype(::Type{<:Diag{T}}) where {T} = eltype(T)

# Deal with uniform Diag conversion
Base.convert(::Type{<:Diag{T}},D::Diag) where {T} = Diag{T}(data(D))

Base.similar(D::Diag{T}) where {T} = Diag{T}(similar(data(D)))

# TODO: write in terms of ::Int, not inds
Base.similar(D::Diag{T},inds) where {T} = Diag{T}(similar(data(D),minimum(dims(inds))))
Base.similar(D::Type{<:NonuniformDiag{T}},inds) where {T} = Diag{T}(similar(T,length(inds)==0 ? 1 : minimum(dims(inds))))

Base.similar(D::UniformDiag{T}) where {T<:Number} = Diag{T}(zero(T))
Base.similar(::Type{<:UniformDiag{T}},inds) where {T<:Number} = Diag{T}(zero(T))

# TODO: make this work for other storage besides Vector
Base.zeros(::Type{Diag{T}},dim::Int64) where
  {T<:AbstractVector} = Diag{T}(zeros(eltype(T),dim))
Base.zeros(::Type{Diag{T}},dim::Int64) where
  {T<:Number} = Diag{T}(zero(T))

function Base.:*(D::Diag{<:Vector{El}},x::S) where {El<:Number,S<:Number}
  return Diag{Vector{promote_type(El,S)}}(x*data(D))
end

function Base.:*(D::UniformDiag{<:El},x::S) where {El<:Number,S<:Number}
  return Diag{promote_type(El,S)}(x*data(D))
end

Base.:*(x::Number,D::Diag) = D*x

#
# Type promotions involving Diag
# Useful for knowing how conversions should work when adding and contracting
#

Base.promote_rule(::Type{<:UniformDiag{T1}},::Type{<:UniformDiag{T2}}) where
  {T1<:Number,T2<:Number} = Diag{promote_type(T1,T2)}

Base.promote_rule(::Type{<:NonuniformDiag{T1}},::Type{<:NonuniformDiag{T2}}) where
  {T1<:AbstractVector,T2<:AbstractVector} = Diag{promote_type(T1,T2)}

# TODO: how do we make this work more generally for T2<:AbstractVector{S2}?
# Make a similar_type(AbstractVector{S2},T1) -> AbstractVector{T1} function?
Base.promote_rule(::Type{<:UniformDiag{T1}},::Type{<:NonuniformDiag{T2}}) where
  {T1<:Number,T2<:Vector{S2}} where S2 = Diag{Vector{promote_type(T1,S2)}}

Base.promote_rule(::Type{<:Dense{T1}},::Type{<:Diag{T2}}) where
  {T1,T2} = Dense{promote_type(T1,eltype(T2))}

# Convert a Diag storage type to the closest Dense storage type
dense(::Type{<:Diag{T}}) where {T} = Dense{eltype(T)}

const DiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Diag}
const NonuniformDiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where
                                               {StoreT<:NonuniformDiag}
const UniformDiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where
                                               {StoreT<:UniformDiag}

Base.IndexStyle(::Type{<:DiagTensor}) = IndexCartesian()

# TODO: this needs to be better (promote element type, check order compatibility,
# etc.
function Base.convert(::Type{<:Tensor{ElT,N,<:Dense}}, T::DiagTensor{ElT,N}) where {ElT,N}
  return dense(T)
end

# These are rules for determining the output of a pairwise contraction of Tensors
# (given the indices of the output tensors)
function contraction_output_type(TensorT1::Type{<:DiagTensor},
                                 TensorT2::Type{<:DenseTensor},
                                 indsR)
  return similar_type(promote_type(TensorT1,TensorT2),indsR)
end
contraction_output_type(TensorT1::Type{<:DenseTensor},
                        TensorT2::Type{<:DiagTensor},
                        indsR) = contraction_output_type(TensorT2,TensorT1,indsR)

# This performs the logic that DiagTensor*DiagTensor -> DiagTensor if it is not an outer
# product but -> DenseTensor if it is
# TODO: if the tensors are both order 2 (or less), or if there is an Index replacement,
# then they remain diagonal. Should we limit DiagTensor*DiagTensor to cases that
# result in a DiagTensor, for efficiency and type stability? What about a general
# SparseTensor result?
function contraction_output_type(TensorT1::Type{<:DiagTensor{<:Number,N1}},
                                 TensorT2::Type{<:DiagTensor{<:Number,N2}},
                                 indsR) where {N1,N2}
  if ValLength(indsR)===Val(N1+N2)
    # Turn into is_outer(inds1,inds2,indsR) function?
    # How does type inference work with arithmatic of compile time values?
    return similar_type(dense(promote_type(TensorT1,TensorT2)),indsR)
  end
  return similar_type(promote_type(TensorT1,TensorT2),indsR)
end

# TODO: move to tensor.jl?
zero_contraction_output(TensorT1::Type{<:Tensor},TensorT2::Type{<:Tensor},indsR) = zeros(contraction_output_type(TensorT1,TensorT2,indsR),indsR)

function contraction_output(TensorT1::Type{<:DiagTensor},
                            TensorT2::Type{<:Tensor},
                            indsR)
  return zero_contraction_output(TensorT1,TensorT2,indsR)
end
contraction_output(TensorT1::Type{<:Tensor},TensorT2::Type{<:DiagTensor},indsR) = contraction_output(TensorT2,TensorT1,indsR)

function contraction_output(TensorT1::Type{<:DiagTensor},
                            TensorT2::Type{<:DiagTensor},
                            indsR)
  return zero_contraction_output(TensorT1,TensorT2,indsR)
end

function array(T::DiagTensor{ElT,N}) where {ElT,N}
  return array(dense(T))
end
matrix(T::DiagTensor{<:Number,2}) = array(T)
vector(T::DiagTensor{<:Number,1}) = array(T)

diag_length(T::DiagTensor) = minimum(dims(T))
diag_length(T::DiagTensor{<:Number,0}) = 1

getdiag(T::DiagTensor,ind::Int) = store(T)[ind]

setdiag!(T::DiagTensor,val,ind::Int) = (store(T)[ind] = val)

setdiag(T::DiagTensor,val::ElR,ind::Int) where {ElR} = Tensor(Diag{ElR}(val),inds(T))

Base.@propagate_inbounds function Base.getindex(T::DiagTensor{ElT,N},inds::Vararg{Int,N}) where {ElT,N}
  if all(==(inds[1]),inds)
    return store(T)[inds[1]]
  else
    return zero(eltype(ElT))
  end
end
Base.@propagate_inbounds Base.getindex(T::DiagTensor{<:Number,1},ind::Int) = store(T)[ind]
Base.@propagate_inbounds Base.getindex(T::DiagTensor{<:Number,0}) = store(T)[1]

# Set diagonal elements
# Throw error for off-diagonal
Base.@propagate_inbounds function Base.setindex!(T::DiagTensor{<:Number,N},
                                            val,inds::Vararg{Int,N}) where {N}
  all(==(inds[1]),inds) || error("Cannot set off-diagonal element of Diag storage")
  return store(T)[inds[1]] = val
end
Base.@propagate_inbounds Base.setindex!(T::DiagTensor{<:Number,1},val,ind::Int) = ( store(T)[ind] = val )
Base.@propagate_inbounds Base.setindex!(T::DiagTensor{<:Number,0},val) = ( store(T)[1] = val )

function Base.setindex!(T::UniformDiagTensor{<:Number,N},val,inds::Vararg{Int,N}) where {N}
  error("Cannot set elements of a uniform Diag storage")
end

# TODO: make a fill!! that works for uniform and non-uniform
Base.fill!(T::DiagTensor,v) = fill!(store(T),v)

function dense(::Type{<:Tensor{ElT,N,StoreT,IndsT}}) where {ElT,N,
                                                            StoreT<:Diag,IndsT}
  return Tensor{ElT,N,dense(StoreT),IndsT}
end

# convert to Dense
function dense(T::TensorT) where {TensorT<:DiagTensor}
  R = zeros(dense(TensorT),inds(T))
  for i = 1:diag_length(T)
    setdiag!(R,getdiag(T,i),i)
  end
  return R
end

function outer!(R::DenseTensor{<:Number,NR},
                T1::DiagTensor{<:Number,N1},
                T2::DiagTensor{<:Number,N2}) where {NR,N1,N2}
  for i1 = 1:diag_length(T1), i2 = 1:diag_length(T2)
    indsR = CartesianIndex{NR}(ntuple(r -> r ≤ N1 ? i1 : i2, Val(NR)))
    R[indsR] = getdiag(T1,i1)*getdiag(T2,i2)
  end
  return R
end

# TODO: write an optimized version of this?
function outer!(R::DenseTensor{ElR},
                T1::DenseTensor,
                T2::DiagTensor) where {ElR}
  R .= zero(ElR)
  outer!(R,T1,dense(T2))
  return R
end

function outer!(R::DenseTensor{ElR},
                T1::DiagTensor,
                T2::DenseTensor) where {ElR}
  R .= zero(ElR)
  outer!(R,dense(T1),T2)
  return R
end

# Right an in-place version
function outer(T1::DiagTensor{ElT1,N1},
               T2::DiagTensor{ElT2,N2}) where {ElT1,ElT2,N1,N2}
  indsR = unioninds(inds(T1),inds(T2))
  R = Tensor(Dense{promote_type(ElT1,ElT2)}(zeros(dim(indsR))),indsR)
  outer!(R,T1,T2)
  return R
end

function Base.permutedims!(R::DiagTensor{<:Number,N},
                           T::DiagTensor{<:Number,N},
                           perm::NTuple{N,Int},f::Function=(r,t)->t) where {N}
  # TODO: check that inds(R)==permute(inds(T),perm)?
  for i=1:diag_length(R)
    @inbounds setdiag!(R,f(getdiag(R,i),getdiag(T,i)),i)
  end
  return R
end

function Base.permutedims(T::DiagTensor{<:Number,N},
                          perm::NTuple{N,Int},f::Function=identity) where {N}
  R = similar(T,permute(inds(T),perm))
  permutedims!(R,T,perm,f)
  return R
end

function Base.permutedims(T::UniformDiagTensor{ElT,N},
                          perm::NTuple{N,Int},
                          f::Function=identity) where {ElR,ElT,N}
  R = Tensor(Diag{promote_type(ElR,ElT)}(f(getdiag(T,1))),
                                         permute(inds(T),perm))
  return R
end

# Version that may overwrite in-place or may return the result
function permutedims!!(R::NonuniformDiagTensor{<:Number,N},
                       T::NonuniformDiagTensor{<:Number,N},
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->t) where {N}
  permutedims!(R,T,perm,f)
  return R
end

function permutedims!!(R::UniformDiagTensor{ElR,N},
                       T::UniformDiagTensor{ElT,N},
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->t) where {ElR,ElT,N}
  R = Tensor(Diag{promote_type(ElR,ElT)}(f(getdiag(R,1),getdiag(T,1))),inds(R))
  return R
end

function Base.permutedims!(R::DenseTensor{ElR,N},
                           T::DiagTensor{ElT,N},
                           perm::NTuple{N,Int},
                           f::Function = (r,t)->t) where {ElR,ElT,N}
  for i = 1:diag_length(T)
    @inbounds setdiag!(R,f(getdiag(R,i),getdiag(T,i)),i)
  end
  return R
end

function permutedims!!(R::DenseTensor{ElR,N},
                       T::DiagTensor{ElT,N},
                       perm::NTuple{N,Int},f::Function=(r,t)->t) where {ElR,ElT,N}
  permutedims!(R,T,perm,f)
  return R
end

function _contract!!(R::UniformDiagTensor{ElR,NR},labelsR,
                     T1::UniformDiagTensor{<:Number,N1},labelsT1,
                     T2::UniformDiagTensor{<:Number,N2},labelsT2) where {ElR,NR,N1,N2}
  if NR==0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    R = setdiag(R,diag_length(T1)*getdiag(T1,1)*getdiag(T2,1),1)
  else
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    R = setdiag(R,getdiag(T1,1)*getdiag(T2,1),1)
  end
  return R
end

function _contract!(R::DiagTensor{ElR,NR},labelsR,
                    T1::DiagTensor{<:Number,N1},labelsT1,
                    T2::DiagTensor{<:Number,N2},labelsT2) where {ElR,NR,N1,N2}
  if NR==0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    Rdiag = zero(ElR)
    for i = 1:diag_length(T1)
      Rdiag += getdiag(T1,i)*getdiag(T2,i)
    end
    setdiag!(R,Rdiag,1)
  else
    min_dim = min(diag_length(T1),diag_length(T2))
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    for i = 1:min_dim
      setdiag!(R,getdiag(T1,i)*getdiag(T2,i),i)
    end
  end
  return R
end

function _contract!(C::DenseTensor{ElC,NC},Clabels,
                    A::DiagTensor{ElA,NA},Alabels,
                    B::DenseTensor{ElB,NB},Blabels) where {ElA,NA,
                                                           ElB,NB,
                                                           ElC,NC}
  if all(i -> i < 0, Blabels)
    # If all of B is contracted
    # TODO: can also check NC+NB==NA
    min_dim = minimum(dims(B))
    if length(Clabels) == 0
      # all indices are summed over, just add the product of the diagonal
      # elements of A and B
      for i = 1:min_dim
        setdiag!(C,getdiag(C,1)+getdiag(A,i)*getdiag(B,i),1)
      end
    else
      # not all indices are summed over, set the diagonals of the result
      # to the product of the diagonals of A and B
      # TODO: should we make this return a Diag storage?
      for i = 1:min_dim
        setdiag!(C,getdiag(A,i)*getdiag(B,i),i)
      end
    end
  else
    astarts = zeros(Int,length(Alabels))
    bstart = 0
    cstart = 0
    b_cstride = 0
    nbu = 0
    for ib = 1:length(Blabels)
      ia = findfirst(==(Blabels[ib]),Alabels)
      if !isnothing(ia)
        b_cstride += stride(B,ib)
        bstart += astarts[ia]*stride(B,ib)
      else
        nbu += 1
      end
    end

    c_cstride = 0
    for ic = 1:length(Clabels)
      ia = findfirst(==(Clabels[ic]),Alabels)
      if !isnothing(ia)
        c_cstride += stride(C,ic)
        cstart += astarts[ia]*stride(C,ic)
      end
    end

    # strides of the uncontracted dimensions of
    # B
    bustride = zeros(Int,nbu)
    custride = zeros(Int,nbu)
    # size of the uncontracted dimensions of
    # B, to be used in CartesianIndices
    busize = zeros(Int,nbu)
    n = 1
    for ib = 1:length(Blabels)
      if Blabels[ib] > 0
        bustride[n] = stride(B,ib)
        busize[n] = size(B,ib)
        ic = findfirst(==(Blabels[ib]),Clabels)
        custride[n] = stride(C,ic)
        n += 1
      end
    end

    boffset_orig = 1-sum(strides(B))
    coffset_orig = 1-sum(strides(C))
    cartesian_inds = CartesianIndices(Tuple(busize))
    for inds in cartesian_inds
      boffset = boffset_orig
      coffset = coffset_orig
      for i in 1:nbu
        ii = inds[i]
        boffset += ii*bustride[i]
        coffset += ii*custride[i]
      end
      for j in 1:diag_length(A)
        C[cstart+j*c_cstride+coffset] += getdiag(A,j)*
                                         B[bstart+j*b_cstride+boffset]
      end
    end
  end
end
_contract!(C::DenseTensor,Clabels,
           A::DenseTensor,Alabels,
           B::DiagTensor,Blabels) = _contract!(C,Clabels,
                                               B,Blabels,
                                               A,Alabels)


text:
function create_network(
	nettopology::String,
	agentcount::Int64,
	networkprops::Dict
)
	if nettopology == "grid"
		network = LightGraphs.grid(agentcount, networkprops)
	elseif nettopology == "erdos_renyi"
		network = LightGraphs.erdos_renyi(agentcount, networkprops["p"])
	elseif nettopology == "watts_strogatz"
		network = LightGraphs.watts_strogatz(agentcount, networkprops["k"], networkprops["beta"])
	elseif nettopology == "barabasi_albert"
		network = LightGraphs.barabasi_albert(agentcount, networkprops["m0"])
	elseif nettopology == "complete"
		network = LightGraphs.complete_graph(agentcount)
	else
		network = default_network(agentcount, networkprops)
	end
	return network
end

function LightGraphs.grid(agentcount::Int64, networkprops::Dict)
	try
		height = Int(networkprops["grid_height"])
		width = Int(agentcount / height)
		network = LightGraphs.grid([height, width])
		return network
	catch
		network = LightGraphs.grid(Int64[agentcount, 1])
		return network
	end
end

function default_network(agentcount::Int64, networkprops::Dict)
	network = LightGraphs.grid(agentcount, networkprops)
	return network
end

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
@testset "phylogenetic_tree(::String)" begin
    tree = phylogenetic_tree("(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5);")
    names = [get_prop(tree.graph, x, :name) for x in 1:6]
    branch_lengths = [get_prop(tree.graph, e, :branch_length) for
                        e in collect(edges(bfs_tree(tree.graph, 1)))[1:end]]
    @test names == ["root", "A", "B", "E", "C", "D"]
    @test branch_lengths == [0.1, 0.2, 0.5, 0.3, 0.4]
end

@testset "phylogenetic_tree(::AbstractHLAData)" begin
    fasta_path = joinpath(@__DIR__, "data", "test.fasta")
    hla_types = rand(HLAType, 5)
    hla_data = HLAData("test", fasta_path, hla_types, missing, missing)

    @test Escape.phylogenetic_tree(hla_data) isa PhylogeneticTree
end

@testset "add_to_graph!(::AbstractMetaGraph, ::Int, ::String)" begin
    graph = MetaDiGraph()
    add_vertex!(graph)

    Escape.add_to_graph!(graph, 1, "(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)")
    names = [get_prop(graph, x, :name) for x in [3,4,5,6,7]]
    branch_lengths = [get_prop(graph, e, :branch_length) for
                        e in collect(edges(graph))[2:end]]

    @test names == ["A", "B", "E", "C", "D"]
    @test branch_lengths == [0.1, 0.2, 0.5, 0.3, 0.4]

end

@testset "extract_name(::String)" begin
    @test ismissing(Escape.extract_name(""))
    @test ismissing(Escape.extract_name("(A:3, B:1)"))
    @test ismissing(Escape.extract_name("(A:3, B:1):3"))
    @test ismissing(Escape.extract_name("(,)"))
    @test Escape.extract_name("B:5") == "B"
    @test Escape.extract_name("(A:3, B:1)E") == "E"
    @test Escape.extract_name("(A:3, B:1)E:5") == "E"
    @test Escape.extract_name("(,)E") == "E"
    @test Escape.extract_name("(,B)E:3") == "E"
end

@testset "extract_length(::String)" begin
    @test ismissing(Escape.extract_length(""))
    @test ismissing(Escape.extract_length("(A:3, B:1)"))
    @test ismissing(Escape.extract_length("(,)"))
    @test ismissing(Escape.extract_length("(A:3, B:1)E"))
    @test ismissing(Escape.extract_length("(,)E"))
    @test Escape.extract_length("(A:3, B:1):3") == 3.0
    @test Escape.extract_length("B:5") == 5.0
    @test Escape.extract_length("(A:3, B:1)E:5") == 5.0
    @test Escape.extract_length("(,B)E:3") == 3.0
end

@testset "newick_component(::String)" begin
    components = ["A:0.1", "B:0.2", "(C:0.3,D:0.4)E:0.5"]
    @test Escape.newick_components("(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5);") == components
    @test Escape.newick_components("(,)") == ["", ""]
    @test Escape.newick_components("(,B)E:3") == ["", "B"]
end

@testset "newick_string(::PhylogeneticTree)" begin
    newick = ["(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5);", "(,);",
              "((B:0.2,(C:0.3,D:0.4)E:0.5)A:0.1);",
              "(B:6.0,(A:5.0,C:3.0,E:4.0):5.0,D:11.0);"]
    for n in newick
        @test Escape.newick_string(phylogenetic_tree(n)) == n
    end
end

@testset "leaves(::PhylogeneticTree)" begin
    fasta_path = joinpath(@__DIR__, "data", "phylogeny.fasta")
    hla_data = HLAData("test", fasta_path, rand(HLAType, 6), missing, missing)
    tree = phylogenetic_tree(hla_data)

    @test leaves(tree) isa Vector{Int}
end

@testset "isleaf(v::Int, tree::Phylogenetic)" begin
    fasta_path = joinpath(@__DIR__, "data", "phylogeny.fasta")
    hla_data = HLAData("test", fasta_path, rand(HLAType, 6), missing, missing)
    tree = phylogenetic_tree(hla_data)

    @test isleaf(leaves(tree)[1], tree)
    @test !isleaf(90, tree)
end

@testset "get_property(::PhylogeneticTree, ::Int, ::Symbol)" begin
    fasta_path = joinpath(@__DIR__, "data", "phylogeny.fasta")
    hla_data = HLAData("test", fasta_path, rand(HLAType, 6), missing, missing)
    tree = phylogenetic_tree(hla_data)
    tree_leaves = Escape.leaves(tree)
    @test get_property(tree, 1, :name) == "root"
    @test get_property(tree, tree_leaves[1], :name) ==
        get_prop(tree.graph, tree_leaves[1], :name)
    @test get_property(tree, tree_leaves[1], :name) ==
        get_prop(tree.graph, tree_leaves[1], :name)
    @test ismissing(get_property(tree, 100, :name))
    @test ismissing(get_property(tree, 1, :nonexisting))
end

@testset "set_property!(::PhylogeneticTree, ::Int, ::Symbol)" begin
    fasta_path = joinpath(@__DIR__, "data", "phylogeny.fasta")
    hla_data = HLAData("test", fasta_path, rand(HLAType, 6), missing, missing)
    tree = phylogenetic_tree(hla_data)

    @test get_property(tree, 1, :name) == "root"
    set_property!(tree, 1, :name, "new name")
    @test get_property(tree, 1, :name) == "new name"
end

@testset "annotate!(::PhylogeneticTree, ::HLAData, ::Replacement)" begin
    fasta_path = joinpath(@__DIR__, "data", "test.fasta")
    r = Replacement("test", 2, 'S', false)
    hla_data = HLAData("test", fasta_path, rand(HLAType, 5), missing, missing)
    tree = phylogenetic_tree(hla_data)

    annotation = Dict("1" => "0", "2" => "0", "3" => "1", "4" => "1", "5" => "1")
    Escape.annotate!(tree, hla_data, r)

    for k in keys(annotation)
        v = filter(x -> get_property(tree, x, :name) == k, leaves(tree))[1]
        @test get_property(tree, v, :state) == annotation[k]
    end
end

@testset "matching(::PhylogeneticTree, ::HLAData)" begin
    fasta_path = joinpath(@__DIR__, "data", "phylogeny.fasta")
    hla_data = HLAData("test", fasta_path, rand(HLAType, 6), missing, missing)
    tree = phylogenetic_tree(hla_data)

    @test Escape.matching(tree, hla_data)

    set_property!(tree, leaves(tree)[2], :name, "test")
    @test Escape.matching(tree, hla_data) isa ErrorException

    short_hla_data = HLAData("short", joinpath(@__DIR__, "data", "test.fasta"),
                             rand(HLAType, 5), missing, missing)
    @test Escape.matching(tree, short_hla_data) isa DimensionMismatch
end

@testset "load PhylogeneticTree from .jls" begin
    tree = phylogenetic_tree("(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5);")
    path, io = mktemp()

    serialize(path * ".jls", tree)
    loaded_tree = deserialize(path * ".jls")
    close(io)
    rm(path)

    @test loaded_tree isa PhylogeneticTree
end