text:
export AbstractSolution,
       ReachSolution,
       CheckSolution,
       solve,
       project

"""
    AbstractSolution

Abstract type representing the solution of a rechability problem.
"""
abstract type AbstractSolution end

"""
    ReachSolution{S<:LazySet} <: AbstractSolution

Type that wraps a the solution of a reachability problem as a sequence of
lazy sets, and a dictionary of options.

### Fields

- `Xk`       -- the list of reachable states
- `options`  -- the dictionary of options

### Notes

If the solution has been projected in 2D, the sequence `Xk` is an array
of polygons in constraint representation. In high-dimensions this is a sequence
of cartesian product arrays of low-dimensional sets.
"""
struct ReachSolution{S<:LazySet} <: AbstractSolution
  Xk::Vector{S}
  options::Options
end
# constructor with no options
ReachSolution(Xk::Vector{S}) where {S<:LazySet} =
    ReachSolution{S}(Xk, Options())

"""
    CheckSolution

Type that wraps a the solution of a property checking problem, which is just the
answer if the property is satisfied.

### Fields

- `satisfied` -- is the property satisfied?
- `violation` -- step at which the property is violated (-1 otherwise)
- `options`   -- the dictionary of options
"""
struct CheckSolution <: AbstractSolution
  satisfied::Bool
  violation::Int
  options::Options
end
# constructor with no options
CheckSolution(satisfied::Bool, violation::Int) =
    CheckSolution(satisfied, violation, Options())

"""
    solve(system, options)  or  solve(system, :key1 => val1, [...], keyK => valK)

Solves a reachability problem s.t. the given options.
If some options are not defined, we may fall back to default values.

### Input

- `system`  -- a (discrete or continuoues) system specification
- `options` -- options for solving the problem

### Output

A solution object whose content depends on the input options.

### Notes

To see all available input options, see
`keys(Reachability.available_keywords.dict)`.
"""
function solve(system::InitialValueProblem,
               options_input::Options)::AbstractSolution

    # ==========
    # Dimensions
    # ==========
    options_input.dict[:n] = statedim(system)

    # =======
    # Options
    # =======
    # solver-specific options (adds default values for unspecified options)
    options = validate_solver_options_and_add_default_values!(options_input)

    # ===================
    # Time discretization
    # ===================
    if system isa InitialValueProblem{<:AbstractContinuousSystem}
        info("Time discretization...")
        tic()
        Δ = discretize(
            system,
            options[:δ],
            approx_model=options[:approx_model],
            pade_expm=options[:pade_expm],
            lazy_expm=options[:lazy_expm_discretize],
            lazy_sih=options[:lazy_sih]
            )
        tocc()
    else
        Δ = system
    end

    # ==============
    # Transformation
    # ==============
    if options[:coordinate_transformation] != ""
        info("Transformation...")
        tic()
        (Δ, transformation_matrix) = transform(Δ, options[:coordinate_transformation])
        tocc()
    else
        transformation_matrix = nothing
    end

    # ===================================
    # Sparse/dense/lazy matrix conversion
    # ===================================
    A = Δ.s.A
    create_new_system = false
    if !options[:lazy_expm] && options[:lazy_expm_discretize]
        # convert SparseMatrixExp to eplicit matrix
        info("Making lazy matrix exponential explicit...")
        tic()
        n = options.dict[:n]
        if options[:assume_sparse]
            B = sparse(Int[], Int[], eltype(A)[], n, n)
        else
            B = Matrix{eltype(A)}(n, n)
        end
        for i in 1:n
            B[i, :] = get_row(A, i)
        end
        A = B
        create_new_system = true
        tocc()
    end
    if options[:assume_sparse]
        if A isa SparseMatrixExp
            # ignore this case
        elseif !method_exists(sparse, Tuple{typeof(A)})
            info("`assume_sparse` option cannot be applied to a matrix of " *
                 "type $(typeof(A)) and will be ignored")
        elseif !(A isa AbstractSparseMatrix)
            # convert to sparse matrix
            A = sparse(A)
            create_new_system = true
        end
    end
    if create_new_system
        # set new matrix
        if method_exists(inputset, Tuple{typeof(Δ.s)})
            Δ = DiscreteSystem(A, Δ.x0, inputset(Δ))
        else
            Δ = DiscreteSystem(A, Δ.x0)
        end
    end

    if options[:mode] == "reach"
        # ================================
        # Input -> Output variable mapping
        # ================================
        options.dict[:inout_map] =
            inout_map_reach(options[:partition], options[:blocks], options[:n])

        # ============================
        # Reachable states computation
        # ============================
        info("Reachable States Computation...")
        tic()
        output_function = options[:project_reachset] ?
            nothing :
            options[:projection_matrix]
        Rsets = reach(
            Δ,
            options[:N];
            algorithm=options[:algorithm],
            ε_init=options[:ε_init],
            set_type_init=options[:set_type_init],
            ε_iter=options[:ε_iter],
            set_type_iter=options[:set_type_iter],
            assume_sparse=options[:assume_sparse],
            assume_homogeneous=options[:assume_homogeneous],
            lazy_X0=options[:lazy_X0],
            blocks=options[:blocks],
            partition=options[:partition],
            block_types_init=options[:block_types_init],
            block_types_iter=options[:block_types_iter],
            template_directions_init=options[:template_directions_init],
            template_directions_iter=options[:template_directions_iter],
            lazy_inputs_interval=options[:lazy_inputs_interval],
            output_function=output_function
            )
        info("- Total")
        tocc()

        # ==========
        # Projection
        # ==========
        if options[:project_reachset] || options[:projection_matrix] != nothing
            info("Projection...")
            tic()
            RsetsProj = project(Rsets, options;
                                transformation_matrix=transformation_matrix)
            tocc()
            return ReachSolution(RsetsProj, options)
        end

        return ReachSolution(Rsets, options)

    elseif options[:mode] == "check"
        # ============================================
        # Input -> Output variable mapping in property
        # ============================================
        options.dict[:property] = inout_map_property(options[:property],
            options[:partition], options[:blocks], options[:n])

        # =================
        # Property checking
        # =================
        info("Property Checking...")
        tic()
        answer = check_property(
            Δ,
            options[:N];
            algorithm=options[:algorithm],
            ε_init=options[:ε_init],
            set_type_init=options[:set_type_init],
            ε_iter=options[:ε_iter],
            set_type_iter=options[:set_type_iter],
            assume_sparse=options[:assume_sparse],
            assume_homogeneous=options[:assume_homogeneous],
            lazy_X0=options[:lazy_X0],
            blocks=options[:blocks],
            partition=options[:partition],
            block_types_init=options[:block_types_init],
            block_types_iter=options[:block_types_iter],
            template_directions_init=options[:template_directions_init],
            template_directions_iter=options[:template_directions_iter],
            eager_checking=options[:eager_checking],
            property=options[:property],
            lazy_inputs_interval=options[:lazy_inputs_interval]
            )
        info("- Total")
        tocc()

        if answer == 0
            info("The property is satisfied!")
            return CheckSolution(true, -1, options)
        else
            info("The property may be violated at index $answer," *
                " (time point $(answer * options[:δ]))!")
            return CheckSolution(false, answer, options)
        end
    else
        error("Unsupported mode.")
    end
end

solve(system::AbstractSystem, options::Pair{Symbol,<:Any}...) =
    solve(system, Options(Dict{Symbol,Any}(options)))

"""
    project(Rsets, options; [transformation_matrix])

Projects a sequence of sets according to the settings defined in the options.

### Input

- `Rsets`   -- solution of a reachability problem
- `options` -- options structure
- `transformation_matrix` -- (optional, default: nothing) matrix implementing
                              the transformation)

### Notes

A projection matrix can be given in the options structure, or passed as a
dictionary entry.
"""
function project(Rsets::Vector{<:LazySet}, options::Options;
                 transformation_matrix=nothing)
    plot_vars = copy(options[:plot_vars])
    for i in 1:length(plot_vars)
        if plot_vars[i] != 0
            plot_vars[i] = options[:inout_map][plot_vars[i]]
        end
    end
    reduced_n = sum(x -> x != 0, options[:inout_map])
    output_function = !options[:project_reachset]
    RsetsProj = project_reach(plot_vars,
                              reduced_n,
                              options[:δ],
                              Rsets,
                              options[:algorithm],
                              ε=options[:ε_proj],
                              set_type=options[:set_type_proj],
                              transformation_matrix=transformation_matrix,
                              projection_matrix=options[:projection_matrix],
                              output_function=output_function
                             )
end

project(reach_sol::AbstractSolution) = project(reach_sol.Xk, reach_sol.options)

project(Rsets::Vector{<:LazySet}, options::Pair{Symbol,<:Any}...) =
    project(Rsets, Options(Dict{Symbol,Any}(options)))


text:
# https://github.com/bokeh/Bokeh.jl

supported_args(::BokehBackend) = merge_with_base_supported([
    # :annotations,
    # :axis,
    # :background_color,
    :linecolor,
    # :color_palette,
    # :fillrange,
    # :fillcolor,
    # :fillalpha,
    # :foreground_color,
    :group,
    # :label,
    # :layout,
    # :legend,
    :seriescolor, :seriesalpha,
    :linestyle,
    :seriestype,
    :linewidth,
    # :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    # :markeralpha,
    # :markerstrokewidth,
    # :markerstrokecolor,
    # :markerstrokestyle,
    # :n,
    # :bins,
    # :nc,
    # :nr,
    # :pos,
    # :smooth,
    # :show,
    :size,
    :title,
    # :window_title,
    :x,
    # :xguide,
    # :xlims,
    # :xticks,
    :y,
    # :yguide,
    # :ylims,
    # :yrightlabel,
    # :yticks,
    # :xscale,
    # :yscale,
    # :xflip,
    # :yflip,
    # :z,
    # :tickfont,
    # :guidefont,
    # :legendfont,
    # :grid,
    # :surface,
    # :levels,
  ])
supported_types(::BokehBackend) = [:path, :scatter]
supported_styles(::BokehBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
supported_markers(::BokehBackend) = [:none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5]
supported_scales(::BokehBackend) = [:identity, :ln]
is_subplot_supported(::BokehBackend) = false

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

function _initialize_backend(::BokehBackend; kw...)
  @eval begin
    warn("Bokeh is no longer supported... many features will likely be broken.")
    import Bokeh
    export Bokeh
  end
end

const _glyphtypes = KW(
    :circle    => :Circle,
    :rect       => :Square,
    :diamond    => :Diamond,
    :utriangle  => :Triangle,
    :dtriangle  => :InvertedTriangle,
    # :pentagon   =>
    # :hexagon    =>
    # :heptagon   =>
    # :octagon    =>
    :cross      => :Cross,
    :xcross     => :X,
    :star5      => :Asterisk,
  )

function bokeh_glyph_type(d::KW)
  st = d[:seriestype]
  mt = d[:markershape]
  if st == :scatter && mt == :none
    mt = :circle
  end

  # if we have a marker, use that
  if st == :scatter || mt != :none
    return _glyphtypes[mt]
  end

  # otherwise return a line
  return :Line
end

function get_stroke_vector(linestyle::Symbol)
  dash = 12
  dot = 3
  gap = 2
  linestyle == :solid && return Int[]
  linestyle == :dash && return Int[dash, gap]
  linestyle == :dot && return Int[dot, gap]
  linestyle == :dashdot && return Int[dash, gap, dot, gap]
  linestyle == :dashdotdot && return Int[dash, gap, dot, gap, dot, gap]
  error("unsupported linestyle: ", linestyle)
end

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

# function _create_plot(pkg::BokehBackend, d::KW)
function _create_backend_figure(plt::Plot{BokehBackend})
  # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
  # TODO: initialize the plot... title, xlabel, bgcolor, etc

  datacolumns = Bokeh.BokehDataSet[]
  tools = Bokeh.tools()
  filename = tempname() * ".html"
  title = plt.attr[:title]
  w, h = plt.attr[:size]
  xaxis_type = plt.attr[:xscale] == :log10 ? :log : :auto
  yaxis_type = plt.attr[:yscale] == :log10 ? :log : :auto
  # legend = plt.attr[:legend] ? xxxx : nothing
  legend = nothing
  extra_args = KW()  # TODO: we'll put extra settings (xlim, etc) here
  Bokeh.Plot(datacolumns, tools, filename, title, w, h, xaxis_type, yaxis_type, legend) #, extra_args)

  # Plot(bplt, pkg, 0, d, KW[])
end

# function _series_added(::BokehBackend, plt::Plot, d::KW)
function _series_added(plt::Plot{BokehBackend}, series::Series)
  bdata = Dict{Symbol, Vector}(:x => collect(series.d[:x]), :y => collect(series.d[:y]))

  glyph = Bokeh.Bokehjs.Glyph(
      glyphtype = bokeh_glyph_type(d),
      linecolor = webcolor(d[:linecolor]),  # shape's stroke or line color
      linewidth = d[:linewidth],          # shape's stroke width or line width
      fillcolor = webcolor(d[:markercolor]),
      size      = ceil(Int, d[:markersize] * 2.5),  # magic number 2.5 to keep in same scale as other backends
      dash      = get_stroke_vector(d[:linestyle])
    )

  legend = nothing  # TODO
  push!(plt.o.datacolumns, Bokeh.BokehDataSet(bdata, glyph, legend))

  # push!(plt.seriesargs, d)
  # plt
end

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

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot_object(plt::Plot{BokehBackend}, d::KW)
end

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

# accessors for x/y data

# function getxy(plt::Plot{BokehBackend}, i::Int)
#   series = plt.o.datacolumns[i].data
#   series[:x], series[:y]
# end
#
# function setxy!(plt::Plot{BokehBackend}, xy::Tuple{X,Y}, i::Integer)
#   series = plt.o.datacolumns[i].data
#   series[:x], series[:y] = xy
#   plt
# end

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

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

function Base.writemime(io::IO, ::MIME"image/png", plt::AbstractPlot{BokehBackend})
  # TODO: write a png to io
  warn("mime png not implemented")
end

function Base.display(::PlotsDisplay, plt::Plot{BokehBackend})
  Bokeh.showplot(plt.o)
end

# function Base.display(::PlotsDisplay, plt::Subplot{BokehBackend})
#   # TODO: display/show the subplot
# end

ext: jl


text:
using Documenter, SimpleTranslations

makedocs(;
    modules=[SimpleTranslations],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "API" => "api.md",
    ],
    repo="https://github.com/dmolina/SimpleTranslations.jl/blob/{commit}{path}#L{line}",
    sitename="SimpleTranslations.jl",
    authors="Daniel Molina",
)

deploydocs(;
    repo="github.com/dmolina/SimpleTranslations.jl",
)


text:
@testset "Trafo" begin
  A = matrix(FlintZZ, [0 1 0 0 0;
                       0 0 2 0 0;
                       0 0 3 0 0;
                       0 0 0 4 0;
                       5 0 0 0 0])
  Asparse = sparse_matrix(A)

  v = fmpz[1, 2, 3, 4, 5]

  T = @inferred sparse_trafo_scale(2, fmpz(-2))
  @inferred Hecke.apply_left!(Asparse, T)
  @test Asparse == sparse_matrix(FlintZZ, [0 1 0 0 0;
                                           0 0 -4 0 0;
                                           0 0 3 0 0;
                                           0 0 0 4 0;
                                           5 0 0 0 0])
  @inferred Hecke.apply_right!(v, T)
  @test v == fmpz[1, -4, 3, 4, 5]

  T = @inferred sparse_trafo_swap(fmpz, 5, 4)
  @inferred Hecke.apply_left!(Asparse, T)
  @test Asparse == sparse_matrix(FlintZZ, [0 1 0 0 0;
                                           0 0 -4 0 0;
                                           0 0 3 0 0;
                                           5 0 0 0 0;
                                           0 0 0 4 0])
  @inferred Hecke.apply_right!(v, T)
  @test v == fmpz[1, -4, 3, 5, 4]

  T = @inferred sparse_trafo_add_scaled(3, 2, fmpz(2))
  @inferred Hecke.apply_left!(Asparse, T)
  @test Asparse == sparse_matrix(FlintZZ, [0 1 0 0 0;
                                           0 0 2 0 0;
                                           0 0 3 0 0;
                                           5 0 0 0 0;
                                           0 0 0 4 0])
  @inferred Hecke.apply_right!(v, T)
  @test v == fmpz[1, -4, -5, 5, 4]

  T = @inferred sparse_trafo_para_add_scaled(2, 3, fmpz(2), fmpz(-1), fmpz(3), fmpz(-2))
  @inferred Hecke.apply_left!(Asparse, T)
  @test Asparse == sparse_matrix(FlintZZ, [0 1 0 0 0;
                                           0 0 1 0 0;
                                           0 0 0 0 0;
                                           5 0 0 0 0;
                                           0 0 0 4 0])
  @inferred Hecke.apply_right!(v, T)
  @test v == fmpz[1, -23, 14, 5, 4]

  Asparse = sparse_matrix(FlintZZ, [1 1 0 0 0;
                                    0 1 2 0 0;
                                    0 0 1 2 1;
                                    0 0 0 4 0;
                                    0 0 0 0 1])

  B = matrix(FlintZZ, [1 2 3; 4 5 6; 7 8 9])
  T = @inferred sparse_trafo_partial_dense(3, 3:5, 3:5, B)
  @inferred Hecke.apply_left!(Asparse, T)
  @test Asparse == sparse_matrix(FlintZZ, [1 1 0 0 0;
                                           0 1 2 0 0;
                                           0 0 1 10 4;
                                           0 0 4 28 10;
                                           0 0 7 46 16])
  @inferred Hecke.apply_right!(v, T)
  @test v == fmpz[1, -23, 62, 85, 108]

  Asparse = sparse_matrix(FlintZZ, [1 1 0 0 0;
                                    0 0 0 0 0;
                                    0 0 1 10 4;
                                    0 0 4 28 10;
                                    0 0 7 46 16])

  T = @inferred sparse_trafo_move_row(fmpz, 2, 5)
  @inferred Hecke.apply_left!(Asparse, T)
  @test Asparse == sparse_matrix(FlintZZ, [1 1 0 0 0;
                                           0 0 1 10 4;
                                           0 0 4 28 10;
                                           0 0 7 46 16;
                                           0 0 0 0 0])
  @inferred Hecke.apply_right!(v, T)
  @test v == fmpz[1, 108, -23, 62, 85]

end


text:
using FermiCG
using Test

@testset "strings" begin
    display(StringCI.string_to_index("110010"))
    @test StringCI.string_to_index("110010") == 19
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ContainerServiceWindowsProfile <: SwaggerModel
    adminUsername::Any # spec type: Union{ Nothing, String } # spec name: adminUsername
    adminPassword::Any # spec type: Union{ Nothing, String } # spec name: adminPassword

    function ContainerServiceWindowsProfile(;adminUsername=nothing, adminPassword=nothing)
        o = new()
        validate_property(ContainerServiceWindowsProfile, Symbol("adminUsername"), adminUsername)
        setfield!(o, Symbol("adminUsername"), adminUsername)
        validate_property(ContainerServiceWindowsProfile, Symbol("adminPassword"), adminPassword)
        setfield!(o, Symbol("adminPassword"), adminPassword)
        o
    end
end # type ContainerServiceWindowsProfile

const _property_map_ContainerServiceWindowsProfile = Dict{Symbol,Symbol}(Symbol("adminUsername")=>Symbol("adminUsername"), Symbol("adminPassword")=>Symbol("adminPassword"))
const _property_types_ContainerServiceWindowsProfile = Dict{Symbol,String}(Symbol("adminUsername")=>"String", Symbol("adminPassword")=>"String")
Base.propertynames(::Type{ ContainerServiceWindowsProfile }) = collect(keys(_property_map_ContainerServiceWindowsProfile))
Swagger.property_type(::Type{ ContainerServiceWindowsProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ContainerServiceWindowsProfile[name]))}
Swagger.field_name(::Type{ ContainerServiceWindowsProfile }, property_name::Symbol) =  _property_map_ContainerServiceWindowsProfile[property_name]

function check_required(o::ContainerServiceWindowsProfile)
    (getproperty(o, Symbol("adminUsername")) === nothing) && (return false)
    (getproperty(o, Symbol("adminPassword")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ContainerServiceWindowsProfile }, name::Symbol, val)
    if name === Symbol("adminUsername")
    end
    if name === Symbol("adminPassword")
    end
end


text:
mutable struct DDSolverData
    max_iterations::Int64
    epsilon::Float64
    time_limit::Float64
    alpha::Float64
    delta::Float64
    max_no_improvement::Int64

    function DDSolverData()
        data = new()
        data.max_iterations = 100
        data.epsilon = 0.001
        data.time_limit = 3600
        data.alpha = 2.0  #2.0
        data.delta = 0.5
        data.max_no_improvement = 10
        return data
    end
end

#IDEA: Create a DDModel using ModelGraph information
mutable struct DDModel
    solver_data::DDSolverData   #solver inputs
    solution::DDSolution        #solution data

    #Callbacks we can expose
    multiplier_update_function::Function                        #function that takes a DDModel argument and returns multiplier updates for equality and inequality multipliers
    upper_bound_function::Union{Nothing,Function}               #optional, can speed up convergence
    subproblem_solver::Union{Nothing,JuMP.OptimizerFactory}     #e.g. GLPK.Optimizer

    subproblems::Vector{JuMP.Model}                             #we assume these are all minimization problems

    linkvar_multipliers::Vector{Float64}
    equality_multipliers::Vector{Float64}                       #link equality constraint multipliers
    inequality_multipliers::Vector{Float64}                     #link inequality constraint multipliers

    link_var_matrix::Union{Nothing,AbstractMatrix}              #nonanticpitivity constraint matrix
    link_eq_matrix::Union{Nothing,AbstractMatrix}               #link equality matrix
    link_ineq_matrix::Union{Nothing,AbstractMatrix}             #link inequality matrix

    b_eq::Vector{Float64}                                       #link equality right hand side
    b_ineq::Vector{Float64}                                     #link inequality right hand side

    link_var_variables::Vector{JuMP.VariableRef}                #Vector of JuMP variables corresponding to nonanticpitivity constraints
    link_eq_variables::Vector{JuMP.VariableRef}                 #Vector of JuMP variables corresponding to link equality constraints
    link_ineq_variables::Vector{JuMP.VariableRef}               #Vector of JuMP variables corresponding to link inequality constraints

    current_lower_bound::Float64
    current_upper_bound::Float64

    residuals_linkvars::Vector           #A_linkvar*x_linkvar
    residuals_equality::Vector           #A_eq*x_eq - b_eq
    residuals_inequality::Vector         #A_ineq*x_ineq - b_ineq
end

# NOTE: We can lift shared variables into equality constraints
# shared_variables::Dict{JuMP.AbstractJuMPScalar,JuMP.AbstractJuMPScalar} #these get dualized into simple linkconstraints
# shared_variable_multipliers::Vector{Float64}
#dd_model.shared_variables = getsharedvariables(graph)  #dictionary

DDModel() = DDModel(DDSolverData(),DDSolution(),subgradient!,nothing,nothing,
                                Vector{JuMP.Model}(),
                                Vector{Float64}(),Vector{Float64}(),Vector{Float64}(),nothing,nothing,nothing,Vector{Float64}(),Vector{Float64}(),
                                Vector{JuMP.VariableRef}(),Vector{JuMP.VariableRef}(),Vector{JuMP.VariableRef}(),-Inf,Inf,
                                Vector{Float64}(),Vector{Float64}(),Vector{Float64}())

#Dual decomposition algorithm
#NOTE: This might get renamed to lagrange_solver depending how we deal with ADMM
function dual_decomposition_solve(graph::ModelGraph,args...;kwargs...)

    #NOTE: figure out which args to pass
    dd_model = DDModel(graph)
    status = solve(dd_model)

    #set solution on graph

    return status
end

function dual_decomposition_pmap_solve(graph::ModelGraph,workers,args...;kwargs...)
end

#Populate a DDModel using a ModelGraph
function DDModel(graph::ModelGraph) #args,kwargs

    #DDMODEL SUPPORTS
    #supports: linkconstraints
    #supports: linkvariables
    #does not support: subgraphs (need to aggregate)
    #does not support: interval constraints
    #does not support: nonlinear constraints

    #TODO Automatically aggregate subgraphs
    if has_subgraphs(graph)
        error("Dual Decomposition Solver does not support subgraphs.  run aggregate_subgraphs(graph) first" )
    end

    #Fill in all the data we need to do dual decomposition
    dd_model = DDModel()

    #Setup the subproblems based on modelnodes and linkvariables
    subproblems,sub_var_map,master_var_map =  build_subproblems(graph)

    dd_model.subproblems = subproblems               #copies model nodes and adds  master constraints to subproblems

    link_eq_constraints = get_link_eq_constraints(graph)          #link equality constraints in graph
    link_ineq_constraints = get_link_ineq_constraints(graph)      #link inequality constraints in graph
    link_variables = [var for var in getlinkvariables(graph)]     #link variables TODO: Make ModelGraphs return concrete type

    #Setup data structures
    link_eq_matrix,b_eq,link_eq_variables = prepare_link_matrix(link_eq_constraints,sub_var_map)
    link_ineq_matrix,b_ineq,link_ineq_variables, = prepare_link_matrix(link_ineq_constraints,sub_var_map)
    link_var_matrix,link_var_variables = prepare_link_matrix(link_variables,master_var_map)                 #This is the big sparse H matrix for nonanticpitivity constraints

    dd_model.link_var_matrix = link_var_matrix
    dd_model.link_var_variables = link_var_variables

    dd_model.link_eq_matrix = link_eq_matrix
    dd_model.b_eq = b_eq
    dd_model.link_eq_variables = link_eq_variables

    dd_model.link_ineq_matrix = link_ineq_matrix
    dd_model.b_ineq = b_ineq
    dd_model.link_ineq_variables = link_ineq_variables

    #Setup multipliers
    dd_model.equality_multipliers = zeros(size(link_eq_matrix,1))
    dd_model.inequality_multipliers = zeros(size(link_ineq_matrix,1))
    dd_model.linkvar_multipliers = zeros(size(link_var_matrix,1))

    #TODO Multiplier warm-start

    #Setup subproblem multiplier vectors
    _prepare_eq_multiplier_map!(link_eq_matrix,link_eq_variables)
    _prepare_ineq_multiplier_map!(link_ineq_matrix,link_ineq_variables)
    _prepare_linkvar_multiplier_map!(link_var_matrix,link_var_variables)

    return dd_model
end

#Solve a DD Model using dual decomposition or ADMM
function optimize!(dd_model::DDModel)

    println("Optimizing with Dual Decomposition:")
    println("# of subproblems: $(length(dd_model.subproblems))")
    println("# of nonanticpitivity constraints: ")
    println("# of link equality constraints:")
    println("# of link inequality constraints:")

    Zk_old = dd_model.current_lower_bound
    max_iterations = dd_model.solver_data.max_iterations
    epsilon = dd_model.solver_data.epsilon
    delta = dd_model.solver_data.delta

    no_improvement = 0
    for iter in 1:max_iterations
        println("Iteration: ",iter)
        # Solve subproblems
        Zk = 0  #objective value
        for subproblem in dd_model.subproblems
           Zkn = solve_subproblem!(dd_model,subproblem)
           Zk += Zkn                       # add objective values
        end
        #steptaken = true

        dd_model.current_lower_bound = Zk  #Update lower bound
        # If no improvement in the lowerbound, increase the no improvement counter
        if Zk < Zk_old
            no_improvement += 1
        end
        Zk_old = Zk

        # If too many iterations have happened without improvement in the lower bound, decrease alpha (step-size)
        if no_improvement >= dd_model.solver_data.max_no_improvement
            no_improvement = 0
            dd_model.solver_data.alpha *= delta
        end

        A_eq = dd_model.link_eq_matrix
        A_ineq = dd_model.link_ineq_matrix
        H_linkvars = dd_model.link_var_matrix

        x_eq = JuMP.value.(dd_model.link_eq_variables)
        x_ineq = JuMP.value.(dd_model.link_ineq_variables)
        x_linkvars = JuMP.value.(dd_model.link_var_variables)

        # Update residuals for multplier calculation
        dd_model.residuals_equality = A_eq*x_eq - dd_model.b_eq         #Ax = b
        dd_model.residuals_inequality = A_ineq*x_ineq - dd_model.b_ineq #Ax <= b
        dd_model.residuals_linkvars = H_linkvars*x_linkvars             #Hx = 0

        # Check convergence
        if norm(dd_model.residuals_equality) + norm(dd_model.residuals_linkvars) < epsilon && all(dd_model.residuals_inequality .>= 0)
            dd_model.solution.termination = :Optimal
            return :Optimal
        end

        #multiplier updates
        eq_multipliers_delta, ineq_multipliers_delta, linkvar_multipliers_delta = dd_model.multiplier_update_function(dd_model)
        dd_model.linkvar_multipliers += linkvar_multipliers_delta
        dd_model.equality_multipliers += eq_multipliers_delta
        dd_model.inequality_multipliers += ineq_multipliers_delta
        dd_model.inequality_multipliers[dd_model.inequality_multipliers .<= 0] .= 0  #NOTE: set to zero, or just don't update?
    end

    dd_model.solution.termination= :MaxIterationsReached
    return :MaxIterationsReached

end

#Use current multiplier values to solve a lagrange sub-problem
function solve_subproblem!(dd_model::DDModel,subproblem::JuMP.Model)

    #IDEA: Use parameter JuMP and make the multipliers into parameters
    m = subproblem

    eq_multipliers = dd_model.equality_multipliers
    ineq_multipliers = dd_model.inequality_multipliers
    linkvar_multipliers = dd_model.linkvar_multipliers

    scale = m.ext[:objective_scale]
    JuMP.set_objective_function(m,scale*m.ext[:original_objective])

    #Add dualized part to objective function for each linkconstraint that hits this node
    obj = JuMP.objective_function(m)

    #TODO: modify multiplier term based on objective sense
    #Add equality multiplier terms
    for term in m.ext[:eq_multipliers]  #could have duplicate terms
        coeff = term[1]
        var = term[2]
        multiplier_index = term[3]
        obj += coeff*var*eq_multipliers[multiplier_index]
    end

    #Add inequality multiplier terms
    for term in m.ext[:ineq_multipliers]
        coeff = term[1]
        var = term[2]
        multiplier_index = term[3]
        obj += coeff*var*ineq_multipliers[multiplier_index]
    end

    #Add nonanticpitivity multiplier terms
    for term in m.ext[:linkvar_multipliers]
        coeff = term[1]
        var = term[2]
        multiplier_index = term[3]
        obj += coeff*var*linkvar_multipliers[multiplier_index]
    end

    JuMP.set_objective_function(m,obj)

    JuMP.optimize!(m)

    m.ext[:current_objective_value] = JuMP.objective_value(m)
    return JuMP.objective_value(m)
end

function _prepare_linkvar_multiplier_map!(link_matrix::SparseMatrixCSC{Float64,Int64},link_variables::Vector)
    for i in 1:size(link_matrix)[1] #loop through rows
        row = link_matrix[i,:]
        for j in row.nzind
            coeff = link_matrix[i,j]
            var = link_variables[j]
            model = var.model
            push!(model.ext[:linkvar_multipliers],(coeff,var,i))
        end
    end
end

function _prepare_eq_multiplier_map!(link_matrix::SparseMatrixCSC{Float64,Int64},link_variables::Vector)
    for i in 1:size(link_matrix)[1] #loop through rows
        row = link_matrix[i,:]
        for j in row.nzind
            coeff = link_matrix[i,j]
            var = link_variables[j]
            model = var.model
            push!(model.ext[:eq_multipliers],(coeff,var,i))
        end
    end
end

function _prepare_ineq_multiplier_map!(link_matrix::SparseMatrixCSC{Float64,Int64},link_variables::Vector)
    for i in 1:size(link_matrix)[1] #loop through rows
        row = link_matrix[i,:]
        for j in row.nzind
            coeff = link_matrix[i,j]
            var = link_variables[j]
            model = var.model
            push!(model.ext[:ineq_multipliers],(coeff,var,i))
        end
    end
end

function build_subproblems(graph::ModelGraph)

    linkvariables = getlinkvariables(graph)
    sub_var_map = Dict()
    master_var_map = Dict((k,[]) for k in linkvariables)

    subproblems = []
    for node in getnodes(graph)
        model = getmodel(node)
        model.ext = Dict()          #need to clear out the model.ext in order to copy

        subproblem = copy(model)
        #Attach optimizer from node

        ref_map = ModelGraphs.AggregationMap(subproblem)

        #ref_map.varmap =  Dict(zip(JuMP.all_variables(model),JuMP.all_variables(subproblem)))
        new_vars = JuMP.all_variables(subproblem)
        for var in JuMP.all_variables(model)
            index = var.index
            new_var = new_vars[index.value]
            ref_map[var] = new_vars[index.value]
            if var in keys(node.linkvariablemap)     #Point linked_variables to master problem variables
                push!(master_var_map[node.linkvariablemap[var]],new_vars[index.value])
                ref_map[node.linkvariablemap[var].vref] = new_vars[index.value]
            end
        end

        subproblem.ext[:modelnode] = node

        #Add extra link variables if necessary to submodels.
        for linkvariable in linkvariables
            if !(linkvariable in values(node.linkvariablemap)) #if the child does not have a master variable
                new_var = @variable(subproblem)
                var_name = JuMP.name(linkvariable.vref)
                new_name = var_name
                JuMP.set_name(new_var,new_name)
                if JuMP.start_value(linkvariable.vref) != nothing
                    JuMP.set_start_value(new_x,JuMP.start_value(linkvariable.vref))
                end
                # ref_map[new_var] = linkvariable.vref
                ref_map[linkvariable.vref] = new_var
            else

            end
        end

        constraint_types = JuMP.list_of_constraint_types(graph.mastermodel)
        for (func,set) in constraint_types
            constraint_refs = JuMP.all_constraints(graph.mastermodel, func, set)

            for constraint_ref in constraint_refs
                constraint = JuMP.constraint_object(constraint_ref)

                # #TODO: Override subproblem child variables with master variable constraints
                # if func == JuMP.VariableRef
                #     #if constraint.func in values(master_var_map)  #If this is a child variable
                #         println("here")
                #         for ref in JuMP.all_constraints(subproblem,func,set)
                #             delete(subproblem,ref)
                #         end
                #     #end
                # end
                new_constraint = ModelGraphs._copy_constraint(constraint,ref_map)
                new_ref = JuMP.add_constraint(subproblem,new_constraint)
                ref_map[constraint_ref] = new_ref
            end
        end

        subproblem.ext[:eq_multipliers] = []
        subproblem.ext[:ineq_multipliers] = []
        subproblem.ext[:linkvar_multipliers] = []
        subproblem.ext[:objective_scale] = 1
        obj = JuMP.objective_function(subproblem)
        if JuMP.objective_sense(subproblem) == MOI.MAX_SENSE
            subproblem.ext[:objective_scale] = -1
            JuMP.set_objective_sense(subproblem,MOI.MIN_SENSE)
        end

        subproblem.ext[:original_objective] = obj
        push!(subproblems,subproblem)

        merge!(sub_var_map,ref_map.varmap)

    end
    return subproblems, sub_var_map, master_var_map
end

# #Flip objective, add dictionary entry for multipliers for this model
# function prepare_subproblem!(graph::ModelGraph,m::JuMP.Model)
#     m.ext[:eq_multipliers] = []
#     m.ext[:ineq_multipliers] = []
#     m.ext[:objective_scale] = 1
#     obj = JuMP.objective_function(m)
#     if JuMP.objective_sense(m) == MOI.MAX_SENSE
#         m.ext[:objective_scale] = -1
#         JuMP.set_objective_sense(m,MOI.MIN_SENSE)
#     end
#     m.ext[:original_objective] = obj
#     #Add master constraints corresponding to linked_variable
# end


text:
export TestSuite

"""
Abstract test suite that can be used for all packages inheriting from GPUArray
"""
module TestSuite

using Base.Test

toarray(T, x::Tuple{X, Vararg{Int}}) where X = fill(first(x), Base.tail(x))
toarray(::Type{T}, x::NTuple{N, Int}) where {T <: Bool, N} = rand(T, x)
toarray(::Type{T}, x::NTuple{N, Int}) where {T <: Integer, N} = rand(T(1):T(10), x)
toarray(T, x::NTuple{N, Int}) where N = rand(T, x)
toarray(T, x) = x
togpu(T, x::AbstractArray) = T(x)
togpu(T, x) = x

"""
Calls function `f` on input arrays generated by `sizes` as Base.Array and converted to `Typ`.
Compares the result of `f` and tests if they agree.
`sizes` can be the shape of the array, a value or a tuple `(val, shape...)` which will create
a `fill(val, shape...)`.
"""
function against_base(f, Typ, sizes...)
    jl_arrays = toarray.(eltype(Typ), sizes)
    gpu_arrays = togpu.(Typ, jl_arrays)
    res_jl = f(jl_arrays...)
    res_gpu = f(gpu_arrays...)
    pass = res_jl ≈ Array(res_gpu) # all the array printing in @test A ≈ B is just too bad
    @test pass
end

include("blas.jl")
include("broadcasting.jl")
include("construction.jl")
include("fft.jl")
include("gpuinterface.jl")
include("linalg.jl")
include("mapreduce.jl")
include("base.jl")
include("indexing.jl")
# include("vector.jl")
include("random.jl")

function supported_eltypes()
    (Float32, Float64, Int32, Int64, Complex64, Complex128)
end

"""
Runs the test suite on array type `Typ`
"""
function run_tests(Typ)
    GPUArrays.allowslow(false)
    run_gpuinterface(Typ)
    run_base(Typ)
    run_blas(Typ)
    run_broadcasting(Typ)
    run_construction(Typ)
    run_fft(Typ)
    run_linalg(Typ)
    run_mapreduce(Typ)
    run_indexing(Typ)
    run_random(Typ)
end

export against_base, run_tests, supported_eltypes

end


text:
abstract type AbstractRate end

struct SimpleRate <: AbstractRate end

struct FZeroRate <: AbstractRate end

"""
Calculate rate formula.
"""
function calc_rate(::FZeroRate, su, rel_reserve, turnover, j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma)
    # Find the type of `rate` so that diff and units work with atol etc
    one_r = oneunit(eltype(turnover))
    # Get rate with a zero finder
    let m=m, turnover=turnover, j_E_mai=j_E_mai, y_E_Ea=y_E_Ea, y_E_Eb=y_E_Eb, y_V_E=y_V_E, κsoma=κsoma
        find_zero((-2one_r, 1one_r), Secant(), one_r*1e-10, 200) do r
            rate_formula(r, su, rel_reserve, turnover, j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma)
        end
    end
end

"""
    rate_formula(r, ureserve::NTuple, turnover::NTuple, j_E_mai, y_V_E, κsoma)
Rate formulas for E, CN or CNE reserves
"""
rate_formula(r, su, rel_reserve::NTuple{1}, turnover::NTuple{1},
             j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    (j_E,) = non_growth_flux.(rel_reserve, turnover, r)
    y_V_E * (κsoma * j_E - j_E_mai) - r

    j_E = reserve * (turnover - y_V_E * (κsoma * j_E - j_E_mai))
end
rate_formula(r, su, rel_reserve::NTuple{2}, turnover::NTuple{2},
             j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    j_Ea, j_Eb = non_growth_flux.(rel_reserve, turnover, r)
    j_E = synthesizing_unit(su, j_Ea * y_E_Ea, j_Eb * y_E_Eb)
    y_V_E * (κsoma * j_E - j_E_mai) - r
end
rate_formula(r, su, rel_reserve::NTuple{3}, turnover::NTuple{3},
             j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    j_Ea, j_Eb, j_E = non_growth_flux.(rel_reserve, turnover, r)
    j_E += synthesizing_unit(su, j_Ea * y_E_Ea, j_Eb * y_E_Eb)
    r = y_V_E * (κsoma * j_E - j_E_mai)
end

calc_rate(::SimpleRate, su, rel_reserve::NTuple{2}, turnover::NTuple{2},
          j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    j_Ea, j_Eb = rel_reserve .* turnover
    j_E = synthesizing_unit(j_Ea * y_E_Ea, j_Eb * y_E_Eb)
    j_E - j_E_mai
end
calc_rate(::SimpleRate, su, rel_reserve::NTuple{3}, turnover::NTuple{3},
          j_E_mai, y_E_Ea, y_E_Eb, y_V_E, κsoma) = begin
    j_Ea, j_Eb, j_E = rel_reserve .* turnover
    j_E += synthesizing_unit(su, j_Ea * y_E_Ea, j_Eb * y_E_Eb)
    j_E - j_E_mai
end

"""
    non_growth_flux(ureserve, turnover, r)
Returns the current non_growth_flux flux at rate r,
or the flux as a proportion of u[V], depending on ureserve values.
"""
non_growth_flux(reserve, turnover, r) = reserve * (turnover - r)


text:
using RigidToricGraph
using Test

@testset "RigidToricGraph.jl" begin
    # Write your tests here.

    using Oscar

    A = [1 0 0 0 0 0 1 0;1 0 0 0 0 0 0 1;0 1 0 0 1 0 0 0; 0 1 0 0 0 1 0 0;0 1 0 0 0 0 1 0;0 1 0 0 0 0 0 1;0 0 1 0 1 0 0 0; 0 0 1 0 0 1 0 0;0 0 1 0 0 0 1 0;0 0 1 0 0 0 0 1;0 0 0 1 1 0 0 0; 0 0 0 1 0 1 0 0;0 0 0 1 0 0 1 0;0 0 0 1 0 0 0 1]
    C = positive_hull(A)
    Cdual = polarize(C)
    HD = Cdual.pm_cone.HASSE_DIAGRAM
    @test length(get_twofaces(HD)) == 34
    @test length(get_threefaces(HD, false)) == 68

    def_degree = [1,0,1,0,-2,1,0,3]
    cs = crosscut_skeleton(Cdual, HD, def_degree)
    @test cs == map(x->Set(x), [[2,4,7],[2,4,8],[2,4,9],[2,7,8],[2,8,9],[4,7,8],[4,8,9]])
end


text:
module ShallowWater

export ShallowWaterModel

using StaticArrays
using ...MPIStateArrays: MPIStateArray
using LinearAlgebra: dot, Diagonal
using CLIMAParameters.Planet: grav

using ..Ocean
using ...VariableTemplates
using ...Mesh.Geometry
using ...DGMethods
using ...DGMethods.NumericalFluxes
using ...BalanceLaws
using ..Ocean: kinematic_stress, coriolis_parameter

import ...DGMethods.NumericalFluxes: update_penalty!
import ...BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_conditions,
    boundary_state!
import ..Ocean:
    ocean_init_state!,
    ocean_init_aux!,
    ocean_boundary_state!,
    _ocean_boundary_state!

using ...Mesh.Geometry: LocalGeometry

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
⊗(a::SVector, b::SVector) = a * b'

abstract type TurbulenceClosure end
struct LinearDrag{L} <: TurbulenceClosure
    λ::L
end
struct ConstantViscosity{L} <: TurbulenceClosure
    ν::L
end

"""
    ShallowWaterModel <: BalanceLaw

A `BalanceLaw` for shallow water modeling.

write out the equations here

# Usage

    ShallowWaterModel(problem)

"""
struct ShallowWaterModel{C, PS, P, T, A, FT} <: BalanceLaw
    param_set::PS
    problem::P
    coupling::C
    turbulence::T
    advection::A
    c::FT
    fₒ::FT
    β::FT
    function ShallowWaterModel{FT}(
        param_set::PS,
        problem::P,
        turbulence::T,
        advection::A;
        coupling::C = Uncoupled(),
        c = FT(0), # m/s
        fₒ = FT(1e-4), # Hz
        β = FT(1e-11), # Hz / m
    ) where {FT <: AbstractFloat, PS, P, T, A, C}
        return new{C, PS, P, T, A, FT}(
            param_set,
            problem,
            coupling,
            turbulence,
            advection,
            c,
            fₒ,
            β,
        )
    end
end
SWModel = ShallowWaterModel

function vars_state(m::SWModel, ::Prognostic, T)
    @vars begin
        η::T
        U::SVector{2, T}
    end
end

function init_state_prognostic!(m::SWModel, state::Vars, aux::Vars, localgeo, t)
    ocean_init_state!(m, m.problem, state, aux, localgeo, t)
end

function vars_state(m::SWModel, ::Auxiliary, T)
    @vars begin
        y::T
        Gᵁ::SVector{2, T} # integral of baroclinic tendency
        Δu::SVector{2, T} # reconciliation Δu = 1/H * (Ū - ∫u)
    end
end

function init_state_auxiliary!(
    m::SWModel,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    init_state_auxiliary!(
        m,
        (m, A, tmp, geom) -> ocean_init_aux!(m, m.problem, A, geom),
        state_auxiliary,
        grid,
        direction,
    )
end

function vars_state(m::SWModel, ::Gradient, T)
    @vars begin
        ∇U::SVector{2, T}
    end
end

function compute_gradient_argument!(
    m::SWModel,
    f::Vars,
    q::Vars,
    α::Vars,
    t::Real,
)
    compute_gradient_argument!(m.turbulence, f, q, α, t)
end

compute_gradient_argument!(::LinearDrag, _...) = nothing

@inline function compute_gradient_argument!(
    T::ConstantViscosity,
    f::Vars,
    q::Vars,
    α::Vars,
    t::Real,
)
    f.∇U = q.U

    return nothing
end

function vars_state(m::SWModel, ::GradientFlux, T)
    @vars begin
        ν∇U::SMatrix{3, 2, T, 6}
    end
end

function compute_gradient_flux!(
    m::SWModel,
    σ::Vars,
    δ::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    compute_gradient_flux!(m, m.turbulence, σ, δ, q, α, t)
end

compute_gradient_flux!(::SWModel, ::LinearDrag, _...) = nothing

@inline function compute_gradient_flux!(
    ::SWModel,
    T::ConstantViscosity,
    σ::Vars,
    δ::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    ν = Diagonal(@SVector [T.ν, T.ν, -0])
    ∇U = δ.∇U

    σ.ν∇U = -ν * ∇U

    return nothing
end

@inline function flux_first_order!(
    m::SWModel,
    F::Grad,
    q::Vars,
    α::Vars,
    t::Real,
    direction,
)
    U = @SVector [q.U[1], q.U[2], -0]
    η = q.η
    H = m.problem.H
    Iʰ = @SMatrix [
        1 -0
        -0 1
        -0 -0
    ]

    F.η += U
    F.U += grav(m.param_set) * H * η * Iʰ

    advective_flux!(m, m.advection, F, q, α, t)

    return nothing
end

advective_flux!(::SWModel, ::Nothing, _...) = nothing

@inline function advective_flux!(
    m::SWModel,
    ::NonLinearAdvectionTerm,
    F::Grad,
    q::Vars,
    α::Vars,
    t::Real,
)
    U = q.U
    H = m.problem.H
    V = @SVector [U[1], U[2], -0]

    F.U += 1 / H * V ⊗ U

    return nothing
end

function flux_second_order!(
    m::SWModel,
    G::Grad,
    q::Vars,
    σ::Vars,
    ::Vars,
    α::Vars,
    t::Real,
)
    flux_second_order!(m, m.turbulence, G, q, σ, α, t)
end

flux_second_order!(::SWModel, ::LinearDrag, _...) = nothing

@inline function flux_second_order!(
    ::SWModel,
    ::ConstantViscosity,
    G::Grad,
    q::Vars,
    σ::Vars,
    α::Vars,
    t::Real,
)
    G.U += σ.ν∇U

    return nothing
end

@inline wavespeed(m::SWModel, n⁻, q::Vars, α::Vars, t::Real, direction) = m.c

@inline function source!(
    m::SWModel{P},
    S::Vars,
    q::Vars,
    d::Vars,
    α::Vars,
    t::Real,
    direction,
) where {P}
    # f × u
    U, V = q.U
    f = coriolis_parameter(m, m.problem, α.y)
    S.U -= @SVector [-f * V, f * U]

    forcing_term!(m, m.coupling, S, q, α, t)
    linear_drag!(m.turbulence, S, q, α, t)

    return nothing
end

@inline function forcing_term!(m::SWModel, ::Uncoupled, S, Q, A, t)
    S.U += kinematic_stress(m.problem, A.y)

    return nothing
end

linear_drag!(::ConstantViscosity, _...) = nothing

@inline function linear_drag!(T::LinearDrag, S::Vars, q::Vars, α::Vars, t::Real)
    λ = T.λ
    U = q.U

    S.U -= λ * U

    return nothing
end

boundary_conditions(shallow::SWModel) = shallow.problem.boundary_conditions

"""
    boundary_state!(nf, ::SWModel, args...)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, bc, shallow::SWModel, args...)
    return _ocean_boundary_state!(nf, bc, shallow, args...)
end

"""
    ocean_boundary_state!(nf, bc::OceanBC, ::SWModel)

splits boundary condition application into velocity
"""
@inline function ocean_boundary_state!(nf, bc::OceanBC, m::SWModel, args...)
    return ocean_boundary_state!(nf, bc.velocity, m, m.turbulence, args...)
end

include("bc_velocity.jl")

end

max_forks_repo_name: jm-c/CLIMA


text:
#=
real symmetric or complex Hermitian (svec scaled triangle) domain 𝕊ᵈ₊
=#

mutable struct MatrixEpiPerSepSpectral{D <: PrimDual, C <: RealCompF} <: Cache
    oa_s::Vector{AE}
    h::SepSpectralFun
    d::Int
    w_temp::Vector{RealF}
    W_temp::Matrix{C}
    function MatrixEpiPerSepSpectral{D, C}() where {D <: PrimDual, C <: RealCompF}
        return new{D, C}()
    end
end

function create_sepspectral_cache(
    ::Type{MatrixCSqr{RealF, C}},
    use_dual::Bool,
    d::Int,
    ::Optimizer,
) where {C <: RealCompF}
    D = primal_or_dual(use_dual)
    cache = MatrixEpiPerSepSpectral{D, C}()
    cache.w_temp = zeros(RealF, svec_length(C, d))
    cache.W_temp = zeros(C, d, d)
    return cache
end

function MOIPajarito.Cones.add_init_cuts(
    cache::MatrixEpiPerSepSpectral{D, C},
    opt::Optimizer,
) where {D, C}
    h = cache.h
    d = cache.d
    (u, v) = swap_epiper(D, cache.oa_s[1:2]...)
    @views w = cache.oa_s[3:end]
    W_diag = [w[svec_idx(C, i, i)] for i in 1:d]

    # add variable bounds
    JuMP.@constraint(opt.oa_model, v >= 0)
    if dom_pos(D, h)
        JuMP.@constraint(opt.oa_model, W_diag .>= 0)
    end
    opt.use_init_fixed_oa || return

    # add cuts at p = 1 and R = r₀ I
    r_vals = init_r_vals(D, h)
    for r0 in r_vals
        R_diag = fill(r0, d)
        q = per_sepspec(conj_or_val(D), h, 1.0, R_diag)
        JuMP.@constraint(opt.oa_model, u + q * v + JuMP.dot(R_diag, W_diag) >= 0)
    end
    return
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::MatrixEpiPerSepSpectral{D},
    opt::Optimizer,
) where {D}
    (p, q) = swap_epiper(D, z[1:2]...)
    R = cache.W_temp
    @views svec_to_smat!(R, z[3:end], rt2)
    F = eigen(Hermitian(R, :U))
    V = F.vectors
    ω = F.values
    cuts = AE[]

    if dom_pos(D, cache.h) && abs(q) < 1e-7
        # TODO decide when to add
        # add eigenvector cuts
        num_pos = count(>(1e-8), ω)
        @views V_neg = V[:, 1:num_pos] * Diagonal(sqrt.(ω[1:num_pos]))
        @views w = cache.oa_s[3:end]
        cuts = _get_psd_cuts(V_neg, w, cache, opt)
    end

    # add epigraph cut
    cut = _get_cut(p, ω, V, cache, opt)
    push!(cuts, cut)
    return cuts
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::MatrixEpiPerSepSpectral{D},
    opt::Optimizer,
) where {D}
    h = cache.h
    Ws = cache.W_temp
    @views svec_to_smat!(Ws, s[3:end], rt2)
    F = eigen(Hermitian(Ws, :U))
    V = F.vectors
    ω = F.values
    cuts = AE[]

    if dom_pos(D, h)
        num_neg = count(<(-opt.tol_feas), ω)
        if !iszero(num_neg)
            # add eigenvector cuts
            V_neg = V[:, 1:num_neg]
            @views w = cache.oa_s[3:end]
            cuts = _get_psd_cuts(V_neg, w, cache, opt)
        end

        ω = max.(ω, 1e-10)
    end

    (us, vs) = swap_epiper(D, s[1:2]...)
    v_pos = max(vs, 1e-9)
    if us - per_sepspec(val_or_conj(D), h, v_pos, ω) > -opt.tol_feas
        return cuts
    end

    # gradient cut is (1, h⋆(rω), V * Diagonal(rω) * V') at rω = -h'(ω / vs)
    rω = -h_grad(D, h, ω / v_pos)
    cut = _get_cut(1.0, rω, V, cache, opt)
    push!(cuts, cut)
    return cuts
end

function _get_cut(
    p::RealF,
    rω::Vector{RealF},
    V::Matrix{C},
    cache::MatrixEpiPerSepSpectral{D, C},
    opt::Optimizer,
) where {D, C}
    # strengthened cut is (p, p * h⋆(rω / p), V * Diagonal(rω) * V')
    q = per_sepspec(conj_or_val(D), cache.h, p, rω)
    (p, q) = swap_epiper(D, p, q)

    R = V * Diagonal(rω) * V'
    r = smat_to_svec!(cache.w_temp, R, rt2)

    z = vcat(p, q, r)
    cut = dot_expr(z, cache.oa_s, opt)
    return cut
end


text:
__precompile__()

module Serie

using Makie, GLVisualize, GeometryTypes, Reactive
import GLVisualize: play_slider, mm

include("utils.jl")
include("populations.jl")

function app(pop::Population, iters::Int64=1000)
    scene = Scene()
    control_scene = Scene(scene, Signal(SimpleRectangle(0, 0, 1800, 50)))
    screen = control_scene[:screen]

    play_gui, play_s = play_slider(screen, 8mm, linspace(0, 4, 10^5), slider_length=1500)
    controls = Pair[
        "fertility" => play_gui,
    ]
    _view(visualize(controls), screen, camera=:fixed_pixel)

    ax = linspace(1, iters, iters) # HACK: As Makie loops endlessly when using regular ranges...
    ay = linspace(0, 10, 11) # HACK: As Makie does not support different scales on x and y axis...
    grid = axis(ax, ay)

    points = map(play_s) do var_fert
        pts = Point2f0[]
        p = pop.initial
        for n = 1:iters
            push!(pts, Point2f0(n, 10p))
            f = pop.fertilities[n % end + 1](var_fert)
            p = pop.law(f, p)
        end
        pts
    end

    # pos = lift_node(getindex.(scene, (:mouseposition, :time))...) do mpos, t
    #     map(linspace(0, 2pi, 60)) do i
    #         circle = Point2f0(sin(i), cos(i))
    #         mouse = to_world(Point2f0(mpos), cam)
    #         secondary = (sin((i * 10f0) + t) * 0.09) * normalize(circle)
    #         (secondary .+ circle) .+ mouse
    #     end
    # end

    lines(points)
    center!(scene)
    # _view(visualize(points, :lines), view_screen)
    # renderloop(window)
end

function __init__()
    pop = Population(logistic, [id, cst(3.8), cst(3.7)], .5)
    app(pop, 100)
end

end


text:
import MPI
using FastTransforms: cheb2leg, leg2cheb
using ApproxFun

#using Plots
Np=Int(1e4)
Nx=Int(32)   #Int(ceil(64*pi/2))
k0=0.5
dt=0.1
tmax=15
L=2*pi/k0

print("Asymptotically identical to ", floor(Nx/pi), " Fourier modes\n\n")
function f(x,v)
  return exp(-0.5.*v.^2)/sqrt(2*pi).*(1+0.1.*cos(k0*x))
end

MPI.Init()

comm = MPI.COMM_WORLD
MPI.Barrier(comm)
mpi_rank = MPI.Comm_rank(comm)
mpi_size = MPI.Comm_size(comm)
#print(mpi_size)

Np_loc=Int(floor(Np/mpi_size))

#Seed the random generator
#We would like to skip ahead values
srand( (Np_loc*4)*(mpi_rank+1) )

xk=rand(Np_loc)*L
vk=randn(Np_loc)
wk=f(xk,vk)*L

# Derivative of a chebyshev series
function chebTdx_series{T}(x::T, c::Array{T,1})
  N=length(c)-1;
  b1::T=0
  b2::T=0
  @inbounds for n=N:-1:2
     @inbounds b0::T=c[n+1]+ 2*x*(n+1)/n*b1 - (n+2)/n*b2
    b2=b1; b1=b0;
  end
  return c[2]+ 4*x.*b1 - 3*b2
end

# Derivative of a chebyshev series
function chebTdx_series{T}(x::Array{T,1}, c::Array{T,1})
  N=length(c)-1;
  alpha=Array(T,N-1)
  beta=Array(T,N-1)
  for n=N:-1:2
    alpha[n-1]=2*(n+1)/n
    beta[n-1]=-(n+2)/n
  end
  E=Array(T,length(x))
   @inbounds for pdx=1:length(x)
   b1::T=0
   b2::T=0
   @inbounds for n=N:-1:2
      @inbounds b0::T=c[n+1]+ x[pdx]*alpha[n-1]*b1 + beta[n-1]*b2
     b2=b1; b1=b0;
    end
    E[pdx]=c[2]+ 4*x[pdx].*b1 - 3*b2
  end
  return E
end

# Legendre polynomials
function LegendreP{T}(x::T,N::Int)
  Pn=Array(T, N)
  Pn[1]=1
  Pn[2]=x
  @inbounds for idx=3:N
    n=idx-2
    Pn[idx]= x.*(2*n+1)/(n+1).*Pn[idx-1]- n/(n+1).*Pn[idx-2]
  end
  return Pn
end

#@time LegendreP(xk[1],32)
function  rhs_particle_legendre{T}(xk::Array{T,1},wk::Array{T,1}, Nx::Int)
    const Np=length(xk)
    rho=zeros(T,Nx)
    # Precalculate coefficients
    alpha=Array(T,Nx-2)
    beta=Array(T,Nx-2)
    @inbounds for n=1:Nx-2
      alpha[n]=(2*n+1)/(n+1)
      beta[n]=-n/(n+1)
    end

    @inbounds for pdx=1:Np
      P0::T=1; P1::T=xk[pdx];
      rho[1]+=P0.*wk[pdx]
      rho[2]+=P1.*wk[pdx]
      @inbounds for idx=1:Nx-2
        P2::T=xk[pdx].*alpha[idx].*P1+ beta[idx].*P0
        P0=P1; P1=P2;
        @inbounds rho[idx+2]+=P2.*wk[pdx]
      end
    end
  return rho
end

function  rhs_particle2{T}(xk::Array{T,1},wk::Array{T,1}, Nx::Int)
    const Np=length(xk)
    rho=zeros(T,Nx)
    Pn=Array(T, Nx)
    Pn[1]=1.0

    @inbounds for pdx=1:Np
      Pn[2]=xk[pdx]
      @inbounds for idx=3:Nx
        n=idx-2
        Pn[idx]= Pn[2].*(2*n+1)/(n+1).*Pn[idx-1]- n/(n+1).*Pn[idx-2]
      end
      rho+=Pn.*wk[pdx]
    end
  return rho
end

function  rhs_particle{T}(xk::Array{T,1},wk::Array{T,1}, Nx::Int)
    Np=length(xk)
    rho=zeros(T,Nx)
    @inbounds for pdx=1:Np
      rho+=LegendreP(xk[pdx],Nx).*wk[pdx]
    end
  return rho
end

#using BenchmarkTools
#t=@benchmark rhs=rhs_particle( xk/L*2-1,wk,Nx)/Np*(2/L)
#t2=@benchmark rhs=rhs_particle2( xk/L*2-1,wk,Nx)/Np*(2/L)
#t=@benchmark rhs2=rhs_particle_legendre( xk/L*2-1,wk,Nx)/Np*(2/L)
#@time rhs2=rhs_particle_legendre( xk/L*2-1,wk,Nx)/Np*(2/L)
#@time rhs=rhs_particle2( xk/L*2-1,wk,Nx)/Np*(2/L)

OM=Interval(0,L) # Domain
S=ApproxFun.Space(OM) #Default chebyshev space
B=periodic(S,0) # Periodic boundary
B2=DefiniteIntegral(S) # Nullspace condition
D2=-Derivative(S,2) #Laplace operator
D1=Derivative(S,1) # First derivative for H1 norm
IN=DefiniteIntegral(S) # First derivative for L2 norm
#Phi = \([B;B2;D2],[0.0;0.0;rho];tolerance=1E-16,maxlength=Nx+2)
QR = qrfact([B;B2;D2])

#OP=[B;B2;D2]
#sparse(IN[1:Nx,1:Nx])
#det(sparse(D2[1:Nx,1:Nx]+B2[1:Nx,1:Nx]+B[1:Nx,1:Nx]))
#AbstractMatrix(D2)

# Coefficients for Runge Kutta
const rksd=[2/3.0, -2/3.0, 1  ]
const rksc=[ 7/24.0, 3/4.0, -1/24.0]

Nt=Int(ceil(tmax/dt))
fieldenergy=zeros(Float64,Nt)
kineticenergy=zeros(Float64,Nt)
momentum=zeros(Float64,Nt)
rho=Fun(S)
Phi=Fun(S)
rhs=Array(Float64,Nx)
rhs=rhs_particle_legendre( xk/L*2-1,wk,Nx)/Np*(2/L)
rho=Fun(S,leg2cheb(rhs))
Phi=QR \[0.0;0.0;rho]
fieldenergy[1]=coefficients(IN*(D1*Phi)^2)[1]/2
Phic=-coefficients(Phi)[1:Nx+2]/L*2

#rho=Fun(S,leg2cheb(rand(Nx*64) ))
#using BenchmarkTools
#@benchmark Phi = \([B;B2;D2],[0.0;0.0;rho];tolerance=1E-16, maxlength=Nx+2)
#@benchmark Phi=QR \[0.0;0.0;rho]

MPI.Barrier(comm)
if (mpi_rank==0)
tic()
end
@inbounds for tdx=1:Nt

  @inbounds for rkdx=1:length(rksd)

    # Bin to right hand side
    rhs=rhs_particle_legendre( xk/L*2-1,wk,Nx)/Np*(2/L)
    rhs=MPI.allreduce(rhs, MPI.SUM, comm)

    #rhs=rhs_particle_legendre( xk/L*2-1,wk,Nx)/Np*(2/L)
    # Legendre L^2 projection
    rhs=rhs./(2./(2*(0:Nx-1)+1.0))
    # Remove ion background (quick'n dirty)
    rhs[1]=0
    #Transform to Chebyshev
    rho=Fun(S,leg2cheb(rhs))

    # Solve periodic Poisson equation with Chebyshev
    # Chebyshev periodic Poisson
    #Phi = \([B;B2;D2],[0.0;0.0;rho];tolerance=1E-16, maxlength=Nx+2)
    Phi=QR \[0.0;0.0;rho]
    #,maxlength=Nx+2

    #E=-rho./(-im*kx)
    if rkdx==1 & mpi_rank==0
      fieldenergy[tdx]=coefficients(IN*(D1*Phi)^2)[1]/2
      kineticenergy[tdx]=0.5.*sum(vk.^2.*wk )/Np
      momentum[tdx]=sum(vk.*wk)/Np
    end
    #Phic=-coefficients(Phi)[1:Nx+2]/L*2
    Phic=Phi.coefficients/L*2
    #@inbounds for pdx=1:Np_loc
    #vk[pdx]+=dt*rksd[rkdx]*chebTdx_series(xk[pdx]./(L*2.)-1.0, Phic)
    #vk+=dt*rksd[rkdx].*map(x->chebTdx_series(x, Phic), xk./(L*2.)-1.0)
  #  end
    vk.+=dt*rksd[rkdx]*chebTdx_series(xk/L*2-1.0,Phic )
    xk.+=dt*rksc[rkdx].*vk
    xk=mod(xk,L)
  end
end
MPI.Barrier(comm)

if (mpi_rank==0)
   toc()

#print("test: ", sum(fieldenergy))
using PyPlot
figure()
ttime=(0:Nt-1)*dt
semilogy(ttime,fieldenergy)
grid()
xlabel("time")
ylabel("electrostatic energy")
end
#MPI.Finalize()

# E=vk
# @benchmark chebTdx_series(xk/L*2-1.0,Phic )
#
# @benchmark rhs_particle_legendre( xk/L*2-1,wk,Nx)/Np*(2/L)

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
module Ex5_10c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.10c - Solution of an EFDIP
println("Example 5.10c")

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

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

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

SFDI = [ 0 1 1; 1 0 1; 1 1 0] .> 0; # enter structure matrix

#
Qt, Rft = efdisyn(sysf, SFDI; atol, rdim = 1);

# normalize Q and Rf to match example
scale = sign.([ Rft.sys[1].D[1,2], Rft.sys[2].D[1,3], Rft.sys[3].D[1,1]])
Q = FDIFilter(scale .* Qt.sys, p, mu)
R = FDIFilterIF(scale .* Rft.sys, 0, 0, mf)

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

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

println("Q = ")
display(Q)
println("R = ")
display(R)

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
"""
    runtimeinfo()

Gather information about runtime.  It returns an object that can be
`show`n using `text/plain` and `text/markdown` MIMEs.
"""
runtimeinfo() = RunTimeInfo()

function safe_lscpu()
    try
        return read(`lscpu`, String)
    catch
        return nothing
    end
end

Base.@kwdef struct RunTimeInfo
    blas_num_threads::Union{Int,Nothing} = blas_num_threads()
    blas_vendor::Symbol = LinearAlgebra.BLAS.vendor()
    lscpu::Union{String,Nothing} = safe_lscpu()
end

function Base.show(io::IO, ::MIME"text/plain", info::RunTimeInfo)
    buffer = IOBuffer()
    show(buffer, MIME"text/markdown"(), info)
    seekstart(buffer)
    show(io, MIME"text/plain"(), Markdown.parse(buffer))
end

function Base.show(io::IO, ::MIME"text/markdown", info::RunTimeInfo)
    println(io, "| Runtime Info | |")
    println(io, "|:--|:--|")
    println(io, "| BLAS #threads | ", something(info.blas_num_threads, "unknown"), " |")
    println(io, "| `BLAS.vendor()` | `", info.blas_vendor, "` |")
    println(io, "| `Sys.CPU_THREADS` | ", Sys.CPU_THREADS, " |")
    # Hiding `nthreads` ATM as it can be misleading when it is set via
    # `BenchmarkConfig`:
    # println(io, "| `Threads.nthreads()` | ", Threads.nthreads(), " |")
    if info.lscpu !== nothing
        println(io)
        println(io, "`lscpu` output:")
        println(io)
        for line in split(info.lscpu, "\n")
            println(io, "    ", line)
        end
    end
end

"""
    blas_num_threads() :: Union{Int, Nothing}

Get the number of threads BLAS is using.

Taken from:
https://github.com/JuliaLang/julia/blob/v1.3.0/stdlib/Distributed/test/distributed_exec.jl#L999-L1019

See also: https://stackoverflow.com/a/37516335
"""
function blas_num_threads()
    blas = LinearAlgebra.BLAS.vendor()
    # Wrap in a try to catch unsupported blas versions
    try
        if blas == :openblas
            return ccall((:openblas_get_num_threads, Base.libblas_name), Cint, ())
        elseif blas == :openblas64
            return ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
        elseif blas == :mkl
            return ccall((:MKL_Get_Max_Num_Threads, Base.libblas_name), Cint, ())
        end

        # OSX BLAS looks at an environment variable
        if Sys.isapple()
            return tryparse(Cint, get(ENV, "VECLIB_MAXIMUM_THREADS", "1"))
        end
    catch
    end

    return nothing
end


text:
"""$(TYPEDEF) Trajectory Optimization Problem.
Contains the full definition of a trajectory optimization problem, including:
* dynamics model (`Model`)
* objective (`Objective`)
* constraints (`ConstraintSet`)
* initial and final states
* Primal variables (state and control trajectories)
* Discretization information: knot points (`N`), time step (`dt`), and total time (`tf`)

# Constructors:
```julia
Problem(model, obj, constraints, x0, xf, Z, N, tf) # defaults to RK3 integration
Problem{Q}(model, obj, constraints, x0, xf, Z, N, tf) where Q<:QuadratureRule
Problem(model, obj, xf, tf; x0, constraints, N, X0, U0, dt, integration)
Problem{Q}(prob::Problem)  # change integration
```
where `Z` is a trajectory (Vector of `KnotPoint`s)

# Arguments
* `model`: Dynamics model. Can be either `Discrete` or `Continuous`
* `obj`: Objective
* `X0`: Initial state trajectory. If omitted it will be initialized with NaNs, to be later overwritten by the solver.
* `U0`: Initial control trajectory. If omitted it will be initialized with zeros.
* `x0`: Initial state. Defaults to zeros.
* `xf`: Final state. Defaults to zeros.
* `dt`: Time step
* `tf`: Final time. Set to zero to specify a time penalized problem.
* `N`: Number of knot points. Defaults to 51, unless specified by `dt` and `tf`.
* `integration`: One of the defined integration types to discretize the continuous dynamics model.
Both `X0` and `U0` can be either a `Matrix` or a `Vector{Vector}`, but must be the same.
At least 2 of `dt`, `tf`, and `N` need to be specified (or just 1 of `dt` and `tf`).
"""
struct Problem{Q<:QuadratureRule,T<:AbstractFloat}
    model::AbstractModel
    obj::AbstractObjective
    constraints::ConstraintList
    x0::MVector
    xf::MVector
    Z::Traj
    N::Int
    t0::T
    tf::T
    function Problem{Q}(model::AbstractModel, obj::AbstractObjective,
            constraints::ConstraintList,
            x0::StaticVector, xf::StaticVector,
            Z::Traj, N::Int, t0::T, tf::T) where {Q,T}
        n,m = size(model)
        @assert length(x0) == length(xf) == n
        @assert length(Z) == N
        @assert tf > t0
        new{Q,T}(model, obj, constraints, x0, xf, Z, N, t0, tf)
    end
end

"Use RK3 as default integration"
Problem(model, obj, constraints, x0, xf, Z, N, t0, tf) =
    Problem{RobotDynamics.RK3}(model, obj, constraints, x0, xf, Z, N, t0, tf)

function Problem(model::L, obj::O, xf::AbstractVector, tf;
        constraints=ConstraintList(size(model)...,length(obj)),
        t0=zero(tf),
        x0=zero(xf), N::Int=length(obj),
        X0=[x0*NaN for k = 1:N],
        U0=[@SVector zeros(size(model)[2]) for k = 1:N-1],
        dt=fill((tf-t0)/(N-1),N-1),
        integration=DEFAULT_Q) where {L,O}
    n,m = size(model)
    if dt isa Real
        dt = fill(dt,N)
    end
	@assert sum(dt[1:N-1]) ≈ tf "Time steps are inconsistent with final time"
    if X0 isa AbstractMatrix
        X0 = [X0[:,k] for k = 1:size(X0,2)]
    end
    if U0 isa AbstractMatrix
        U0 = [U0[:,k] for k = 1:size(U0,2)]
    end
    t = pushfirst!(cumsum(dt), 0)
    Z = Traj(X0,U0,dt,t)

    Problem{integration}(model, obj, constraints, SVector{n}(x0), SVector{n}(xf),
        Z, N, t0, tf)
end

"$(TYPEDSIGNATURES)
Get number of states, controls, and knot points"
Base.size(prob::Problem) = size(prob.model)..., prob.N

"""```julia
integration(::Problem)
integration(::DynamicsConstraint)
```
Get the integration rule"""
integration(prob::Problem{Q}) where Q = Q

"```julia
controls(::Problem)
controls(::Traj)
```
Get the control trajectory
"
controls(prob::Problem) = controls(prob.Z)

"```julia
states(::Problem)
states(::Traj)
```
Get the state trajectory
"
states(prob::Problem) = states(prob.Z)

"""
	get_times(::Problem)

Get the times for all the knot points in the problem.
"""
@inline RobotDynamics.get_times(prob::Problem) = get_times(get_trajectory(prob))

"""
	initial_trajectory!(prob::Problem, Z)

Copy the trajectory
"""
function initial_trajectory!(prob, Z0::AbstractTrajectory)
	Z = get_trajectory(prob)
    for k = 1:prob.N
        Z[k].z = Z0[k].z
    end
end

"""
	initial_states!(::Problem, X0::Vector{<:AbstractVector})
	initial_states!(::Problem, X0::AbstractMatrix)

Copy the state trajectory
"""
@inline initial_states!(prob, X0) = RobotDynamics.set_states!(get_trajectory(prob), X0)

"""
	set_initial_state!(prob::Problem, x0::AbstractVector)

Set the initial state in `prob` to `x0`
"""
function set_initial_state!(prob::Problem, x0::AbstractVector)
    prob.x0 .= x0
end

"""
    set_initial_time!(prob, t0)

Set the initial time of the optimization problem, shifting the time of all points in the trajectory.
Returns the updated final time.
"""
function set_initial_time!(prob, t0::Real)
    Z = get_trajectory(prob)
    Δt = t0 - Z[1].t
    for k in eachindex(Z)
        Z[k].t += Δt
    end
    return Z[end].t
end

function set_goal_state!(prob::Problem, xf::AbstractVector; objective=true, constraint=true)
    if objective
        obj = get_objective(prob)
        for k in eachindex(obj.cost)
            set_LQR_goal!(obj[k], xf)
        end
    end
    if constraint
        for con in get_constraints(prob)
            if con isa GoalConstraint
                set_goal_state!(con, xf)
            end
        end
    end
    copyto!(prob.xf, xf)
    return nothing
end

"""
	initial_controls!(::Problem, U0::Vector{<:AbstractVector})
	initial_controls!(::Problem, U0::AbstractMatrx)

Copy the control trajectory
"""
@inline initial_controls!(prob, U0) = RobotDynamics.set_controls!(get_trajectory(prob), U0)

"```julia
cost(::Problem)
cost(::AbstractSolver)
```
Compute the cost for the current trajectory"
@inline cost(prob::Problem) = cost(prob.obj, prob.Z)

"Copy the problem"
function copy(prob::Problem{Q}) where Q
    Problem{Q}(prob.model, copy(prob.obj), copy(prob.constraints), prob.x0, prob.xf,
        copy(prob.Z), prob.N, prob.t0, prob.tf)
end

function max_violation(prob::Problem, Z::Traj=prob.Z)
    conSet = get_constraints(prob)
    evaluate!(conSet, Z)
    max_violation!(conSet)
    return maximum(conSet.c_max)
end

num_constraints(prob::Problem) = get_constraints(prob).p

@inline get_constraints(prob::Problem) = prob.constraints
@inline get_model(prob::Problem) = prob.model
@inline get_objective(prob::Problem) = prob.obj
@inline get_trajectory(prob::Problem) = prob.Z
@inline is_constrained(prob) = isempty(get_constraints(prob))
@inline get_initial_state(prob::Problem) = prob.x0

states(x) = states(get_trajectory(x))
controls(x) = controls(get_trajectory(x))

"```julia
change_integration(prob::Problem, Q<:QuadratureRule)
```
Change dynamics integration for the problem"
change_integration(prob::Problem, ::Type{Q}) where Q<:QuadratureRule =
    Problem{Q}(prob)

function Problem{Q}(p::Problem) where Q
    Problem{Q}(p.model, p.obj, p.constraints, p.x0, p.xf, p.Z, p.N, p.t0, p.tf)
end

"""
	rollout!(::Problem)

Simulate the dynamics forward from the initial condition `x0` using the controls in the
trajectory `Z`.
If a problem is passed in, `Z = prob.Z`, `model = prob.model`, and `x0 = prob.x0`.
"""
@inline rollout!(prob::Problem{Q}) where {Q} = rollout!(Q, get_model(prob), get_trajectory(prob), get_initial_state(prob))

function Problem(p::Problem; model=p.model, obj=p.obj, constraints=p.constraints,
    x0=p.x0, xf=p.xf, t0=p.t0, tf=p.tf)
    Problem(model, obj, constraints, x0, xf, p.Z, p.N, t0, tf)
end

"```julia
add_dynamics_constraints!(prob::Problem)
```
Add dynamics constraints to the constraint set"
function add_dynamics_constraints!(prob::Problem{Q}, integration=Q, idx=-1) where Q
	n,m = size(prob)
    conSet = prob.constraints

    # Implicit dynamics
    dyn_con = DynamicsConstraint{integration}(prob.model, prob.N)
    add_constraint!(conSet, dyn_con, 1:prob.N-1, idx) # add it at the end

    # Initial condition
    init_con = GoalConstraint(prob.x0)
    add_constraint!(conSet, init_con, 1, 1)  # add it at the top

    return nothing
end


text:
# Accessing a.L directly might involve an extra copy();
# instead, always use the stored Cholesky factor:
chol_lower(a::Cholesky) = a.uplo === 'L' ? a.L : a.U'
chol_upper(a::Cholesky) = a.uplo === 'U' ? a.U : a.L'

# For a dense Matrix, the following allows us to avoid the Adjoint wrapper:
chol_lower(a::Matrix) = cholesky(Hermitian(a, :L)).L

if HAVE_CHOLMOD
    CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}

    chol_lower(cf::CholTypeSparse) = cf.L
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Compute the Greenwich Mean Sideral Time (GMST).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] http://www.navipedia.net/index.php/CEP_to_ITRF, accessed 2015-12-01.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export j2000_to_gmst, jd_to_gmst

"""
    j2000_to_gmst(J2000_UT1::Number)

Compute the Greenwich Mean Sideral Time (GMST) \\[rad] given the instant
`J2000_UT1` in J2000.0 reference [UT1].

!!! info
    The algorithm is based in **[2]**.

# References

- **[2]** http://www.navipedia.net/index.php/CEP_to_ITRF, accessed 2015-12-01.
"""
function j2000_to_gmst(J2000_UT1::Number)
	# Julian centuries elapsed from the epoch J2000.0.
	T_UT1 = J2000_UT1/36525

	# Greenwich Mean Sideral Time at T_UT1 [s].
    θ_GMST = @evalpoly(
        T_UT1,
        + 67310.54841,
        + (876600.0 * 3600 + 8640184.812866),
        + 0.093104,
        - 6.2e-6
    )

    # Reduce to the interval [0, 86400]s.
    θ_GMST = mod(θ_GMST, 86400)

    # Convert to radian and return.
    return θ_GMST * π / 43200
end

"""
    jd_to_gmst(JD_UT1::Number)

Compute the Greenwich Mean Sideral Time (GMST) \\[rad] for the Julian Day
`JD_UT1` [UT1].

!!! info
    The algorithm is based in **[1]**(p. 188).

# References

- **[1]** Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
jd_to_gmst(JD_UT1::Number) = j2000_to_gmst(JD_UT1 - JD_J2000)


text:
"""
    FairTensor{C}

Fairness Tensor with C classes. It consists of C 2 x 2 matrices stacked up to form a Matrix
of size C x 2 x 2. Each 2 x 2 matrix contains values [[TP, FP], [FN, TN]].
"""
struct FairTensor{C}
    mat # No type has been specified here due to the LinProgWrapper postprocessing algorithm
    labels::Vector{String}
end

"""
    FairTensor(m, labels)

Instantiates a FairTensor using the matrix m.
"""
function FairTensor(m, labels::Vector{String})
    s = size(m)
    (s[3] == s[2] && s[3] == 2) || throw(ArgumentError("Expected a C*2*2 type Matrix."))
    length(labels) == s[1] ||
        throw(ArgumentError("As many labels as classes must be provided."))
    FairTensor{s[1]}(m, labels)
end

# allow to access ft[i,j]
Base.getindex(ft::FairTensor, inds...) = getindex(ft.mat, inds...)
Base.setindex!(ft::FairTensor, value, inds...) = setindex!(ft.mat, value, inds...)

"""
    fair_tensor(ŷ, y, grp)

Computes the fairness tensor, where ŷ are the predicted classes,
y are the ground truth values, grp are the group values.
The ordering follows that of `levels(y)`.

Note that ŷ, y and grp are all categorical arrays
"""
function fair_tensor(ŷ::Vec{<:CategoricalElement}, y::Vec{<:CategoricalElement},
                          grp::Vec{<:CategoricalElement})

    check_dimensions(ŷ, y)
    check_dimensions(ŷ, grp)
    length(levels(y))==2 || throw(ArgumentError("Binary Targets are only supported"))
    labels = levels(y)
    favLabel = labels[2]
    unfavLabel = labels[1]

    levels_ = levels(grp)
    c = length(levels_)
    # Dictionary data-structure is used now to map group labels and the corresponding index.
    # Other alternative could be binary search on levels_ everytime. But it would be slow by log(length(levels_)).
    grp_idx = Dict()
    for i in 1:c
        grp_idx[levels_[i]] = i
    end

    # Coverting Categorical Vector to Bool Vector.
    # TODO: Can think of adding another dispatch where user directly passes Bool Vec
    y = y.==favLabel
    ŷ = ŷ.==favLabel
    n = length(y)

    fact = zeros(Int, c, 2, 2)
    @inbounds for i in 1:n
        if ŷ[i] && y[i]
            fact[grp_idx[grp[i]], 1, 1] += 1
        elseif ŷ[i] && !y[i]
            fact[grp_idx[grp[i]], 1, 2] += 1
        elseif !ŷ[i] && y[i]
            fact[grp_idx[grp[i]], 2, 1] += 1
        elseif !y[i] && !y[i]
            fact[grp_idx[grp[i]], 2, 2] += 1
        end
    end
    return FairTensor(fact, string.(levels_))

end

# synonym
fact = fair_tensor

# aggregation:
Base.round(m::FairTensor; kws...) = m
function Base.:+(t1::FairTensor, t2::FairTensor)
    if t1.labels != t2.labels
        throw(ArgumentError("Tensor labels must agree"))
    end
    FairTensor(t1.mat + t2.mat, t1.labels)
end


text:
@testset "errors" begin
# Unexpected character in array
@test_throws ArgumentError JSON3.read("[1,2,3/4,5,6,7]")
# Unexpected character in object
@test_throws ArgumentError JSON3.read("{\"1\":2, \"2\":3 _ \"4\":5}")
# Invalid escaped character
@test_throws ArgumentError repr(JSON3.read("[\"alpha\\α\"]"))
# Invalid 'simple' and 'unknown value'
@test_throws ArgumentError JSON3.read("[tXXe]")
@test_throws ArgumentError JSON3.read("[fail]")
@test_throws ArgumentError JSON3.read("∞")
# Invalid number
@test_throws ArgumentError JSON3.read("[5,2,-]")
@test_throws ArgumentError JSON3.read("[5,2,+β]")
# Incomplete escape
@test_throws ArgumentError JSON3.read("\"\\")
@test_throws ArgumentError JSON3.read("[\"🍕\"_\"🍕\"")

end # @testset "errors"

@enum Animal zebra aardvark horse
@test JSON3.write(zebra) == "\"zebra\""
@test JSON3.write([aardvark, horse, Dict("z" => zebra)]) == "[\"aardvark\",\"horse\",{\"z\":\"zebra\"}]"

@testset "Symbol" begin
    symtest = Dict(:symbolarray => [:apple, :pear], :symbolsingleton => :hello)
    @test (JSON3.write(symtest) == "{\"symbolarray\":[\"apple\",\"pear\"],\"symbolsingleton\":\"hello\"}"
             || JSON3.write(symtest) == "{\"symbolsingleton\":\"hello\",\"symbolarray\":[\"apple\",\"pear\"]}")
end

@testset "Floats" begin
    @test_throws ErrorException JSON3.write([NaN])
    @test_throws ErrorException JSON3.write([Inf])
    @test JSON3.write([NaN], allow_inf=true) == "[NaN]"
    @test JSON3.write([Inf], allow_inf=true) == "[Infinity]"
    @test JSON3.write([-Inf], allow_inf=true) == "[-Infinity]"
    @test JSON3.read("[Inf]", Vector{Float64}) == [Inf]
    @test JSON3.read("[Infinity]", Vector{Float64}) == [Inf]
    @test JSON3.read("[NaN]", Vector{Float64})[1] === NaN
    @test_throws ArgumentError JSON3.read("NaN")
    @test_throws ArgumentError JSON3.read("Inf")
    @test_throws ArgumentError JSON3.read("Infinity")
    @test_skip(@test_throws ArgumentError JSON3.read("-Infinity"))
    @test JSON3.read("NaN"; allow_inf=true) === NaN
    @test JSON3.read("Inf"; allow_inf=true) === Inf
    @test JSON3.read("Infinity"; allow_inf=true) === Inf
    @test JSON3.read("-Infinity"; allow_inf=true) === -Inf
end

@testset "Char" begin
    @test JSON3.write('a') == "\"a\""
    @test JSON3.write('\\') == "\"\\\\\""
    @test JSON3.write('\n') == "\"\\n\""
    @test JSON3.write('🍩') =="\"🍩\""
end

@test sprint(JSON3.write, Float64) == "\"Float64\""

@testset "Null bytes" begin
    zeros = Dict("\0" => "\0")
    json_zeros = JSON3.write(zeros)
    @test occursin("\\u0000", json_zeros)
    @test !occursin("\\0", json_zeros)
    @test JSON3.read(json_zeros, Dict) == zeros
end

@testset "Arrays" begin
# Printing an empty array or Dict shouldn't cause a BoundsError
@test JSON3.write(String[]) == "[]"
@test JSON3.write(Dict()) == "{}"

#Multidimensional arrays
@test_broken JSON3.write([0 1; 2 0]) == "[[0,2],[1,0]]" #TODO
end

@testset "Pairs" begin
@test JSON3.write(1 => 2) == "{\"1\":2}"
@test JSON3.write(:foo => 2) == "{\"foo\":2}"
@test JSON3.write([1, 2] => [3, 4]) == "{\"$([1, 2])\":[3,4]}"
@test JSON3.write([1 => 2]) == "[{\"1\":2}]"
end

@testset "Sets" begin
@test JSON3.write(Set()) == "[]"
@test JSON3.write(Set([1, 2])) in ["[1,2]", "[2,1]"]
end

#Examples from http://json.org/example.html
JSON3.read("{\"menu\": {
    \"id\": \"file\",
    \"value\": \"File\",
    \"popup\": {
      \"menuitem\": [
        {\"value\": \"New\", \"onclick\": \"CreateNewDoc()\"},
        {\"value\": \"Open\", \"onclick\": \"OpenDoc()\"},
        {\"value\": \"Close\", \"onclick\": \"CloseDoc()\"}
      ]
    }
  }}
  ")

JSON3.read("{
\"glossary\": {
   \"title\": \"example glossary\",
\"GlossDiv\": {
       \"title\": \"S\",
 \"GlossList\": {
           \"GlossEntry\": {
               \"ID\": \"SGML\",
     \"SortAs\": \"SGML\",
     \"GlossTerm\": \"Standard Generalized Markup Language\",
     \"Acronym\": \"SGML\",
     \"Abbrev\": \"ISO 8879:1986\",
     \"GlossDef\": {
                   \"para\": \"A meta-markup language, used to create markup languages such as DocBook.\",
       \"GlossSeeAlso\": [\"GML\", \"XML\"]
               },
     \"GlossSee\": \"markup\"
           }
       }
   }
}
}
")

JSON3.read("""
{"widget": {
"debug": "on",
"window": {
   "title": "Sample Konfabulator Widget",
   "name": "main_window",
   "width": 500,
   "height": 500
},
"image": {
   "src": "Images/Sun.png",
   "name": "sun1",
   "hOffset": 250,
   "vOffset": 250,
   "alignment": "center"
},
"text": {
   "data": "Click Here",
   "size": 36.5,
   "style": "bold",
   "name": "text1",
   "hOffset": 250,
   "vOffset": 100,
   "alignment": "center",
   "onMouseUp": "sun1.opacity = (sun1.opacity / 100) * 90;"
}
}}""")

JSON3.read("{\"web-app\": {
\"servlet\": [
{
 \"servlet-name\": \"cofaxCDS\",
 \"servlet-class\": \"org.cofax.cds.CDSServlet\",
 \"init-param\": {
   \"configGlossary:installationAt\": \"Philadelphia, PA\",
   \"configGlossary:adminEmail\": \"ksm@pobox.com\",
   \"configGlossary:poweredBy\": \"Cofax\",
   \"configGlossary:poweredByIcon\": \"/images/cofax.gif\",
   \"configGlossary:staticPath\": \"/content/static\",
   \"templateProcessorClass\": \"org.cofax.WysiwygTemplate\",
   \"templateLoaderClass\": \"org.cofax.FilesTemplateLoader\",
   \"templatePath\": \"templates\",
   \"templateOverridePath\": \"\",
   \"defaultListTemplate\": \"listTemplate.htm\",
   \"defaultFileTemplate\": \"articleTemplate.htm\",
   \"useJSP\": false,
   \"jspListTemplate\": \"listTemplate.jsp\",
   \"jspFileTemplate\": \"articleTemplate.jsp\",
   \"cachePackageTagsTrack\": 200,
   \"cachePackageTagsStore\": 200,
   \"cachePackageTagsRefresh\": 60,
   \"cacheTemplatesTrack\": 100,
   \"cacheTemplatesStore\": 50,
   \"cacheTemplatesRefresh\": 15,
   \"cachePagesTrack\": 200,
   \"cachePagesStore\": 100,
   \"cachePagesRefresh\": 10,
   \"cachePagesDirtyRead\": 10,
   \"searchEngineListTemplate\": \"forSearchEnginesList.htm\",
   \"searchEngineFileTemplate\": \"forSearchEngines.htm\",
   \"searchEngineRobotsDb\": \"WEB-INF/robots.db\",
   \"useDataStore\": true,
   \"dataStoreClass\": \"org.cofax.SqlDataStore\",
   \"redirectionClass\": \"org.cofax.SqlRedirection\",
   \"dataStoreName\": \"cofax\",
   \"dataStoreDriver\": \"com.microsoft.jdbc.sqlserver.SQLServerDriver\",
   \"dataStoreUrl\": \"jdbc:microsoft:sqlserver://LOCALHOST:1433;DatabaseName=goon\",
   \"dataStoreUser\": \"sa\",
   \"dataStorePassword\": \"dataStoreTestQuery\",
   \"dataStoreTestQuery\": \"SET NOCOUNT ON;select test='test';\",
   \"dataStoreLogFile\": \"/usr/local/tomcat/logs/datastore.log\",
   \"dataStoreInitConns\": 10,
   \"dataStoreMaxConns\": 100,
   \"dataStoreConnUsageLimit\": 100,
   \"dataStoreLogLevel\": \"debug\",
   \"maxUrlLength\": 500}},
{
 \"servlet-name\": \"cofaxEmail\",
 \"servlet-class\": \"org.cofax.cds.EmailServlet\",
 \"init-param\": {
 \"mailHost\": \"mail1\",
 \"mailHostOverride\": \"mail2\"}},
{
 \"servlet-name\": \"cofaxAdmin\",
 \"servlet-class\": \"org.cofax.cds.AdminServlet\"},

{
 \"servlet-name\": \"fileServlet\",
 \"servlet-class\": \"org.cofax.cds.FileServlet\"},
{
 \"servlet-name\": \"cofaxTools\",
 \"servlet-class\": \"org.cofax.cms.CofaxToolsServlet\",
 \"init-param\": {
   \"templatePath\": \"toolstemplates/\",
   \"log\": 1,
   \"logLocation\": \"/usr/local/tomcat/logs/CofaxTools.log\",
   \"logMaxSize\": \"\",
   \"dataLog\": 1,
   \"dataLogLocation\": \"/usr/local/tomcat/logs/dataLog.log\",
   \"dataLogMaxSize\": \"\",
   \"removePageCache\": \"/content/admin/remove?cache=pages&id=\",
   \"removeTemplateCache\": \"/content/admin/remove?cache=templates&id=\",
   \"fileTransferFolder\": \"/usr/local/tomcat/webapps/content/fileTransferFolder\",
   \"lookInContext\": 1,
   \"adminGroupID\": 4,
   \"betaServer\": true}}],
\"servlet-mapping\": {
\"cofaxCDS\": \"/\",
\"cofaxEmail\": \"/cofaxutil/aemail/*\",
\"cofaxAdmin\": \"/admin/*\",
\"fileServlet\": \"/static/*\",
\"cofaxTools\": \"/tools/*\"},

\"taglib\": {
\"taglib-uri\": \"cofax.tld\",
\"taglib-location\": \"/WEB-INF/tlds/cofax.tld\"}}}")

JSON3.read("""
{"menu": {
"header": "SVG\\tViewer\\u03b1",
"items": [
   {"id": "Open"},
   {"id": "OpenNew", "label": "Open New"},
   null,
   {"id": "ZoomIn", "label": "Zoom In"},
   {"id": "ZoomOut", "label": "Zoom Out"},
   {"id": "OriginalView", "label": "Original View"},
   null,
   {"id": "Quality"},
   {"id": "Pause"},
   {"id": "Mute"},
   null,
   {"id": "Find", "label": "Find..."},
   {"id": "FindAgain", "label": "Find Again"},
   {"id": "Copy"},
   {"id": "CopyAgain", "label": "Copy Again"},
   {"id": "CopySVG", "label": "Copy SVG"},
   {"id": "ViewSVG", "label": "View SVG"},
   {"id": "ViewSource", "label": "View Source"},
   {"id": "SaveAs", "label": "Save As"},
   null,
   {"id": "Help"},
   {"id": "About", "label": "About Adobe SVG Viewer..."}
]
}}""")

#Example JSON strings from http://www.jquery4u.com/json/10-example-json-files/
JSON3.read("{\"markers\": [
   {
       \"point\":\"new GLatLng(40.266044,-74.718479)\",
       \"homeTeam\":\"Lawrence Library\",
       \"awayTeam\":\"LUGip\",
       \"markerImage\":\"images/red.png\",
       \"information\": \"Linux users group meets second Wednesday of each month.\",
       \"fixture\":\"Wednesday 7pm\",
       \"capacity\":\"\",
       \"previousScore\":\"\"
   },
   {
       \"point\":\"new GLatLng(40.211600,-74.695702)\",
       \"homeTeam\":\"Hamilton Library\",
       \"awayTeam\":\"LUGip HW SIG\",
       \"markerImage\":\"images/white.png\",
       \"information\": \"Linux users can meet the first Tuesday of the month to work out harward and configuration issues.\",
       \"fixture\":\"Tuesday 7pm\",
       \"capacity\":\"\",
       \"tv\":\"\"
   },
   {
       \"point\":\"new GLatLng(40.294535,-74.682012)\",
       \"homeTeam\":\"Applebees\",
       \"awayTeam\":\"After LUPip Mtg Spot\",
       \"markerImage\":\"images/newcastle.png\",
       \"information\": \"Some of us go there after the main LUGip meeting, drink brews, and talk.\",
       \"fixture\":\"Wednesday whenever\",
       \"capacity\":\"2 to 4 pints\",
       \"tv\":\"\"
   }
] }")

JSON3.read("{
\"colorsArray\":[{
       \"colorName\":\"red\",
       \"hexValue\":\"#f00\"
   },
   {
       \"colorName\":\"green\",
       \"hexValue\":\"#0f0\"
   },
   {
       \"colorName\":\"blue\",
       \"hexValue\":\"#00f\"
   },
   {
       \"colorName\":\"cyan\",
       \"hexValue\":\"#0ff\"
   },
   {
       \"colorName\":\"magenta\",
       \"hexValue\":\"#f0f\"
   },
   {
       \"colorName\":\"yellow\",
       \"hexValue\":\"#ff0\"
   },
   {
       \"colorName\":\"black\",
       \"hexValue\":\"#000\"
   }
]
}")

JSON3.read("{
\"colorsArray\":[{
       \"red\":\"#f00\",
       \"green\":\"#0f0\",
       \"blue\":\"#00f\",
       \"cyan\":\"#0ff\",
       \"magenta\":\"#f0f\",
       \"yellow\":\"#ff0\",
       \"black\":\"#000\"
   }
]
}")

JSON3.read("{
\"red\":\"#f00\",
\"green\":\"#0f0\",
\"blue\":\"#00f\",
\"cyan\":\"#0ff\",
\"magenta\":\"#f0f\",
\"yellow\":\"#ff0\",
\"black\":\"#000\"
}")

JSON3.read("{\"results\":[

{\"text\":\"@twitterapi  http://tinyurl.com/ctrefg\",
\"to_user_id\":396524,
\"to_user\":\"TwitterAPI\",
\"from_user\":\"jkoum\",
\"metadata\":
    {
     \"result_type\":\"popular\",
     \"recent_retweets\": 109
    },
\"id\":1478555574,
\"from_user_id\":1833773,
\"iso_language_code\":\"nl\",
\"source\":\"twitter\",
\"profile_image_url\":\"http://s3.amazonaws.com/twitter_production/profile_images/118412707/2522215727_a5f07da155_b_normal.jpg\",
\"created_at\":\"Wed, 08 Apr 2009 19:22:10 +0000\"}],
\"since_id\":0,
\"max_id\":1480307926,
\"refresh_url\":\"?since_id=1480307926&q=%40twitterapi\",
\"results_per_page\":15,
\"next_page\":\"?page=2&max_id=1480307926&q=%40twitterapi\",
\"completed_in\":0.031704,
\"page\":1,
\"query\":\"%40twitterapi\"}")

JSON3.read("{
\"data\": [
 {
    \"id\": \"X999_Y999\",
    \"from\": {
       \"name\": \"Tom Brady\", \"id\": \"X12\"
    },
    \"message\": \"Looking forward to 2010!\",
    \"actions\": [
       {
          \"name\": \"Comment\",
          \"link\": \"http://www.facebook.com/X999/posts/Y999\"
       },
       {
          \"name\": \"Like\",
          \"link\": \"http://www.facebook.com/X999/posts/Y999\"
       }
    ],
    \"type\": \"status\",
    \"created_time\": \"2010-08-02T21:27:44+0000\",
    \"updated_time\": \"2010-08-02T21:27:44+0000\"
 },
 {
    \"id\": \"X998_Y998\",
    \"from\": {
       \"name\": \"Peyton Manning\", \"id\": \"X18\"
    },
    \"message\": \"Where's my contract?\",
    \"actions\": [
       {
          \"name\": \"Comment\",
          \"link\": \"http://www.facebook.com/X998/posts/Y998\"
       },
       {
          \"name\": \"Like\",
          \"link\": \"http://www.facebook.com/X998/posts/Y998\"
       }
    ],
    \"type\": \"status\",
    \"created_time\": \"2010-08-02T21:27:44+0000\",
    \"updated_time\": \"2010-08-02T21:27:44+0000\"
 }
]
}")

JSON3.read("""{
"title": "Talk On Travel Pool",
"link": "http://www.flickr.com/groups/talkontravel/pool/",
"description": "Travel and vacation photos from around the world.",
"modified": "2009-02-02T11:10:27Z",
"generator": "http://www.flickr.com/",
"totalItems":222,
"items": [
       {
       "title": "View from the hotel",
       "link": "http://www.flickr.com/photos/33112458@N08/3081564649/in/pool-998875@N22",
       "media": {"m":"http://farm4.static.flickr.com/3037/3081564649_4a6569750c_m.jpg"},
       "date_taken": "2008-12-04T04:43:03-08:00",
       "description": "
 Talk On Travel has added a photo to the pool:

 ",
       "published": "2008-12-04T12:43:03Z",
       "author": "nobody@flickr.com (Talk On Travel)",
       "author_id": "33112458@N08",
       "tags": "spain dolphins tenerife canaries lagomera aqualand playadelasamericas junglepark losgigantos loscristines talkontravel"
       }
]
}""")

JSON3.read("{\"apiVersion\":\"2.0\",
\"data\":{
\"updated\":\"2010-01-07T19:58:42.949Z\",
\"totalItems\":800,
\"startIndex\":1,
\"itemsPerPage\":1,
\"items\":[
   {\"id\":\"hYB0mn5zh2c\",
    \"uploaded\":\"2007-06-05T22:07:03.000Z\",
    \"updated\":\"2010-01-07T13:26:50.000Z\",
    \"uploader\":\"GoogleDeveloperDay\",
    \"category\":\"News\",
    \"title\":\"Google Developers Day US - Maps API Introduction\",
    \"description\":\"Google Maps API Introduction ...\",
    \"tags\":[
       \"GDD07\",\"GDD07US\",\"Maps\"
    ],
    \"thumbnail\":{
       \"default\":\"http://i.ytimg.com/vi/hYB0mn5zh2c/default.jpg\",
       \"hqDefault\":\"http://i.ytimg.com/vi/hYB0mn5zh2c/hqdefault.jpg\"
    },
    \"player\":{
       \"default\":\"http://www.youtube.com/watch?v\u003dhYB0mn5zh2c\"
    },
    \"content\":{
       \"1\":\"rtsp://v5.cache3.c.youtube.com/CiILENy.../0/0/0/video.3gp\",
       \"5\":\"http://www.youtube.com/v/hYB0mn5zh2c?f...\",
       \"6\":\"rtsp://v1.cache1.c.youtube.com/CiILENy.../0/0/0/video.3gp\"
    },
    \"duration\":2840,
    \"aspectRatio\":\"widescreen\",
    \"rating\":4.63,
    \"ratingCount\":68,
    \"viewCount\":220101,
    \"favoriteCount\":201,
    \"commentCount\":22,
    \"status\":{
       \"value\":\"restricted\",
       \"reason\":\"limitedSyndication\"
    },
    \"accessControl\":{
       \"syndicate\":\"allowed\",
       \"commentVote\":\"allowed\",
       \"rate\":\"allowed\",
       \"list\":\"allowed\",
       \"comment\":\"allowed\",
       \"embed\":\"allowed\",
       \"videoRespond\":\"moderated\"
    }
   }
]
}
}")

JSON3.read("{
\"menu\": {
   \"header\": \"xProgress SVG Viewer\",
   \"items\": [
       {
           \"id\": \"Open\"
       },
       {
           \"id\": \"OpenNew\",
           \"label\": \"Open New\"
       },
       null,
       {
           \"id\": \"ZoomIn\",
           \"label\": \"Zoom In\"
       },
       {
           \"id\": \"ZoomOut\",
           \"label\": \"Zoom Out\"
       },
       {
           \"id\": \"OriginalView\",
           \"label\": \"Original View\"
       },
       null,
       {
           \"id\": \"Quality\"
       },
       {
           \"id\": \"Pause\"
       },
       {
           \"id\": \"Mute\"
       },
       null,
       {
           \"id\": \"Find\",
           \"label\": \"Find...\"
       },
       {
           \"id\": \"FindAgain\",
           \"label\": \"Find Again\"
       },
       {
           \"id\": \"Copy\"
       },
       {
           \"id\": \"CopyAgain\",
           \"label\": \"Copy Again\"
       },
       {
           \"id\": \"CopySVG\",
           \"label\": \"Copy SVG\"
       },
       {
           \"id\": \"ViewSVG\",
           \"label\": \"View SVG\"
       },
       {
           \"id\": \"ViewSource\",
           \"label\": \"View Source\"
       },
       {
           \"id\": \"SaveAs\",
           \"label\": \"Save As\"
       },
       null,
       {
           \"id\": \"Help\"
       },
       {
           \"id\": \"About\",
           \"label\": \"About xProgress CVG Viewer...\"
       }
   ]
}
}")

JSON3.read("{
\"firstName\": \"John\",
\"lastName\": \"Smith\",
\"age\": 25,
\"address\":
{
    \"streetAddress\": \"21 2nd Street\",
    \"city\": \"New York\",
    \"state\": \"NY\",
    \"postalCode\": \"10021\"
},
\"phoneNumber\":
[
    {
      \"type\": \"home\",
      \"number\": \"212 555-1234\"
    },
    {
      \"type\": \"fax\",
      \"number\": \"646 555-4567\"
    }
]
}")

JSON3.read("{
   \"name\":\"Product\",
   \"properties\":
   {
           \"id\":
           {
                   \"type\":\"number\",
                   \"description\":\"Product identifier\",
                   \"required\":true
           },
           \"name\":
           {
                   \"description\":\"Name of the product\",
                   \"type\":\"string\",
                   \"required\":true
           },
           \"price\":
           {
                   \"type\":\"number\",
                   \"minimum\":0,
                   \"required\":true
           },
           \"tags\":
           {
                   \"type\":\"array\",
                   \"items\":
                   {
                           \"type\":\"string\"
                   }
           }
   }
}")

JSON3.read("{
\"ResultSet\": {
   \"totalResultsAvailable\": \"1827221\",
   \"totalResultsReturned\": 2,
   \"firstResultPosition\": 1,
   \"Result\": [
       {
           \"Title\": \"potato jpg\",
           \"Summary\": \"Kentang Si bungsu dari keluarga Solanum tuberosum L ini ternyata memiliki khasiat untuk mengurangi kerutan  jerawat  bintik hitam dan kemerahan pada kulit  Gunakan seminggu sekali sebagai\",
           \"Url\": \"http://www.mediaindonesia.com/spaw/uploads/images/potato.jpg\",
           \"ClickUrl\": \"http://www.mediaindonesia.com/spaw/uploads/images/potato.jpg\",
           \"RefererUrl\": \"http://www.mediaindonesia.com/mediaperempuan/index.php?ar_id=Nzkw\",
           \"FileSize\": 22630,
           \"FileFormat\": \"jpeg\",
           \"Height\": \"362\",
           \"Width\": \"532\",
           \"Thumbnail\": {
               \"Url\": \"http://thm-a01.yimg.com/nimage/557094559c18f16a\",
               \"Height\": \"98\",
               \"Width\": \"145\"
           }
       },
       {
           \"Title\": \"potato jpg\",
           \"Summary\": \"Introduction of puneri aloo This is a traditional potato preparation flavoured with curry leaves and peanuts and can be eaten on fasting day  Preparation time   10 min\",
           \"Url\": \"http://www.infovisual.info/01/photo/potato.jpg\",
           \"ClickUrl\": \"http://www.infovisual.info/01/photo/potato.jpg\",
           \"RefererUrl\": \"http://sundayfood.com/puneri-aloo-indian-%20recipe\",
           \"FileSize\": 119398,
           \"FileFormat\": \"jpeg\",
           \"Height\": \"685\",
           \"Width\": \"1024\",
           \"Thumbnail\": {
               \"Url\": \"http://thm-a01.yimg.com/nimage/7fa23212efe84b64\",
               \"Height\": \"107\",
               \"Width\": \"160\"
           }
       }
   ]
}
}")

JSON3.read("""
{"অলিম্পিকস": {
"অ্যাথলেট": "২২টি দেশ থেকে ২,০৩৫ জন প্রতিযোগী",
"ইভেন্ট": "২২টি ইভেন্টের মধ্যে ছিল দড়ি টানাটানি",
"রেকর্ড": [
   {"১০০মি. স্প্রিন্ট": "রেজি ওয়াকার, দক্ষিণ আফ্রিকা"},
   {"Marathon": "জনি হেইস"},
   {" ফ্রি-স্টাইল সাঁতার": "Henry Taylor, Britain"}
]
}}
""")

# issue #5
issue5 = "[\"A\",\"B\",\"C\\n\"]"
@test JSON3.read(issue5) == ["A", "B", "C\n"]

# $ escaping issue
dollars = ["all of the \$s", "µniçø∂\$"]
@test JSON3.read(JSON3.write(dollars)) == dollars

# unmatched brackets
brackets = Dict("foo"=>"ba}r", "be}e]p"=>"boo{p")
@test JSON3.read(JSON3.write(brackets), Dict{String, String}) == brackets

test21 = "[\r\n{\r\n\"a\": 1,\r\n\"b\": 2\r\n},\r\n{\r\n\"a\": 3,\r\n\"b\": 4\r\n}\r\n]"
a = JSON3.read(test21, Vector{Any})
@test isa(a, Vector{Any})
@test length(a) == 2

obj = JSON3.read("{\"a\":2e10}")
@test obj.a == 2e10

obj = JSON3.read("{\"\U0001d712\":\"\\ud835\\udf12\"}")
@test obj.𝜒 == "𝜒"

struct t109
    i::Int
end
JSON3.StructType(::Type{t109}) = JSON3.Struct()

let iob = IOBuffer()
    JSON3.write(iob, t109(1))
    @test (JSON3.read(String(take!(iob)))).i == 1
end

# issue #127 support AbstractCmd
if !Sys.iswindows()
@test (JSON3.read(`echo \{\"i\":1\}`)).i == 1
@test (JSON3.read(`echo` & `echo \{\"i\":1\}`)).i == 1
end

@test_broken JSON3.write([Int64[] Int64[]]) == "[[],[]]" #TODO
@test JSON3.write([Int64[] Int64[]]') == "[]"

@test Float32(JSON3.read(JSON3.write(2.1f-8))) == 2.1f-8

# PR 60
@test JSON3.write(String([0xff])) == "\"\xff\""

# PR 116
@test mktemp() do path, io
    JSON3.write(path, [1, 2, 3])
    open(JSON3.read, path) == [1, 2, 3]
end


text:
@testset "utils" begin
    import BayesianOptimization: IterationCounter, DurationCounter, step!,
    init!, isdone
    it = IterationCounter(0, 0, 10)
    for _ in 1:10
        step!(it)
    end
    @test isdone(it) == true
    maxiterations!(it, 40)
    @test it.N == 40
    init!(it)
    @test it.c == 0
    @test it.i == 10

    now = time()
    it = DurationCounter(now, .2, now, now + .2)
    sleep(.21)
    @test isdone(it) == true
    maxduration!(it, .5)
    init!(it)
    @test isdone(it) == false
    sleep(.51)
    @test isdone(it) == true
end