text:
using Documenter, Makie
cd(@__DIR__)
include("../examples/library.jl")
include("documenter_extension.jl")
import AbstractPlotting: _help, to_string, to_func, to_type

pathroot = Pkg.dir("Makie")
docspath = Pkg.dir("Makie", "docs")
srcpath = joinpath(pathroot, "docs", "src")
srcmediapath = joinpath(pathroot, "docs", "media")
buildpath = joinpath(pathroot, "docs", "build")
mediapath = joinpath(pathroot, "docs", "build", "media")
expdbpath = joinpath(buildpath, "examples-database.html")
# TODO can we teach this to documenter somehow?
ispath(mediapath) || mkpath(mediapath)
plotting_functions = (atomics..., contour, arrows, barplot, poly)

function output_path(entry, ending; subdir = nothing)
    if subdir == nothing
        joinpath(mediapath, string(entry.unique_name, ending))
    else
        joinpath(mediapath, subdir, string(entry.unique_name, ending))
    end
end

function save_example(entry, x::Scene)
    path = output_path(entry, ".jpg")
    Makie.save(path, x)
    path
end

save_example(entry, x::String) = x # nothing to do

function save_example(entry, x::Makie.Stepper) #TODO: this breaks thumbnail generation
    # return a list of all file names
    path = [output_path(entry, "-$i.jpg"; subdir = string(entry.unique_name)) for i = 1:x.step - 1]
    return path
end

AbstractPlotting.set_theme!(resolution = (500, 500))
eval_examples(outputfile = output_path) do example, value
    AbstractPlotting.set_theme!(resolution = (500, 500))
    Random.seed!(42)
    path = save_example(example, value)
    if isa(value, Makie.Stepper)
        name = [string.("thumb-", example.unique_name, "-$i", ".jpg") for i = 1:value.step - 1]
    else
        name = string("thumb-", example.unique_name, ".jpg")
    end
    try
        generate_thumbnail.(path, joinpath.(dirname.(path), name))
    catch e
        @warn("generate_thumbnail failed with path $path, entry $(example.unique_name), and filename $name")
        Base.showerror(stderr, e)
        println(stderr)
        Base.show_backtrace(stderr, Base.catch_backtrace())
        println(stderr)
    end
end

# =============================================
# automatically generate an overview of the atomic functions, using a source md file
@info("Generating functions overview")
path = joinpath(srcpath, "functions-overview.md")
srcdocpath = joinpath(srcpath, "src-functions.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    medialist = readdir(mediapath)
    isempty(medialist) && error("media folder is empty -- perhaps you forgot to generate the plots? :)")
    println(io, "# Atomic functions overview")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    for func in plotting_functions
        fname = to_string(func)
        expdbpath = joinpath(buildpath, "examples-$fname.html")
        println(io, "## `$fname`\n")
        try
            println(io, "```@docs")
            println(io, "$fname")
            println(io, "```\n")
            help_attributes(io, func; extended = true)
            embed_thumbnail_link(io, func, buildpath, expdbpath)
        catch e
            println("ERROR: Didn't work with $fname\n")
            Base.showerror(stderr, e)
        end
        println(io, "\n")
    end
end

# =============================================
# automatically generate gallery based on looping through the database
# using pre-generated plots from generate_plots.jl
cd(docspath)
example_pages = nothing
example_list = String[]
for func in plotting_functions
    fname = to_string(func)
    info("Generating examples gallery for $fname")
    path = joinpath(srcpath, "examples-$fname.md")
    indices = find_indices(func)
    open(path, "w") do io
        println(io, "# `$fname`")
        examples2source(fname, scope_start = "", scope_end = "", indent = "") do entry, source
            # print bibliographic stuff
            println(io, "## $(entry.title)")
            print(io, "Tags: ")
            tags = sort(collect(entry.tags))
            for j = 1:length(tags) - 1; print(io, "`$(tags[j])`, "); end
            println(io, "`$(tags[end])`.\n")
            uname = string(entry.unique_name)
            src_lines = entry.file_range
            println(io, """
                ```julia
                $source
                ```
                """
            )
            embed_plot(io, uname, mediapath, buildpath; src_lines = src_lines)
        end
    end
    push!(example_list, "examples-$fname.md")
end

# =============================================
# automatically generate an overview of the plot attributes (keyword arguments), using a source md file
@info("Generating attributes page")
include("../src/plot_attr_desc.jl")
path = joinpath(srcpath, "plot-attributes.md")
srcdocpath = joinpath(srcpath, "src-plot-attributes.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plot attributes")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    println(io, "## List of attributes")
    print_table(io, plot_attr_desc)
end

# =============================================
# automatically generate an overview of the axis attributes, using a source md file
@info("Generating axis page")
path = joinpath(srcpath, "axis.md")
srcdocpath = joinpath(srcpath, "src-axis.md")
include("../src/Axis2D_attr_desc.jl")
include("../src/Axis3D_attr_desc.jl")

open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    src = read(srcdocpath, String)
    println(io, src)
    print(io)
    # Axis2D section
    println(io, "## `Axis2D`")
    println(io, "### `Axis2D` attributes groups")
    print_table(io, Axis2D_attr_desc)
    print(io)
    for (k, v) in Axis2D_attr_groups
        println(io, "#### `:$k`\n")
        print_table(io, v)
        println(io)
    end
    # Axis3D section
    println(io, "## `Axis3D`")
    println(io, "### `Axis3D` attributes groups")
    print_table(io, Axis3D_attr_desc)
    print(io)
    for (k, v) in Axis3D_attr_groups
        println(io, "#### `:$k`\n")
        print_table(io, v)
        println(io)
    end
end

# automatically generate an overview of the function signatures, using a source md file
@info("Generating signatures page")
path = joinpath(srcpath, "signatures.md")
srcdocpath = joinpath(srcpath, "src-signatures.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plot function signatures")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    println(io, "```@docs")
    println(io, "convert_arguments")
    println(io, "```\n")

    println(io, "See [Plot attributes](@ref) for the available plot attributes.")
end

# build docs with Documenter
@info("Running `makedocs` with Documenter.")
makedocs(
    modules = [Makie, AbstractPlotting],
    doctest = false, clean = false,
    format = :html,
    sitename = "Makie.jl",
    pages = Any[
        "Home" => "index.md",
        "Basics" => [
            "basic-tutorials.md",
            "help_functions.md",
            "functions-overview.md",
            "signatures.md",
            "plot-attributes.md",
            # "extending.md",
            "axis.md",
            "interaction.md",
            "output.md",
            # "layout.md"
        ],
        # atomics_pages,
        "Examples" => [
            "index-examples.md",
            example_list...
        ],
        "Developer Documentation" => [
            "why-makie.md",
        #     "devdocs.md",
        ],
    ]
)

# deploy
ENV["DOCUMENTER_DEBUG"] = "true"
if !haskey(ENV, "DOCUMENTER_KEY")
    # Workaround for when deploying locally and silly Windows truncating the env variable
    # on the CI these should be set!
    ENV["TRAVIS_BRANCH"] = "latest"
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_REPO_SLUG"] = "github.com/JuliaPlots/Makie.jl.git"
    ENV["TRAVIS_TAG"] = "v1.0.0"
    ENV["TRAVIS_OS_NAME"] = ""
    ENV["TRAVIS_JULIA_VERSION"] = ""
    ENV["PATH"] = string(ENV["PATH"], raw";C:\Python27\Scripts")
    ENV["DOCUMENTER_KEY"] = open(x->String(read(x)), joinpath(homedir(), "documenter.key"))
end

deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-cinder"),
    repo = "github.com/JuliaPlots/Makie.jl.git",
    julia = "",
    osname = "",
    latest = "v1.0.0",
    target = "build",
    make = nothing
)


text:
mutable struct Term
    courses::Array{Course}      # Array of courses for the term
    totalEnrolled::Int          # Total number of students enrolled
    failures::Int               # Number of failures within the term

    function Term(courses::Array{Course})
        this = new()

        this.courses = courses
        this.totalEnrolled = 0
        this.failures = 0

        return this
    end
end

ext: jl

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
export eigs, entropy, polar, random_orthog, random_unitary, Spectrum, svd, truncerror

#
# Linear Algebra of order 2 NDTensors
#
# Even though DenseTensor{_,2} is strided
# and passable to BLAS/LAPACK, it cannot
# be made <: StridedArray

function (
  T1::Tensor{ElT1,2,StoreT1} * T2::Tensor{ElT2,2,StoreT2}
) where {ElT1,StoreT1<:Dense,IndsT1,ElT2,StoreT2<:Dense,IndsT2}
  RM = matrix(T1) * matrix(T2)
  indsR = (ind(T1, 1), ind(T2, 2))
  return tensor(Dense(vec(RM)), indsR)
end

function LinearAlgebra.exp(T::DenseTensor{ElT,2}) where {ElT<:Union{Real,Complex}}
  expTM = exp(matrix(T))
  return tensor(Dense(vec(expTM)), inds(T))
end

function LinearAlgebra.exp(
  T::Hermitian{ElT,<:DenseTensor{ElT,2}}
) where {ElT<:Union{Real,Complex}}
  # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric,
  # so extract the parent matrix
  expTM = parent(exp(matrix(T)))
  return tensor(Dense(vec(expTM)), inds(T))
end

"""
  Spectrum
contains the (truncated) density matrix eigenvalue spectrum which is computed during a
decomposition done by `svd` or `eigen`. In addition stores the truncation error.
"""
struct Spectrum{VecT<:Union{AbstractVector,Nothing}}
  eigs::VecT
  truncerr::Float64
end

eigs(s::Spectrum) = s.eigs
truncerror(s::Spectrum) = s.truncerr

function entropy(s::Spectrum)
  S = 0.0
  eigs_s = eigs(s)
  isnothing(eigs_s) &&
    error("Spectrum does not contain any eigenvalues, cannot compute the entropy")
  for p in eigs_s
    p > 1e-13 && (S -= p * log(p))
  end
  return S
end

function svd_catch_error(A; kwargs...)
  USV = try
    svd(A; kwargs...)
  catch
    return nothing
  end
  return USV
end

function lapack_svd_error_message(alg)
  return "The SVD algorithm `\"$alg\"` has thrown an error,\n" *
         "likely because of a convergance failure. You can try\n" *
         "other SVD algorithms that may converge better using the\n" *
         "`alg` (or `svd_alg` if called through `factorize` or MPS/MPO functionality) keyword argument:\n\n" *
         " - \"divide_and_conquer\" is a divide-and-conquer algorithm\n" *
         "   (LAPACK's `gesdd`). It is fast, but may lead to some innacurate\n" *
         "   singular values for very ill-conditioned matrices.\n" *
         "   It also may sometimes fail to converge, leading to errors\n" *
         "   (in which case `\"qr_iteration\"` or `\"recursive\"` can be tried).\n\n" *
         " - `\"qr_iteration\"` (LAPACK's `gesvd`) is typically slower \n" *
         "   than \"divide_and_conquer\", especially for large matrices,\n" *
         "   but is more accurate for very ill-conditioned matrices \n" *
         "   compared to `\"divide_and_conquer\"`.\n\n" *
         " - `\"recursive\"` is ITensor's custom SVD algorithm. It is very\n" *
         "   reliable, but may be slow if high precision is needed.\n" *
         "   To get an `svd` of a matrix `A`, an eigendecomposition of\n" *
         "   ``A^{\\dagger} A`` is used to compute `U` and then a `qr` of\n" *
         "   ``A^{\\dagger} U`` is used to compute `V`. This is performed\n" *
         "   recursively to compute small singular values.\n\n" *
         "Returning `nothing`. For an output `F = svd(A, ...)` you can check if\n" *
         "`isnothing(F)` in your code and try a different algorithm.\n\n" *
         "To suppress this message in the future, you can wrap the `svd` call in the\n" *
         "`@suppress` macro from the `Suppressor` package.\n"
end

"""
    svd(T::DenseTensor{<:Number,2}; kwargs...)

svd of an order-2 DenseTensor
"""
function LinearAlgebra.svd(T::DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)

  #
  # Keyword argument deprecations
  #
  use_absolute_cutoff = false
  if haskey(kwargs, :absoluteCutoff)
    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
  end

  use_relative_cutoff = true
  if haskey(kwargs, :doRelCutoff)
    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
  end

  if haskey(kwargs, :fastsvd) || haskey(kwargs, :fastSVD)
    error(
      "In svd, fastsvd/fastSVD keyword arguments are removed in favor of alg, see documentation for more details.",
    )
  end

  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
  mindim::Int = get(kwargs, :mindim, 1)
  cutoff::Float64 = get(kwargs, :cutoff, 0.0)
  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)
  alg::String = get(kwargs, :alg, "divide_and_conquer")

  #@timeit_debug timer "dense svd" begin
  if alg == "divide_and_conquer"
    MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.DivideAndConquer())
  elseif alg == "qr_iteration"
    MUSV = svd_catch_error(matrix(T); alg=LinearAlgebra.QRIteration())
  elseif alg == "recursive"
    MUSV = svd_recursive(matrix(T))
  else
    error(
      "svd algorithm $alg is not currently supported. Please see the documentation for currently supported algorithms.",
    )
  end
  if isnothing(MUSV)
    if any(isnan, T)
      println("SVD failed, the matrix you were trying to SVD contains NaNs.")
    else
      println(lapack_svd_error_message(alg))
    end
    return nothing
  end
  MU, MS, MV = MUSV
  conj!(MV)
  #end # @timeit_debug

  P = MS .^ 2
  if truncate
    truncerr, _ = truncate!(
      P;
      mindim=mindim,
      maxdim=maxdim,
      cutoff=cutoff,
      use_absolute_cutoff=use_absolute_cutoff,
      use_relative_cutoff=use_relative_cutoff,
    )
  else
    truncerr = 0.0
  end
  spec = Spectrum(P, truncerr)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:, 1:dS]
    resize!(MS, dS)
    MV = MV[:, 1:dS]
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dS)
  v = eltype(IndsT)(dS)
  Uinds = IndsT((ind(T, 1), u))
  Sinds = IndsT((u, v))
  Vinds = IndsT((ind(T, 2), v))
  U = tensor(Dense(vec(MU)), Uinds)
  S = tensor(Diag(MS), Sinds)
  V = tensor(Dense(vec(MV)), Vinds)
  return U, S, V, spec
end

function LinearAlgebra.eigen(
  T::Hermitian{ElT,<:DenseTensor{ElT,2,IndsT}}; kwargs...
) where {ElT<:Union{Real,Complex},IndsT}
  # Keyword argument deprecations
  use_absolute_cutoff = false
  if haskey(kwargs, :absoluteCutoff)
    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
  end
  use_relative_cutoff = true
  if haskey(kwargs, :doRelCutoff)
    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
  end

  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
  mindim::Int = get(kwargs, :mindim, 1)
  cutoff::Float64 = get(kwargs, :cutoff, 0.0)
  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)

  DM, VM = eigen(matrix(T))

  # Sort by largest to smallest eigenvalues
  p = sortperm(DM; rev=true, by=abs)
  DM = DM[p]
  VM = VM[:, p]

  if truncate
    truncerr, _ = truncate!(
      DM;
      maxdim=maxdim,
      cutoff=cutoff,
      use_absolute_cutoff=use_absolute_cutoff,
      use_relative_cutoff=use_relative_cutoff,
    )
    dD = length(DM)
    if dD < size(VM, 2)
      VM = VM[:, 1:dD]
    end
  else
    dD = length(DM)
    truncerr = 0.0
  end
  spec = Spectrum(DM, truncerr)

  # Make the new indices to go onto V
  l = eltype(IndsT)(dD)
  r = eltype(IndsT)(dD)
  Vinds = IndsT((dag(ind(T, 2)), dag(r)))
  Dinds = IndsT((l, dag(r)))
  V = tensor(Dense(vec(VM)), Vinds)
  D = tensor(Diag(DM), Dinds)
  return D, V, spec
end

"""
    random_unitary(n::Int,m::Int)::Matrix{ComplexF64}
    random_unitary(::Type{ElT},n::Int,m::Int)::Matrix{ElT}

Return a random matrix U of dimensions (n,m)
such that if n >= m, U'*U is the identity, or if
m > n U*U' is the identity. Optionally can pass a numeric
type as the first argument to obtain a matrix of that type.

Sampling is based on https://arxiv.org/abs/math-ph/0609050
such that in the case `n==m`, the unitary matrix will be sampled
according to the Haar measure.
"""
function random_unitary(::Type{ElT}, n::Int, m::Int) where {ElT<:Number}
  if n < m
    return Matrix(random_unitary(ElT, m, n)')
  end
  F = qr(randn(ElT, n, m))
  Q = Matrix(F.Q)
  # The upper triangle of F.factors
  # are the elements of R.
  # Multiply cols of Q by the signs
  # that would make diagonal of R
  # non-negative:
  for c in 1:size(Q, 2)
    Q[:, c] .*= sign(F.factors[c, c])
  end
  return Q
end

random_unitary(n::Int, m::Int) = random_unitary(ComplexF64, n, m)

"""
    random_orthog(n::Int,m::Int)::Matrix{Float64}
    random_orthog(::Type{ElT},n::Int,m::Int)::Matrix{ElT}

Return a random, real matrix O of dimensions (n,m)
such that if n >= m, transpose(O)*O is the
identity, or if m > n O*transpose(O) is the
identity. Optionally can pass a real number type
as the first argument to obtain a matrix of that type.
"""
random_orthog(::Type{ElT}, n::Int, m::Int) where {ElT<:Real} = random_unitary(ElT, n, m)

random_orthog(n::Int, m::Int) = random_orthog(Float64, n, m)

"""
    qr_positive(M::AbstractMatrix)

Compute the QR decomposition of a matrix M
such that the diagonal elements of R are
non-negative. Such a QR decomposition of a
matrix is unique. Returns a tuple (Q,R).
"""
function qr_positive(M::AbstractMatrix)
  sparseQ, R = qr(M)
  Q = convert(Matrix, sparseQ)
  nc = size(Q, 2)
  for c in 1:nc
    if real(R[c, c]) < 0.0
      R[c, c:end] *= -1
      Q[:, c] *= -1
    end
  end
  return (Q, R)
end

function LinearAlgebra.eigen(
  T::DenseTensor{ElT,2,IndsT}; kwargs...
) where {ElT<:Union{Real,Complex},IndsT}
  # Keyword argument deprecations
  use_absolute_cutoff = false
  if haskey(kwargs, :absoluteCutoff)
    @warn "In svd, keyword argument absoluteCutoff is deprecated in favor of use_absolute_cutoff"
    use_absolute_cutoff = get(kwargs, :absoluteCutoff, use_absolute_cutoff)
  end
  use_relative_cutoff = true
  if haskey(kwargs, :doRelCutoff)
    @warn "In svd, keyword argument doRelCutoff is deprecated in favor of use_relative_cutoff"
    use_relative_cutoff = get(kwargs, :doRelCutoff, use_relative_cutoff)
  end

  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
  maxdim::Int = get(kwargs, :maxdim, minimum(dims(T)))
  mindim::Int = get(kwargs, :mindim, 1)
  cutoff::Float64 = get(kwargs, :cutoff, 0.0)
  use_absolute_cutoff::Bool = get(kwargs, :use_absolute_cutoff, use_absolute_cutoff)
  use_relative_cutoff::Bool = get(kwargs, :use_relative_cutoff, use_relative_cutoff)

  DM, VM = eigen(matrix(T))

  # Sort by largest to smallest eigenvalues
  #p = sortperm(DM; rev = true)
  #DM = DM[p]
  #VM = VM[:,p]

  if truncate
    truncerr, _ = truncate!(
      DM;
      maxdim=maxdim,
      cutoff=cutoff,
      use_absolute_cutoff=use_absolute_cutoff,
      use_relative_cutoff=use_relative_cutoff,
    )
    dD = length(DM)
    if dD < size(VM, 2)
      VM = VM[:, 1:dD]
    end
  else
    dD = length(DM)
    truncerr = 0.0
  end
  spec = Spectrum(abs.(DM), truncerr)

  i1, i2 = inds(T)

  # Make the new indices to go onto D and V
  l = typeof(i1)(dD)
  r = dag(sim(l))
  Dinds = (l, r)
  Vinds = (dag(i2), r)
  D = complex(tensor(Diag(DM), Dinds))
  V = complex(tensor(Dense(vec(VM)), Vinds))
  return D, V, spec
end

function LinearAlgebra.qr(T::DenseTensor{ElT,2,IndsT}; kwargs...) where {ElT,IndsT}
  positive = get(kwargs, :positive, false)
  # TODO: just call qr on T directly (make sure
  # that is fast)
  if positive
    QM, RM = qr_positive(matrix(T))
  else
    QM, RM = qr(matrix(T))
  end
  # Make the new indices to go onto Q and R
  q, r = inds(T)
  q = dim(q) < dim(r) ? sim(q) : sim(r)
  Qinds = IndsT((ind(T, 1), q))
  Rinds = IndsT((q, ind(T, 2)))
  Q = tensor(Dense(vec(Matrix(QM))), Qinds)
  R = tensor(Dense(vec(RM)), Rinds)
  return Q, R
end

# TODO: support alg keyword argument to choose the svd algorithm
function polar(T::DenseTensor{ElT,2,IndsT}) where {ElT,IndsT}
  QM, RM = polar(matrix(T))
  dim = size(QM, 2)
  # Make the new indices to go onto Q and R
  q = eltype(IndsT)(dim)
  # TODO: use push/pushfirst instead of a constructor
  # call here
  Qinds = IndsT((ind(T, 1), q))
  Rinds = IndsT((q, ind(T, 2)))
  Q = tensor(Dense(vec(QM)), Qinds)
  R = tensor(Dense(vec(RM)), Rinds)
  return Q, R
end


text:
# TODO remove this type piracy when we require MOI v0.9 as we need this fix: https://github.com/JuliaOpt/MathOptInterface.jl/pull/800
function Base.:*(f::MOI.Utilities.ScalarLike, p::AbstractPolynomialLike)
    return MultivariatePolynomials.multconstant(f, p)
end
function Base.:*(p::AbstractPolynomialLike, f::MOI.Utilities.ScalarLike)
    return MultivariatePolynomials.multconstant(p, f)
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

if !isdefined(@__MODULE__(), :NamespaceA) @__MODULE__().eval(:(module NamespaceA __precompile__(false); import FlatBuffers end)) end
if !isdefined(NamespaceA, :NamespaceB) NamespaceA.eval(:(module NamespaceB __precompile__(false); import FlatBuffers end)) end
include("NamespaceA/NamespaceB/EnumInNestedNS.jl")
include("NamespaceA/NamespaceB/TableInNestedNS.jl")
include("NamespaceA/NamespaceB/StructInNestedNS.jl")

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
using Test, OrdinaryDiffEq, DiffEqProblemLibrary, DiffEqCallbacks, LinearAlgebra
import LinearAlgebra: norm
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

# save_everystep, scalar problem
prob = prob_ode_linear
saved_values = SavedValues(Float64, Float64)
cb = SavingCallback((u,t,integrator)->u, saved_values)
sol = solve(prob, Tsit5(), callback=cb)
print("\n", saved_values, "\n")
@test all(idx -> sol.t[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> sol.u[idx] == saved_values.saveval[idx], eachindex(saved_values.t))

# save_everystep, inplace problem
prob2D = prob_ode_2Dlinear
saved_values = SavedValues(eltype(prob2D.tspan), typeof(prob2D.u0))
cb = SavingCallback((u,t,integrator)->copy(u), saved_values)
sol = solve(prob2D, Tsit5(), callback=cb)
@test all(idx -> sol.t[idx] .== saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> all(sol.u[idx] .== saved_values.saveval[idx]), eachindex(saved_values.t))

saved_values = SavedValues(eltype(prob2D.tspan), eltype(prob2D.u0))
cb = SavingCallback((u,t,integrator)->u[1], saved_values)
sol = solve(prob2D, Tsit5(), callback=cb)
@test all(idx -> sol.t[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> sol.u[idx][1] == saved_values.saveval[idx], eachindex(saved_values.t))

# saveat, scalar problem
saved_values = SavedValues(Float64, Float64)
saveat = range(prob.tspan[1], stop=prob.tspan[2], length=50)
cb = SavingCallback((u,t,integrator)->u, saved_values, saveat=saveat)
sol = solve(prob, Tsit5(), callback=cb)
@test all(idx -> saveat[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> abs(sol(saveat[idx]) - saved_values.saveval[idx]) < 8.e-15, eachindex(saved_values.t))

# saveat, inplace problem
saved_values = SavedValues(eltype(prob2D.tspan), typeof(prob2D.u0))
saveat = range(prob2D.tspan[1], stop=prob.tspan[2], length=50)
cb = SavingCallback((u,t,integrator)->copy(u), saved_values, saveat=saveat)
sol = solve(prob2D, Tsit5(), callback=cb)
@test all(idx -> saveat[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> norm(sol(saveat[idx]) - saved_values.saveval[idx]) < 8.e-15, eachindex(saved_values.t))

saved_values = SavedValues(eltype(prob2D.tspan), eltype(prob2D.u0))
saveat = range(prob2D.tspan[1], stop=prob.tspan[2], length=50)
cb = SavingCallback((u,t,integrator)->u[1], saved_values, saveat=saveat)
sol = solve(prob2D, Tsit5(), callback=cb)
@test all(idx -> saveat[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> abs(sol(saveat[idx])[1] - saved_values.saveval[idx]) < 8.e-15, eachindex(saved_values.t))

# saveat, tdir<0, scalar problem
prob_inverse = ODEProblem(prob.f, prob.u0, (prob.tspan[end], prob.tspan[1]),1.01)
saved_values = SavedValues(Float64, Float64)
saveat = range(prob_inverse.tspan[1], stop=prob.tspan[2], length=50)
cb = SavingCallback((u,t,integrator)->u, saved_values, saveat=saveat, tdir=-1)
sol = solve(prob_inverse, Tsit5(), callback=cb)
@test all(idx -> saveat[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> abs(sol(saveat[idx]) - saved_values.saveval[idx]) < 8.e-15, eachindex(saved_values.t))

# saveat, tdir<0, inplace problem
prob2D_inverse = ODEProblem(prob2D.f, prob2D.u0, (prob2D.tspan[end], prob2D.tspan[1]),1.01)
saved_values = SavedValues(eltype(prob2D_inverse.tspan), typeof(prob2D_inverse.u0))
saveat = range(prob2D_inverse.tspan[1], stop=prob2D_inverse.tspan[2], length=50)
cb = SavingCallback((u,t,integrator)->copy(u), saved_values, saveat=saveat, tdir=-1)
sol = solve(prob2D_inverse, Tsit5(), callback=cb)
@test all(idx -> saveat[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> norm(sol(saveat[idx]) - saved_values.saveval[idx]) < 8.e-15, eachindex(saved_values.t))

saved_values = SavedValues(eltype(prob2D_inverse.tspan), eltype(prob2D_inverse.u0))
saveat = range(prob2D_inverse.tspan[1], stop=prob2D_inverse.tspan[2], length=50)
cb = SavingCallback((u,t,integrator)->u[1], saved_values, saveat=saveat, tdir=-1)
sol = solve(prob2D_inverse, Tsit5(), callback=cb)
@test all(idx -> saveat[idx] == saved_values.t[idx], eachindex(saved_values.t))
@test all(idx -> abs(sol(saveat[idx])[1] - saved_values.saveval[idx]) < 8.e-15, eachindex(saved_values.t))

# Make sure it doesn't error with mutable in oop
prob = ODEProblem((u,p,t) -> u, rand(4,4), (0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,t,integrator)->(tr(u),norm(u)), saved_values, saveat=0.0:0.1:1.0)
sol = solve(prob, Tsit5(), callback=cb)
println(saved_values.saveval)

# Save only end
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4,4), (0.0,1.0))
saved_values = SavedValues(Float64, Tuple{Float64,Float64})
cb = SavingCallback((u,t,integrator)->(tr(u),norm(u)), saved_values,
                    save_everystep = false, save_start = false)
sol = solve(prob, Tsit5(), callback=cb)
print(saved_values.saveval)
@test length(saved_values.t) == 1
@test saved_values.t[1] == 1.0


text:
using StatsBase
using Test

a = [1, 2, 3, 4, 5, 6, 7]
b = [1, 3, 3, 4, 6, 7, 8]

@test counteq(a, b) == 3
@test countne(a, b) == 4

a = rand(5, 6)
b = rand(5, 6)

@test sqL2dist(a, b) ≈ sum(abs2.(a - b))
@test L2dist(a, b)   ≈ sqrt(sqL2dist(a, b))
@test L1dist(a, b)   ≈ sum(abs.(a - b))
@test Linfdist(a, b) ≈ maximum(abs.(a - b))

@test gkldiv(a, b) ≈ sum(a .* log.(a ./ b) - a + b)

@test meanad(a, b)               ≈ mean(abs.(a - b))
@test maxad(a, b)                ≈ maximum(abs.(a - b))
@test msd(a, b)                  ≈ mean(abs2.(a - b))
@test rmsd(a, b)                 ≈ sqrt(msd(a, b))
@test rmsd(a, b; normalize=true) ≈ rmsd(a, b) / (maximum(a) - minimum(a))
@test psnr(a, b, 2)              ≈ 10 * log10(4 / msd(a, b))


text:
# # Non-linear effects

using Unfold
using CairoMakie
using DataFrames
using Random
using Colors

# ## Generating a non-linea signal
# We start with generating
rng = MersenneTwister(2) # make repeatable
n = 20 # datapoints
evts = DataFrame(:x=>rand(rng,n))
signal = -(3*(evts.x .-0.5)).^2 .+ 0.5 .* rand(rng,n)

plot(evts.x,signal)
#
# Looks perfectly non-linear - great!
#
# # Compare linear & non-linear fit
# First we have to reshape it to fit to Unfold format, a 3d array,  1 channel x 1 timepoint x 20 datapoints
signal = reshape(signal,length(signal),1,1)
signal = permutedims(signal,[3,2,1])
size(signal)

# Next we define three different models. **linear** **3 splines** and **10 splines**
# note the different formulas, one `x` the other `spl(x,3)`
design_linear = Dict(Any=>	(@formula(0~1+x),[0]))
design_spl3   = Dict(Any=>	(@formula(0~1+spl(x,3)),[0]))
design_spl10  = Dict(Any=>	(@formula(0~1+spl(x,10)),[0])) #hide

# and we fit the parameters
uf_linear = fit(UnfoldModel,design_linear,evts,signal);
uf_spl3 = fit(UnfoldModel,design_spl3,evts,signal);
uf_spl10 = fit(UnfoldModel,design_spl10,evts,signal); #hide

# extract the fitted values using Unfold.effects

p_linear = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_linear);
p_spl3 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl3);
p_spl10 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl10);
first(p_linear,5)

# And plot them
pl = plot(evts.x,signal[1,1,:])
	lines!(p_linear.x,p_linear.yhat)
	lines!(p_spl3.x,p_spl3.yhat)
	lines!(p_spl10.x,p_spl10.yhat)
pl

# What we can clearly see here, that the linear effect (blue) underfits the data, the `spl(x,10)` overfits it, but the `spl(x,3)` fits it perfectly.

# ## Looking under the hood
# Let's have a brief look how the splines manage what they are managing.
#
# The most important bit to understand is, that we are replacing `x` by a set of coefficients `spl(x)`.
# These new coefficients each tile the range of `x` (in our case from [0-1]) in overlapping areas, each which will be fit by one coefficient.
# Because the ranges are overlapping, we get a smooth function.
#
# Maybe this becomes clear after looking at such a basisfunction:

term_spl = Unfold.formula(uf_spl10).rhs.terms[2]

# This is the spline term.
typeof(term_spl)

# a special type generated in Unfold.jl
#
# it has a field .fun which is the spline function. We can evaluate it at a point
term_spl.fun([0.2])

# each column of this 1 row matrix is a coefficient for our regression model.
lines(disallowmissing(term_spl.fun([0.2]))[1,:])

# Note: We have to use disallowmissing, because our splines return a `missing` whenever we ask it to return a value outside it's defined range, e.g.:
term_spl.fun([-0.2])

# because it never has seen any data outside and can't extrapolate!

# Okay back to our main issue: Let's plot the whole basis set
basisSet = term_spl.fun(0.:0.01:1)
basisSet = disallowmissing(basisSet[.!any(ismissing.(basisSet),dims=2)[:,1],:]) # remove missings
ax = Axis(Figure()[1,1])
[lines!(ax,basisSet[:,k]) for k in 1:size(basisSet,2)]
current_figure()

# notice how we flipped the plot around, i.e. now on the x-axis we do not plot the coefficients, but the `x`-values
# Now each line is one basis-function of the spline.
#
# Unfold returns us one coefficient per basis-function
β=coef(uf_spl10)[1,1,:]
β = Float64.(disallowmissing(β))

# but because we used an intercept, we have to do some remodelling in the basisSet
X = hcat(ones(size(basisSet,1)), basisSet[:,1:5], basisSet[:,7:end])

# now we can weight the spline by the basisfunction
weighted = (β .*X')

# Plotting them makes for a nice looking plot!
ax = Axis(Figure()[1,1])
[lines!(weighted[k,:]) for k = 1:11]
current_figure()

# and sum them up
lines(sum(weighted,dims=1)[1,:])
plot!(X*β,color="gray") #(same as matrixproduct X*β directly!)
current_figure()

# And this is how you can think about splines

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
#
# Run the simulation steps by position incrementation
#

using Plots
using Printf

function runsim(nsteps, dt, m, data)
    p = initial_point(data)
    f = [ Array{Float64}(undef, 2) for i = 1:data.N ]
    v = initial_velocity(data.T, m, data.N)
    e = runsim!(p, f, v, nsteps, dt, m, data)
    return p, f, v, e
end

function runsim!(p, f, v, nsteps, dt, m, data)
    iprint = 50_000
    np = trunc(Int64, nsteps/iprint) + 1
    e = Array{Float64}(undef, np, 2)
    @printf("        u        |        k        |       et        \n")
    @printf("-----------------+-----------------+-----------------\n")
    for i = 1:nsteps
        u = utotal3!(f, p, data)
        k = 0.
        if i%iprint == 0 || i == 1
            x=Array{Float64}(undef, data.N)
            y=Array{Float64}(undef, data.N)
        end
        for ip = 1:data.N
            k += m*(v[ip][1]^2 + v[ip][2]^2)/2
            @. p[ip] = p[ip] + v[ip]*dt + f[ip]*dt^2/(2*m)
            @. v[ip] = v[ip] + f[ip]*dt/m
            if i%iprint == 0 || i == 1
                x[ip] = p[ip][1]
                y[ip] = p[ip][2]
            end
        end
        if i%iprint == 0 || i == 1
            ipr = trunc(Int64, i/iprint) + 1
            e[ipr,:] = [u; k]
            @printf("  % e  |  % e  |  % e  \n", u, k, u+k)
            display(scatter(x,y,xlims=(0,100),ylims=(0,100), markersize=4, label=false))
        end
    end
    return e
end


text:
module CUSPARSE

import CUDAdrv: CUDAdrv, CuContext, CuStream_t

using ..CuArrays
using ..CuArrays: libcusparse, active_context

using SparseArrays
using LinearAlgebra

import Base.one
import Base.zero

const SparseChar = Char
import Base.one
import Base.zero

export CuSparseMatrixCSC, CuSparseMatrixCSR,
       CuSparseMatrixHYB, CuSparseMatrixBSR,
       CuSparseMatrix, AbstractCuSparseMatrix,
       CuSparseVector

include("libcusparse_types.jl")
include("error.jl")
include("libcusparse.jl")

const _handles = Dict{CuContext,cusparseHandle_t}()
const _handle = Ref{cusparseHandle_t}()

function handle()
    if _handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cusparseCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusparseDestroy(handle))
            handle
        end
    end

    return _handle[]
end

include("array.jl")
include("util.jl")
include("wrappers.jl")
include("highlevel.jl")

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

abstract type AbstractSOCtoNonConvexQuadBridge{T} <: AbstractBridge end

"""
    SOCtoNonConvexQuadBridge{T}

Constraints of the form `VectorOfVariables`-in-`SecondOrderCone` can be
transformed into a `ScalarQuadraticFunction`-in-`LessThan` and a
`ScalarAffineFunction`-in-`GreaterThan`. Indeed, the definition of the
second-order cone
```math
t \\ge \\lVert x \\rVert_2 \\  (1)
```
is equivalent to
```math
\\sum x_i^2 \\le t^2  (2)
```
with ``t \\ge 0``.  (3)

!!! warning
    This transformation starts from a convex constraint (1) and creates a
    non-convex constraint (2), because the Q matrix associated with the constraint (2)
    has one negative eigenvalue. This might be wrongly interpreted by a solver.
    Some solvers can look at (2) and understand that it is a second order cone, but
    this is not a general rule.
    For these reasons this bridge is not automatically added by [`MOI.Bridges.full_bridge_optimizer`](@ref).
    Care is recommended when adding this bridge to a optimizer.
"""
struct SOCtoNonConvexQuadBridge{T} <: AbstractSOCtoNonConvexQuadBridge{T}
    quad::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{T},MOI.LessThan{T}}
    var_pos::Vector{
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}},
    }
    vars::Vector{MOI.VariableIndex}
end

const SOCtoNonConvexQuad{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoNonConvexQuadBridge{T},OT}

function bridge_constraint(
    ::Type{SOCtoNonConvexQuadBridge{T}},
    model,
    func::MOI.VectorOfVariables,
    set::MOI.SecondOrderCone,
) where {T}
    vis = func.variables

    t = vis[1]
    x = vis[2:end]
    a_terms = MOI.ScalarAffineTerm{T}[]
    q_terms = MOI.ScalarQuadraticTerm{T}[]
    push!(q_terms, MOI.ScalarQuadraticTerm(-T(2), t, t))
    for var in x
        push!(q_terms, MOI.ScalarQuadraticTerm(T(2), var, var))
    end

    fq = MOI.ScalarQuadraticFunction(q_terms, a_terms, zero(T))
    quad = MOI.add_constraint(model, fq, MOI.LessThan(zero(T)))
    # ScalarAffineFunction's are added instead of VariableIndex's
    # because models can only have one VariableIndex per variable.
    # Hence, adding a VariableIndex constraint here could conflict with
    # a user defined VariableIndex
    fp = convert(MOI.ScalarAffineFunction{T}, t)
    var_pos = MOI.add_constraint(model, fp, MOI.GreaterThan(zero(T)))

    return SOCtoNonConvexQuadBridge(quad, [var_pos], vis)
end

"""
    RSOCtoNonConvexQuadBridge{T}

Constraints of the form `VectorOfVariables`-in-`SecondOrderCone` can be
transformed into a `ScalarQuadraticFunction`-in-`LessThan` and a
`ScalarAffineFunction`-in-`GreaterThan`. Indeed, the definition of the
second-order cone
```math
2tu \\ge \\lVert x \\rVert_2^2, t,u \\ge 0  (1)
```
is equivalent to
```math
\\sum x_i^2 \\le 2tu  (2)
```
with ``t,u \\ge 0``.  (3)

*WARNING* This transformation starts from a convex constraint (1) and creates a
non-convex constraint (2), because the Q matrix associated with the constraint 2
has two negative eigenvalues. This might be wrongly interpreted by a solver.
Some solvers can look at (2) and understand that it is a rotated second order cone, but
this is not a general rule.
For these reasons, this bridge is not automatically added by [`MOI.Bridges.full_bridge_optimizer`](@ref).
Care is recommended when adding this bridge to an optimizer.
"""
struct RSOCtoNonConvexQuadBridge{T} <: AbstractSOCtoNonConvexQuadBridge{T}
    quad::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{T},MOI.LessThan{T}}
    var_pos::Vector{
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}},
    }
    vars::Vector{MOI.VariableIndex}
end

const RSOCtoNonConvexQuad{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoNonConvexQuadBridge{T},OT}

function bridge_constraint(
    ::Type{RSOCtoNonConvexQuadBridge{T}},
    model,
    func::MOI.VectorOfVariables,
    set::MOI.RotatedSecondOrderCone,
) where {T}
    vis = func.variables

    t = vis[1]
    u = vis[2]
    x = vis[3:end]
    a_terms = MOI.ScalarAffineTerm{T}[]
    q_terms = MOI.ScalarQuadraticTerm{T}[]
    push!(q_terms, MOI.ScalarQuadraticTerm(-T(2), t, u))
    for var in x
        push!(q_terms, MOI.ScalarQuadraticTerm(T(2), var, var))
    end

    fq = MOI.ScalarQuadraticFunction(q_terms, a_terms, zero(T))
    quad = MOI.add_constraint(model, fq, MOI.LessThan(zero(T)))
    # ScalarAffineFunction's are added instead of VariableIndex's
    # because models can only have one VariableIndex per variable.
    # Hence, adding a VariableIndex constraint here could conflict with
    # a user defined VariableIndex
    fp1 = convert(MOI.ScalarAffineFunction{T}, t)
    var_pos1 = MOI.add_constraint(model, fp1, MOI.GreaterThan(zero(T)))
    fp2 = convert(MOI.ScalarAffineFunction{T}, u)
    var_pos2 = MOI.add_constraint(model, fp2, MOI.GreaterThan(zero(T)))

    return RSOCtoNonConvexQuadBridge(quad, [var_pos1, var_pos2], vis)
end

function MOI.supports_constraint(
    ::Type{SOCtoNonConvexQuadBridge{T}},
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.SecondOrderCone},
) where {T}
    return true
end

function MOI.supports_constraint(
    ::Type{RSOCtoNonConvexQuadBridge{T}},
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.RotatedSecondOrderCone},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(
    ::Type{<:AbstractSOCtoNonConvexQuadBridge},
)
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(
    ::Type{<:AbstractSOCtoNonConvexQuadBridge{T}},
) where {T}
    return Tuple{Type,Type}[
        (MOI.ScalarQuadraticFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
    ]
end

function concrete_bridge_type(
    ::Type{SOCtoNonConvexQuadBridge{T}},
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.SecondOrderCone},
) where {T}
    return SOCtoNonConvexQuadBridge{T}
end

function concrete_bridge_type(
    ::Type{RSOCtoNonConvexQuadBridge{T}},
    ::Type{MOI.VectorOfVariables},
    ::Type{MOI.RotatedSecondOrderCone},
) where {T}
    return RSOCtoNonConvexQuadBridge{T}
end

# Attributes, Bridge acting as a model
function MOI.get(
    ::AbstractSOCtoNonConvexQuadBridge{T},
    ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{T},MOI.LessThan{T}},
)::Int64 where {T}
    return 1
end

function MOI.get(
    bridge::AbstractSOCtoNonConvexQuadBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarQuadraticFunction{T},
        MOI.LessThan{T},
    },
) where {T}
    return [bridge.quad]
end

function MOI.get(
    bridge::AbstractSOCtoNonConvexQuadBridge{T},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}},
)::Int64 where {T}
    return length(bridge.var_pos)
end

function MOI.get(
    bridge::AbstractSOCtoNonConvexQuadBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T},
        MOI.GreaterThan{T},
    },
) where {T}
    return copy(bridge.var_pos)
end

# References
function MOI.delete(
    model::MOI.ModelLike,
    bridge::AbstractSOCtoNonConvexQuadBridge,
)
    MOI.delete(model, bridge.quad)
    MOI.delete.(model, bridge.var_pos)
    return
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::AbstractSOCtoNonConvexQuadBridge,
)
    vals = MOI.get.(model, MOI.VariablePrimal(attr.result_index), bridge.vars)
    return vals
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    b::SOCtoNonConvexQuadBridge{T},
) where {T}
    return MOI.SecondOrderCone(length(b.vars))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    b::RSOCtoNonConvexQuadBridge{T},
) where {T}
    return MOI.RotatedSecondOrderCone(length(b.vars))
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    b::AbstractSOCtoNonConvexQuadBridge{T},
) where {T}
    return MOI.VectorOfVariables(b.vars)
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using StatsBase

"""
https://stats.stackexchange.com/questions/111759/exponential-family-observed-vs-expected-sufficient-statistics
"""

mutable struct MvWatson
    u::Vector{Float64}
    k::Float64
end
const SufficientStatistic = Vector

mutable struct MvWatsonExponential <: ExponentialFamily
    dim::Integer
    concentration::Float64
    n::Vector{Float64}
    theta::Vector{Float64}
    DLNF::Float64
end

function MvWatsonExponential(n::Vector{Float64})
    dim = round(Int64,(-1+sqrt(1+(8*length(n))))/2)
    w   = MvWatsonExponential(dim, 1, n, Float64[],0)
    update_parameters!(w)
    return w
end
"""
Apply Newton-Raphson method root finder  to approximate ``κ`` from
 ``\\vert \\vert n \\vert \\vert_2`` using the following iterative update iteration:
 ``k_{l+1} = k_l - \\dfrac{q(a,c;k_l) - \\vert \\vert n \\vert \\vert_2}{q'(a,c,k_l)} ``

"""
function estimate_concentration(dim::Integer, n::Vector{<:AbstractFloat}; tol::Float64 = 0.000001)
    a     = 0.5
    c     = (dim)/2
    r = norm(n)
    kappa_upper_bound = 500;
    k = 1.0
    diff = 100;
    ii=2;
    prevk = Inf
    while  abs(k - prevk) > tol && (k < kappa_upper_bound)
        prevk=k
        q = kummer_ratio(a,c,k, tol = eps()*2)
        dq = (1 - (c/k))*q + (a/k) -q^2
        k = k - ( (q - r) / (dq))
    end
    return k
end
function update_parameters!(watson::MvWatsonExponential)
    d = dimension(watson)
    watson.concentration = estimate_concentration(d,watson.n);
    g_norm_theta = kummer_ratio(0.5,d/2, watson.concentration)
    R_norm_theta = g_norm_theta/ watson.concentration;
    watson.theta = watson.n ./ R_norm_theta; # natural parameter
    LNF = log(M(0.5, d/2, watson.concentration));# The log normalizing function
    watson.DLNF = dot(watson.theta, watson.n) - LNF; # Dual (expectation parameter) of the log normalizing function

end
function bregman_divergence(wmm::MvWatsonExponential, x::Vector)
    n     = wmm.n
    theta = wmm.theta
    return  wmm.DLNF +  dot(x-n,theta)
end
"""
Computes the bregman divergence of every pattern located in the columns of x
"""
function bregman_divergence(wmm::MvWatsonExponential, x::Matrix)
    n = wmm.n
    theta = wmm.theta
    res =  wmm.DLNF .+ theta'*(x.-n)
    return res
end
function bregman_divergence(wmm1::MvWatsonExponential, wmm2::MvWatsonExponential)
    n1 = wmm1.n
    n2 = wmm2.n
    theta2 = wmm2.theta
    return  wmm1.DLNF - wmm2.DLNF - dot(n1-n2,theta2)
end

function bregman_divergence(wmm1::MvWatson, wmm2::MvWatson)

    a     = 0.5
    c     = dimension(wmm1)/2
    theta1     = wmm1.k* sufficient_statistic(MvWatsonExponential, wmm1.u)
    theta2     = wmm2.k* sufficient_statistic(MvWatsonExponential, wmm2.u)
    return  log(M(a,c,wmm1.k)) - log(M(a,c,wmm2.k)) - dot(theta1-theta2, kummer_ratio(a,c,wmm2.k)*(theta2/wmm2.k))
end
function sufficient_statistic_dimension(d::Integer)
    return round(Integer,d+ (d*(d-1))/2)
end
function MvWatsonExponential(dim::Integer)
    p = sufficient_statistic_dimension(dim)
    return MvWatsonExponential(dim,0.0, zeros(p), zeros(p), 0.0)
end
function dimension(w::MvWatsonExponential)
    return w.dim
end
function sufficient_statistic(::Type{MvWatsonExponential}, d::Integer, n::Integer=1)
    p = sufficient_statistic_dimension(d)
    if n==1
	return zeros(p)
    else
	return zeros(p,n)
    end
end
function sufficient_statistic(t::Type{MvWatsonExponential}, x::Vector{<:AbstractFloat})
    y = sufficient_statistic(t,length(x))
    sufficient_statistic!(t,y,x)
    return y
end
function sufficient_statistic(t::Type{MvWatsonExponential}, x::AbstractArray{<:Number,2})
    y = zeros(
        sufficient_statistic_dimension(size(x,1)),
        size(x,2))
    for i=1:size(x,2)
        y[:, i] = sufficient_statistic(t,x[:, i])
    end
    return y
end

function sufficient_statistic(t::Type{MvWatsonExponential},  x::Array{<:AbstractFloat,3})
    y = zeros(
        sufficient_statistic_dimension(size(x,1)),
        size(x,2)*size(x,3))
    idx = 1
    for i=1:size(x,2), j=1:size(x,3)
        y[:,idx] = sufficient_statistic(t,x[:, i ,j])
        idx += 1
    end
    return y
end
function sufficient_statistic!(::Type{MvWatsonExponential}, y::Vector{<:AbstractFloat}, x::Vector{<:AbstractFloat})
    y[1:length(x)] .= x.^2
    r = length(x)+1
    for j=1:length(x)-1
        for i=j+1:length(x)
            y[r] = x[j]*x[i]*sqrt(2)
            r=r+1
        end
    end
end
function dimension(d::MvWatson)
    return length(d.u)
end
"""

Kummer’s confluent hypergeometric function
https://www.mathworks.com/matlabcentral/fileexchange/29766-confluent-hypergeometric-function?focused=5192808&tab=function

 Estimates the value by summing powers of the generalized hypergeometric
 series:

 ``sum(n=0-->Inf)[(a)_n*x^n/{(b)_n*n!}``

 until the specified tolerance is acheived.

"""
function M1(a::Float64, b::Float64, z::T; tol=0.00001) where T<:AbstractFloat
    term = BigFloat((z*a)/b)
    f = BigFloat(1 + term)
    n = 1;
    an = a
    bn = b;
    nmin = 50;
    while (n < nmin) || abs(term) > tol
        n = n + 1;
        an = an + 1;
        bn = bn + 1;
        term = (term*an*z)/(bn*n);
        prevf = f
        f = f + term;
        if isinf(f)
            return Inf
        end
    end
    return Float64(f)
end
function kummer_ratio(a::Float64, c::Float64, x::T; tol::Float64=0.00001) where T<:AbstractFloat
    val =  (a/c)*(M(a+1,c+1,x)/(M(a,c,x)))
    if isnan(val)
        return 1
    else
        return val
    end
end
import Distributions.pdf
import Distributions.logpdf!
function pdf(w::MvWatson,x::Vector{<:AbstractFloat})::Float64
    p = dimension(x)
    return ((gamma(p/2)/(2*pi^(p/2)*M(0.5,p)/2,w.k)))*exp(w.k*(dot(w.u,x)^2))
end
function pdf(w::MvWatson, m::Matrix{<:AbstractFloat})
    p = dimension(w)
    c = ((gamma(p/2)/(2*pi^(p/2)*M(0.5,p)/2,w.k)))
    return prod(c*exp.(w.k.*(w.u'*v).^2))
end
function logpdf!(output, w::MvWatson, data::Matrix)
    dim = size(data,2)
    coeff = log(1/M(0.5,dim/2,w.k))
    output .=vec(w.k.*(data'*w.u).^2 .+ coeff)
    ind = output.>700;
    output[ind] = output[ind] .- 1000;
end

function M(a,b,x)
    nl = 0
    ta=[]
    tb=[]
    xg=[]
    tba=[]
    a0=a;
    a1=a;
    x0=x
    hg=0.0
    if (b == 0.0 || b == -abs(floor(b)))
	hg=1.0+300;
    elseif (a == 0.0 || x == 0.0)
	hg=1.0;
    elseif (a == -1.0)
    	hg=1.0-x./b
    elseif (a == b)
    	hg=exp(x)
    elseif (a-b == 1.0)
    	hg=(1.0+x./b).*exp(x)
    elseif (a == 1.0&&b == 2.0)
    	hg=(exp(x)-1.0)./x
    elseif (a == floor(a)&& a < 0.0)
    	m=floor(-a)
	r=1.0
    	hg=1.0
	for  k=1:m
	    r=r.*(a+k-1.0)./k./(b+k-1.0).*x
	    hg=hg+r
	end
    end
    if (hg != 0.0)
	return hg
    end
    if (x < 0.0)
	a=b-a
	a0=a
    	x=abs(x)
    end
    if (a < 2.0)
	nl=0
    end
    if (a >= 2.0)
	nl=1
	la=floor(a)
	a=a-la-1.0
    end
    y0 = 0
    y1 = 1

    for n=0:nl
	if (a0 >= 2.0)
	    a=a+1.0
	end
	if (x <= 30.0+abs(b)||a < 0.0)
    	    hg=1.0
            rg=1.0
	    for  j=1:500
    	        rg=rg.*(a+j-1.0)./(j.*(b+j-1.0)).*x
        	hg=hg+rg;
            	if (abs(rg./hg) < 1.0e-15)
		    break
		end
    	    end
	else
    	    (a,ta)=gamma(a,ta)
	    (b,tb)=gamma(b,tb)
    	    xg=b-a
	    (xg,tba)=gamma(xg,tba)
	    sum1=1.0
	    sum2=1.0
	    r1=1.0
	    r2=1.0
	    for  i=1:8
	        r1=-r1.*(a+i-1.0).*(a-b+i)./(x.*i)
	        r2=-r2.*(b-a+i-1.0).*(a-i)./(x.*i)
	        sum1=sum1+r1
	        sum2=sum2+r2
	    end
	    hg1=tb./tba.*x.^(-a).*cos(pi.*a).*sum1;
	    hg2=tb./ta.*exp(x).*x.^(a-b).*sum2;
    	    hg=hg1+hg2;
	end
    	if (n == 0)
	    y0=hg
	end
    	if (n == 1)
	    y1=hg
	end
    end
    if (a0 >= 2.0)
	for  i=1:la-1
	    hg=((2.0.*a-b+x).*y1+(b-a).*y0)./a;
    	    y0=y1;
            y1=hg;
	    a=a+1.0;
    	end
    end
    if (x0 < 0.0)
	hg=hg.*exp(x0)
    end
    a=a1
    x=x0
    if isinf(hg)
        return realmax()
    else
        return hg
    end
end

function gamma(x,ga)
    g=zeros(26);
    if (x == floor(x))
    	if (x > 0.0)
            ga=1.0
	    m1=x-1
    	    for  k=2:m1
        	ga=ga.*k;
	    end
    	else
            ga=1.0e300
	end
    else
    	if (abs(x) > 1.0)
	    z=abs(x)
    	    m=floor(z);
	    r=1.0;
	    for  k=1:m;
    	        r=r.*(z-k);
	    end;
    	    z=z-m;
	else
    	    z=x;
	end
	g=[1.0,0.5772156649015329,-0.6558780715202538,-0.420026350340952e-1,0.1665386113822915,-.421977345555443e-1,-.96219715278770e-2,.72189432466630e-2,-.11651675918591e-2,-.2152416741149e-3,.1280502823882e-3,-.201348547807e-4,-.12504934821e-5,.11330272320e-5,-.2056338417e-6,.61160950e-8,.50020075e-8,-.11812746e-8,.1043427e-9,.77823e-11,-.36968e-11,.51e-12,-.206e-13,-.54e-14,.14e-14,.1e-15];
        gr=g[26]
	for  k=25:-1:1
            gr=gr.*z+g[k];
	end
    	ga=1.0./(gr.*z)
	if (abs(x) > 1.0)
    	    ga=ga.*r;
            if (x < 0.0)
		ga=-pi./(x.*ga.*sin(pi.*x))
	    end
    	end
    end
    return (x,ga)
end
"""
Return the same dsitribution using the natural parameters
"""
function  source_parameters(m::MvWatsonExponential)
    μ = source_vector(MvWatsonExponential, m.n)
    return MvWatson(μ/norm(μ), m.concentration)
end
function source_vector(m::Type{MvWatsonExponential}, W::SufficientStatistic)
    dim = round(Int64,(-1+sqrt(1+(8*length(W))))/2)
    x_init = sqrt.(W[1:dim])
    W = W.>0
    x = x_init
    if (W[4] & !W[5] & !W[6])
        x[3] = -x_init[3]
    elseif !W[4] && !W[5] & W[6]
        x[1] = -x_init[1]
    elseif !W[4] && W[5] && !W[6]
        x[2] = -x_init[2]
    end
    return x
end

# The Multivariate Watson Distribution: Maximum-Likelihood Estimation and other Aspects
#Suvrit Sra and Dmitrii Karp
import Distributions.fit_mle
function L(r::T, c::T, a::T) where T<:Number
    return ((r*c-a)/(r*(1-r)))*(1+((1-r)/(c-a)))
end
function B(r::T, c::T, a::T) where T<:Number
    return ((r*c-a)/(2*r*(1-r)))*(1+sqrt(1+ ((r*(c+1)*r*(1-r))/(a*(c-a)))))
end
function U(r::T, c::T, a::T) where T<:Number
    return ((r*c-a)/(r*(1-r)))*(1+(r/a))
end
function concentration(r::T1, c::T2, a::T3) where T1<:Number where T2<:Number where T3<:Number
    if (r>0) && (r<a/(2*c))
        return U(r,c,a)
    elseif (r>=(a/(2*c))) && (r<((2*a)/(2*sqrt(c))))
        return B(r,c,a)
    elseif (r>=((2*a)/(2*sqrt(c)))) && (r<1)
        return L(r,c,a)
    end
end
function concentration(v::Vector, S::Matrix, c::Float64, a::Float64)
    r = v'*S*v
    if (r<=0)
        r=eps()
    end
    if r>=1
        r=1-5*eps()
    end
    return concentration(float64(r), float64(c), float64(a))

end

import Distributions.loglikelihood
function  loglikelihood(d::MvWatson, X::AbstractMatrix)
    n = size(X,2)
    p = size(X,1)
    S = StatsBase.scattermat_zm(X,2)/n
    return loglikelihood_from_scatter(d, n, S)
end
function loglikelihood_from_scatter(d::MvWatson, n::Integer, S::AbstractMatrix)
    r = d.u'*S*d.u
    p = size(S,1)
    return n*(d.k*r - log(M(1/2, p/2, d.k)))
end
function fit_mle(m::Type{MvWatson}, X::Matrix)
    n = size(X,2)
    p = size(X,1)
    S = StatsBase.scattermat_zm(X,2)/n
    a = 1/2
    c = p/2

    (eval,evec) = eig(S)
    w_1 = MvWatson(evec[:,1], concentration(evec[:,1], S, c, a))
    w_2 = MvWatson(evec[:,3], concentration(evec[:,3], S, c, a))

    llh1 = loglikelihood_from_scatter(w_1, n, S)
    llh2 = loglikelihood_from_scatter(w_2, n, S)

    if (llh1 > llh2)
        return w_1
    else
        return w_2
    end

end
#=facts("Watson Distribution") do
we = MvWatsonExponential(3)
@fact sufficientStatistic( [3.0;1.0;5.0]) --> roughly([9.0; 1.0; 25.0; 4.24264; 21.2132; 7.07107]; atol=0.0001)
@fact M(0.5, 3/2, 0.5) --> roughly(1.19496, atol=0.001)
@fact kummer_ratio(0.5, 3/2, 0.5) --> roughly(0.37973, atol=0.001)
end=#


text:
mutable struct GeneratedAttributes
    length::Int64
    counts::Dict{ParticleType, Int64}
    energy_min::Dict{ParticleType, Float32}
    energy_max::Dict{ParticleType, Float32}
    energy_sum::Dict{ParticleType, Float64}
    x_min::Float32
    x_max::Float32
    y_min::Float32
    y_max::Float32
    z_min::Float32
    z_max::Float32
end

function GeneratedAttributes()
    GeneratedAttributes(0,
        Dict(),
        Dict(),
        Dict(),
        Dict(),
        Inf, -Inf,
        Inf, -Inf,
        Inf, -Inf
    )
end

mutable struct IAEAPhspWriter{R, I <: IO}
    record_contents::R
    generated_attributes::GeneratedAttributes
    io_header::I
    io_phsp::I
    function IAEAPhspWriter(record_contents::R, generated_attributes,
                            io_header::I, io_phsp::I) where {R, I}
        w = new{R,I}(record_contents, generated_attributes,
                io_header, io_phsp)
        finalizer(close, w)
        w
    end
end

function increment(d, key, val=1)
    T = eltype(values(d))
    d[key] = get!(d,key,zero(T)) + val
end

function Base.write(io::IAEAPhspWriter, p::IAEAParticle)
    ret = write_particle(io.io_phsp, p, io.record_contents)

    ga = io.generated_attributes
    typ = p.typ
    increment(ga.counts, typ, 1)
    ga.length += 1
    increment(ga.energy_sum, typ, p.E * p.weight)
    ga.energy_min[typ] = min(get!(ga.energy_min, typ, Inf), p.E)
    ga.energy_max[typ] = max(get!(ga.energy_max, typ, -Inf), p.E)

    ga.x_min = min(ga.x_min, p.x)
    ga.y_min = min(ga.y_min, p.y)
    ga.z_min = min(ga.z_min, p.z)
    ga.x_max = max(ga.x_max, p.x)
    ga.y_max = max(ga.y_max, p.y)
    ga.z_max = max(ga.z_max, p.z)

    ret
end

function Base.write(w::IAEAPhspWriter, ps)
    ret = 0
    for p in ps
        ret += write(w,p)
    end
    ret
end

function print_key(io::IO, k)
    println(io,'$',k,":")
end
function print_val(io, v)
    if v != ""
        println(io,v)
    end
end
function println_kv(io::IO, k, v)
    print_key(io,k)
    print_val(io,v)
    println(io)
end

function println_kv_get(io::IO, attr::AbstractDict, k, default)
    v = get(attr, k, default)
    println_kv(io, k, v)
end

function println_particle_count(io, d, p::ParticleType)
    key = if p == photon
            :PHOTONS
        elseif p == electron
            :ELECTRONS
        elseif p == positron
            :POSITRONS
        elseif p == proton
            :PROTONS
        elseif p == neutron
            :NEUTRONS
        else
            error("Unknown particle type $p")
        end
    println_kv(io, key, d[p])
end

function print_header(io::IO, w::IAEAPhspWriter)
    attr = w.record_contents.attributes
    println_kv_get(io, attr, :IAEA_INDEX, "0 // test header")
    println_kv_get(io, attr, :TITLE, "")
    println_kv(io, :FILE_TYPE, 0)
    r = w.record_contents
    ga = w.generated_attributes
    record_length = ptype_disksize(r)
    checksum = ga.length * record_length
    println_kv(io, :CHECKSUM, checksum)
    print_record_contents(io, r)
    println_kv(io, :RECORD_LENGTH, record_length)

    println_kv(io, :BYTE_ORDER, "1234")
    println_kv_get(io, attr, :ORIG_HISTORIES, "$(typemax(Int32))")
    println_kv(io, :PARTICLES, ga.length)

    for p in keys(ga.counts)
        println_particle_count(io, ga.counts, p)
    end

    println_kv_get(io, attr, :TRANSPORT_PARAMETERS, "")
    println_kv_get(io, attr, :MACHINE_TYPE, "")
    println_kv_get(io, attr, :MONTE_CARLO_CODE_VERSION, "")
    #
    println_kv_get(io, attr, :GLOBAL_PHOTON_ENERGY_CUTOFF, 0.0)
    println_kv_get(io, attr, :GLOBAL_PARTICLE_ENERGY_CUTOFF, 0.0)
    println_kv_get(io, attr, :COORDINATE_SYSTEM_DESCRIPTION, "")

    println(io)
    println(io, "// OPTIONAL INFORMATION")
    println(io)
    println_kv_get(io, attr, :BEAM_NAME, "")
    println_kv_get(io, attr, :FIELD_SIZE, "")
    println_kv_get(io, attr, :NOMINAL_SSD, "")
    println_kv_get(io, attr, :MC_INPUT_FILENAME, "")
    println_kv_get(io, attr, :VARIANCE_REDUCTION_TECHNIQUES, "")
    println_kv_get(io, attr, :INITIAL_SOURCE_DESCRIPTION, "")
    println_kv_get(io, attr, :PUBLISHED_REFERENCE, "")
    println_kv_get(io, attr, :AUTHORS, "")
    println_kv_get(io, attr, :INSTITUTION, "")
    println_kv_get(io, attr, :LINK_VALIDATION, "")
    println_kv_get(io, attr, :ADDITIONAL_NOTES, "Generated via PhaseSpaceIO.jl")
    # # TODO:
    # println_kv(io, :STATISTICAL_INFORMATION_PARTICLES, "")
    stats = """
    // cm
       $(ga.x_min)  $(ga.x_max)
       $(ga.y_min)  $(ga.y_max)
       $(ga.z_min)  $(ga.z_max)"""
    println_kv(io, :STATISTICAL_INFORMATION_GEOMETRY, stats)
end

function extra_float_count(r::IAEAHeader{Nf, Ni, Nt}) where {Nf, Ni, Nt}
    Nf
end

function extra_long_count(r::IAEAHeader{Nf, Ni, Nt}) where {Nf, Ni, Nt}
    Ni
end

function print_record_contents(io::IO, r::IAEAHeader)
    # $RECORD_CONTENTS:
    # 1     // X is stored ?
    # 1     // Y is stored ?
    # 1     // Z is stored ?
    # 1     // U is stored ?
    # 1     // V is stored ?
    # 1     // W is stored ?
    # 1     // Weight is stored ?
    # 0     // Extra floats stored ?
    # 1     // Extra longs stored ?
    # 0     // Generic integer variable stored in the extralong array [ 0]

    # $RECORD_CONSTANT:
    #
    # $RECORD_LENGTH:
    # 33
    t = r.record_contents

    print_key(io, "RECORD_CONTENTS")
    record_constants = Float32[]
    for propstr in ["X","Y","Z", "U","V","W", "Weight"]
        field = Symbol(lowercase(string(propstr)))
        if field in propertynames(t)
            val = 0
            push!(record_constants, getproperty(t, field))
        else
            val = 1
        end
        println(io, val, " "^5, "// ", propstr, " is stored ?")
    end
    Nf = extra_float_count(r)
    Ni = extra_long_count(r)
    println(io, Nf, " "^5, "// Extra floats stored")
    println(io, Ni, " "^5, "// Extra longs stored")
    println(io, "0     // Generic integer variable stored in the extralong array [ 0]")
    println(io)

    print_key(io, "RECORD_CONSTANT")
    for c in record_constants
        println(io, c)
    end
    println(io)
end

function iaea_writer(path::IAEAPath, r::IAEAHeader)
    io_header = open(path.header, "w")
    io_phsp   = open(path.phsp  , "w")
    ga = GeneratedAttributes()
    writer = IAEAPhspWriter(r, ga, io_header, io_phsp)
end
iaea_writer(path) = iaea_writer(IAEAPath(path))

function Base.flush(w::IAEAPhspWriter)
    if isopen(w.io_header)
        print_header(w.io_header, w)
    end
    flush(w.io_header)
    flush(w.io_phsp)
end

function Base.close(w::IAEAPhspWriter)
    flush(w)
    close(w.io_header)
    close(w.io_phsp)
end

"""
    iaea_writer(f, path, r::IAEAHeader)

Write particles in IAEA format to `path`:
```jldoctest
julia> using PhaseSpaceIO

julia> h = IAEAHeader{0,0}();

julia> path = IAEAPath(tempname());

julia> iaea_writer(path, h) do w
           p = IAEAParticle(x=1,y=2,z=3,u=0,v=1,w=0,weight=5, typ=photon, E=6);
           write(w, p)
       end
29

julia> phsp_iterator(collect, path)
1-element Array{IAEAParticle{0,0},1}:
 IAEAParticle(typ=photon, E=6.0, weight=5.0, x=1.0, y=2.0, z=3.0, u=0.0, v=1.0, w=0.0, new_history
=true, extra_floats=(), extra_ints=())
```
"""
function iaea_writer(f, path, r::IAEAHeader)
    w = iaea_writer(IAEAPath(path), r)
    ret = call_fenced(f, w)
    close(w)
    ret
end


text:
# General function for merge, not used for the moment
#function merge(cliques::AbstractVector, cliquetree::AbstractVector, nb_lc::Int;
#               treshold::Union{Int,Function}=treshold_percent(cliques, 0.1),
#               minimizes::Vector{Function}=[minimize_molzhan],
#               choose_heuritic::Union{Bool,Function}=false,
#               heuristics_args...)
#    if typeof(treshold) <: Int
#        treshold = (;cliques, kwargs...) -> length(cliques) >= treshold
#    end
#    stop = false
#    it = 0
#    nb_cliques_old = length(cliques)
#    while !stop
#        h_index = choose_heuristic(;cliques=cliques,
#                                   cliquetree=cliquetree,
#                                   it=it,
#                                   nb_cliques_old=nb_cliques_old,
#                                   heuristics_args...)
#        edge = minimizes[h_index](cliques, cliquetree, heuristics_args...)
#        edges = vcat(edges..., get_added_edges(cliques[edge[1]], cliques[edge[2]])...)
#        cliques, cliquetree = merge_cliques(cliques, cliquetree, edge)
#        stop = treshold(;cliques=cliques,
#                        cliquetree=cliquetree,
#                        it=it,
#                        nb_cliques_old=nb_cliques_old,
#                        heuristics_args...)
#    end
#    return cliques, cliquetree, edges
#end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


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

function import_data(dir::String, source::String, revision::String, year::Integer)

    if source == "OECD"

        if revision == "2021"

            file = source * "/" * revision * "/"
            path = ifelse(contains(dir[end-5:end], '.'), dir, dir * file)

            time = 1995:2018
            index = findfirst(time .== year) # get index of specified year

            Z = RData.load(path * "ICIO2021econZ" * ".RData")["ICIO2021econZ"][index,:,:] # intermediate demand
            Y = RData.load(path * "ICIO2021econFD" * ".RData")["ICIO2021econFD"][index,:,:] # final demand
            W = RData.load(path * "ICIO2021econVA" * ".RData")["ICIO2021econVA"][index,:,:] # value added
            X = RData.load(path * "ICIO2021econX" * ".RData")["ICIO2021econX"][index,:,:] # gross output

            println(" ✓ Raw data from $source (rev. $revision) for $year was successfully imported!")

        else
            println(" × The revision - $revision - of $source IO tables is not available!")
        end

    end

    return Z, Y, W, X

end

function remove_reorder(Z::Matrix{Float64}, Y::Matrix{Float64}, W::Matrix{Float64}, X::Matrix{Float64}, N::Integer, S::Integer)

    # merge MX1 and MX2 and substitute for MEX
    # rows
    Z[(N-43)*S+1:(N-42)*S, :] .= Z[(N-43)*S+1:(N-42)*S, :] .+ Z[N*S+1:(N+1)*S, :] .+ Z[(N+1)*S+1:(N+2)*S, :]
    Y[(N-43)*S+1:(N-42)*S, :] .= Y[(N-43)*S+1:(N-42)*S, :] .+ Y[N*S+1:(N+1)*S, :] .+ Y[(N+1)*S+1:(N+2)*S, :]
    W[(N-43)*S+1:(N-42)*S] .= W[(N-43)*S+1:(N-42)*S] .+ W[N*S+1:(N+1)*S] .+ W[(N+1)*S+1:(N+2)*S]
    X[(N-43)*S+1:(N-42)*S] .= X[(N-43)*S+1:(N-42)*S] .+ X[N*S+1:(N+1)*S] .+ X[(N+1)*S+1:(N+2)*S]
    # columns
    Z[:, (N-43)*S+1:(N-42)*S] .= Z[:, (N-43)*S+1:(N-42)*S] .+ Z[:, N*S+1:(N+1)*S] .+ Z[:, (N+1)*S+1:(N+2)*S]
    Y[:, (N-42)] .= Y[:, (N-42)] .+ Y[:, (N+1)] .+ Y[:, (N+2)]
    # added domest trade twice! (actually no double counting because columns/rows are designed to avoid)
    # Z[(N-43)*S+1:(N-42)*S, (N-43)*S+1:(N-42)*S] .= Z[(N-43)*S+1:(N-42)*S, (N-43)*S+1:(N-42)*S] ./ 2
    # Y[(N-43)*S+1:(N-42)*S, (N-42)] .= Y[(N-43)*S+1:(N-42)*S, (N-42)] ./ 2

    # merge CH1 and CH2 and substitute for CHN
    # rows
    Z[(N-24)*S+1:(N-23)*S, :] .= Z[(N-24)*S+1:(N-23)*S, :] .+ Z[(N+2)*S+1:(N+3)*S, :] .+ Z[(N+3)*S+1:(N+4)*S, :]
    Y[(N-24)*S+1:(N-23)*S, :] .= Y[(N-24)*S+1:(N-23)*S, :] .+ Y[(N+2)*S+1:(N+3)*S, :] .+ Y[(N+3)*S+1:(N+4)*S, :]
    W[(N-24)*S+1:(N-23)*S] .= W[(N-24)*S+1:(N-23)*S] .+ W[(N+2)*S+1:(N+3)*S] .+ W[(N+3)*S+1:(N+4)*S]
    X[(N-24)*S+1:(N-23)*S] .= X[(N-24)*S+1:(N-23)*S] .+ X[(N+2)*S+1:(N+3)*S] .+ X[(N+3)*S+1:(N+4)*S]
    # columns
    Z[:, (N-24)*S+1:(N-23)*S] .= Z[:, (N-24)*S+1:(N-23)*S] .+ Z[:, (N+2)*S+1:(N+3)*S] .+ Z[:, (N+3)*S+1:(N+4)*S]
    Y[:, (N-23)] .= Y[:, (N-23)] .+ Y[:, (N+3)] .+ Y[:, (N+4)]
    # added domestic trade twice! (actually no double counting because columns/rows are designed to avoid)
    # Z[(N-24)*S+1:(N-23)*S, (N-24)*S+1:(N-23)*S] .= Z[(N-24)*S+1:(N-23)*S, (N-24)*S+1:(N-23)*S] ./ 2
    # Y[(N-24)*S+1:(N-23)*S, (N-23)] .= Y[(N-24)*S+1:(N-23)*S, (N-23)] ./ 2

    # delete MX1/CH1 and MX2/CH2
    Z = Z[1:N*S, 1:N*S]
    Y = Y[1:N*S, 1:N]
    W = W[1:N*S]
    X = X[1:N*S]

    return Z, Y, W, X

end

function trade(Z::Matrix{Float64}, Y::Matrix{Float64}, N::Integer, S::Integer, dimension::String, flow::String)

    if dimension == "industry"

        # bilateral trade by exporting industry/country and importing country
        #   - remove domestic trade
        GRTR_INT = [sum(Z[i, j:j+S-1]) for i in 1:N*S, j in 1:S:N*S] # NS×N
        GRTR_FNL = copy(Y) # NS×N

        for i in 1:N*S
            GRTR_INT[i, ceil(Int, i/S)] = 0.0
            GRTR_FNL[i, ceil(Int, i/S)] = 0.0
        end

    else

        # bilateral trade by exporting country and importing country
        #   - remove domestic trade
        GRTR_INT = [sum(Z[i:i+S-1, j:j+S-1]) for i in 1:S:N*S, j in 1:S:N*S] # N×N
        GRTR_FNL = [sum(Y[i:i+S-1, j]) for i in 1:S:N*S, j in 1:N] # N×N

        for i in 1:N
            GRTR_INT[i, i] = 0.0
            GRTR_FNL[i, i] = 0.0
        end

    end

    GRTR = GRTR_INT .+ GRTR_FNL

    if flow == "exports"
        return GRTR, GRTR_INT, GRTR_FNL
    else
        return GRTR', GRTR_INT', GRTR_FNL'
    end

end

function exports(Z::Matrix{Float64}, Y::Matrix{Float64}, N::Integer, S::Integer, dimension::String)

    if dimension == "industry"

        # bilateral trade by exporting industry/country and importing country
        #   - remove domestic trade
        EXGR_INT = [sum(Z[i, j:j+S-1]) for i in 1:N*S, j in 1:S:N*S] # NS×N
        EXGR_FNL = copy(Y) # NS×N

        for i in 1:N*S
            EXGR_INT[i, ceil(Int, i/S)] = 0.0
            EXGR_FNL[i, ceil(Int, i/S)] = 0.0
        end

    else

        # bilateral trade by exporting country and importing country
        #   - remove domestic trade
        EXGR_INT = [sum(Z[i:i+S-1, j:j+S-1]) for i in 1:S:N*S, j in 1:S:N*S] # N×N
        EXGR_FNL = [sum(Y[i:i+S-1, j]) for i in 1:S:N*S, j in 1:N] # N×N

        for i in 1:N
            EXGR_INT[i, i] = 0.0
            EXGR_FNL[i, i] = 0.0
        end

    end

    EXGR = EXGR_INT .+ EXGR_FNL

    return EXGR, EXGR_INT, EXGR_FNL

end

function imports(Z::Matrix{Float64}, Y::Matrix{Float64}, N::Integer, S::Integer, dimension::String)

    if dimension == "industry"

        # bilateral trade by exporting industry/country and importing country
        #   - remove domestic trade
        IMGR_INT = [sum(Z[i:i+S-1, j]) for i in 1:S:N*S, j in 1:N*S] # N×NS
        IMGR_FNL = [sum(Y[i:i+S-1, j]) for i in 1:S:N*S, j in 1:N] # N×N

        for i in 1:N*S
            IMGR_INT[ceil(Int, i/S), i] = 0.0
            IMGR_FNL[ceil(Int, i/S), ceil(Int, i/S)] = 0.0
        end

        # cannot add up since final demand is not differentiated along industries

        return IMGR_INT, IMGR_FNL

    else

        # bilateral trade by exporting country and importing country
        #   - remove domestic trade
        IMGR_INT = [sum(Z[i:i+S-1, j:j+S-1]) for i in 1:S:N*S, j in 1:S:N*S] # N×N
        IMGR_FNL = [sum(Y[i:i+S-1, j]) for i in 1:S:N*S, j in 1:N] # N×N

        for i in 1:N
            IMGR_INT[i, i] = 0.0
            IMGR_FNL[i, i] = 0.0
        end

        IMGR = IMGR_INT .+ IMGR_FNL # N×N

        return IMGR, IMGR_INT, IMGR_FNL

    end

end


text:
abstract type AbstractApproximator end
include("PolynomialBasis.jl")
include("LinearApproximator.jl")


text:
### Tests all functionality provided by the oackage.

#Fetches the packages required to run the tests.
using BRNEquilibrate
using DifferentialEquations
using Random
using Test

#Runs the tests.
@time begin
  @time @testset "Model Macro" begin include("reaction_networks.jl") end
  @time @testset "Model Macro" begin include("test_equilibrate.jl") end
  @time @testset "Model Macro" begin include("test_bifurcations.jl") end
  @time @testset "Model Macro" begin include("test_plots.jl") end
end

max_stars_repo_name: TorkelE/BRNEquilibrate

max_issues_repo_name: TorkelE/BRNEquilibrate

max_forks_repo_name: TorkelE/BRNEquilibrate


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Example taken from
# https://github.com/JuliaIO/ImageMagick.jl/blob/master/deps/build.jl
dependencies = [
    "build_SQLite.v3.28.0.jl",
]

for elem in dependencies
    # it's a bit faster to run the build in an anonymous module instead of
    # starting a new julia process
    m = Module(:__anon__)
    Core.include(m, (joinpath(@__DIR__, elem)))
end

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libproj"], :libproj),
    ExecutableProduct(prefix, "cct", :cct_path),
    ExecutableProduct(prefix, "cs2cs", :cs2cs_path),
    ExecutableProduct(prefix, "geod", :geod_path),
    ExecutableProduct(prefix, "gie", :gie_path),
    ExecutableProduct(prefix, "proj", :proj_path),
    ExecutableProduct(prefix, "projinfo", :projinfo_path),
    FileProduct(prefix, "share/proj/CH", :ch_path),
    FileProduct(prefix, "share/proj/GL27", :gl27_path),
    FileProduct(prefix, "share/proj/ITRF2000", :itrf2000_path),
    FileProduct(prefix, "share/proj/ITRF2008", :itrf2008_path),
    FileProduct(prefix, "share/proj/ITRF2014", :itrf2014_path),
    FileProduct(prefix, "share/proj/nad.lst", :nad_lst_path),
    FileProduct(prefix, "share/proj/nad27", :nad27_path),
    FileProduct(prefix, "share/proj/nad83", :nad83_path),
    FileProduct(prefix, "share/proj/null", :null_path),
    FileProduct(prefix, "share/proj/other.extra", :other_extra_path),
    FileProduct(prefix, "share/proj/proj.db", :proj_db_path),
    FileProduct(prefix, "share/proj/world", :world_path),
    FileProduct(prefix, "share/proj/alaska", :alaska_path),
    FileProduct(prefix, "share/proj/conus", :conus_path),
    FileProduct(prefix, "share/proj/egm96_15.gtx", :egm96_15_path),
    FileProduct(prefix, "share/proj/ntv1_can.dat", :ntv1_can_path),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaGeo/PROJBuilder/releases/download/v6.1.0-1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/PROJ.v6.1.0.aarch64-linux-gnu.tar.gz", "424e8a5e50ad0ea3458c4fa4245ee67570bb641887a3b7cab4c9fa0809bd0f22"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/PROJ.v6.1.0.aarch64-linux-musl.tar.gz", "0f4ca2cb518cb6d40a30ca47ac40447bfb6dd7496cc05944e704510d5ddcc9ef"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/PROJ.v6.1.0.arm-linux-gnueabihf.tar.gz", "7d92b09644616e8ae354312beb8e33a751017806f38d5a1c98aee8df598c8304"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/PROJ.v6.1.0.arm-linux-musleabihf.tar.gz", "867ceb4c28f5fd13761c01f15b6def155a45e181f5fa974827b8aca041284c37"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/PROJ.v6.1.0.i686-linux-gnu.tar.gz", "aeb5e670b395e2f60e5f4f35c2a95c997bdbe1e4c4520e7e76dc10f2cc1072d3"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/PROJ.v6.1.0.i686-linux-musl.tar.gz", "5c8f003938f5ac0785afd41106f7b4e0d2500a1fa2ddf7369887498fd47a3cad"),
    # removed compiler_abi as suggested in https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/407#issuecomment-473688254
    # such that GCC4 platforms will also pick up this GCC7 build, ref https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/407
    Windows(:i686) => ("$bin_prefix/PROJ.v6.1.0.i686-w64-mingw32-gcc7.tar.gz", "9e1a042a74861263295a8be06630a7554bbc3de944bd73ec6bd15b3f00c107b9"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/PROJ.v6.1.0.powerpc64le-linux-gnu.tar.gz", "912998fa0f459bf43ace29eca6efe30997136cb0b21cb84a3e42e62f9e2ac9eb"),
    MacOS(:x86_64) => ("$bin_prefix/PROJ.v6.1.0.x86_64-apple-darwin14.tar.gz", "eb3252074d7513d203ce9b7d86d9370c5f6e2f183db5eb060418545147a07488"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/PROJ.v6.1.0.x86_64-linux-gnu.tar.gz", "ecd78e16e613bd888e7b74371688f52da62140795a7e23088e8bc8009a53c90d"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/PROJ.v6.1.0.x86_64-linux-musl.tar.gz", "d725d227c8636a3a4ea389f5ae385f2fdf31d4cf966dc2cef077feb047f75bfa"),
    FreeBSD(:x86_64) => ("$bin_prefix/PROJ.v6.1.0.x86_64-unknown-freebsd11.1.tar.gz", "146a3f7d91115d45fe00d79d53e701c27c5e7c241ba8a691b709274c39fdb9ab"),
    # removed compiler_abi as suggested in https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/407#issuecomment-473688254
    # such that GCC4 platforms will also pick up this GCC7 build, ref https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/407
    Windows(:x86_64) => ("$bin_prefix/PROJ.v6.1.0.x86_64-w64-mingw32-gcc7.tar.gz", "92b330c66647ac63870865546d717e095c34bb09dd11b81eb2ac0915aa49bc95"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
import RData
using RNGPool
using CoupledConditionalSMC

# data generated with A = 0.95 (not 0.9)
dataJLS = RData.load("demo/ar1data.RData")
ys = vec(dataJLS["observations"])

include("lgModel.jl")

# they start at time 0 with observations from time 1,2,...
# so we have an initial distribution of N(0,1+0.9^2)
theta = LGTheta(0.9, 1.0, 1.0, 1.0, 0.0, 1.81)
model = makeLGModel(theta, ys)
lM = LinearGaussian.makelM(theta)
ko = kalman(theta, ys)

import Statistics: mean, std

setRNGs(12345)

function runDemo(model, N, n, m, AT, AS, BS)
  println("LG Model with N = ", N, ", n = ", n)

## JLS use independent initialization

  if AT
    println("Ancestral Tracing:")
    vs = CoupledConditionalSMC.couplingTimes(model, N, n, m, true, true)
    println(mean(vs), ", ", std(vs))
    println()

    println("JLS Ancestral Tracing:")
    vs = CoupledConditionalSMC.couplingTimesJLS(model, N, n, m, true)
    println(mean(vs), ", ", std(vs))
    println()
  end

  if AS
    println("Ancestor Sampling:")
    vs = CoupledConditionalSMC.couplingTimes(model, lM, N, n, m, :AS, true, true)
    println(mean(vs), ", ", std(vs))
    println()

    println("JLS Ancestor Sampling:")
    vs = CoupledConditionalSMC.couplingTimesJLS(model, lM, N, n, m, :AS, true)
    println(mean(vs), ", ", std(vs))
    println()
  end

  if BS
    println("Backward Sampling:")
    vs = CoupledConditionalSMC.couplingTimes(model, lM, N, n, m, :BS, true, true)
    println(mean(vs), ", ", std(vs))
    println()
  end

end

m = 1000

runDemo(model, 64, 50, m, true, true, true)
runDemo(model, 128, 50, m, true, true, true)

runDemo(model, 128, 100, m, true, true, true)
runDemo(model, 256, 100, m, true, true, true)

runDemo(model, 256, 200, m, true, true, true)
runDemo(model, 512, 200, m, true, true, true)

runDemo(model, 512, 400, m, true, true, true)
runDemo(model, 1024, 400, m, true, true, true)

runDemo(model, 64, 800, m, false, true, true)
runDemo(model, 128, 800, m, false, true, true)
runDemo(model, 256, 800, m, false, true, true)
runDemo(model, 512, 800, m, false, true, true)
runDemo(model, 1024, 800, m, true, true, true)
runDemo(model, 2048, 800, m, true, true, true)


text:
module Test_bouncingBall_no_additional_values_at_events

#=
SavedValues not stored at events #557

https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/557
=#

using DifferentialEquations, PyPlot

function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -9.81
end

function outputs(u, t, integrator) # called at communication points
    return 0.5*u[1]
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
 z = u[1] + 1e-12
 return z
end

function affect_neg(integrator)
  integrator.u[2] = -0.8*integrator.u[2]
end

u0 = [1.0,0.0]
tspan = (0.0,3.0)

t_start = 0.0
t_inc   = 0.01
t_end   = 3.0
tspan   = (t_start, t_end)
tspan2  = t_start:t_inc:t_end

saved_values=SavedValues(Float64, Float64)
cb1 = SavingCallback(outputs, saved_values, saveat=tspan2)
cb2 = ContinuousCallback(condition,nothing,affect_neg! = affect_neg)

prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),saveat=t_inc,callback=CallbackSet(cb1,cb2))

figure(1)
clf()
plot(sol.t, getindex.(sol.u, 1), label="h")
plot(saved_values.t, saved_values.saveval, label="outputs")
grid(true)
legend()

println("length(sol.t)          = ", length(sol.t))
println("length(saved_values.t) = ", length(saved_values.t))
end


text:
module MadNLPMumps

import StaticArrays: SVector, setindex
import MUMPS_seq_jll
import MadNLP:
    parsefile, dlopen,
    @kwdef, Logger, @debug, @warn, @error,
    SparseMatrixCSC, SubVector, StrideOneVector,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    AbstractOptions, AbstractLinearSolver, set_options!,
    introduce, factorize!, solve!, improve!, is_inertia, inertia, findIJ, nnz

const INPUT_MATRIX_TYPE = :csc
const version = parsefile(joinpath(dirname(pathof(MUMPS_seq_jll)),"..","Project.toml"))["version"]

@kwdef mutable struct Options <: AbstractOptions
    mumps_dep_tol::Float64 = 0.
    mumps_mem_percent::Int = 1000
    mumps_permuting_scaling::Int = 7
    mumps_pivot_order::Int = 7
    mumps_pivtol::Float64 = 1e-6
    mumps_pivtolmax::Float64 = .1
    mumps_scaling::Int = 77
end

if version == "5.3.5+0"
    @kwdef mutable struct Struc
        sym::Cint = 0
        par::Cint = 0
        job::Cint = 0

        comm_fortran::Cint = 0

        icntl::SVector{60,Cint} = zeros(60)
        keep::SVector{500,Cint} = zeros(500)
        cntl::SVector{15,Cdouble} = zeros(15)
        dkeep::SVector{230,Cdouble} = zeros(230)
        keep8::SVector{150,Int64} = zeros(150)
        n::Cint = 0
        nblk::Cint = 0

        nz_alloc::Cint = 0

        nz::Cint = 0
        nnz::Int64 = 0
        irn::Ptr{Cint} = C_NULL
        jcn::Ptr{Cint} = C_NULL
        a::Ptr{Cdouble} = C_NULL

        nz_loc::Cint = 0
        nnz_loc::Int64 = 0
        irn_loc::Ptr{Cint} = C_NULL
        jcn_loc::Ptr{Cint} = C_NULL
        a_loc::Ptr{Cdouble} = C_NULL ###

        nelt::Cint = 0
        eltptr::Ptr{Cint} = C_NULL
        eltvar::Ptr{Cint} = C_NULL
        a_elt::Ptr{Cdouble} = C_NULL

        blkptr::Ptr{Cint} = C_NULL
        blkvar::Ptr{Cint} = C_NULL

        perm_in::Ptr{Cint} = C_NULL

        sym_perm::Ptr{Cint} = C_NULL
        uns_perm::Ptr{Cint} = C_NULL

        colsca::Ptr{Cdouble} = C_NULL
        rowsca::Ptr{Cdouble} = C_NULL
        colsca_from_mumps::Cint = 0
        rowsca_from_mumps::Cint = 0

        rhs::Ptr{Cdouble} = C_NULL
        redrhs::Ptr{Cdouble} = C_NULL
        rhs_sparse::Ptr{Cdouble} = C_NULL
        sol_loc::Ptr{Cdouble} = C_NULL
        rhs_loc::Ptr{Cdouble} = C_NULL

        irhs_sparse::Ptr{Cint} = C_NULL
        irhs_ptr::Ptr{Cint} = C_NULL
        isol_loc::Ptr{Cint} = C_NULL
        irhs_loc::Ptr{Cint} = C_NULL

        nrhs::Cint = 0
        lrhs::Cint = 0
        lredrhs::Cint = 0
        nz_rhs::Cint = 0
        lsol_loc::Cint = 0
        nloc_rhs::Cint = 0
        lrhs_loc::Cint = 0

        schur_mloc::Cint = 0
        schur_nloc::Cint = 0
        schur_lld::Cint = 0

        mblock::Cint = 0
        nblock::Cint = 0
        nprow::Cint = 0
        npcol::Cint = 0

        info::SVector{80,Cint} = zeros(80)
        infog::SVector{80,Cint} = zeros(80)
        rinfo::SVector{40,Cdouble} = zeros(40)
        rinfog::SVector{40,Cdouble} = zeros(40)

        deficiency::Cint = 0
        pivnul_list::Ptr{Cint} = C_NULL
        mapping::Ptr{Cint} = C_NULL

        size_schur::Cint = 0
        listvar_schur::Ptr{Cint} = C_NULL
        schur::Ptr{Cdouble} = C_NULL ##

        instance_number::Cint = 0
        wk_user::Ptr{Cdouble} = C_NULL

        version_number::SVector{32,Cchar} = zeros(32)

        ooc_tmpdir::SVector{256,Cchar} = zeros(256)
        ooc_prefix::SVector{64,Cchar} = zeros(64)

        write_problem::SVector{256,Cchar} = zeros(256)
        lwk_user::Cint = 0

        save_dir::SVector{256,Cchar} = zeros(256)
        save_prefix::SVector{256,Cchar} = zeros(256)

        metis_options::SVector{40,Cint} = zeros(40)
    end
elseif version == "5.2.1+4"
    @kwdef mutable struct Struc
        sym::Cint = 0
        par::Cint = 0
        job::Cint = 0

        comm_fortran::Cint = 0

        icntl::SVector{60,Cint} = zeros(60)
        keep::SVector{500,Cint} = zeros(500)
        cntl::SVector{15,Cdouble} = zeros(15)
        dkeep::SVector{230,Cdouble} = zeros(230)
        keep8::SVector{150,Int64} = zeros(150)
        n::Cint = 0

        nz_alloc::Cint = 0

        nz::Cint = 0
        nnz::Int64 = 0
        irn::Ptr{Cint} = C_NULL
        jcn::Ptr{Cint} = C_NULL
        a::Ptr{Cdouble} = C_NULL

        nz_loc::Cint = 0
        nnz_loc::Int64 = 0
        irn_loc::Ptr{Cint} = C_NULL
        jcn_loc::Ptr{Cint} = C_NULL
        a_loc::Ptr{Cdouble} = C_NULL ###

        nelt::Cint = 0
        eltptr::Ptr{Cint} = C_NULL
        eltvar::Ptr{Cint} = C_NULL
        a_elt::Ptr{Cdouble} = C_NULL

        perm_in::Ptr{Cint} = C_NULL

        sym_perm::Ptr{Cint} = C_NULL
        uns_perm::Ptr{Cint} = C_NULL

        colsca::Ptr{Cdouble} = C_NULL
        rowsca::Ptr{Cdouble} = C_NULL
        colsca_from_mumps::Cint = 0
        rowsca_from_mumps::Cint = 0

        rhs::Ptr{Cdouble} = C_NULL
        redrhs::Ptr{Cdouble} = C_NULL
        rhs_sparse::Ptr{Cdouble} = C_NULL
        sol_loc::Ptr{Cdouble} = C_NULL
        rhs_loc::Ptr{Cdouble} = C_NULL

        irhs_sparse::Ptr{Cint} = C_NULL
        irhs_ptr::Ptr{Cint} = C_NULL
        isol_loc::Ptr{Cint} = C_NULL
        irhs_loc::Ptr{Cint} = C_NULL

        nrhs::Cint = 0
        lrhs::Cint = 0
        lredrhs::Cint = 0
        nz_rhs::Cint = 0
        lsol_loc::Cint = 0
        nloc_rhs::Cint = 0
        lrhs_loc::Cint = 0

        schur_mloc::Cint = 0
        schur_nloc::Cint = 0
        schur_lld::Cint = 0

        mblock::Cint = 0
        nblock::Cint = 0
        nprow::Cint = 0
        npcol::Cint = 0

        info::SVector{80,Cint} = zeros(80)
        infog::SVector{80,Cint} = zeros(80)
        rinfo::SVector{40,Cdouble} = zeros(40)
        rinfog::SVector{40,Cdouble} = zeros(40)

        deficiency::Cint = 0
        pivnul_list::Ptr{Cint} = C_NULL
        mapping::Ptr{Cint} = C_NULL

        size_schur::Cint = 0
        listvar_schur::Ptr{Cint} = C_NULL
        schur::Ptr{Cdouble} = C_NULL ##

        instance_number::Cint = 0
        wk_user::Ptr{Cdouble} = C_NULL

        version_number::SVector{32,Cchar} = zeros(32)

        ooc_tmpdir::SVector{256,Cchar} = zeros(256)
        ooc_prefix::SVector{64,Cchar} = zeros(64)

        write_problem::SVector{256,Cchar} = zeros(256)
        lwk_user::Cint = 0

        save_dir::SVector{256,Cchar} = zeros(256)
        save_prefix::SVector{256,Cchar} = zeros(256)

        metis_options::SVector{40,Cint} = zeros(40)
    end
# elseif version == "4.10.0+0"
#     @kwdef mutable struct Struc
#         sym::Cint = 0
#         par::Cint = 0
#         job::Cint = 0

#         comm_fortran::Cint = 0

#         icntl::SVector{40,Cint} = zeros(40)
#         cntl::SVector{15,Cdouble} = zeros(15)

#         n::Cint = 0
#         nz_alloc::Cint = 0

#         nz::Cint = 0
#         irn::Ptr{Cint} = C_NULL
#         jcn::Ptr{Cint} = C_NULL
#         a::Ptr{Cdouble} = C_NULL

#         nz_loc::Cint = 0
#         irn_loc::Ptr{Cint} = C_NULL
#         jcn_loc::Ptr{Cint} = C_NULL
#         a_loc::Ptr{Cdouble} = C_NULL

#         nelt::Cint = 0
#         eltptr::Ptr{Cint} = C_NULL
#         eltvar::Ptr{Cint} = C_NULL
#         a_elt::Ptr{Cdouble} = C_NULL

#         perm_in::Ptr{Cint} = C_NULL

#         sym_perm::Ptr{Cint} = C_NULL
#         uns_perm::Ptr{Cint} = C_NULL

#         colsca::Ptr{Cdouble} = C_NULL
#         rowsca::Ptr{Cdouble} = C_NULL

#         rhs::Ptr{Cdouble} = C_NULL
#         redrhs::Ptr{Cdouble} = C_NULL
#         rhs_sparse::Ptr{Cdouble} = C_NULL
#         sol_loc::Ptr{Cdouble} = C_NULL

#         irhs_sparse::Ptr{Cint} = C_NULL
#         irhs_ptr::Ptr{Cint} = C_NULL
#         isol_loc::Ptr{Cint} = C_NULL

#         nrhs::Cint = 0
#         lrhs::Cint = 0
#         lredrhs::Cint = 0
#         nz_rhs::Cint = 0
#         lsol_loc::Cint = 0

#         schur_mloc::Cint = 0
#         schur_nloc::Cint = 0
#         schur_lld::Cint = 0

#         mblock::Cint = 0
#         nblock::Cint = 0
#         nprow::Cint = 0
#         npcol::Cint = 0

#         info::SVector{40,Cint} = zeros(40)
#         infog::SVector{40,Cint} = zeros(40)
#         rinfo::SVector{40,Cdouble} = zeros(40)
#         rinfog::SVector{40,Cdouble} = zeros(40)

#         deficiency::Cint = 0
#         pivnul_list::Ptr{Cint} = C_NULL
#         mapping::Ptr{Cint} = C_NULL

#         size_schur::Cint = 0
#         listvar_schur::Ptr{Cint} = C_NULL
#         schur::Ptr{Cdouble} = C_NULL

#         instance_number::Cint = 0
#         wk_user::Ptr{Cdouble} = C_NULL

#         version_number::SVector{16,Cchar} = zeros(16)

#         ooc_tmpdir::SVector{256,Cchar} = zeros(256)
#         ooc_prefix::SVector{64,Cchar} = zeros(64)

#         write_problem::SVector{256,Cchar} = zeros(256)
#         lwk_user::Cint = 0
#     end
else
    error("MUMPS_seq_jll version not supported")
end

mutable struct Solver <: AbstractLinearSolver
    csc::SparseMatrixCSC{Float64,Int32}
    I::Vector{Int32}
    J::Vector{Int32}
    sym_perm::Vector{Int32}
    pivnul_list::Vector{Int32}
    mumps_struc::Struc
    is_singular::Bool
    opt::Options
    logger::Logger
end

dmumps_c(mumps_struc::Struc)=ccall(
    (:dmumps_c,MUMPS_seq_jll.libdmumps),
    Cvoid,
    (Ref{Struc},),
    mumps_struc)

# this is necessary, when multi-threaded calls are made with Mumps, not to clash with MPI
mumps_lock = Threads.SpinLock()
function locked_dmumps_c(mumps_struc::Struc)
    lock(mumps_lock)
    try
        dmumps_c(mumps_struc)
    finally
        unlock(mumps_lock)
    end
end
# ---------------------------------------------------------------------------------------

function Solver(csc::SparseMatrixCSC{Float64,Int32};
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(),logger=Logger(),
                kwargs...)

    set_options!(opt,option_dict,kwargs)

    I,J = findIJ(csc)
    sym_perm = zeros(Int32,csc.n)
    pivnul_list = zeros(Int32,csc.n)

    mumps_struc = Struc()

    mumps_struc.sym =  2
    mumps_struc.par =  1
    mumps_struc.job = -1
    mumps_struc.comm_fortran = -987654 # MPI.COMM_WORLD.val

    locked_dmumps_c(mumps_struc)
    mumps_struc.n = csc.n;
    mumps_struc.nz= nnz(csc);
    mumps_struc.a = pointer(csc.nzval)
    mumps_struc.irn = pointer(I)
    mumps_struc.jcn = pointer(J)
    mumps_struc.sym_perm = pointer(sym_perm)
    mumps_struc.pivnul_list = pointer(pivnul_list)

    # symbolic factorization
    mumps_struc.job = 1;

    mumps_struc.icntl = setindex(mumps_struc.icntl,0,2)
    mumps_struc.icntl = setindex(mumps_struc.icntl,0,3)
    mumps_struc.icntl = setindex(mumps_struc.icntl,0,4)
    mumps_struc.icntl = setindex(mumps_struc.icntl,opt.mumps_permuting_scaling,6)
    mumps_struc.icntl = setindex(mumps_struc.icntl,opt.mumps_pivot_order,7)
    mumps_struc.icntl = setindex(mumps_struc.icntl,opt.mumps_scaling,8)
    mumps_struc.icntl = setindex(mumps_struc.icntl,0,10)
    mumps_struc.icntl = setindex(mumps_struc.icntl,1,13)
    mumps_struc.icntl = setindex(mumps_struc.icntl,opt.mumps_mem_percent,14)

    mumps_struc.cntl = setindex(mumps_struc.cntl,opt.mumps_pivtol,1)

    a = copy(csc.nzval) # would there be a better way?
    csc.nzval.=1

    locked_dmumps_c(mumps_struc);
    mumps_struc.info[1] < 0 && throw(SymbolicException())

    csc.nzval.=a

    M = Solver(csc,I,J,sym_perm,pivnul_list,mumps_struc,false,opt,logger)
    finalizer(finalize,M)

    return M
end

function factorize!(M::Solver)
    M.is_singular = false
    M.mumps_struc.job = 2;
    cnt = 0
    while true
        locked_dmumps_c(M.mumps_struc)
        if M.mumps_struc.info[1] in [-8,-9]
            cnt >= 10 && throw(FactorizationException())
            M.mumps_struc.icntl = setindex(M.mumps_struc.icntl,M.mumps_struc.icntl[14]*2.,14)
            cnt += 1
        elseif M.mumps_struc.info[1] == -10
            M.is_singular = true
            break
        elseif M.mumps_struc.info[1] < 0
            throw(FactorizationException())
        else
            break
        end
    end
    return M
end

function solve!(M::Solver,rhs::StrideOneVector{Float64})
    M.is_singular && return rhs
    M.mumps_struc.rhs = pointer(rhs)
    M.mumps_struc.job = 3
    locked_dmumps_c(M.mumps_struc)
    M.mumps_struc.info[1] < 0 && throw(SolveException())
    return rhs
end

is_inertia(::Solver) = true
function inertia(M::Solver)
    return (M.csc.n-M.is_singular-M.mumps_struc.infog[12],
            M.is_singular,
            M.mumps_struc.infog[12])
end

function improve!(M::Solver)
    if M.mumps_struc.cntl[1] == M.opt.mumps_pivtolmax
        @debug(M.logger,"improve quality failed.")
        return false
    end
    M.mumps_struc.cntl = setindex(M.mumps_struc.cntl,min(M.opt.mumps_pivtolmax,M.mumps_struc.cntl[1]^.5),1)
    @debug(M.logger,"improved quality: pivtol = $(M.mumps_struc.cntl[1])")
    return true
end

function finalize(M::Solver)
    M.mumps_struc.job = -2
    locked_dmumps_c(M.mumps_struc);
end

introduce(::Solver)="mumps"

end # module


text:
"""
六边型相关的功能
"""

struct AbstractHexagon{T} <: AbstractPolygon{T}
    vertex :: Array{Point2D, 1}
    edges :: Array{Segment, 1}
    center :: Point2D
end

"""
创建一个等边的六角型，平行与坐标轴
"""
function EqHexagon(center::Point2D, edgelen)
    cnx, cny = center.x, center.y
    xlim = edgelen / 2.0
    ylim = sqrt(3) * xlim
    verts = [
        Point2D(cnx+2xlim, cny),
        Point2D(cnx+xlim, cny+ylim),
        Point2D(cnx-xlim, cny+ylim),
        Point2D(cnx-2xlim, cny),
        Point2D(cnx-xlim, cny-ylim),
        Point2D(cnx+xlim, cny-ylim)
    ]
    edges = [Segment(verts[idx], verts[idx+1]) for idx in 1:1:5]
    push!(edges, Segment(verts[6], verts[1]))
    return AbstractHexagon{:EQ}(verts, edges, center)
end

"""
六边型的面积
"""
function area(hexa::AbstractHexagon{N}) where N
    if N == :EQ
        return (hexa.edges[1].length)^2 * 3sqrt(3) / 2
    end
    throw(error("not implement"))
end

#display函数在中间加::MIME"text/plain"
function Base.show(io::IO, hexa::AbstractHexagon{N}) where N
    name = Dict(
        :EQ => "EqHexagon",
    )
    @printf(io, "%s\n center: %s\n edgelen: %.6f",
    name[N], hexa.center, hexa.edges[1].length)
end


text:
#=================================
Methods using Interval Arithmetic
=================================#
using IntervalArithmetic

"""
    enclose_IntervalArithmetic(f::Function, dom::Interval_or_IntervalBox)

Compute a range enclosure using interval arithmetic substitution.

### Input

- `f`   -- function
- `dom` -- hyperrectangular domain, either a unidimensional  `Interval` or
           a multidimensional `IntervalBox`

### Output

An interval representing the range enclosure (minimum and maximum) of `f` over
its domain `dom`.

### Algorithm

The result is obtained by substitution of `dom` on the function `f` and
application of the rules of interval arithmetic.

We refer to the documentation and source code of `IntervalArithmetic` for details.
"""
function enclose_IntervalArithmetic(f::Function, dom::Interval_or_IntervalBox)
    return f(dom...)
end

#===================================
Methods using Interval Optimisation
===================================#
using IntervalOptimisation

"""
    range_and_optimisers(f::Function, dom::Interval_or_IntervalBox;
                         structure=HeapedVector, tol=1e-3)

Return an interval that contains the global minimum (resp. an interval that
contains the global maximum) of a univariate or multivariate function `f`,
as well as a list of boxes containing the minimisers (resp. minimisers).

### Input

- `f`         -- function
- `dom`       -- hyperrectangular domain, either a unidimensional  `Interval` or
                 a multidimensional `IntervalBox`
- `structure` -- (optional, default: `HeapedVector`) the way in which vector elements
                 are kept arranged; possible options are `HeapedVector` and `SortedVector`
- `tol`       -- tolerance to which the optima are computed

### Output

The 4-tuple `(global_min, minimisers, global_max, maximisers)`, where:

- `global_min`: interval that is guaranteed to contain the global minimum of `f`
                over `dom`
- `minimisers`: list of boxes containing the minimisers of `f` over `dom`
- `global_max`: interval that is guaranteed to contain the global maximum of `f`
                over `dom`
- `maximisers`: list of boxes containing the maximisers of `f` over `dom`

### Algorithm

The solver finds the global minimum of the function `f` over the `Interval` or
`IntervalBox` `dom` using the Moore-Skelboe algorithm.

We refer to the documentation and source code of `IntervalOptimisation` for details.
"""
function range_and_optimisers(f::Function, dom::Interval_or_IntervalBox;
                              structure=HeapedVector, tol=1e-3)
    _f = wrap_f(f, dom)
    return _range_and_optimisers(_f, dom, structure, tol)
end

@inline wrap_f(f::Function, dom::Interval) = f

# for multidimensional domains, it is required that f has a single argument
@inline wrap_f(f::Function, dom::IntervalBox) = X -> f(X...)

function _range_and_optimisers(f, dom::Interval_or_IntervalBox, structure, tol)
    global_min, minimisers = minimise(f, dom, structure=structure, tol=tol)
    global_max, maximisers = maximise(f, dom, structure=structure, tol=tol)
    return global_min, minimisers, global_max, maximisers
end

"""
    enclose_IntervalOptimisation(f::Function, dom::Interval_or_IntervalBox;
                                 structure=HeapedVector, tol=1e-3)

Compute a range enclosure using interval optimisation.

### Input

- `f`         -- function
- `dom`       -- hyperrectangular domain, either a unidimensional  `Interval` or
                 a multidimensional `IntervalBox`
- `structure` -- (optional, default: `HeapedVector`) the way in which vector elements
                 are kept arranged; possible options are `HeapedVector` and `SortedVector`
- `tol`       -- tolerance to which the optima are computed

### Output

An interval representing the range enclosure (minimum and maximum) of `f` over
its domain `dom`.

### Algorithm

The solver finds the global minimum of the function `f` over the `Interval` or
`IntervalBox` `dom` using the Moore-Skelboe algorithm. The method actually returns
an *interval* containing the global minimum (resp. an *interval* containing the
global maximum), as well as a list of boxes containing the minimisers
(resp. maximisers); they can be obtained calling [`range_and_optimisers`](@ref).

This function conservatively chooses the infimum (resp. supremum) over each interval.

We refer to the documentation and source code of `IntervalOptimisation` for details
on the implemented methods.
"""
function enclose_IntervalOptimisation(f::Function, dom::Interval_or_IntervalBox;
                                      structure=HeapedVector, tol=1e-3)

    global_min, _, global_max, _ = range_and_optimisers(f, dom,
                                                        structure=structure, tol=tol)

    return Interval(inf(global_min), sup(global_max))
end