text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

const FACTOR_LABEL_MAP = OrderedDict(
    1e8 => "hundreds of millions",
    1e7 => "tens of millions",
    1e6 => "millions",
    1e5 => "hundreds of thousands",
    1e4 => "tens of thousands",
    1e3 => "thousands",
    1e2 => "hundreds",
    1e1 => "tens",
    1.0 => "",
)

function _get_yscale(max_yvalue)
    for (yscale, ylabel) ∈ FACTOR_LABEL_MAP
        if max_yvalue / yscale >= 2.0
            return (yscale, ylabel)
        end
    end
    (1.0, "")
end

function _get_plot_title(df)
    city = "city" ∈ names(df) ? df[1, :city] : missing
    state = "state" ∈ names(df) ? df[1, :state] : missing
    country = "country" ∈ names(df) ? df[1, :country] : missing
    location = skipmissing([city, state, country])
    join(location, ", ")
end

plot(args...; kwargs...) = Plots.plot(args...; kwargs...)

function plot(
    df::DataFrame;
    keycolumn = :date,
    columns = intersect(string.(option(:plot_columns)), names(df)),
    plot_begin = 1,
    plot_end = nrow(df),
    labels = nothing,
    title = _get_plot_title(df),
    date_format = option(:plot_date_format),
    new_window::Bool = true,
    ylabel = "People",
    yscale = nothing,
    legend = false,
    left_margin = 0mm,
    seriestypes = :path,
    seriescolors = :auto,
    kwargs...,
)
    if plot_begin isa Integer
        if plot_begin <= 0
            ini = max(1, nrow(df) + plot_begin)
        else
            ini = min(plot_begin, nrow(df))
        end
    elseif plot_begin isa Period
        ini = Dates.days(plot_begin) + 1
    else
        column = plot_begin isa Date ? df.date : df[!, keycolumn]
        ini = findfirst(isequal(plot_begin), column)
        if isnothing(ini)
            throw(ArgumentError("plot_begin not found in DataFrame: $(plot_begin)"))
        end
    end

    if plot_end isa Integer
        if plot_end <= 0
            iend = max(1, nrow(df) + plot_end)
        else
            iend = min(plot_end, nrow(df))
        end
    elseif plot_end isa Period
        iend = nrow(df) - Dates.days(plot_end)
    else
        column = plot_end isa Date ? df.date : df[!, keycolumn]
        iend = findfirst(isequal(plot_end), column)
        if isnothing(ini)
            throw(ArgumentError("plot_end not found in DataFrame: $(plot_end)"))
        end
    end
    @argcheck ini <= iend

    if !(isa(columns, AbstractVector) || isa(columns, Tuple))
        columns = [columns]
        #fillcolors = [fillcolors]
    elseif (legend == false) && (length(columns) > 1)
        legend = :topleft
    end

    if !(isa(seriestypes, AbstractVector) || isa(seriestypes, Tuple))
        seriestypes = fill(seriestypes, length(columns))
    end
    if !(isa(seriescolors, AbstractVector) || isa(seriescolors, Tuple))
        seriescolors = fill(seriescolors, length(columns))
    end

    df = df[ini:iend, :]
    columns = string.(columns)
    max_yvalue = 0
    for cname ∈ columns
        max_yvalue = max(maximum(skipmissing(df[!, cname])), max_yvalue)
    end
    numpeople = if hasproperty(df, :estimated_population)
        df.estimated_population[1]
    else
        2 * max_yvalue
    end
    if isnothing(yscale)
        (yscale, scalestr) = _get_yscale(max_yvalue)
        !isempty(scalestr) && (ylabel *= " ($(scalestr))")
    end

    left_margin isa Real && (left_margin *= mm)

    if keycolumn == :date
        X = Dates.format.(df.date, date_format)
    else
        X = string.(df[!, keycolumn])
        isnothing(labels) && seriestype == :pie && (labels = [X])
    end
    if isnothing(labels)
        labels = string.(prettify.(Symbol.(columns)))
    end
    Y = df[!, columns[1]] ./ yscale

    @debug("Plotting series", column = columns[1], X = Tuple(X), Y = Tuple(Y))

    win = Plots.plot(
        X,
        Y;
        xrotation = 90,
        xticks = min(30, length(X)),
        label = labels[1],
        ylabel = ylabel,
        title = title,
        legend = legend,
        left_margin = left_margin,
        seriestype = seriestypes[1],
        seriescolors = seriescolors[1],
        kwargs...,
    )
    #for (yn, label) ∈ zip(columns[2:end], labels[2:end])
    for (column, label, seriestype, seriescolor) ∈
        zip(columns[2:end], labels[2:end], seriestypes[2:end], seriescolors[2:end])

        @debug("Plotting series", column = column, X = Tuple(X), Y = Tuple(Y))
        Y = df[!, column] ./ yscale
        Plots.plot!(
            win,
            X,
            Y;
            label = label,
            seriestype = seriestype,
            seriescolor = seriescolor,
        )
    end
    new_window && Plots.gui(win)
    win
end

function model_combined_data(
    model::AbstractEndemicModel;
    columns = option(:plot_columns),
    plot_period = Day(14),
    kwargs...,
)
    if !(isa(columns, AbstractVector) || isa(columns, Tuple))
        columns = [columns]
    end
    columns = string.(columns)

    realdf = realdata(model)
    model_initial_date = default_kwarg(model, :initial_date, today() - Day(15))
    plot_initial_date = realdf.date[end] - plot_period
    ndays = Dates.days(model_initial_date - plot_initial_date)

    @debug "Combining data" model_initial_date plot_initial_date ndays

    realdf = realdf[(end-Dates.days(plot_period)):end, :]
    if ndays > 1
        model = model_step(model, -ndays; initial_date = plot_initial_date)
    end
    modeldf = @where(
        modeldata(model),
        (:date .<= realdf.date[end]) .& (:date .>= plot_initial_date)
    )

    @debug "Getting columns" columns Tuple(names(modeldf)) Tuple(names(realdf))
    columns = intersect(columns, names(modeldf), names(realdf))
    @debug "Got" columns
    rcolumns = [Symbol(sym, " (real)") for sym ∈ columns]
    mcolumns = [Symbol(sym, " (model)") for sym ∈ columns]

    realdf = select(realdf, :date, (columns .=> rcolumns)...)
    modeldf = select(modeldf, :date, (columns .=> mcolumns)...)

    @debug "Model dataframes computed" _debuginfo(realdf) _debuginfo(modeldf)
    (leftjoin(realdf, modeldf; on = :date), [rcolumns; mcolumns])
end

function plot(
    model::AbstractEndemicModel,
    kind = :combined_data;
    columns = option(:plot_columns),
    new_window::Bool = true,
    title = nothing,
    kwargs...,
)
    kind = Symbol(kind)
    if kind == :combined_data
        (df, columns) = model_combined_data(model; columns = columns, kwargs...)
    elseif kind == :realdata
        df = realdata(model)
    elseif kind == :modeldata
        df = modeldata(model)
    else
        throw(ArgumentError("Unrecognized plot kind: $(kind)"))
    end
    isnothing(title) && (title = _get_plot_title(realdata(model)))

    plot(df; title = title, new_window = new_window, columns = columns, kwargs...)
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the Multivector type.
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# External packages
import DataStructures.SDMKeyIteration, DataStructures.SortedDict

# GeometricAlgebra.jl
using GeometricAlgebra

#=
# --- Constructor tests

@testset "Multivector: inner constructor" begin
    # Notes
    # -----
    # * Test value of constructed instance

    # --- Preparations

    vectors = [3 3; -4 -4; 0 1]
    one_vector = [3; 4; 0]

    test_value = rand() + 1  # add 1 to avoid 0
    test_value = rand() > 0.5 ? test_value : -test_value

    test_dim = length(one_vector)

    # --- Multivector{T}(blades::Vector{AbstractBlade{T}};
    #                    reduced::Bool=false) where {T<:AbstractFloat}

    for precision_type in subtypes(AbstractFloat)
        # Preparations
        scalar = Scalar{precision_type}(test_value)
        one_blade = Blade{precision_type}(one_vector)
        two_blade = Blade{precision_type}(vectors)
        pseudoscalar = Pseudoscalar{precision_type}(test_dim, test_value)

        blades = Vector([scalar, one_blade, two_blade, pseudoscalar])

        # default value for `reduced`
        M = Multivector{precision_type}(blades)
        @test length(M.parts) == length(blades)
        for B in blades
            @test haskey(M.parts, grade(B))
            @test length(M.parts[grade(B)]) == 1
        end
    end
end

@testset "Multivector: outer constructor" begin
    # Notes
    # -----
    # * Test type of constructed instances. Correct construction of instances
    #   is tested by the inner constructor tests.

    # --- Preparations

    vectors = [3 3; -4 -4; 0 1]
    one_vector = [3; 4; 0]

    test_value = rand() + 1  # add 1 to avoid 0
    test_value = rand() > 0.5 ? test_value : -test_value

    test_dim = length(one_vector)

    # --- Multivector{T}(blades::Vector{AbstractBlade{T}};
    #                    reduced::Bool=false) where {T<:AbstractFloat}

    for precision_type in subtypes(AbstractFloat)
        # Preparations
        scalar = Scalar{precision_type}(test_value)
        two_blade = Blade{precision_type}(vectors)
        pseudoscalar = Pseudoscalar{precision_type}(test_dim, test_value)

        blades = Vector([scalar, two_blade, pseudoscalar])

        # default value for `reduced`
        M = Multivector(blades)
        @test M isa Multivector{precision_type}
    end
end

# --- Test attribute methods

@testset "Multivector: AbstractMultivector attribute functions" begin
    # --- Preparations

    vectors = [3 3; -4 -4; 0 1]
    one_vector = [3; 4; 0]

    test_value = rand() + 1  # add 1 to avoid 0
    test_value = rand() > 0.5 ? test_value : -test_value

    test_dim = length(one_vector)

    # --- Tests

    for precision_type in subtypes(AbstractFloat)
        # Preparations
        scalar = Scalar{precision_type}(test_value)
        one_blade = Blade{precision_type}(one_vector)
        pseudoscalar = Pseudoscalar{precision_type}(test_dim, test_value)
        expected_blades = Vector([scalar, one_blade, pseudoscalar])
        M = Multivector{precision_type}(expected_blades)

        expected_parts = SortedDict(0=>Vector([scalar]),
                                    1=>Vector([one_blade]),
                                    test_dim=>Vector([pseudoscalar]))
        expected_blades = reduce(vcat, values(expected_parts))

        # dim()
        @test dim(M) == test_dim

        # grades()
        @test grades(M) == [0, 1, 3]
        @test grades(M, collect=false) isa SDMKeyIteration

        # blades()
        @test blades(M) isa Vector{AbstractBlade}
        @test blades(M) == expected_blades

        # norm()
        @test_skip norm(M) == 0

        # getindex()
        for k in 1:test_dim
            if k in grades(M)
                k_vectors = M[k]
                @test k_vectors isa Vector{AbstractBlade}
                @test k_vectors == expected_parts[k]
            end
        end
    end
end

@testset "Multivector: convert(M)" begin
    # Preparations
    vectors = [3 3; -4 -4; 0 1]
    one_vector = [3; 4; 0]

    test_value = rand() + 1  # add 1 to avoid 0
    test_value = rand() > 0.5 ? test_value : -test_value

    test_dim = length(one_vector)

    # Tests
    for precision_type_converted in subtypes(AbstractFloat)
        for precision_type_src in subtypes(AbstractFloat)
            # Preparations
            scalar = Scalar{precision_type_src}(test_value)
            one_blade = Blade{precision_type_src}(one_vector)
            pseudoscalar = Pseudoscalar{precision_type_src}(test_dim,
                                                            test_value)
            blades = Vector([scalar, one_blade, pseudoscalar])
            M = Multivector{precision_type_src}(blades)

            # Exercise functionality and check results
            M_converted = convert(Multivector{precision_type_converted}, M)

            @test M_converted isa Multivector{precision_type_converted}

            if precision_type_src == precision_type_converted
                @test M_converted === M
            else
                @test M_converted !== M
                @test_skip M_converted ≈ M
            end
        end
    end
end
=#

ext: jl


text:
# create a Dict
d1 = Dict("name" => "Richard", "age" => 21)
println(typeof(d1))
println(d1)
println(d1["name"])

# to explicitly state the types of 'keys' and 'values'
d2 = Dict{String,String}("name" => "Dawn", "city" => "New York")
println(typeof(d2))
println(d2)
println(d2["name"])

# and to add an entry ...
d2["employee"] = "Google"
println(d2)

# to check if a Dict has a 'key'
println(haskey(d2, "employee"))
println(haskey(d2, "status"))

# to remove an entry from a Dict
delete!(d2, "employee")
println(d2)

# to retrieve the keys / values of a Dict as arrays
println(keys(d2))
println(values(d2))

# it is possible to iterate over a Dict
for (k, v) in d2
    println("key : $k / value : $v")
end

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
"""
	SeisToSegy(filename_in,filename_out;<keyword arguments>)

Convert seis data to SU or SEGY format. The function needs input and output filenames.

# Arguments
- `su=true` : If the flag equals true, converts tu SU format, otherwise to SEGY format

*Credits: AS, 2015*

"""
function SeisToSegy(in,out;su=true)
    if (su==true)
	file_hsize = 0
    else
	file_hsize = 900
	# add commands here to read text and binary headers out.thead and
        # out.bhead and write them to out
    end

    filename_headers = ParseHeaderName(in)
    extent = ReadTextHeader(in)
    n1 = extent.n1
    nx = extent.n2*extent.n3*extent.n4*extent.n5

    stream = open(out,"w")
    total = 60 + n1
    println("ok")
    h_segy = Array{SeisMain.SegyHeader}(undef,1)
    h_seis = Array{Header}(undef,1)
    h1 = Header[]
    push!(h1,SeisMain.InitSeisHeader())
    h1[1].o1 = extent.o1
    h1[1].d1 = extent.d1
    h1[1].n1 = extent.n1

    for j = 1 : nx
	if filename_headers != "NULL"
	    d,h1,e = SeisRead(in,itrace=j,ntrace=1,group="some")
	else
	    println("j=",j)
	    d,e = SeisRead(in,itrace=j,ntrace=1,group="some")
	    println("size(d)=",size(d))
	    h1[1].tracenum = j
    end
    println("test")
	#h_segy[1] = MapHeaders(h1,j,"SeisToSegy")
    h_segy = MapHeaders(h1,j,"SeisToSegy")

    positn = file_hsize + total*(j-1)*4 + segy_count["trace"]
	seek(stream,positn)
	write(stream,convert(Array{Float32,1},d[:]))
	PutSegyHeader(stream,h_segy,n1,file_hsize,j)
    end
    close(stream)

end


text:
function BruteForceOpt_MaxFreq(Proc_DF::ProcessedDF,FreqList,Inputs)
# inputs should be an Nx1 array of pairs where first item in pair is the name of the variable, and the second is a 1xM vector of possible values

#     VarLengths = zeros(Int8,1,length(Inputs))
#     VarIndex = zeros(Int8,size(VarLengths))
#     for i in 1:length(Inputs)
#         VarLengths[i]=length(Inputs[i][2])
#         VarIndex[i]=findfirst(x -> x==Inputs[i][1],convert(Array{String,1},Proc_DF.DF.Name) )

#     end
#     #FinalData = [[zeros(size(VarLengths))] [zeros(length(Proc_DF.InputList)-1,length(FreqList))]]

#     for ii in 1:VarLengths[1], j in 1:VarLengths[2],k in 1:VarLengths[3]

#         Proc_DF.DF.Value[VarIndex[1]]=Inputs[1][2][ii]
#         Proc_DF.DF.Value[VarIndex[2]]=Inputs[2][2][j]
#         Proc_DF.DF.Value[VarIndex[3]]=Inputs[3][2][k]
#         TrialRes = RunDF(Proc_DF,FreqList)
#         if (ii==1)&(j==1)&(k==1)
#             FinalData=TrialRes
#             println(FinalData)
#         else
#             FinalData = [FinalData;TrialRes]
#         end
#         #FinalData = vcat(FinalData,TrialRes[3,:])

#     end
#     #plot(log10.(abs.(TrialRes[3,:])))

# return FinalData
# end

# function RunDF(DF::ProcessedDF,FreqList)

#     CM = SPICE_DF2Matrix(DF)
#     CondMat(w) =
#         CM.RMat .+ CM.invLMat ./ (im .* w) .+ im .* w .* CM.CMat .+ CM.SrcMat
#     inputs = [0; 0; 0; 0; 0; 1]

#     Results = zeros(length(inputs), length(FreqList))
#     for i = 1:length(FreqList)
#         Freq = FreqList[i]
#         Results[:, i] = abs.(inv(CondMat(2 * pi * Freq)) * inputs)
#     end
#     return Results
end


text:
include("../utils.jl")

@testset "Torus" begin
    M = Torus(2)
    @testset "Torus Basics" begin
        @test Circle()^3 === Torus(3)
        @test Circle()^(3,) === PowerManifold(Circle(), 3)
        @test ^(Circle(), 2) === Torus(2)
        @test typeof(^(Circle(), 2)) == Torus{2}
        @test repr(M) == "Torus(2)"
        @test representation_size(M) == (2,)
        @test manifold_dimension(M) == 2
        @test !is_point(M, 9.0)
        @test_throws DomainError is_point(M, 9.0, true)
        @test !is_point(M, [9.0; 9.0])
        @test_throws CompositeManifoldError is_point(M, [9.0 9.0], true)
        @test !is_vector(M, [9.0; 9.0], 0.0)
        @test_throws DomainError is_vector(M, 9.0, 0.0, true) # point false and checked
        @test !is_vector(M, [9.0; 9.0], [0.0; 0.0])
        @test_throws DomainError is_vector(M, [0.0, 0.0], 0.0, true)
        @test injectivity_radius(M) ≈ π
        x = [1.0, 2.0]
        y = [-1.0, 2.0]
        z = [0.0, 0.0]
        basis_types = (DefaultOrthonormalBasis(),)
        test_manifold(
            M,
            [x, y, z],
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=true,
            test_project_tangent=false,
            test_musical_isomorphisms=true,
            test_default_vector_transport=false,
            basis_types_to_from=basis_types,
            is_tangent_atol_multiplier=1,
            test_inplace=true,
        )
    end
end


text:
"""
    module ArgumentModes
\n The module provides type `Mode` which could be seen as an extension of `Val`.
`Mode` is indended to be used as a type for a function method argument. A
specialization of `Mode` type contains a list of accepted symbols (flags). The
dispatch would choose the method only if the argument value (a `Mode` instance)
containes only symbols from the list of accepted symbols as declared in the type
of the argument.

Presumed uses of `Mode`:
1) Replacement of using ordinary `Symbol` as a function argument as a flag
   parameter. Here `Mode` allows to explicitly declare a function argument as a
   set of symbols (flags) with a distinct list of accepted symbols for each of
   a function's methods. The dispatch process then chooses a method with
   respect to those lists. That way it also indirectly performs the typo
   control. For example, `open(f, Mode(:read))`, `open(f, Mode(:write))`,
   `open(f, Mode(:write, :sync))` might correspond (depending on the design)
   to 3 different methods.

2) The way to explicitly show in which meaning a value to a function argument
   is provided. This might be useful when it is not possible to distinct the
   meaning of an argument only by the type.

   For example, suppose a function **f** processes array-typed objects with
   arbitrary dimensions number. We might want to declare methods for both
   processing a single object and an iteratable collection of objects. It
   would be difficult to distinct the methods using only the type of the
   argument since `f([x,y])` could both mean to process a single object `[x,y]`
   or to process two objects `x` and `y`. However, using `Mode` in the
   declaration of method for a collection, the user would be allowed to
   explicitly indicate what is passed in the call: `f([x,y])` for a single
   object and `f(Mode(:collection)=> [x, y])` for a collection of objects.

Series of tests showed that the current implementation of `Mode` fully compiles
out when it used for function arguments and method dispatch, so it seems that
there is no runtime overhead for using it (at least for the use cases
considered in the tests).

# Types
- `Mode`: similar to a set of symbols (possibly, with parameters) with ability
  for dispatch process to choose the function method based on a value of a `Mode` instance.

# Type constructors
- `Mode[ s₁ [=> t₁] [, s₂ [=> t₂]]... ]`: declare a type specialization which
  would accept the symbols `s₁`, `s₂`, …. Intended to be used in a function
  method declaration as a type of an argument.

# Instance constructors
- `Mode( s₁[=> v₁] [, v₂ [=> v₂]]... )`: create an instance with symbols `s₁`,
  `s₂`, … and corresponding values `v₁`, `v₂`, … in it. Intended to be used
  in a function call as an argument value.

# Operations on a `Mode` instance
- `m => v` given `m::Mode{[s => Nothing]}`: return a new instance with single
  symbol in it.
- `m₁ ~ m₂`: join `m₁` and `m₂`.
- `keys(m), values(m), pairs(m)`: return symbols / values / pairs of symbols
  and values.
- `m[s₁ [, s₂]... ]`: return the value of `s₁` / the tuple of
  values of `s₁, s₂, ...`.
- `m[]`: if `m` contains only one value, return it; throw an ArgumentError
  otherwise.
- `checkmode(m, ...)`: check if `m` is a `Mode` and contains the prescribed
  symbols, and, possibly, do an action.
"""
module ArgumentModes
  if isdefined(Base, :Experimental) &&
      isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1 # 0 disables inlining in user-code
  end

  include("Impl.jl")
  using .Impl
  Mode.body.body.name.module = ArgumentModes
  export checkmode
  export Mode

  # doesn't seem to help
  # include("precompile.jl")
end # module


text:
module ApplicationBuilder

using Compat
Compat.@warn """ApplicationBuilder has changed since JuliaCon 2018.
  - `module BuildApp` has been removed. You should remove `using BuildApp`
      from your build scripts, and just use `using ApplicationBuilder`.
  - You should also _remove_ `using ApplicationBuilder` from the source-code of
      programs being built, since the `change_dir_if_bundle` behavior will
      now come default for all applications. `using ApplicationBuilder` in
      your program source may cause it to fail to build.
 This message will be removed in the next version of ApplicationBuilder.
"""

using Glob, PackageCompiler

export build_app_bundle

@static if Compat.Sys.isapple()

include("sign_mac_app.jl")
include("mac_commandline_app.jl")

"""
    build_app_bundle(juliaprog_main;
        appname, builddir, resources, libraries, verbose, bundle_identifier,
        app_version, icns_file, certificate, entitlements_file, commandline_app)

Compile `juliaprog_main` into an executable, and bundle it together with all its
`resources` and `libraries` into an App called `appname`.

juliaprog_main: Path to a ".jl" file that defines the function `julia_main()`

# Examples
```julia-repl
julia> build_app_bundle("main.jl", appname="MyApp", resources=["img.jpg"],
           libraries=[MyPackage._libp])
```
"""
function build_app_bundle(juliaprog_main;
        appname=splitext(basename(juliaprog_main))[1], builddir = "builddir",
        resources = String[], libraries = String[], verbose = false,
        bundle_identifier = nothing, app_version = "0.1", icns_file = nothing,
        certificate = nothing, entitlements_file = nothing,
        snoopfile = nothing, autosnoop = false, commandline_app = false,
    )

    builddir = abspath(builddir)
    if bundle_identifier == nothing
        bundle_identifier = make_bundle_identifier(appname)
        println("  Using calculated bundle_identifier: '$bundle_identifier'")
    end
    binary_name = splitext(basename(juliaprog_main))[1]

    # ----------- Input sanity checking --------------

    if !isfile(juliaprog_main) throw(ArgumentError("Cannot build application. No such file '$juliaprog_main'")) end
    # Bundle identifier requirements: https://apple.stackexchange.com/a/238381/52530
    if occursin(r"\s", bundle_identifier) throw(ArgumentError("Bundle identifier must not contain whitespace.")) end
    if occursin(r"[^A-Za-z0-9-.]", bundle_identifier) throw(ArgumentError("Bundle identifier must contain only alphanumeric characters (A-Z,a-z,0-9), hyphen (-), and period (.).")) end

    # ----------- Initialize App ---------------------
    println("~~~~~~ Creating mac app in \"$builddir/$appname.app\" ~~~~~~~")

    appDir="$builddir/$appname.app/Contents"

    launcherDir="$appDir/MacOS"
    resourcesDir="$appDir/Resources"
    libsDir="$appDir/Libraries"

    applet_name = nothing
    if commandline_app
        # TODO: What if the user specifies Resources that could overwrite
        #  applet resources? (ie Scripts/ or applet.rsrc)
        applet_name = build_commandline_app_bundle(builddir, binary_name, appname, verbose)
    end

    mkpath(launcherDir)

    function has_files(files)
        return !isempty(files) && !all(isempty, [strip(f) for f in files])
    end

    mkpath(resourcesDir)  # Always create resources directory (For .icns file)
    has_files(libraries) && mkpath(libsDir)

    # ----------- Copy user libs & assets -------------
    run_verbose(verbose, cmd) = (verbose && println("    julia> run($cmd)") ; run(cmd))

    if has_files(libraries) || has_files(resources)
        # Copy assets and libs early so they're available during compilation.
        #  This is mostly relevant if you have .dylibs in your assets, which the compiler needs to look at.

        println("~~~~~~ Copying user-specified libraries & resources to bundle... ~~~~~~~")

        function copy_file_dir_or_glob(pattern, dest)
            if isfile(pattern)
                run_verbose(verbose, `cp -f $pattern $dest/`) # Copy the file to dest
            elseif isdir(pattern)
              run_verbose(verbose, `cp -rf $pattern $dest/`) # Copy the entire *dir* (not its contents) to dest.
            elseif pattern[1] == '/'
                files = glob(pattern[2:end], pattern[1:1])
                run_verbose(verbose, `cp -rf $files $dest/`) # Copy the specified glob pattern to dest.
            elseif !isempty(glob(pattern))
                run_verbose(verbose, `cp -rf $(glob(pattern)) $dest/`) # Copy the specified glob pattern to dest.
            else
                Compat.@warn "Skipping unknown file '$pattern'!"
            end
        end

        function clean_file_pattern(pattern, errorMsg_fileCmd)
            clean_pattern = strip(pattern)
            if isempty(clean_pattern)
                throw(ArgumentError("ERROR: $errorMsg_fileCmd '$pattern' must not be empty."))
            elseif clean_pattern[1] == '~'
                clean_pattern = homedir() * pattern[2:end]
            end
            return clean_pattern
        end
        println("  Resources:")
        for res in resources
            res = clean_file_pattern(res, "-R")
            print("    - $res ...")
            copy_file_dir_or_glob(res, resourcesDir)
            println("............ done")
        end
        println("  Libraries:")
        for lib in libraries
            lib = clean_file_pattern(lib, "-L")
            print("    - $lib ...")
            copy_file_dir_or_glob(lib, libsDir)
            println("............ done")
        end
    end

    # ----------- Compile a binary ---------------------
    # Compile the binary right into the app.
    println("~~~~~~ Compiling a binary from '$juliaprog_main'... ~~~~~~~")

    if autosnoop
        if snoopfile != nothing
            println("WARNING: autosnoop is overwriting user-specified snoopfile.")
        end
        snoopfile = "$builddir/$appname-snoopfile.jl"
        write(snoopfile, """include("$(abspath("$juliaprog_main"))");  julia_main([""]); """)
    end

    # When Apple launches an app, it sets the current-working-directory to homedir.
    # Therefore, we inject this function definition into the app, and call it from
    # the C driver program.
    utils_injection_file = "$launcherDir/applicationbuilderutils.jl"
    write(utils_injection_file,
        """
            include("$(abspath(juliaprog_main))")
        """*raw"""
            using Compat

            Base.@ccallable function cd_to_bundle_resources()::Compat.Nothing
                full_binary_name = PROGRAM_FILE  # PROGRAM_FILE is set manually in program.c
                if Compat.Sys.isapple()
                    m = match(r".app/Contents/MacOS/[^/]+$", full_binary_name)
                    if m != nothing
                        resources_dir = joinpath(dirname(dirname(full_binary_name)), "Resources")
                        cd(resources_dir)
                    end
                    println("cd_to_bundle_resources(): Changed to new pwd: $(pwd())")
                end
            end
            precompile(cd_to_bundle_resources, ())  # Compile it for the binary.
        """)
    custom_program_c = "$(@__DIR__)/program.c"
    # Provide an environment variable telling the code it's being compiled into a mac bundle.
    withenv("LD_LIBRARY_PATH"=>"$libsDir:$libsDir/julia",
            "COMPILING_APPLE_BUNDLE"=>"true") do
        verbose && println("  PackageCompiler.static_julia(...)")
        # Compile executable and copy julia libs to $launcherDir.
        PackageCompiler.build_executable(utils_injection_file, binary_name, custom_program_c;
                builddir=launcherDir, verbose=verbose, optimize="3",
                snoopfile=snoopfile, debug="0", cpu_target="x86-64",
                cc_flags=`-mmacosx-version-min=10.10 -headerpad_max_install_names`)
    end

    for b in ["$launcherDir/$binary_name", "$launcherDir/$binary_name.dylib"]
        run_verbose(verbose, `install_name_tool -add_rpath "@executable_path/../Frameworks/" $b`)
        run_verbose(verbose, `install_name_tool -add_rpath "@executable_path/../Libraries/" $b`)
    end

    # In order to pass Apple's GateKeeper, the App must not reference any external Julia libs.
    for binary_file in glob("*", launcherDir)
        try
            #  an example output line from otool: "         path /Applications/Dev Apps/Julia-0.6.app/Contents/Resources/julia/lib (offset 12)"
            external_julia_deps = readlines(pipeline(`otool -l $binary_file`,
                 `grep $(dirname(Compat.Sys.BINDIR))`,  # filter julia lib deps
                 `sed 's/\s*path//'`, # remove leading "  path"
                 `sed 's/(.*)$//'`)) # remove trailing parens
            for line in external_julia_deps
                path = strip(line)
                run_verbose(verbose, `install_name_tool -delete_rpath "$path" $binary_file`)
            end
        catch
        end
        # Also need to strip any non-x86 architectures to make Apple happy.
        # (It looks like this only affects libgcc_s.1.dylib.)
        try
            if success(pipeline(`file $binary_file`, `grep 'i386'`))
                run_verbose(verbose, `lipo $binary_file -thin x86_64 -output $binary_file`)
            end
        catch
        end
    end

    # ---------- Create Info.plist to tell it where to find stuff! ---------
    # This lets you have a different app name from your juliaprog_main.
    println("~~~~~~ Generating 'Info.plist' for '$bundle_identifier'... ~~~~~~~")

    info_plist() = """
        <?xml version="1.0" encoding="UTF-8"?>
        <!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
        <plist version="1.0">
        <dict>
        	<key>CFBundleAllowMixedLocalizations</key>
        	<true/>
        	<key>CFBundleDevelopmentRegion</key>
        	<string>en</string>
        	<key>CFBundleDisplayName</key>
        	<string>$appname</string>
        	<key>CFBundleExecutable</key>
        	<string>$(commandline_app ? applet_name : binary_name)</string>
        	<key>CFBundleIconFile</key>
        	<string>$appname.icns</string>
        	<key>CFBundleIdentifier</key>
        	<string>$bundle_identifier</string>
            <key>CFBundleInfoDictionaryVersion</key>
            <string>6.0</string>
        	<key>CFBundleName</key>
        	<string>$appname</string>
        	<key>CFBundlePackageType</key>
        	<string>APPL</string>
        	<key>CFBundleShortVersionString</key>
        	<string>$app_version</string>
        	<key>CFBundleVersion</key>
        	<string>$app_version</string>
        	<key>NSHighResolutionCapable</key>
            <string>YES</string>
        	<key>LSMinimumSystemVersionByArchitecture</key>
        	<dict>
        		<key>x86_64</key>
        		<string>10.10</string>
        	</dict>
        	<key>LSRequiresCarbon</key>
        	<true/>
        	<key>NSHumanReadableCopyright</key>
        	<string>© 2018 $bundle_identifier </string>
        </dict>
        </plist>
        """

    verbose && println(info_plist())
    write("$appDir/Info.plist", info_plist());

    # Copy Julia icons
    julia_app_resources_dir() = joinpath(Compat.Sys.BINDIR, "..","..")
    if (icns_file == nothing)
        icns_file = joinpath(julia_app_resources_dir(),"julia.icns")
        verbose && println("Attempting to copy default icons from Julia.app: $icns_file")
    end
    if isfile(icns_file)
        Compat.cp(icns_file, "$resourcesDir/$appname.icns", force=true);
    else
        Compat.@warn "Skipping nonexistent icons file: '$icns_file'"
    end

    # --------------- CLEAN UP before distributing ---------------
    println("~~~~~~ Cleaning up temporary files... ~~~~~~~")

    # Delete the file we added to inject code.
    rm(utils_injection_file)

    # Delete the tmp build files
    function delete_if_present(file, path)
        files = glob(file, path)
        if !isempty(files) run(`rm -r $(files)`) end
    end
    delete_if_present("*.ji",launcherDir)
    delete_if_present("*.o",launcherDir)

    # Remove debug .dylib libraries and any precompiled .ji's
    delete_if_present("*.dSYM",libsDir)
    delete_if_present("*.dSYM","$libsDir/julia")
    delete_if_present("*.backup","$libsDir/julia")
    delete_if_present("*.ji","$libsDir/julia")
    delete_if_present("*.o","$libsDir/julia")

    if certificate != nothing
        println("~~~~~~ Signing the binary and all libraries ~~~~~~~")
        sign_application_libs(launcherDir, certificate)
        if entitlements_file != nothing
            set_entitlements("$launcherDir/$binary_name", certificate, entitlements_file)
        end
    end

    println("~~~~~~ Done building '$builddir/$appname.app'! ~~~~~~~")
    return 0
end

function make_bundle_identifier(appname)
    cleanregex = r"[^a-zA-Z0-9]"
    cleanuser = replace(lowercase(ENV["USER"]), cleanregex => "")
    cleanapp = replace(lowercase(appname), cleanregex => "")
    "com.$cleanuser.$cleanapp"
end

end  # isapple

@static if Compat.Sys.islinux() || Compat.Sys.iswindows()
	include("bundle.jl")
end

@static if Compat.Sys.iswindows()
    include("installer.jl")
end

end  # module


text:
using Test

using ARCHModels
using Random
using GLM
using DataFrames
using StableRNGs

rng = StableRNG(1)
T = 10^4;

@testset "lgamma" begin
    @test ARCHModels.lgamma(1.0f0) == 0.0f0
end

@testset "TGARCH" begin
    @test ARCHModels.nparams(TGARCH{1, 2, 3}) == 7
    @test ARCHModels.presample(TGARCH{1, 2, 3}) == 3
    Random.seed!(1)
    spec = TGARCH{1,1,1}([1., .05, .9, .01]);
    str = sprint(show, spec)
    if VERSION < v"1.5.5"
        @test startswith(str, "TGARCH{1,1,1} specification.\n\n─────────────────────────────────\n               ω    γ₁   β₁    α₁\n─────────────────────────────────\nParameters:  1.0  0.05  0.9  0.01\n─────────────────────────────────\n")
    else
        @test startswith(str, "TGARCH{1, 1, 1} specification.\n\n─────────────────────────────────\n               ω    γ₁   β₁    α₁\n─────────────────────────────────\nParameters:  1.0  0.05  0.9  0.01\n─────────────────────────────────\n")
    end
    am = simulate(spec, T);
    am = selectmodel(TGARCH, am.data; meanspec=NoIntercept(), show_trace=true, maxlags=2)
    @test all(isapprox.(coef(am), [0.9439667311150648
                                   0.04573706835008625
                                   0.9043902283152758
                                   0.012555948398277313], rtol=1e-4))
   #everything below is just pure GARCH, in fact
    Random.seed!(1)
    spec = GARCH{1, 1}([1., .9, .05])
    am0 = simulate(spec, T);
    am00 = deepcopy(am0)
    Random.seed!(1)
    am00.data .= 0.
    simulate!(am00)
    @test all(am00.data .== am0.data)
    Random.seed!(1)
    am00 = simulate(am0)
    @test all(am00.data .== am0.data)
    Random.seed!(1)
    am000 = simulate(am0, nobs(am0))
    @test all(am000.data .== am0.data)
    am = selectmodel(GARCH, am0.data; meanspec=NoIntercept(), show_trace=true)
    @test isfitted(am) == true
    #with unconditional as presample:
    #@test all(isapprox.(coef(am), [0.9086632896184081,
    #                               0.9055268468427705,
    #                               0.050367854809777915], rtol=1e-4))
    @test all(isapprox.(coef(am), [0.9086479266110243,
                                   0.905531642067773,
                                   0.050324600594884535], rtol=1e-4))
    #with unconditional as presample:
    #@test all(isapprox.(stderror(am), [0.14582381264705224,
    #                                   0.010354562480367474,
    #                                   0.005222817398477784], rtol=1e-4))
    @test all(isapprox.(stderror(am), [0.14578736059501485,
                                       0.010356284482676704,
                                       0.005228247833454602], rtol=1e-3))
    @test sum(volatilities(am0)) ≈ 44768.17421580251
    @test sum(abs, residuals(am0)) ≈ 8022.163087384836
    @test sum(abs, residuals(am0, standardized=false)) ≈ 35939.07066637026
    am2 = UnivariateARCHModel(spec, am0.data)
    @test isfitted(am2) == false
    io = IOBuffer()
    str = sprint(io -> show(io, am2))
    if VERSION < v"1.5.5"
        @test startswith(str, "\nTGARCH{0,1,1}")
    else
        @test startswith(str, "\nGARCH{1, 1}")
    end
    fit!(am2)
    @test isfitted(am2) == true
    io = IOBuffer()
    str = sprint(io -> show(io, am2))
    if VERSION < v"1.5.5"
        @test startswith(str, "\nTGARCH{0,1,1}")
    else
        @test startswith(str, "\nGARCH{1, 1}")
    end
    am3 = fit(am2)
    @test isfitted(am3) == true
    @test all(am2.spec.coefs .== am.spec.coefs)
    @test all(am3.spec.coefs .== am2.spec.coefs)
end
@testset "ARCH" begin
    Random.seed!(1);
    spec = ARCH{2}([1., .3, .4]);
    am = simulate(spec, T);
    @test selectmodel(ARCH, am.data).spec.coefs == fit(ARCH{2}, am.data).spec.coefs
    Random.seed!(1);
    spec = ARCH{0}([1.]);
    am = simulate(spec, T);
    fit!(am)
    @test all(isapprox.(coef(am), 1.013031276122647, rtol=1e-4))

end

@testset "EGARCH" begin
    Random.seed!(1)
    @test ARCHModels.nparams(EGARCH{1, 2, 3}) == 7
    @test ARCHModels.presample(EGARCH{1, 2, 3}) == 3
    am = simulate(EGARCH{1, 1, 1}([.1, 0., .9, .1]), T; meanspec=Intercept(3))
    am7 = selectmodel(EGARCH, am.data; maxlags=2, show_trace=true)
    #with unconditional as presample:
    #@test all(isapprox(coef(am7), [0.08502955535533116,
    #                               0.004709708474515596,
    #                               0.9164935566284109,
    #                               0.09325947325535855,
    #                               3.0137461089470308], rtol=1e-4))
    @test all(isapprox(coef(am7), [0.08504883253172882,
                                   0.0047015720582706125,
                                   0.9164488571272553,
                                   0.09323297680588628,
                                   3.013732273404755], rtol=1e-4))

    @test coefnames(EGARCH{2, 2, 2}) == ["ω", "γ₁", "γ₂", "β₁", "β₂", "α₁", "α₂"]
end
@testset "StatisticalModel" begin
    #not implemented: adjr2, deviance, mss, nulldeviance, r2, rss, weights
    Random.seed!(1);
    spec = GARCH{1, 1}([1., .9, .05])
    am = simulate(spec, T)
    fit!(am)
    @test loglikelihood(am) ==  ARCHModels.loglik!(Float64[],
                                                                Float64[],
                                                                Float64[],
                                                                Float64[],
                                                                typeof(spec),
                                                                StdNormal{Float64},
                                                                NoIntercept{Float64}(),
                                                                am.data,
                                                                spec.coefs
                                                                )
    @test nobs(am) == T
    @test dof(am) == 3
    @test coefnames(GARCH{1, 1}) == ["ω", "β₁", "α₁"]
    #with unconditional as presample:
    #@test aic(am) ≈ 58369.082969298106 rtol=1e-4
    #@test bic(am) ≈ 58390.713990414035 rtol=1e-4
    #@test aicc(am) ≈ 58369.085370258494 rtol=1e-4
    @test aic(am) ≈ 58369.08203915834 rtol=1e-4
    @test bic(am) ≈ 58390.71306027427 rtol=1e-4
    @test aicc(am) ≈ 58369.08444011873 rtol=1e-4

    @test all(coef(am) .== am.spec.coefs)
    #with unconditional as presample:
    #@test all(isapprox(confint(am), [0.6228537382166024 1.1944723245606814;
    #                                0.8852323068993577 0.9258214298904262;
    #                                0.040131313548448275 0.060604377407810675],
    #                   rtol=1e-4)
    #                   )
    @test all(isapprox(confint(am), [0.6229099504436408 1.1943859027784078;
                                     0.8852336974680757 0.9258295866674704;
                                     0.04007742313906393 0.06057177805070514],
                       rtol=1e-4)
                       )
    #with unconditional as presample:
    #@test all(isapprox(informationmatrix(am; expected=false), [0.15326216336912968 2.9536982257433135 2.618124940552642;
    #                                                           2.9536982257433135 58.956837321202826 53.74888605159925;
    #                                                           2.618124940552642 53.74888605159925 53.29656483617587],
    #                   rtol=1e-4)
    #                   )
    @test all(isapprox(informationmatrix(am; expected=false), [0.15267905577531846 2.9430054037767617 2.607992890001237;
                                                               2.9430054037767617 58.76705213893014 53.57462949757515;
                                                               2.607992890001237 53.57462949757517 53.14231838982629],
                       rtol=1e-4)
                       )
    @test_throws ErrorException informationmatrix(am)
    #with unconditional as presample:
    #@test all(isapprox(score(am), [-4.091261171623728e-6 3.524550271549742e-5 -6.989366926291041e-5], rtol=1e-4))
    @test all(isapprox(score(am), [0. 0. 0.], atol=1e-3))
    @test islinear(am::UnivariateARCHModel) == false
    @test predict(am) ≈ 4.361606730361275
    @test predict(am, :variance) ≈ 19.023613270332767
    @test predict(am, :return) == 0.0
    @test predict(am, :VaR) ≈ 10.146614544578197
    @test predict.(am, :variance, 1:3) == [predict(am, :variance, h) for h in 1:3]
    @test_throws Base.ErrorException predict.(am, :VaR, 1:3)
end

@testset "MeanSpecs" begin
    Random.seed!(1);
    spec = GARCH{1, 1}([1., .9, .05])
    am = simulate(spec, T; meanspec=Intercept(0.))
    fit!(am)
    #with unconditional as presample:
    #@test all(isapprox(coef(am), [0.910496430719689,
    #                               0.9054120402733519,
    #                               0.05039127076312942,
    #                               0.027705636765390795], rtol=1e-4))
    @test all(isapprox(coef(am), [0.9104793904880137,
                                  0.9054169985282575,
                                  0.05034724930058784,
                                  0.027707836268720806], rtol=1e-4))
    @test ARCHModels.coefnames(Intercept(0.)) == ["μ"]
    @test ARCHModels.nparams(Intercept) == 1
    @test ARCHModels.presample(Intercept(0.)) == 0
    @test ARCHModels.constraints(Intercept{Float64}, Float64) == (-Float64[Inf], Float64[Inf])
    @test typeof(NoIntercept()) == NoIntercept{Float64}
    @test ARCHModels.coefnames(NoIntercept()) == []
    @test ARCHModels.constraints(NoIntercept{Float64}, Float64) == (Float64[], Float64[])
    @test ARCHModels.nparams(NoIntercept) == 0
    @test ARCHModels.presample(NoIntercept()) == 0
    @test ARCHModels.uncond(NoIntercept()) == 0
    @test mean(zeros(5), zeros(5), zeros(5), zeros(5), NoIntercept(), zeros(5), 4) == 0.
    ms = ARMA{2, 2}([1., .5, .2, -.1, .3])
    @test ARCHModels.nparams(typeof(ms)) == length(ms.coefs)
    @test ARCHModels.presample(ms) == 2
    @test ARCHModels.coefnames(ms) == ["c", "φ₁", "φ₂", "θ₁", "θ₂"]
    Random.seed!(1)
    spec = GARCH{1, 1}([1., .9, .05])
    am = simulate(spec, T; meanspec=ms)
    fit!(am)
    @test all(isapprox(coef(am), [0.9063506916409171,
                                  0.905682443482137,
                                  0.05021834228447521,
                                  1.079454022288992,
                                  0.45097800554911116,
                                  0.2357782619617334,
                                  -0.05909354030019596,
                                  0.2878312346045116], rtol=1e-4))
    @test predict(am, :return) ≈ -1.4572460532296017 rtol = 1e-6
    am = selectmodel(ARCH, BG96;  meanspec=AR, maxlags=2);
    @test all(isapprox(coef(am), [0.11916340875306261,
                                  0.3156862868133263,
                                  0.18331803992622006,
                                  -0.00685700871019875,
                                  0.0358362785070197], rtol=1e-4))
    @test typeof(Regression([1 2; 3 4])) == Regression{2, Float64}
    @test typeof(Regression([1. 2.; 3. 4.])) == Regression{2, Float64}
    @test typeof(Regression{Float32}([1 2; 3 4])) == Regression{2, Float32}
    @test typeof(Regression([1 2; 3 4])) == Regression{2, Float64}
    @test typeof(Regression([1, 2], [1 2; 3 4.0f0])) ==  Regression{2, Float32}
    @test typeof(Regression([1, 2.], [1 2; 3 4.0f0])) ==  Regression{2, Float64}
    @test typeof(Regression([1], [1, 2, 3, 4.0f0])) ==  Regression{1, Float32}
    @test typeof(Regression([1, 2, 3, 4.0f0])) ==  Regression{1, Float32}
    @test ARCHModels.nparams(Regression{2, Float64}) == 2

    rng = StableRNG(1)
    beta = [1, 2]
    reg = Regression(beta, rand(rng, 2000, 2))
    u = randn(rng, 2000)*.1
    y = reg.X*reg.coefs+u
    @test ARCHModels.coefnames(reg) == ["β₀", "β₁"]
    @test ARCHModels.presample(reg) == 0
    @test ARCHModels.constraints(typeof(reg), Float64) == ([-Inf, -Inf], [Inf, Inf])
    @test all(isapprox(ARCHModels.startingvals(reg, y),
        [0.9923610899808355, 2.0036469645073303], rtol=1e-4))
    @test ARCHModels.uncond(reg) === 0.
    Random.seed!(1)
    am = simulate(GARCH{1, 1}([1., .9, .05]), 2000; meanspec=reg, warmup=0)
    fit!(am)
    @test_throws Base.ErrorException predict(am, :return)

    @test all(isapprox(coef(am), [1.522147151343697,
                                  0.869732130237804,
                                  0.060598935466617924,
                                  0.7992856525779547,
                                  2.0369823982928796], rtol=1e-4))
    Random.seed!(1)
    am = simulate(GARCH{1, 1}([1., .9, .05]), 1999; meanspec=reg, warmup=0)
    @test predict(am, :return) ≈ 2.1435094119979254
    data = DataFrame(X=ones(1974), Y=BG96)
    model = lm(@formula(Y ~ -1 + X), data)
    am = fit(GARCH{1, 1}, model)
    @test all(isapprox(coef(am), coef(fit(GARCH{1, 1}, BG96, meanspec=Intercept)), rtol=1e-4))
    @test coefnames(am)[end] == "X"
    @test all(isapprox(coef(am), coef(fit(GARCH{1, 1}, model.model)), rtol=1e-4))
end

@testset "VaR" begin
    am = fit(GARCH{1, 1}, BG96)
    @test sum(VaRs(am)) ≈ 2077.0976454790807
end
@testset "Errors" begin
    #with unconditional as presample:
    #@test_warn "Fisher" stderror(UnivariateARCHModel(GARCH{3, 0}([1., .1, .2, .3]), [.1, .2, .3, .4, .5, .6, .7]))
    #@test_warn "non-positive" stderror(UnivariateARCHModel(GARCH{3, 0}([1., .1, .2, .3]), -5*[.1, .2, .3, .4, .5, .6, .7]))

    # the following are temporarily disabled while we use FiniteDiff for Hessians:
    #@test_logs (:warn, "Fisher information is singular; vcov matrix is inaccurate.") stderror(UnivariateARCHModel(GARCH{1, 0}( [1.0, .1]), [0., 1.]))
    #@test_logs (:warn, "non-positive variance encountered; vcov matrix is inaccurate.") stderror(UnivariateARCHModel(GARCH{1, 0}( [1.0, .1]), [1., 1.]))
    e = @test_throws ARCHModels.NumParamError ARCHModels.loglik!(Float64[], Float64[], Float64[], Float64[], GARCH{1, 1}, StdNormal{Float64},
                                                     NoIntercept{Float64}(), zeros(T),
                                                     [0., 0., 0., 0.]
                                                     )
    str = sprint(showerror, e.value)
    @test startswith(str, "incorrect number of parameters")
    @test_throws ARCHModels.NumParamError GARCH{1, 1}([.1])
    e = @test_throws ErrorException predict(UnivariateARCHModel(GARCH{0, 0}([1.]), zeros(10)), :blah)
    str = sprint(showerror, e.value)
    @test startswith(str, "Prediction target blah unknown")
    @test_throws ARCHModels.NumParamError ARMA{1, 1}([1.])
    @test_throws ARCHModels.NumParamError Intercept([1., 2.])
    @test_throws ARCHModels.NumParamError NoIntercept([1.])
    @test_throws ARCHModels.NumParamError StdNormal([1.])
    @test_throws ARCHModels.NumParamError StdT([1., 2.])
    @test_throws ARCHModels.NumParamError StdGED([1., 2.])
    @test_throws ARCHModels.NumParamError Regression([1], [1 2; 3 4])
    at = zeros(10)
    data = rand(10)
    reg = Regression(data[1:5])
    @test_throws ErrorException mean(at, at, at, data, reg, [0.], 6)
end

@testset "Distributions" begin
    Random.seed!(1)
    a=rand(StdT(3))
    Random.seed!(1)
    b=rand(StdT(3), 1)[1]
    @test a==b # https://github.com/JuliaStats/Distributions.jl/issues/846
    Random.seed!(1)
    @test rand(MersenneTwister(), StdNormal()) ≈ 0.2972879845354616
    @testset "Gaussian" begin
        Random.seed!(1)
        data = rand(T)
        @test typeof(StdNormal())==typeof(StdNormal(Float64[]))
        @test fit(StdNormal, data).coefs == Float64[]
        @test coefnames(StdNormal) == String[]
        @test ARCHModels.distname(StdNormal) == "Gaussian"
        @test quantile(StdNormal(), .05) ≈ -1.6448536269514724
        @test ARCHModels.constraints(StdNormal{Float64}, Float64) == (Float64[], Float64[])
    end
    @testset "Student" begin
        Random.seed!(1)
        data = rand(StdT(4), T)
        spec = GARCH{1, 1}([1., .9, .05])
        @test fit(StdT, data).coefs[1] ≈ 3.972437329588246 rtol=1e-4
        @test coefnames(StdT) == ["ν"]
        @test ARCHModels.distname(StdT) == "Student's t"
        @test quantile(StdT(3), .05) ≈ -1.3587150125838563
        Random.seed!(1);
        datat = simulate(spec, T; dist=StdT(4)).data
        Random.seed!(1);
        datam = simulate(spec, T; dist=StdT(4), meanspec=Intercept(3)).data
        am4 = selectmodel(GARCH, datat; dist=StdT, meanspec=NoIntercept{Float64}(), show_trace=true)
        am5 = selectmodel(GARCH, datam; dist=StdT, show_trace=true)
        @test coefnames(am5) == ["ω", "β₁", "α₁", "ν", "μ"]
        @test all(coeftable(am4).cols[2] .== stderror(am4))
        #with unconditional as presample:
        #@test all(isapprox(coef(am4), [0.8307014299672306,
        #                               0.9189503152734588,
        #                               0.042080807758329355,
        #                               3.835646488238764], rtol=1e-4))
        @test all(isapprox(coef(am4), [0.831327902922751,
                                       0.9189082073187017,
                                       0.04211450420515991,
                                       3.834812845564172], rtol=1e-4))
        #with unconditional as presample:
        #@test all(isapprox(coef(am5), [0.8306175556436268,
        #                               0.9189538270625667,
        #                               0.04208964132482301,
        #                               3.8348509665880797,
        #                               2.9918445831618024], rtol=1e-4))
        @test all(isapprox(coef(am5), [0.8312482931785029,
                                       0.9189112458594789,
                                       0.04212329204579697,
                                       3.834033118707376,
                                       2.9918481871096083], rtol=1e-4))
    end
    @testset "GED" begin
        Random.seed!(1)
        @test typeof(StdGED(3)) == typeof(StdGED(3.)) == typeof(StdGED([3.]))
        data = rand(StdGED(1), T)
        @test fit(StdGED, data).coefs[1] ≈ 1.0193004687300224 rtol=1e-4
        @test coefnames(StdGED) == ["p"]
        @test ARCHModels.nparams(StdGED) == 1
        @test ARCHModels.distname(StdGED) == "GED"
        @test quantile(StdGED(1), .05) ≈ -1.6281735335151468
    end
    @testset "Standardized" begin
        using Distributions
        @test eltype(StdNormal{Float64}()) == Float64
        MyStdT=Standardized{TDist}
        @test typeof(MyStdT([1.])) == typeof(MyStdT(1.))
        @test ARCHModels.logconst(MyStdT, [0]) == 0.
        @test coefnames(MyStdT{Float64}) == ["ν"]
        @test ARCHModels.distname(MyStdT{Float64}) == "TDist"
        @test all(isapprox.(ARCHModels.startingvals(MyStdT, [0.]), eps()))
        @test quantile(MyStdT(3.), .1) ≈ quantile(StdT(3.), .1)
        ARCHModels.startingvals(::Type{<:MyStdT}, data::Vector{T}) where T = T[3.]
        Random.seed!(1)
        am = simulate(GARCH{1, 1}([1, 0.9, .05]), 1000, dist=MyStdT(3.))
        @test  loglikelihood(fit(am)) ≈ -2700.9089012063323
    end
end
@testset "tests" begin
    am = fit(GARCH{1, 1}, BG96)
    LM = ARCHLMTest(am)
    @test pvalue(LM) ≈ 0.1139758664282619
    str = sprint(show, LM)
    @test startswith(str, "ARCH LM test for conditional heteroskedasticity")
    @test ARCHModels.testname(LM) == "ARCH LM test for conditional heteroskedasticity"

    vars = VaRs(am, 0.01)
    DQ = DQTest(BG96, VaRs(am), 0.01)
    @test pvalue(DQ) ≈ 2.3891461144184955e-11
    str = sprint(show, DQ)
    @test startswith(str, "Engle and Manganelli's (2004) DQ test (out of sample)")
    @test ARCHModels.testname(DQ) == "Engle and Manganelli's (2004) DQ test (out of sample)"
end
@testset "multivariate" begin
    am1 = fit(DCC, DOW29[:, 1:2])
    am2 = fit(DCC, DOW29[:, 1:2]; method=:twostep)
    am3 = MultivariateARCHModel(DCC{1, 1}([1. 0.; 0. 1.], [0., 0.], [GARCH{1, 1}([1., 0., 0.]), GARCH{1, 1}([1., 0., 0.])]), DOW29[:, 1:2]) # not fitted
    am4 = fit(DCC, DOW29[1:20, 1:29]) # shrinkage n<p
    @test all(fit(am1).spec.coefs .== am1.spec.coefs)
    @test all(isapprox(am1.spec.coefs, [0.8912884521017908, 0.05515419379547665], rtol=1e-3))
    @test all(isapprox(am2.spec.coefs,    [0.8912161306136979, 0.055139392936998946], rtol=1e-3))
    @test all(isapprox(am4.spec.coefs, [0.8935938309400944, 6.938893903907228e-18], atol=1e-3))
    @test all(isapprox(stderror(am1)[1:2], [0.0434344187103969, 0.020778846682313102], rtol=1e-3))
    @test all(isapprox(stderror(am2)[1:2], [0.030405542205923865, 0.014782869078355866], rtol=1e-4))
    @test all(isapprox(predict(am1; what=:correlation)[:], [1.0, 0.4365129466277069, 0.4365129466277069, 1.0], rtol=1e-4))
    @test all(isapprox(predict(am1; what=:covariance)[:], [6.916591739333349, 1.329392154000225, 1.329392154000225,  1.340972349032465], rtol=1e-4))
    @test_throws ErrorException predict(am1; what=:bla)
    @test residuals(am1)[1, 1] ≈ 0.5107042609407892
    @test_throws ErrorException fit(DCC, DOW29; method=:bla)
    @test_throws ARCHModels.NumParamError DCC{1, 1}([1. 0.; 0. 1.], [1., 0., 0.], [GARCH{1, 1}([1., 0., 0.]), GARCH{1, 1}([1., 0., 0.])])
    @test_throws AssertionError DCC{1, 1}([1. 0.; 0. 1.], [0., 0.], [GARCH{1, 1}([1., 0., 0.]), GARCH{1, 1}([1., 0., 0.])]; method=:bla)
    @test coefnames(am1) == ["β₁", "α₁", "ω₁", "β₁₁", "α₁₁", "μ₁", "ω₂", "β₁₂", "α₁₂", "μ₂"]
    @test ARCHModels.nparams(DCC{1, 1}) == 2
    ARCHModels.nparams(DCC{1, 1, GARCH{1, 1}, Float64, 2}) == 8
    @test ARCHModels.presample(DCC{1, 2, GARCH{3, 4}}) == 4
    @test ARCHModels.presample(DCC{1, 2, GARCH{3, 4, Float64}, Float64, 2}) == 4
    io = IOBuffer()
    str = sprint(io -> show(io, am1))
    @test startswith(str, "\n2-dim")
    str = sprint(io -> show(io, am3))
    @test startswith(str, "\n2-dim")
    str = sprint(io -> show(io, am3.spec))
    @test startswith(str, "DCC{1, 1")
    str = sprint(io -> show(IOContext(io, :se=>true), am1))
    @test occursin("Std.Error", str)
    @test_throws ErrorException fit(DCC, DOW29[1:11, :]) # shrinkage requires n>=12
    @test loglikelihood(am1) ≈ -9810.905799585276

    @test ARCHModels.nparams(MultivariateStdNormal) == 0
    @test typeof(MultivariateStdNormal{Float64, 3}()) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal(Float64, 3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal(Float64[], 3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal{Float64}(3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    @test typeof(MultivariateStdNormal(3)) == typeof(MultivariateStdNormal{Float64, 3}(Float64[]))
    Random.seed!(1)
    @test all(isapprox(rand(MultivariateStdNormal(2)), [0.2972879845354616, 0.3823959677906078], rtol=1e-6))
    @test coefnames(MultivariateStdNormal) == String[]
    @test ARCHModels.distname(MultivariateStdNormal) == "Multivariate Normal"

    Random.seed!(1)
    am = am1
    am.spec.coefs .= [.7, .2]
    ams  = simulate(am)
    @test isfitted(ams) == false
    fit!(ams)
    @test isfitted(ams) == true
    @test all(isapprox(ams.spec.coefs, [0.7122483516102956, 0.19156028421431875], rtol=1e-4))
    Random.seed!(2)
    simulate!(ams)
    @test ams.fitted == false
    fit!(ams)
    @test all(isapprox(ams.spec.coefs, [0.6630049669013613, 0.22885770926598498], rtol=1e-4))
    Random.seed!(1)
    amc = fit(DCC{1, 2, GARCH{3, 2}}, DOW29[:, 1:4]; meanspec=AR{3})
    ams = simulate(amc, T)
    fit!(ams)
    @test all(isapprox(ams.meanspec[1].coefs, [-0.09394176323811071, 0.05159711352207107, 0.011348428433666473, -0.020175913191330077], rtol=1e-4))
    Random.seed!(1)
    ame = fit(DCC{1, 2, EGARCH{1, 1, 1}}, DOW29[:, 1:4])
    ams = simulate(ame, T)
    fit!(ams)
    @test all(isapprox(ams.spec.univariatespecs[1].coefs, [0.046700648921491957, -0.07258062140595305, 0.9664860227494515, 0.22051944930571496], rtol=1e-4))

    Random.seed!(1)
    ccc = fit(CCC, DOW29[:, 1:4])
    @test dof(ccc) == 16
    @test ccc.spec.R[1, 2] ≈ 0.37095654552885643
    @test isapprox(stderror(ccc)[1], 0.06298215515406534, rtol=1e-3)
    cccs = simulate(ccc, T)
    @test  cccs.data[end, 1] ≈ -1.5061782364569236
    @test coefnames(ccc) == ["ω₁", "β₁₁", "α₁₁", "μ₁", "ω₂", "β₁₂", "α₁₂", "μ₂", "ω₃", "β₁₃", "α₁₃", "μ₃", "ω₄", "β₁₄", "α₁₄", "μ₄"]
    io = IOBuffer()
    str = sprint(io -> show(io, ccc))
    @test startswith(str, "\n4-dim")
    io = IOBuffer()
    str = sprint(io -> show(io, ccc.spec))
    @test startswith(str, "DCC{0, 0")
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent
using Qaintessent.MaxKColSubgraphQAOA

##==----------------------------------------------------------------------------------------------------------------------

# adapted from https://github.com/FluxML/Zygote.jl/blob/master/test/gradcheck.jl
function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2-y1)/δ
        if eltype(x) <: Complex
            # derivative with respect to imaginary part
            x[i] = tmp - im*δ/2
            y1 = f(xs...)
            x[i] = tmp + im*δ/2
            y2 = f(xs...)
            x[i] = tmp
            Δ[i] += im*(y2-y1)/δ
        end
    end
    return grads
end

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "max-k-col. subgraph QAOA gate gradients" begin
    @testset "MaxKColGates gates" begin
        # fictitious gradients of cost function with respect to quantum gate
        Δ = randn(ComplexF64, 16, 16)
        κ = 2
        graph = Graph(2, [(1,2)])
        g  = MaxKColSubgraphPhaseSeparationGate
        f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(g(θ[], κ, graph))))
        θ = 2π*rand()
        ngrad = ngradient(f, [θ])
        dg = Qaintessent.backward(g(θ, κ, graph), conj(Δ))
        @test isapprox(dg.γ, ngrad[1], rtol=1e-4, atol=1e-4)
    end

    @testset "RNearbyValuesMixerGate gates" begin
        κs = 2:6
        rs = [1, 1, 2, 1, 4, 5]
        for (κ, r) ∈ zip(κs, rs)
            Δ = randn(ComplexF64, Base.:^(2, κ), Base.:^(2, κ))
            g = RNearbyValuesMixerGate
            f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(g(θ[], r, κ))))
            θ = 2π*rand()
            ngrad = ngradient(f, [θ])
            dg = Qaintessent.backward(g(θ, r, κ), conj(Δ))
            # sadly, with less tolerance it doesn't pass (could be improved by gradients that
            # are not just numerically estimated)
            @test isapprox(dg.β, ngrad[1], rtol=1e-4, atol=1e-4)
        end
    end

    @testset "ParityRingMixerGate gates" begin
        for κ in 2:6
            Δ = randn(ComplexF64, Base.:^(2, κ), Base.:^(2, κ))
            g = ParityRingMixerGate
            f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(g(θ[], κ))))
            θ = 2π*rand()
            ngrad = ngradient(f, [θ])
            dg = Qaintessent.backward(g(θ, κ), conj(Δ))
            # sadly, with less tolerance it doesn't pass (could be improved by gradients that
            # are not just numerically estimated)
            @test isapprox(dg.β, ngrad[1], rtol=1e-4, atol=1e-4)
        end
    end

    @testset "PartitionMixerGate gates" begin
        κ = 4
        Δ = randn(ComplexF64, Base.:^(2, κ), Base.:^(2, κ))
        partition = [[(1, 2), (3, 4)], [(2, 3), (4, 1)], [(1, 2), (3, 4)]]
        g = PartitionMixerGate
        f(θ) = 2*real(sum(Δ .* Qaintessent.sparse_matrix(g(θ[], κ, partition))))
        θ = 2π*rand()
        ngrad = ngradient(f, [θ])
        dg = Qaintessent.backward(g(θ, κ, partition), conj(Δ))
        @test isapprox(dg.β, ngrad[1], rtol=1e-4, atol=1e-4)
    end
end


text:
using BenchmarkTools
using FinancialMonteCarlo
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;
Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;
rfCurve = ZeroRate(r);
FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));
@btime EuPrice = pricer(Model, rfCurve, mc, EUData);


text:
function kmp(p,ATMp,Wip1)
    return p[C.kmp0]*(p[C.Dp0]+p[C.Dwp]*Wip1)/(p[C.Sp0]+p[C.Dp0]+p[C.Sap]*ATMp+p[C.Dwp]*Wip1)
end

function Kpm(p,ATMp,Wip1)
    return p[C.Kpm0]*(p[C.Sp0]+p[C.Dp0]+p[C.Sap]*ATMp+p[C.Dwp]*Wip1)/(p[C.Sp0]+p[C.Sap]*ATMp)
end

function Kpw(p,ATMp,Wip1)
    return p[C.Kpw0]*(p[C.Sp0]+p[C.Dp0]+p[C.Sap]*ATMp+p[C.Dwp]*Wip1)/(p[C.Sp0]+p[C.Sap]*ATMp)
end

function γm(p,ATMp,Wip1,γm0,γm1)
    return (γm0*(p[C.Dm0]+p[C.Dwm]*Wip1)+γm1*(p[C.Sm0]+p[C.Sam]*ATMp))/(p[C.Sm0]+p[C.Dm0]+p[C.Sam]*ATMp+p[C.Dwm]*Wip1)
end

function diffeq(du,u,h,p,t)

    du[V.p53] = p[C.kp0] - kmp(p,u[V.ATMp],u[V.Wip1])*u[V.p53]*u[V.Mdm2] - p[C.γp]*u[V.p53];

    du[V.Mdm2] = p[C.km0] + p[C.kpm]*(h(p,t-p[C.τm])[V.p53]^p[C.nm])/(Kpm(p,h(p,t-p[C.τm])[V.ATMp],h(p,t-p[C.τm])[V.Wip1])^p[C.nm]+h(p,t-p[C.τm])[V.p53]^p[C.nm]) - γm(p,u[V.ATMp],u[V.Wip1],p[C.γm0],p[C.γm1])*u[V.Mdm2];

    du[V.Wip1] = p[C.kw0] + p[C.kpw]*(h(p,t-p[C.τw])[V.p53]^p[C.nw])/(Kpw(p,h(p,t-p[C.τw])[V.ATMp],h(p,t-p[C.τw])[V.Wip1])^p[C.nw]+h(p,t-p[C.τw])[V.p53]^p[C.nw]) - p[C.γw]*u[V.Wip1];

    du[V.ATMp] = p[C.Sea]*p[C.Eto]*(p[C.ATMt] - u[V.ATMp]) - p[C.Da0]*u[V.ATMp];

end

function param_values_U2OS()
    p = zeros(C.NUM)

    p[C.kp0] = 1.3396
    p[C.km0] = 0.08
    p[C.kw0] = 0.105
    p[C.kpm] = 0.9172
    p[C.kpw] = 0.1496
    p[C.nm] = 4
    p[C.nw] = 4
    p[C.Kpm0] = 0.4025
    p[C.Kpw0] = 0.4025
    p[C.τm] = 2.1
    p[C.τw] = 2.1
    p[C.kmp0] = 4.038
    p[C.γp] = 0.1
    p[C.γm0] = 0.2579
    p[C.γm1] = 7.7362
    p[C.γw] = 0.4
    p[C.Sea] = 1
    p[C.Da0] = 80
    p[C.Sp0] = 0.25
    p[C.Sap] = 0.8335
    p[C.Dp0] = 0.75
    p[C.Dwp] = 0.25
    p[C.Sm0] = 0.0153
    p[C.Sam] = 0.283
    p[C.Dm0] = 0.5
    p[C.Dwm] = 0.5
    p[C.ATMt] = 16

    return p
end

function param_values_MCF7()
    p = zeros(C.NUM)

    p[C.kp0] = 2.4138
    p[C.km0] = 0.2
    p[C.kw0] = 0.35
    p[C.kpm] = 2.3854
    p[C.kpw] = 0.8702
    p[C.nm] = 4
    p[C.nw] = 4
    p[C.Kpm0] = 0.4025
    p[C.Kpw0] = 0.4025
    p[C.τm] = 2.1
    p[C.τw] = 2.1
    p[C.kmp0] = 2.7048
    p[C.γp] = 0.25
    p[C.γm0] = 0.1804
    p[C.γm1] = 3.6088
    p[C.γw] = 0.4
    p[C.Sea] = 1
    p[C.Da0] = 30
    p[C.Sp0] = 0.25
    p[C.Sap] = 0.4335
    p[C.Dp0] = 0.75
    p[C.Dwp] = 0.25
    p[C.Sm0] = 0.0479
    p[C.Sam] = 0.2830
    p[C.Dm0] = 0.5
    p[C.Dwm] = 0.5
    p[C.ATMt] = 3.5

    return p
end

function initialValues_U2OS()
    u0 = zeros(V.NUM)

    u0[V.p53] = 1
    u0[V.Mdm2] = 0.4
    u0[V.Wip1] = 0.3
    u0[V.ATMp] = 0

    return u0
end

function initialValues_MCF7()
    u0 = zeros(V.NUM)

    u0[V.p53] = 1
    u0[V.Mdm2] = 1
    u0[V.Wip1] = 1
    u0[V.ATMp] = 0

    return u0
end


text:
module Test_rmtools

using LinearAlgebra
using MatrixPencils
using Polynomials
using Test

@testset "Rational Matrix Tools" begin

@testset "Tests polgcdvw, pollcm, conv, poldivrem, poldiv" begin

# test example for Polynomials functions
num = Polynomial([0.8581454436924945, 0.249671302254737, 0.8048498901050951, 0.1922713965697087])
den = Polynomial([0.9261520696359462, 0.07141031902098072, 0.378071465860349])

@test poldeg(coeffs(num)) == 3 && poldeg1(coeffs(num)) == 4

q1, r1 = poldivrem(coeffs(num),coeffs(den))
@test coeffs(num) ≈ conv(coeffs(den),q1) + [r1;zeros(2) ]

d, = polgcdvw(coeffs(num),coeffs(den))
@test d == [1]

m = pollcm(coeffs(num),coeffs(den))
p = conv(coeffs(num),coeffs(den))
@test m*p[end] ≈ m[end]*p

# test example which fails for polgcd but not for polgcdvw
num = [ 0.34137237259273556
0.1815967953870914
0.005829782404846551
0.7934778818753858
0.5308749468599383
0.14488329929437715
0.488795838965558
0.806095325402598
0.29135959590683647
0.22853582789218474];
den = [ 0.7371989151075657
0.3525356298341029
0.7811847379591219
0.6173718776940151
0.32725779892403306
];
d = [ 0.015852680715187084
0.3602221067929683
0.7487570819767699
0.5882036947609526
0.011233093587261544
];

d2, = polgcdvw(conv(num,d),conv(den,d),atol=1.e-10);  # OK
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7

# Example Table 7: Zeng Math. Comp. 2005

num = reverse([1.00000000
23.35360257
29.89831582
10.75803809
15.57240922
18.76038493
13.73079603
30.45600101
46.21275197
44.89871211
30.17981700
8.33455813]);
den = reverse([1.00000000
23.01829201
22.05776405]);
sol=reverse([1.00000000
0.33531056
0.12227539
0.54726624
0.27815340
0.28629915
1.00523653
1.00205392
0.97391204
0.37785145 ]);

# this test fails
q, r = poldivrem(num,den)
@test !(q ≈ sol)

# long division has small residual
@test conv(q,den) ≈ num - [r;zeros(length(num)-length(r))] # OK

@test conv(sol,den)≈ num # OK

q = poldiv(num,den)  # accurate deconvolution for approximate quotient
@test q ≈ sol

# long division has small residual
q, r = poldivrem(conv(sol,den),den)
@test conv(q,den) ≈ conv(sol,den) - [r;zeros(length(num)-length(r))] # OK

q = poldiv(conv(sol,den),den)  # accurate deconvolution for exact quotient
@test q ≈ sol

# testing conv and poldiv for random data
num = rand(100);
d = rand(25); #d = Polynomial([100])
q = poldiv(conv(num,d),d);
@test q ≈ num

@test conv(num,d) ≈ convmtx(d,length(num))*num
@test conv(d,num) ≈ convmtx(d,length(num))*num

@test conv(d,num) ≈ convmtx(num,length(d))*d
@test conv(num,d) ≈ convmtx(num,length(d))*d

Ty = Complex{Float64}
num = rand(Ty,100);
d = rand(Ty,25); #d = Polynomial([100])
q = poldiv(conv(num,d),d)
@test q ≈ num

@test conv(num,d) ≈ convmtx(d,length(num))*num
@test conv(d,num) ≈ convmtx(d,length(num))*num

@test conv(d,num) ≈ convmtx(num,length(d))*d
@test conv(num,d) ≈ convmtx(num,length(d))*d

Ty = Complex{Float64}
num = rand(Ty,100);
d = rand(25); #d = Polynomial([100])
q = poldiv(conv(num,d),d)
@test q ≈ num

@test conv(num,d) ≈ convmtx(d,length(num))*num
@test conv(d,num) ≈ convmtx(d,length(num))*num

@test conv(d,num) ≈ convmtx(num,length(d))*d
@test conv(num,d) ≈ convmtx(num,length(d))*d

# tests polgcdvw, pollcm, conv, poldiv
p1 = [-6, 11, -6, 1]
p2 = [4, -5, 1]
d1, v1, w1, δ1 = polgcdvw(p1,p2)
@test poldiv(p1,d1) ≈ v1*norm(p1) && poldiv(p2,d1) ≈ w1*norm(p2)
m = pollcm(p1,p2)
mc = poldiv(conv(p1,p2),d1)
@test  mc*m[end]≈ m*mc[end]

p1 = [-6, 11, -6, 1]*rand(Complex{Float64})
p2 = [4, -5, 1]*rand(Complex{Float64})
d1, v1, w1, δ1 = polgcdvw(p1,p2)
@test poldiv(p1,d1) ≈ v1*norm(p1) && poldiv(p2,d1) ≈ w1*norm(p2)
m = pollcm(p1,p2)
mc = poldiv(conv(p1,p2),d1)
@test  mc*m[end]≈ m*mc[end]

# Example 2.1 Zeng 2011
p1 = [10, 1, 0, 0,  0, 0, 0, 0, 10/3, 31/3, 1];
p2 = [-60/7, -6/7, 0, 0,  0, 0, 0, 0, 10/7, 71/7, 1];
d = [10, 1];

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-5);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test p1/norm(p1) ≈ conv(v2,d2) && p2/norm(p2) ≈ conv(w2,d2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# Example 2.2 Zeng 2011 "exact"
p2 = [10, 1, 0, 0,  0, 0, 0, 0, 10/3, 31/3, 1]
p1 = [10, 1]
d = [10, 1]

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-5);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test p1/norm(p1) ≈ conv(v2,d2) && p2/norm(p2) ≈ conv(w2,d2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# Example 2.2 Zeng 2011 "approximate"
p1 = [10, 1, 0, 0,  0, 0, 0, 0,  3.333333333,  10.33333333, 1]
p2 = [10, 1]
d = [10, 1]

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-5);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test p1/norm(p1) ≈ conv(v2,d2) && p2/norm(p2) ≈ conv(w2,d2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# q, r = poldivrem(p1,p3) # inexact remainder of order 3!!

# Example of Zeng 2009, with corrected p1
p1 = [-.999999999999, 1, -.999999999999, 1, -.999999999999, 1, -.999999999999, 1];
p2 = [2,-2, 2, -3, 1, -1, 1];
d = [1, -1, 1, -1];   # approximate gcd

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-9);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test poldiv(p1,d2) ≈ v2*norm(p1) && poldiv(p2,d2) ≈ w2*norm(p2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# Example of Zeng 2009, with corrected p1
p1 = [-1, 1, -1, 1, -1, 1, -1, 1]
p2 = [2,-2, 2, -3, 1, -1, 1]
d = [1, -1, 1, -1];   # exact gcd

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-9);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test poldiv(p1,d2) ≈ v2*norm(p1) && poldiv(p2,d2) ≈ w2*norm(p2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# Example in Test 3 of Zeng 2011
v = [1,1,1,1];
w = [1,-1,1,-1,1];
d = rand(20);
p1 = conv(v,d);
p2 = conv(w,d);

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-9);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test poldiv(p1,d2) ≈ v2*norm(p1) && poldiv(p2,d2) ≈ w2*norm(p2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# testing polgcdvw for random data and fixed d
v = rand(5); w = rand(13); d = ones(20);
#d = Polynomial([100])
p1 = conv(v,d); p2 = conv(w,d);

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-9);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
@test poldiv(p1,d2) ≈ v2*norm(p1) && poldiv(p2,d2) ≈ w2*norm(p2)

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# testing polgcdvw for random data
v = rand(5); w = rand(13); d = rand(20);
#d = Polynomial([100])
p1 = conv(v,d); p2 = conv(w,d);

d2, v2, w2 = polgcdvw(p1,p2,atol=1.e-9);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
q1 = poldiv(p1,d2); q2 = poldiv(p2,d2);
@test q1*v2[end] ≈ v2*q1[end] && q2*w2[end] ≈ w2*q2[end]

#
p=.01
d3, v3, w3, δ3 = gcdvwupd(p1,p2,d2+p*rand(length(d2)),v2+p*rand(length(v2)),w2+p*rand(length(w2)),maxnit=100);
@test norm(d*d3[end] - d3*d[end]) / norm(d) < 1.e-7
q1 = poldiv(p1,d3); q2 = poldiv(p2,d3);
@test q1*v3[end] ≈ v3*q1[end] && q2*w3[end] ≈ w3*q2[end]

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# testing polgcdvw for complex random data
Ty = Complex{Float64}
v = rand(Ty,5); w = rand(Ty,13); d = rand(Ty,20);
#d = Polynomial([100])
p1 = conv(v,d); p2 = conv(w,d);

d2, v2, w2, δ2 = polgcdvw(p1,p2,atol=1.e-9);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
q1 = poldiv(p1,d2); q2 = poldiv(p2,d2);
@test q1*v2[end] ≈ v2*q1[end] && q2*w2[end] ≈ w2*q2[end]

d2, v2, w2, δ2 = polgcdvw(p1,p2,atol=1.e-9,maxnit=10);
@test norm(d*d2[end] - d2*d[end]) / norm(d) < 1.e-7
q1 = poldiv(p1,d2); q2 = poldiv(p2,d2);
@test q1*v2[end] ≈ v2*q1[end] && q2*w2[end] ≈ w2*q2[end]

p = 0.01
d3, v3, w3, δ3 = gcdvwupd(p1,p2,d2+p*rand(Ty,length(d2)),v2+p*rand(Ty,length(v2)),w2+p*rand(Ty,length(w2)),maxnit=10);
@test norm(d*d3[end] - d3*d[end]) / norm(d) < 1.e-7
q1 = poldiv(p1,d3); q2 = poldiv(p2,d3);
@test q1*v3[end] ≈ v3*q1[end] && q2*w3[end] ≈ w3*q2[end]

m = pollcm(p1,p2);
mc = poldiv(conv(p1,p2),d);
@test  mc*m[end]≈ m*mc[end]

# some particular cases

num = Polynomial([2, 1]);
d, v, w, δ = polgcdvw(coeffs(num),[0])
@test d == [1.] && v == [2.0, 1.0] && w == [0.] && δ == 0.

d, v, w, δ = polgcdvw([0],coeffs(num))
@test d == [1.] && w == [2.0, 1.0] && v == [0.] && δ == 0.

d, v, w, δ = polgcdvw(coeffs(num),[5])
@test d == [1.] && v*norm(num) ≈ [2.0, 1.0] && w*5 == [5.] && δ == 0.

d, v, w, δ = polgcdvw([5],coeffs(num))
@test d == [1.] && w*norm(num) ≈ [2.0, 1.0] && v*5 == [5.] && δ == 0.

end

@testset "Rational Matrix Tools" begin

# Example Dopico et al. 2020
e1 = Polynomial(rand(6))
e2 = Polynomial(rand(2))
P = [e1 0; 0 e2]
vals, info = pmeigvals(P,atol=1.e-7)
@test info.id  == [4] # infinite eigenvalue of multiplicity 4

zer, iz, info = pmzeros(P,atol=1.e-7)
@test iz == []

sys = rm2ls(P,atol = 1.e-7,minimal=true)
zer, iz, info = spzeros(sys[1:5]...,atol1=1.e-7,atol2 = 1.e-7)
@test iz == []

c1 = [0.06155268911372547, 0.006545378521362721, 0.40445990039119284, 0.8829892254580274, 0.7573496766341161, 0.4094804382958148];
c2 = [0.4059208742974014, 0.48094667571705574];
e1 = Polynomial(c1)
e2 = Polynomial(c2)

NUM = [e1 0; 1 e2]
DEN = [1 1; Polynomial([0, 1]) 1]

A, E, B, C, D, bl = rm2ls(NUM,DEN,atol = 1.e-7,minimal=true)

z1, iz1, info1 = spzeros(A,E,B,C,D)
@test iz1 == [] && any(abs.(z1) .< 1.e-10)

p1, ip1, info = pzeros(A,E)
@test ip1 == [1, 5]  && any(abs.(p1) .< 1.e-10)

# Example Vlad Ionescu
s = Polynomial([0, 1],:s);
NUM = [s^2+3*s+3 1; -1 2*s^2+7*s+4];
DEN = [(s+1)^2 s+2; (s+1)^3 (s+1)*(s+2)]

sys = rm2ls(NUM,DEN,atol = 1.e-7, minimal=true)
val = float(pi)
@test pmeval(NUM,val) ./ pmeval(DEN,val) ≈ lseval(sys[1:5]...,val)

num1, den1 = ls2rm(sys[1:5]...,atol1 = 1.e-7)
@test all(NUM .* pm2poly(den1,:s) .≈ DEN .* pm2poly(num1,:s))

num1, den1 = ls2rm(sys[1],sys[3:5]...,atol1 = 1.e-7)
@test all(NUM .* pm2poly(den1,:s) .≈ DEN .* pm2poly(num1,:s))

psys = rm2lps(NUM,DEN,atol = 1.e-3, minimal=true)
val = float(pi)
@test pmeval(NUM,val) ./ pmeval(DEN,val) ≈ lpseval(psys[1:8]...,val)

λ = Polynomial([0,1],:λ)
syst=rm2lps(λ)
info = pkstruct(syst[7],syst[8])
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [], 1, 1)

s = Polynomial([0, 1],:s);
NUM1 = [(s+1)*(s+2) s+1; -1 s];
DEN1 = [(s+1)^3 (s+1)*(s+2); (s+1)^3 (s+1)*(s+2)]
sys1 = rm2ls(NUM1,DEN1,atol = 1.e-5, minimal=true);
val = float(pi)
@test pmeval(NUM1,val) ./ pmeval(DEN1,val) ≈ lseval(sys1[1:5]...,val)

@test lsequal(sys[1:5]...,sys1[1:4]...,sys[5])

s = Polynomial([0, 1],:s);
N = [(s+1)*(s+2) s+1; -1 s];
D = [(s+1)^3 0; 0 (s+1)*(s+2)];

sys2 = rpmfd2ls(D,N,atol=1.e-7,minimal = true);
@test lsequal(sys[1:5]...,sys2[1:4]...,sys[5],atol1=1.e-5)
val = float(pi)
@test pmeval(N,val) / pmeval(D,val) ≈ lseval(sys2...,val)

psys2 = rpmfd2lps(D,N,atol=1.e-7,minimal = true);
@test lpsequal(psys[1:8]...,psys2[1:6]...,psys2[7]+psys[7],psys[8],atol1=1.e-5)
val = float(pi)
@test pmeval(N,val) / pmeval(D,val) ≈ lpseval(psys2...,val)

s = Polynomial([0, 1],:s);
N = [(s+1)*(s+2) s+1; -1 s];
D = [(s+1)^3 0; 0 (s+1)*(s+2)];

sys2 = lpmfd2ls(poly2pm(D),poly2pm(N),atol=1.e-7,minimal = true);
val = float(pi)
@test  pmeval(D,val) \ pmeval(N,val)  ≈ lseval(sys2...,val)

sys2 = lpmfd2ls(D,N,atol=1.e-7,minimal = true);
val = float(pi)
@test  pmeval(D,val) \ pmeval(N,val)  ≈ lseval(sys2...,val)

psys2 = lpmfd2lps(poly2pm(D),poly2pm(N),atol=1.e-7,minimal = true);
val = float(pi)
@test  pmeval(D,val) \ pmeval(N,val)  ≈ lpseval(psys2...,val)

psys2 = lpmfd2lps(D,N,atol=1.e-7,minimal = true);
val = float(pi)
@test  pmeval(D,val) \ pmeval(N,val)  ≈ lpseval(psys2...,val)

s = Polynomial([0, 1],:s);
D = [(s+1)^3 0; 0 (s+1)*(s+2)];
sys = pminv2ls(D,atol=1.e-7,minimal = true);
val = float(pi)
@test  pmeval(D,val) * lseval(sys...,val)  ≈ I

NUM, DEN = ls2rm(sys...)
@test  rmeval(NUM,DEN,val) * pmeval(D,val)  ≈ I

sys1 = pminv2lps(poly2pm(D),atol=1.e-7,minimal = true)
val = float(pi)
@test  pmeval(D,val) * lpseval(sys1...,val)  ≈ I

sys1 = pminv2lps(D,atol=1.e-7,minimal = true)
val = float(pi)
@test  pmeval(D,val) * lpseval(sys1...,val)  ≈ I

NUM1, DEN1 = lps2rm(sys1...)
@test  rmeval(NUM1,DEN1,val) * pmeval(D,val)  ≈ I

s = Polynomial([0, 1],:x)
U = [-s  0   s^2 + 2*s + 1
s^2 - s + 1  0 -s^3 - s^2 - 1
0 -1 s];
V = [3*s + 1 3*s^3 + s^2 - 3*s - 4  6
-s -s^3 + s + 1  -2
0  0  1];
@test ispmunimodular(U)
@test ispmunimodular(V)

sys1 = pminv2ls(U,atol=1.e-7,minimal = true)
val = float(pi)
@test  pmeval(U,val) * lseval(sys1...,val)  ≈ I

NUM1, DEN1 = ls2rm(sys1...)
@test  rmeval(NUM1,DEN1,val) * pmeval(U,val)  ≈ I

@test ispmunimodular(NUM1,atol=1.e-7) && DEN1[:,:,1] ≈ ones(3,3)

sys2 = pminv2ls(poly2pm(V),atol=1.e-7,minimal = true)
val = float(pi)
@test  pmeval(V,val) * lseval(sys2...,val)  ≈ I

NUM2, DEN2 = ls2rm(sys2...)
@test  rmeval(NUM2,DEN2,val) * pmeval(V,val)  ≈ I

@test ispmunimodular(NUM2,atol=1.e-7) && DEN2[:,:,1] ≈ ones(3,3)

sys1 = pminv2lps(poly2pm(U),atol=1.e-7,minimal = true)
val = float(pi)
@test  pmeval(U,val) * lpseval(sys1...,val)  ≈ I

NUM1, DEN1 = lps2rm(sys1...,atol1=1.e-7,atol2=1.e-7)
@test  rmeval(NUM1,DEN1,val) * pmeval(U,val)  ≈ I

@test ispmunimodular(NUM1,atol=1.e-7) && DEN1[:,:,1] ≈ ones(3,3)

sys2 = pminv2lps(poly2pm(V),atol=1.e-7,minimal = true)
val = float(pi)
@test  pmeval(V,val) * lpseval(sys2...,val)  ≈ I

NUM2, DEN2 = lps2rm(sys2...,atol1=1.e-7,atol2=1.e-7)
@test  rmeval(NUM2,DEN2,val) * pmeval(V,val)  ≈ I

@test ispmunimodular(NUM2,atol=1.e-7) && DEN2[:,:,1] ≈ ones(3,3)

end

@testset "Rational matrix realizations" begin
s = Polynomial([0, 1],:s)
num = Polynomial([4:-1:1...],:s)
den = Polynomial([7:-1:4...,1],:s)
A, B, C, D, blkdims = rm2lspm(num,den)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,D,1)

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,D,1)

num1, den1 = ls2rm(A,Matrix{Int}(I,4,4),B,C,D)
@test num*pm2poly(den1,:s)[1,1] ≈ pm2poly(num1,:s)[1,1]*den

A1, B1, C1, D1 = rm2lspm(den,num)
@test rmeval(den,num,1) ≈ lseval(A1,I,B1,C1,zeros(1,1),1) + pmeval(D1,1)

num1, den1 = ls2rm(A1,Matrix{Int}(I,3,3),B1,C1,zeros(1,1))
@test den*pm2poly(den1,:s)[1,1] ≈ (pm2poly(num1,:s)[1,1]+pm2poly(D1,:s)[1,1]*pm2poly(den1,:s)[1,1])*num

A, E, B, C, D = rm2ls(num,den,minimal = true)
@test rmeval(num,den,1) ≈ lseval(A,E,B,C,D,1)

num1, den1 = ls2rm(A, E, B, C, D)
@test num*pm2poly(den1,:s)[1,1] ≈ pm2poly(num1,:s)[1,1]*den

A, E, B, C, D = rm2ls(den,num,minimal = true)
@test rmeval(den,num,1) ≈ lseval(A,E,B,C,D,1)

den1, num1 = ls2rm(A, E, B, C, D)
@test num*pm2poly(den1,:s)[1,1] ≈ pm2poly(num1,:s)[1,1]*den

A, E, B, C, D = rm2ls(num;minimal = true)
@test rmeval(num,1) ≈ lseval(A,E,B,C,D,1)

num1, den1 = ls2rm(A, E, B, C, D)
@test num*pm2poly(den1,:s)[1,1] ≈ pm2poly(num1,:s)[1,1]

num = rand(2,3,4);
den = rand(2,3,2) .+ 0.5;
# num = rand(1,1,4)
# den = rand(1,1,2)
A, E, B, C, D = rm2ls(num,den,minimal = true,atol = 1.e-7);
@test rmeval(num,den,1) ≈ lseval(A,E,B,C,D,1)

num1, den1 = ls2rm(A, E, B, C, D, atol1 = 1.e-7, atol2 = 1.e-7);
@test rmeval(num1,den1,1) ≈ lseval(A,E,B,C,D,1)

@test all(pm2poly(num) .* pm2poly(den1) .≈ pm2poly(num1) .* pm2poly(den))

num = rand(2,3,2)
den = rand(2,3,4)
# num = rand(1,1,4)
# den = rand(1,1,2)
A, E, B, C, D = rm2ls(num,den,minimal = true,atol = 1.e-7);
@test rmeval(num,den,1) ≈ lseval(A,E,B,C,D,1)

num1, den1 = ls2rm(A, E, B, C, D, atol1 = 1.e-7, atol2 = 1.e-7);
@test rmeval(num1,den1,1) ≈ lseval(A,E,B,C,D,1)

@test all(pm2poly(num) .* pm2poly(den1) .≈ pm2poly(num1) .* pm2poly(den))

# Example Varga, Sima 1997
s = Polynomial([0, 1],:s)
num = [s 2; 1 s]
den = [s+1 (s+1)*(s+3); s+4 (s+2)*(s+4)]

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,D,1)

A, B, C, D, blkdims = rm2lspm(num,den,obs=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,D,1)

A, E, B, C, D, blkdims = rm2ls(num,den,minimal = true, atol = 1.e-7);
@test rmeval(num,den,1) ≈ lseval(A,E,B,C,D,1)

A, E, B, C, D, blkdims = rm2ls(num,den,minimal = true, contr = true, atol = 1.e-7);
@test rmeval(num,den,1) ≈ lseval(A,E,B,C,D,1)

A, E, B, C, D, blkdims = rm2ls(num,den,minimal = true, obs = true, atol = 1.e-7);
@test rmeval(num,den,1) ≈ lseval(A,E,B,C,D,1)

num1, den1 = ls2rm(A, E, B, C, D, atol1 = 1.e-7, atol2 = 1.e-7);
@test rmeval(num1,den1,1) ≈ lseval(A,E,B,C,D,1)

@test all(num .* pm2poly(den1,:s) .≈ pm2poly(num1,:s) .* den)

# Example 4.3 Antsaklis, Michel 2006
s = Polynomial([0, 1],:s)
num = [s^2+1 s+1]
den = [s^2 s^3]

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,D,1)

A, B, C, D, blkdims = rm2lspm(num,den,obs=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,D,1)

# Example 4.3 (modified) Antsaklis, Michel 2006
s = Polynomial([0, 1],:s)
num = [s^3+1 s+1]
den = [s^2 s^3]

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(1,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,obs=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(1,2),1) + pmeval(D,1)

# Example 4.4 Antsaklis, Michel 2006
s = Polynomial([0, 1],:s)
num = [2 1; 1 0];
den = [s+1 1; s 1];

A, B, C, D, blkdims = rm2lspm(num,den)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,obs=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

s = Polynomial([0, 1],:s)
num = [2+s 1; 1 0];
den = [0.5 1; 0.5 1];

A, B, C, D, blkdims = rm2lspm(num,den)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,obs=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,obs=true)
@test rmeval(num,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)
@test rmeval(poly2pm(num),1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

# Example 4.4 (transposed) Antsaklis, Michel 2006
s = Polynomial([0, 1],:s)*im
num = [2 1; 1 0]
den = [s+1 s; 1 1]

A, B, C, D, blkdims = rm2lspm(num,den)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,contr=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

A, B, C, D, blkdims = rm2lspm(num,den,obs=true)
@test rmeval(num,den,1) ≈ lseval(A,I,B,C,zeros(2,2),1) + pmeval(D,1)

P = zeros(Int,3,3,3);
P[:,:,1] = [1 2 -2; 0 -1 -2; 0 0 0];
P[:,:,2] = [1 3 0; 1 4 2; 0 -1 -2];
P[:,:,3] = [1 4 2; 0 0 0; 1 4 2];

A, B, C, D, blkdims = rm2lspm(P)
@test D == P

A, E, B, C, D, blkdims = rm2ls(P,contr=true)
@test pmeval(P,1) ≈ lseval(A,E,B,C,D,1)

@test P ≈ ls2pm(rm2ls(P)[1:5]...)
@test P ≈ ls2rm(rm2ls(P)[1:5]...)[1]
N, D = ls2rm(rm2ls(P)[1:5]...);
@test all(pm2poly(P) .* pm2poly(D) .≈ pm2poly(N))
@test P ≈ lps2pm(rm2lps(P)[1:8]...;atol1=1.e-7,atol2=1.e-7)
@test P ≈ lps2rm(rm2lps(P)[1:8]...;atol1=1.e-7,atol2=1.e-7)[1]
N, D = lps2rm(rm2lps(P)[1:8]...);
@test all(pm2poly(P) .* pm2poly(D) .≈ pm2poly(N))

end

end  # testset

end # module


text:
# ------------------------------------------------------------

Base.string(subplt::Subplot) = "Subplot{$(subplt.backend) p=$(subplt.p) n=$(subplt.n)}"
Base.print(io::IO, subplt::Subplot) = print(io, string(subplt))
Base.show(io::IO, subplt::Subplot) = print(io, string(subplt))

function Base.copy(subplt::Subplot)
  subplot(subplt.plts, subplt.layout, subplt.plotargs)
end

Base.getindex(subplt::Subplot, args...) = subplt.plts[subplt.layout[args...]]

# --------------------------------------------------------------------

getplot(subplt::Subplot, idx::Int = subplt.n) = subplt.plts[mod1(idx, subplt.p)]
getplotargs(subplt::Subplot, idx::Int) = getplot(subplt, idx).plotargs
convertSeriesIndex(subplt::Subplot, n::Int) = ceil(Int, n / subplt.p)

# ------------------------------------------------------------

function validateSubplotSupported()
  if !subplotSupported()
    error(CURRENT_BACKEND.sym, " does not support the subplot/subplot! commands at this time.  Try one of: ", join(filter(pkg->subplotSupported(_backend_instance(pkg)), backends()),", "))
  end
end

"""
Create a series of plots:
```
  y = rand(100,3)
  subplot(y; n = 3)             # create an automatic grid, and let it figure out the nr/nc... will put plots 1 and 2 on the first row, and plot 3 by itself on the 2nd row
  subplot(y; n = 3, nr = 1)     # create an automatic grid, but fix the number of rows to 1 (so there are n columns)
  subplot(y; n = 3, nc = 1)     # create an automatic grid, but fix the number of columns to 1 (so there are n rows)
  subplot(y; layout = [1, 2])   # explicit layout by row... plot #1 goes by itself in the first row, plots 2 and 3 split the 2nd row (note the n kw is unnecessary)
  subplot(plts, n; nr = -1, nc = -1)  # build a layout from existing plots
  subplot(plts, layout)               # build a layout from existing plots
```
"""
function subplot(args...; kw...)
  validateSubplotSupported()
  d = KW(kw)
  preprocessArgs!(d)

  # for plotting recipes, swap out the args and update the parameter dictionary
  args = RecipesBase.apply_recipe(d, KW(kw), args...; issubplot=true)
  _add_markershape(d)

  # figure out the layout
  layoutarg = get(d, :layout, nothing)
  if layoutarg != nothing
    layout = subplotlayout(layoutarg)
  else
    n = get(d, :n, -1)
    if n < 0
      error("You must specify either layout or n when creating a subplot: ", d)
    end
    layout = subplotlayout(n, get(d, :nr, -1), get(d, :nc, -1))
  end

  # initialize the individual plots
  pkg = backend()
  plts = Plot{typeof(pkg)}[]
  for i in 1:length(layout)
    di = getPlotArgs(pkg, d, i)
    di[:subplot] = true
    dumpdict(di, "Plot args (subplot $i)")
    push!(plts, _create_plot(pkg, di))
  end

  # create the object and do the plotting
  subplt = Subplot(nothing, plts, pkg, length(layout), 0, layout, d, false, false, false, (r,c) -> (nothing,nothing))
  subplot!(subplt, args...; kw...)

  subplt
end

# ------------------------------------------------------------------------------------------------

# NOTE: for the subplot calls building from existing plots, we need the first plot to be separate to ensure dispatch calls this instead of the more general subplot(args...; kw...)

# grid layout
function subplot{P}(plt1::Plot{P}, plts::Plot{P}...; kw...)
  d = KW(kw)
  layout = if haskey(d, :layout)
      subplotlayout(d[:layout])
  else
      subplotlayout(length(plts)+1, get(d, :nr, -1), get(d, :nc, -1))
  end
  # layout = subplotlayout(length(plts)+1, get(d, :nr, -1), get(d, :nc, -1))
  subplot(vcat(plt1, plts...), layout, d)
end

# explicit layout
function subplot{P,I<:Integer}(pltsPerRow::AVec{I}, plt1::Plot{P}, plts::Plot{P}...; kw...)
  layout = subplotlayout(pltsPerRow)
  subplot(vcat(plt1, plts...), layout, KW(kw))
end

# this will be called internally
function subplot{P<:AbstractBackend}(plts::AVec{Plot{P}}, layout::SubplotLayout, d::KW)
  validateSubplotSupported()
  p = length(layout)
  n = sum([plt.n for plt in plts])
  subplt = Subplot(nothing, collect(plts), P(), p, n, layout, KW(), false, false, false, (r,c) -> (nothing,nothing))

  _preprocess_subplot(subplt, d)
  _postprocess_subplot(subplt, d)

  subplt
end

# TODO: hcat/vcat subplots and plots together arbitrarily

# ------------------------------------------------------------------------------------------------

function _preprocess_subplot(subplt::Subplot, d::KW, args = ())
  validateSubplotSupported()
  userkw = preprocessArgs!(d)

  # for plotting recipes, swap out the args and update the parameter dictionary
  args = RecipesBase.apply_recipe(d, userkw, args...; issubplot=true)
  _add_markershape(d)

  dumpdict(d, "After subplot! preprocessing")

  # get the full plotargs, overriding any new settings
  # TODO: subplt.plotargs should probably be merged sooner and actually used
  #       for color selection, etc.  (i.e. if we overwrite the subplot palettes to [:heat :rainbow])
  #       then we need to overwrite plt[1].plotargs[:color_palette] to :heat before it's actually used
  #       for color selection!

  # first merge the new args into the subplot's plotargs. then process the plot args and merge
  # those into the plot's plotargs.  (example... `palette = [:blues :reds]` goes into subplt.plotargs,
  # then the ColorGradient for :blues/:reds is merged into plot 1/2 plotargs, which is then used for color selection)
  for i in 1:length(subplt.layout)
    subplt.plts[i].plotargs = getPlotArgs(backend(), merge(subplt.plts[i].plotargs, d), i)
  end
  merge!(subplt.plotargs, d)

  # process links.  TODO: extract to separate function
  for s in (:linkx, :linky, :linkfunc)
    if haskey(d, s)
      setfield!(subplt, s, d[s])
      delete!(d, s)
    end
  end

  args
end

function _postprocess_subplot(subplt::Subplot, d::KW)
  # init (after plot creation)
  if !subplt.initialized
    subplt.initialized = _create_subplot(subplt, false)
  end

  # add title, axis labels, ticks, etc
  for (i,plt) in enumerate(subplt.plts)
    di = plt.plotargs
    dumpdict(di, "Updating sp $i")
    _update_plot(plt, di)
  end

  _update_plot_pos_size(subplt, d)

  # handle links
  subplt.linkx && link_axis(subplt, true)
  subplt.linky && link_axis(subplt, false)

  # set this to be current
  current(subplt)
end

# ------------------------------------------------------------------------------------------------

"""
Adds to a subplot.
"""

# current subplot
function subplot!(args...; kw...)
  validateSubplotSupported()
  subplot!(current(), args...; kw...)
end

# not allowed:
function subplot!(plt::Plot, args...; kw...)
  error("Can't call subplot! on a Plot!")
end

# # this adds to a specific subplot... most plot commands will flow through here
function subplot!(subplt::Subplot, args...; kw...)
  # validateSubplotSupported()

  d = KW(kw)
  args = _preprocess_subplot(subplt, d, args)

  # create the underlying object (each backend will do this differently)
  # note: we call it once before doing the individual plots, and once after
  #       this is because some backends need to set up the subplots and then plot,
  #       and others need to do it the other way around
  if !subplt.initialized
    subplt.initialized = _create_subplot(subplt, true)
  end

  groupby = if haskey(d, :group)
      extractGroupArgs(d[:group], args...)
  else
      nothing
  end

  _add_series_subplot(subplt, d, groupby, args...)
  _postprocess_subplot(subplt, d)

  # show it automatically?
  if haskey(d, :show) && d[:show]
    gui()
  end

  subplt
end

# not allowed:
function plot!(subplt::Subplot, args...; kw...)
  error("Can't call plot! on a Subplot!")
end

# given a fully processed KW, add the series to the Plot
function _add_series_subplot(plt::Plot, d::KW)
  # setTicksFromStringVector(d, d, :x, :xticks)
  # setTicksFromStringVector(d, d, :y, :yticks)
  setTicksFromStringVector(plt, d, d, "x")
  setTicksFromStringVector(plt, d, d, "y")
  setTicksFromStringVector(plt, d, d, "z")

  # this is the actual call to the backend
  _add_series(plt.backend, plt, d)

  _add_annotations(plt, d)
  warnOnUnsupportedScales(plt.backend, d)
end

# handle the grouping... add a series for each group
function _add_series_subplot(subplt::Subplot, d::KW, groupby::GroupBy, args...)
    starting_n = subplt.n
    for (i, glab) in enumerate(groupby.groupLabels)
        tmpd = copy(d)
        tmpd[:numUncounted] = subplt.n - starting_n
        _add_series_subplot(subplt, tmpd, nothing, args...;
                    idxfilter = groupby.groupIds[i],
                    grouplabel = string(glab))
    end
end

# process, filter, and add to the correct plot
function _add_series_subplot(subplt::Subplot, d::KW, ::Void, args...;
                             idxfilter = nothing,
                             grouplabel = "")
    process_inputs(subplt, d, args...)

    if idxfilter != nothing
        # add the group name as the label if there isn't one passed in
        get!(d, :label, grouplabel)
        # filter the data
        filter_data!(d, idxfilter)
    end

    kwList, xmeta, ymeta = build_series_args(subplt, d)

    # TODO: something useful with meta info?

    for (i,di) in enumerate(kwList)

      subplt.n += 1
      plt = getplot(subplt)
      plt.n += 1

      # cleanup the dictionary that we pass into the plot! command
      di[:show] = false
      di[:subplot] = true
      for k in (:title, :xlabel, :xticks, :xlims, :xscale, :xflip,
                        :ylabel, :yticks, :ylims, :yscale, :yflip)
        delete!(di, k)
      end
      dumpdict(di, "subplot! kwList $i")
      dumpdict(plt.plotargs, "plt.plotargs before plotting")

      _replace_linewidth(di)
      _add_series_subplot(plt, di)
    end
end

# --------------------------------------------------------------------

# handle "linking" the subplot axes together
# each backend should implement the _remove_axis and _expand_limits methods
function link_axis(subplt::Subplot, isx::Bool)

  # collect the list of plots and the expanded limits for those plots that should be linked on this axis
  includedPlots = Any[]
  # lims = [Inf, -Inf]
  lims = Dict{Int,Any}()  # maps column to xlim
  for (i,(r,c)) in enumerate(subplt.layout)

    # shouldlink will be a bool or nothing.  if nothing, then use linkx/y (which is true if we get to this code)
    shouldlink = subplt.linkfunc(r,c)[isx ? 1 : 2]
    if shouldlink == nothing || shouldlink
      plt = subplt.plts[i]

      # if we don't have this
      k = isx ? c : r
      if (firstone = !haskey(lims, k))
        lims[k] = [Inf, -Inf]
      end

      isinner = (isx && r < nrows(subplt.layout)) || (!isx && !firstone)
      push!(includedPlots, (plt, isinner, k))

      _expand_limits(lims[k], plt, isx)
    end

  end

  # do the axis adjustments
  for (plt, isinner, k) in includedPlots
    if isinner
      _remove_axis(plt, isx)
    end
    (isx ? xlims! : ylims!)(plt, lims[k]...)
  end
end


text:
using Test, Astro

# earth tests from meeus chapter 11

print("Geodesic distance")
# geodesic distance between Paris and Washington observatories
# should be 6181.628423726299 km

d = geodesic_distance(deg2rad(dms_to_d(-2, 20, 14)),deg2rad(dms_to_d(48, 50, 11)), deg2rad(dms_to_d(77,3,56)), deg2rad(dms_to_d(38,55,17)))
@test isapprox(d, 6181.6, atol=.1)
println(" passed")

print("Geocentric latitude")
# p sin phi and p cos phi

r = geographical_to_geocentric_lat(deg2rad(dms_to_d(33, 21, 22)), 1706)
# TODO earth-tests.jl Is this accurate enough?
@test rad2deg(r[1]) ≈  33.356111 atol= .2# <<<----------------------?
@test isapprox(r[2], 0.546, atol=.1)
@test isapprox(r[3], 0.83, atol=.1)
println(" passed")


text:
module COVIDResourceAllocation

include("models/PatientAllocation.jl")
include("models/ReusableResourceAllocation.jl")
include("models/DisposableResourceAllocation.jl")

include("processing/BedsData.jl")
include("processing/ForecastData.jl")
include("processing/GeographicData.jl")
include("processing/NurseData.jl")

include("util/PatientAllocationResults.jl")
include("util/NurseAllocationResults.jl")

import .PatientAllocation: patient_allocation
import .ReusableResourceAllocation: reusable_resource_allocation
import .ForecastData: forecast
import .BedsData: n_beds
import .GeographicData: adjacencies
import .NurseData: n_nurses
import .PatientAllocationResults
import .NurseAllocationResults

export patient_allocation, reusable_resource_allocation
export forecast, adjacencies, n_beds, n_nurses
export PatientAllocationResults, NurseAllocationResults

end

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation