text:
module StencilBuilder

export @fast3d, @fd

## Stencil builder
function make_index(a, arrays, i, j, k)
    if a in arrays
        i += (a in [ Symbol("u"), Symbol("ut")]) ? 0.5 : 0
        j += (a in [ Symbol("v"), Symbol("vt")]) ? 0.5 : 0
        k += (a in [ Symbol("w"), Symbol("wt")]) ? 0.5 : 0

        if i < 0
            i_int = convert(Int, abs(i))
            ex_i = :( i-$i_int )
        elseif i > 0
            i_int = convert(Int, i)
            ex_i = :( i+$i_int )
        else
            ex_i = :( i )
        end

        if j < 0
            j_int = convert(Int, abs(j))
            ex_j = :( j-$j_int )
        elseif j > 0
            j_int = convert(Int, j)
            ex_j = :( j+$j_int )
        else
            ex_j = :( j )
        end

        if k < 0
            k_int = convert(Int, abs(k))
            ex_k = :( k-$k_int )
        elseif k > 0
            k_int = convert(Int, k)
            ex_k = :( k+$k_int )
        else
            ex_k = :( k )
        end

        return :( $a[ $ex_i, $ex_j, $ex_k] )
    else
        return :( $a )
    end
end

function process_expr(ex, arrays, i, j, k)
    ci1 = 1//2; ci2 = 1//2;
    cg1 = -1; cg2 = 1;

    n = 1

    if (isa(ex.args[1], Symbol) && ex.args[1] == Symbol("gradx"))
        ex.args[1] = Symbol("gradx_")
        ex = :( $ex * dxi )
    elseif (isa(ex.args[1], Symbol) && ex.args[1] == Symbol("grady"))
        ex.args[1] = Symbol("grady_")
        ex = :( $ex * dyi )
    elseif (isa(ex.args[1], Symbol) && ex.args[1] == Symbol("gradz"))
        ex.args[1] = Symbol("gradz_")
        if isinteger(k)
            k_int = convert(Int, abs(k))
            if k > 0
                ex = :( $ex * dzi[k+$k_int] )
            elseif k < 0
                ex = :( $ex * dzi[k-$k_int] )
            else
                ex = :( $ex * dzi[k] )
            end
        else
            k_int = convert(Int, abs(k + 1/2))
            if k > -1/2
                ex = :( $ex * dzhi[k+$k_int] )
            elseif k < -1/2
                ex = :( $ex * dzhi[k-$k_int] )
            else
                ex = :( $ex * dzhi[k] )
            end
        end
    end

    args = ex.args
    while n <= length(args)
        if isa(args[n], Expr)
            args[n] = process_expr(args[n], arrays, i, j, k)
            n += 1
        elseif isa(args[n], Symbol)
            if args[n] == Symbol("gradx_")
                if isa(args[n+1], Expr)
                    args[n] = copy(args[n+1])
                    args[n  ] = process_expr(args[n  ], arrays, i-0.5, j, k)
                    args[n+1] = process_expr(args[n+1], arrays, i+0.5, j, k)
                elseif isa(args[n+1], Symbol)
                    args[n] = args[n+1]
                    args[n  ] = make_index(args[n  ], arrays, i-0.5, j, k)
                    args[n+1] = make_index(args[n+1], arrays, i+0.5, j, k)
                end
                args[n  ] = :( $cg1 * $(args[n  ])  )
                args[n+1] = :( $cg2 * $(args[n+1])  )
                insert!(args, n, Symbol("+"))
                n += 3
            elseif args[n] == Symbol("grady_")
                if isa(args[n+1], Expr)
                    args[n] = copy(args[n+1])
                    args[n  ] = process_expr(args[n  ], arrays, i, j-0.5, k)
                    args[n+1] = process_expr(args[n+1], arrays, i, j+0.5, k)
                elseif isa(args[n+1], Symbol)
                    args[n] = args[n+1]
                    args[n  ] = make_index(args[n  ], arrays, i, j-0.5, k)
                    args[n+1] = make_index(args[n+1], arrays, i, j+0.5, k)
                end
                args[n  ] = :( $cg1 * $(args[n  ])  )
                args[n+1] = :( $cg2 * $(args[n+1])  )
                insert!(args, n, Symbol("+"))
                n += 3
            elseif args[n] == Symbol("gradz_")
                if isa(args[n+1], Expr)
                    args[n] = copy(args[n+1])
                    args[n  ] = process_expr(args[n  ], arrays, i, j, k-0.5)
                    args[n+1] = process_expr(args[n+1], arrays, i, j, k+0.5)
                elseif isa(args[n+1], Symbol)
                    args[n] = args[n+1]
                    args[n  ] = make_index(args[n  ], arrays, i, j, k-0.5)
                    args[n+1] = make_index(args[n+1], arrays, i, j, k+0.5)
                end
                args[n  ] = :( $cg1 * $(args[n  ])  )
                args[n+1] = :( $cg2 * $(args[n+1])  )
                insert!(args, n, Symbol("+"))
                n += 3
            elseif args[n] == Symbol("interpx")
                if isa(args[n+1], Expr)
                    args[n] = copy(args[n+1])
                    args[n  ] = process_expr(args[n  ], arrays, i-0.5, j, k)
                    args[n+1] = process_expr(args[n+1], arrays, i+0.5, j, k)
                elseif isa(args[n+1], Symbol)
                    args[n] = args[n+1]
                    args[n  ] = make_index(args[n  ], arrays, i-0.5, j, k)
                    args[n+1] = make_index(args[n+1], arrays, i+0.5, j, k)
                end
                args[n  ] = :( $ci1 * $(args[n  ])  )
                args[n+1] = :( $ci2 * $(args[n+1])  )
                insert!(args, n, Symbol("+"))
                n += 3
            elseif args[n] == Symbol("interpy")
                if isa(args[n+1], Expr)
                    args[n] = copy(args[n+1])
                    args[n  ] = process_expr(args[n  ], arrays, i, j-0.5, k)
                    args[n+1] = process_expr(args[n+1], arrays, i, j+0.5, k)
                elseif isa(args[n+1], Symbol)
                    args[n] = args[n+1]
                    args[n  ] = make_index(args[n  ], arrays, i, j-0.5, k)
                    args[n+1] = make_index(args[n+1], arrays, i, j+0.5, k)
                end
                args[n  ] = :( $ci1 * $(args[n  ])  )
                args[n+1] = :( $ci2 * $(args[n+1])  )
                insert!(args, n, Symbol("+"))
                n += 3
            elseif args[n] == Symbol("interpz")
                if isa(args[n+1], Expr)
                    args[n] = copy(args[n+1])
                    args[n  ] = process_expr(args[n  ], arrays, i, j, k-0.5)
                    args[n+1] = process_expr(args[n+1], arrays, i, j, k+0.5)
                elseif isa(args[n+1], Symbol)
                    args[n] = args[n+1]
                    args[n  ] = make_index(args[n  ], arrays, i, j, k-0.5)
                    args[n+1] = make_index(args[n+1], arrays, i, j, k+0.5)
                end
                args[n  ] = :( $ci1 * $(args[n  ])  )
                args[n+1] = :( $ci2 * $(args[n+1])  )
                insert!(args, n, Symbol("+"))
                n += 3
            else
                args[n] = make_index(args[n], arrays, i, j, k)
                n += 1
            end
        else
            n += 1
        end
    end
    return ex
end

macro fd(arrays, ex)
    i = (ex.args[1] in [ Symbol("u"), Symbol("ut")]) ? -0.5 : 0
    j = (ex.args[1] in [ Symbol("v"), Symbol("vt")]) ? -0.5 : 0
    k = (ex.args[1] in [ Symbol("w"), Symbol("wt")]) ? -0.5 : 0

    if isa(arrays, Symbol)
        arrays = :( [ $arrays ] )
    end

    ex = process_expr(ex, arrays.args, i, j, k)

    println("Generated stencil: ")
    println(ex)
    println("")

    return esc(ex)
end

macro fast3d(ex)
    ex_loop = quote
        @tturbo unroll=8 for k in ks:ke
            for j in js:je
                for i in is:ie
                    $ex
                end
            end
        end
    end
    return esc(ex_loop)
end

end


text:
module CodeMathOverrides
  #
  # Julia requires these imports
  #
  import Base.+
  import Base.-
  import Base.*
  import Base./
  import Base.&
  import Base.|
  import Base.$
  import Base.>>
  import Base.<<
  import Base.%

  export +
  export -
  export *
  export /
  export &
  export |
  export $
  export >>
  export <<
  export %
  #
  # Shorthand for Float16 type
  #
  const f16 = Float16::DataType
  #
  # Max positive Float16 value
  #
  const F16_MAX = f16(65504.0)
  #
  # Max negative Float16 value
  #
  const F16_MIN = f16(-65504.0)
  #
  # Max positive Float64 value
  #
  const F16_MAX_FLOAT = 65504.0::Float64
  #
  # Max negative Float64 value
  #
  const F16_MIN_FLOAT = -65504.0::Float64
  #
  # Max positive Int value
  #
  const F16_MAX_INT = 65504::Int
  #
  # Max negative Int value
  #
  const F16_MIN_INT = -65504::Int
  #
  # Override for + operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function +(left::f16, right::f16)
    local res::Float64 = Float64(left) + Float64(right)

    if res > F16_MAX_FLOAT return F16_MAX
    elseif res < F16_MIN_FLOAT return F16_MIN end
    #
    # We have to use + operator for Float64, because
    # + for Float16 we have already overridden. This is
    # how we prevent StackOverflow error
    #
    f16(res)
  end
  #
  # Override for - operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function -(left::f16, right::f16)
    local res::Float64 = Float64(left) - Float64(right)

    if res > F16_MAX_FLOAT return F16_MAX
    elseif res < F16_MIN_FLOAT return F16_MIN end
    #
    # We have to use - operator for Float64, because
    # - for Float16 we have already overridden. This is
    # how we prevent StackOverflow error
    #
    f16(res)
  end
  #
  # Override for * operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function *(left::f16, right::f16)
    local res::Float64 = Float64(left) * Float64(right)

    if res > F16_MAX_FLOAT return F16_MAX
    elseif res < F16_MIN_FLOAT return F16_MIN end
    #
    # We have to use * operator for Float64, because
    # * for Float16 we have already overridden. This is
    # how we prevent StackOverflow error
    #
    f16(res)
  end
  #
  # Override for / operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function /(left::f16, right::f16)
    local res::Float64 = Float64(left) / Float64(right)

    if res > F16_MAX_FLOAT return F16_MAX
    elseif res < F16_MIN_FLOAT return F16_MIN end
    #
    # We have to use / operator for Float64, because
    # / for Float16 we have already overridden. This is
    # how we prevent StackOverflow error
    #
    f16(res)
  end
  #
  # Override for & operator for Float16 type. We have to
  # convert Float16 values into Int, because & operator
  # works only with integers
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function (&)(left::f16, right::f16)
    f16(round(Int, left) & round(Int, right))
  end
  #
  # Override for | operator for Float16 type. We have to
  # convert Float16 values into Int, because | operator
  # works only with integers
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function (|)(left::f16, right::f16)
    f16(round(Int, left) | round(Int, right))
  end
  #
  # Override for $ operator for Float16 type. We have to
  # convert Float16 values into Int, because $ operator
  # works only with integers
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function ($)(left::f16, right::f16)
    f16(round(Int, left) $ round(Int, right))
  end
  #
  # Override for >> operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function (>>)(left::f16, right::f16)
    local res::Int = round(Int, left) >> round(Int, right)

    if res > F16_MAX_INT return F16_MAX
    elseif res < F16_MIN_INT return F16_MIN end

    f16(res)
  end
  #
  # Override for << operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function (<<)(left::f16, right::f16)
    local res::Int = round(Int, left) << round(Int, right)

    if res > F16_MAX_INT return F16_MAX
    elseif res < F16_MIN_INT return F16_MIN end

    f16(res)
  end
  #
  # Override for << operator for Float16 type. Checks type
  # overflow and returns max or min value in this case.
  # @param left Left operand
  # @param right Right operand
  # @return {Float16}
  #
  function (%)(left::f16, right::f16)
    local r::Int = round(Int, right)
    f16(round(Int, left) % (r === 0 ? 1 : r))
  end
end

ext: jl

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"""
    package_devpath(pkg::Module)
Returns the path to the top level directory of a devved out package source tree. The package
is identified by its top level module `pkg`.
"""
function package_devpath(pkg::Module)
    pkg == parentmodule(pkg) || throw(ArgumentError("$(pkg) is a submodule. Use the package top-level module."))
    path = pathof(pkg)
    path === nothing && throw(ArgumentError("could not find path to $(pkg)."))
    name = String(nameof(pkg))

    # check that pkg is not originating from a standard installation directory
    # since those are supposed to be immutable.
    for depot in DEPOT_PATH
        sep = Sys.iswindows() ? "\\\\" : "/"
        if startswith(path, joinpath(depot, "packages", name)) &&
            occursin(Regex(name * sep * "\\w{4,5}" * sep * "src" * sep * name * ".jl"), path)
            throw(ArgumentError(string(
                "module $(name) was found in a standard installation directory. ",
                "Please make sure that $(name) is ready for development by running ",
                "`pkg> develop $(name)` from the Pkg REPL, or ",
                "`Pkg.develop(\"$(name)\")` from the Julia REPL, and try again.")))
        end
    end
    # We assume that the path to source file of pkg is ../Package/src/Package.jl, but we
    # return simply the top level directory of the package (i.e. ../Package)
    return normpath(joinpath(path, "..", ".."))
end


text:
abstract type AbstractProjected{O,Sp,Sa} <: AbstractSampled{O,Sp,Sa} end

"""
    Projected(order::Order, span, sampling, crs, usercrs)
    Projected(; order=Ordered(), span=AutoSpan(), sampling=Points(), crs, usercrs=nothing)

An [`AbstractSampled`]($DDabssampleddocs) `IndexMode` with projections attached.

Fields and behaviours are identical to [`Sampled`]($DDsampleddocs)
with the addition of `crs` and `usercrs` fields.

If both `crs` and `usercrs` fields contain CRS data (in a `GeoFormat` wrapper
from GeoFormatTypes.jl) the selector inputs and plot axes will be converted
from and to the specified `usercrs` projection automatically. A common use case
would be to pass `usercrs=EPSG(4326)` to the constructor when loading eg. a GDALarray:

```julia
GDALarray(filename; usercrs=EPSG(4326))
```

The underlying `crs` will be detected by GDAL.

If `usercrs` is not supplied (ie. `usercrs=nothing`), the base index will be
shown on plots, and selectors will need to use whatever format it is in.
"""
struct Projected{O<:Order,Sp<:Regular,Sa,C,IC} <: AbstractProjected{O,Sp,Sa}
    order::O
    span::Sp
    sampling::Sa
    crs::C
    usercrs::IC
end
Projected(; order=Ordered(), span=Regular(),
          sampling=Points(), crs, usercrs=nothing) =
    Projected(order, span, sampling, crs, usercrs)

crs(mode::Projected, dim) = crs(mode)
crs(mode::Projected) = mode.crs

usercrs(mode::Projected, dim) = usercrs(mode)
usercrs(mode::Projected) = mode.usercrs

rebuild(g::Projected, order=order(g), span=span(g),
        sampling=sampling(g), crs=crs(g), usercrs=usercrs(g)) =
    Projected(order, span, sampling, crs, usercrs)

"""
    Converted(order::Order, span, sampling, crs, dimcrs)
    Converted(; order=Ordered(), span=AutoSpan(), sampling=Points(), crs, dimcrs)

An [`AbstractSampled`]($DDabssampleddocs) `IndexMode` with projections, where the
dimension index has already been converted to another projection as a vector,
usually `EPSG(4326)`.

Fields and behaviours are identical to [`Sampled`]($DDsampleddocs) with the addition of
`crs` and `dimcrs` fields.

The dimension will be indexed as for [`Sampled`]($DDsampleddocs), but to save in another
format the underlying projection will be used.  ```julia GDALarray(filename; usercrs=EPSG(4326))
```

The underlying `crs` will be detected by GDAL.

If `usercrs` is not supplied (ie. `usercrs=nothing`), the base index will be
shown on plots, and selectors will need to use whatever format it is in.
"""
struct Converted{O<:Order,Sp,Sa,C,DC} <: AbstractProjected{O,Sp,Sa}
    order::O
    span::Sp
    sampling::Sa
    crs::C
    dimcrs::DC
end
Converted(; order=Ordered(), span=AutoSpan(),
          sampling=Points(), crs, dimcrs) =
    Converted(order, span, sampling, crs, dimcrs)

crs(mode::Converted, dim) = crs(mode)
crs(mode::Converted) = mode.crs

dimcrs(mode::Converted, dim) = dimcrs(mode)
dimcrs(mode::Converted) = mode.dimcrs

rebuild(g::Converted, order=order(g), span=span(g),
        sampling=sampling(g), crs=crs(g), dimcrs=dimcrs(g)) =
    Converted(order, span, sampling, crs, dimcrs)

"""
    LatLon(order, span, sampling)
    LatLon(; order=Ordered(), span=AutoSpan(), sampling=Points())

An [`AbstractSampled`]($DDabssampleddocs) mode for standard latitude/longitude dimensions.
"""
struct LatLon{O<:Order,Sp,Sa} <: AbstractSampled{O,Sp,Sa}
    order::O
    span::Sp
    sampling::Sa
end
LatLon(; order=Ordered(), span=AutoSpan(), sampling=Points()) =
    LatLon(order, span, sampling)

crs(mode::LatLon, args...) = EPSG(4326)

rebuild(g::LatLon, order=order(g), span=span(g), sampling=sampling(g)) =
    LatLon(order, span, sampling)

"""
    convertmode(dstmode::Type{<:IndexMode}, x)

Convert the dimension mode between `Projected` and `Converted`.
Other dimension modes pass through unchanged.

This is used to e.g. save a netcdf file to GeoTiff.
"""
convertmode(dstmode::Type{<:IndexMode}, A::AbstractArray) =
    rebuild(A, data(A), convertmode(dstmode, dims(A)))
convertmode(dstmode::Type{<:IndexMode}, dims::Tuple) =
    map(d -> convertmode(dstmode, d), dims)
convertmode(dstmode::Type{<:IndexMode}, dim::Dimension) =
    convertmode(dstmode, basetypeof(mode(dim)), dim)
# Regular modes pass through
convertmode(dstmode::Type, srcmode::Type{<:IndexMode}, dim::Dimension) = dim
# AbstractProjected passes through if it's the same as dstmode
convertmode(dstmode::Type{M}, srcmode::Type{M}, dim::Dimension) where M<:AbstractProjected = dim
# Otherwise AbstractProjected needs ArchGDAL
convertmode(dstmode::Type, srcmode::Type{<:AbstractProjected}, dim::Dimension) =
    error("Load ArchGDAL.jl to convert projected dimensions")


text:
"A grid corresponding to an unstructured collection of points."
struct ScatteredGrid{T} <: AbstractGrid{T,1}
    points     ::  Vector{T}
    domain     ::  Domain
    ScatteredGrid(points::Vector{T}, domain=DomainSets.FullSpace{T}()) where T =
        new{T}(points, domain)
end

size(g::ScatteredGrid) = (length(g.points),)

unsafe_grid_getindex(g::ScatteredGrid, idx) = (@inbounds g.points[idx])

covering(g::ScatteredGrid) = g.domain

## Some functionality for random grids

"""
Compute a scattered grid of M points randomly distributed in `Ω`, using
the uniform probability measure on `Ω`.
"""
function randomgrid(Ω::Domain, M::Int)
    points = [randompoint(Ω, boundingbox(Ω)) for m in 1:M]
    ScatteredGrid(sort(points))
end

"Generate a single random point inside the given box, with `eltype` `T`."
randompoint(dom::ProductDomain) =
    DomainSets.toexternalpoint(dom, map(randompoint,components(dom)))

# Don't return an SVector in 1d, just a value
function randompoint(dom::AbstractInterval)
    T = prectype(dom)
    val = rand(T)
    convert(T,val * infimum(dom) + (1-val) * supremum(dom))
end

"""
Generate a single random point inside the given domain, with `eltype` `T`.
Random points are generated inside the given box, until one is inside the domain.
"""
function randompoint(Ω::Domain, box = boundingbox(Ω))
    local point
    indomain = false
    while ~indomain
        point = randompoint(box)
        indomain = point ∈ Ω
    end
    point
end


text:
@testset "extractor" begin

    xml = parsexml("""
    <?xml version="1.0" encoding="UTF-8"?>
    <university university-name="Julia University">
      <person id="1">
        <age>24</age>
        <study-field>Mechanical Engineering</study-field>
        <GPA>4.5</GPA>
        <taken-courses>Artificial Intelligence</taken-courses>
        <taken-courses>Robotics</taken-courses>
        He is a genius
      </person>
      <person id="2">
        <age>18</age>
        <study-field>Computer Engineering</study-field>
        <GPA>4.2</GPA>
        <taken-courses>Julia</taken-courses>
      </person>
    </university>
    """)

    # extract University (root of the document)
    U = University(xml)

    @test U.name == "Julia University"

    # extract Person

    P1 = U.people[1]

    @test stripall(P1.comment) == stripall("He is a genius") # TODO
    @test P1.age == 24
    @test P1.field == "Mechanical Engineering"
    @test P1.GPA == 4.5
    @test P1.courses == ["Artificial Intelligence", "Robotics"]
    @test P1.id == 1

    P2 = U.people[2]

    @test stripall(@capture_out(pprint(P2))) == stripall("""
    <person id="2">
        <age>18</age>
        <study-field>Computer Engineering</study-field>
        <GPA>4.2</GPA>
        <taken-courses>Julia</taken-courses>
      </person>
    """)

end


text:
"""
$(TYPEDEF)

Real symmetric or complex Hermitian positive semidefinite cone of squares.
"""
struct MatrixCSqr{T <: Real, R <: RealOrComplex{T}} <: ConeOfSquares{T} end

"""
$(TYPEDSIGNATURES)

The rank of the matrix cone of squares, equal to the side dimension of the matrix.
"""
vector_dim(::Type{<:MatrixCSqr{<:Real, R}}, d::Int) where R = svec_length(R, d)

mutable struct MatrixCSqrCache{T <: Real, R <: RealOrComplex{T}} <: CSqrCache{T}
    is_complex::Bool
    rt2::T
    viw_X::Matrix{R}
    viw_λ::Vector{T}
    w_λ::Vector{T}
    w_λi::Vector{T}
    ϕ::T
    ζ::T
    ζi::T
    σ::T
    ∇h::Vector{T}
    ∇2h::Vector{T}
    ∇3h::Vector{T}
    viw_λ_Δ::Vector{T}
    Δh::Matrix{T}
    Δ2h::Matrix{T}
    θ::Matrix{T}
    # TODO try delete some aux fields
    wd::Vector{T}
    wT::Matrix{T}
    w1::Matrix{R}
    w2::Matrix{R}
    w3::Matrix{R}
    w4::Matrix{R}
    α::Vector{T}
    γ::Vector{T}
    c0::T
    c4::T
    c5::T

    MatrixCSqrCache{T, R}() where {T <: Real, R <: RealOrComplex{T}} = new{T, R}()
end

function setup_csqr_cache(cone::EpiPerSepSpectral{MatrixCSqr{T, R}}) where {T, R}
    cone.cache = cache = MatrixCSqrCache{T, R}()
    cache.is_complex = (R <: Complex{T})
    cache.rt2 = sqrt(T(2))
    d = cone.d
    cache.viw_X = zeros(R, d, d)
    cache.w_λ = zeros(T, d)
    cache.w_λi = zeros(T, d)
    cache.∇h = zeros(T, d)
    cache.∇2h = zeros(T, d)
    cache.∇3h = zeros(T, d)
    cache.viw_λ_Δ = zeros(T, svec_length(d))
    cache.Δh = zeros(T, d, d)
    cache.Δ2h = zeros(T, d, svec_length(d))
    cache.θ = zeros(T, d, d)
    cache.wd = zeros(T, d)
    cache.wT = zeros(T, d, d)
    cache.w1 = zeros(R, d, d)
    cache.w2 = zeros(R, d, d)
    cache.w3 = zeros(R, d, d)
    cache.w4 = zeros(R, d, d)
    cache.α = zeros(T, d)
    cache.γ = zeros(T, d)
    return
end

function set_initial_point!(
    arr::AbstractVector,
    cone::EpiPerSepSpectral{<:MatrixCSqr},
    )
    (arr[1], arr[2], w0) = get_initial_point(cone.d, cone.h)
    @views fill!(arr[3:end], 0)
    incr = (cone.cache.is_complex ? 2 : 1)
    idx = 3
    @inbounds for i in 1:cone.d
        arr[idx] = 1
        idx += incr * i + 1
    end
    return arr
end

# TODO check whether it is faster to do chol before eigdecomp
function update_feas(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    @assert !cone.feas_updated
    cache = cone.cache
    v = cone.point[2]

    cone.is_feas = false
    if v > eps(T)
        w = viw_X = cache.viw_X
        svec_to_smat!(w, cone.w_view, cache.rt2)
        w_chol = cholesky!(Hermitian(w, :U), check = false)
        if isposdef(w_chol)
            svec_to_smat!(w, cone.w_view, cache.rt2)
            w ./= v
            viw_λ = cache.viw_λ = update_eigen!(viw_X)
            if all(>(eps(T)), viw_λ)
                cache.ϕ = h_val(viw_λ, cone.h)
                cache.ζ = cone.point[1] - v * cache.ϕ
                cone.is_feas = (cache.ζ > eps(T))
            end
        end
    end

    cone.feas_updated = true
    return cone.is_feas
end

# TODO check whether it is faster to do chol before eigdecomp
# TODO check if this is faster or slower than only using prox check
function is_dual_feas(cone::EpiPerSepSpectral{MatrixCSqr{T, R}}) where {T, R}
    u = cone.dual_point[1]
    (u < eps(T)) && return false
    @views w = cone.dual_point[3:end]

    uiw = cone.cache.w1
    if h_conj_dom_pos(cone.h)
        # use cholesky to check conjugate domain feasibility
        svec_to_smat!(uiw, w, cone.cache.rt2)
        w_chol = cholesky!(Hermitian(uiw, :U), check = false)
        isposdef(w_chol) || return false
    end

    svec_to_smat!(uiw, w, cone.cache.rt2)
    uiw ./= u
    # TODO in-place and dispatch to GLA or LAPACK.geevx! directly for efficiency
    uiw_eigen = eigen(Hermitian(uiw, :U))
    return (cone.dual_point[2] - u * h_conj(uiw_eigen.values, cone.h) > eps(T))
end

function update_grad(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    @assert !cone.grad_updated && cone.is_feas
    v = cone.point[2]
    grad = cone.grad
    cache = cone.cache
    ζi = cache.ζi = inv(cache.ζ)
    viw_λ = cache.viw_λ
    viw_X = cache.viw_X
    ∇h = cache.∇h
    h_der1(∇h, viw_λ, cone.h)
    cache.σ = cache.ϕ - dot(viw_λ, ∇h)
    @. cache.w_λ = v * viw_λ
    @. cache.w_λi = inv(cache.w_λ)

    grad[1] = -ζi
    grad[2] = -inv(v) + ζi * cache.σ
    @. cache.wd = ζi * ∇h - cache.w_λi
    mul!(cache.w1, Diagonal(cache.wd), viw_X') # TODO check efficient
    gw = mul!(cache.w2, viw_X, cache.w1)
    @views smat_to_svec!(cone.grad[3:end], gw, cache.rt2)

    cone.grad_updated = true
    return grad
end

function update_hess_aux(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    @assert !cone.hess_aux_updated
    @assert cone.grad_updated
    cache = cone.cache
    viw_λ = cache.viw_λ
    w_λi = cache.w_λi
    ∇h = cache.∇h
    ∇2h = cache.∇2h
    viw_λ_Δ = cache.viw_λ_Δ
    Δh = cache.Δh

    h_der2(∇2h, viw_λ, cone.h)

    # setup viw_λ_Δ
    rteps = sqrt(eps(T))
    idx = 1
    @inbounds for j in 1:cone.d
        λ_j = viw_λ[j]
        for i in 1:(j - 1)
            t = viw_λ[i] - λ_j
            if abs(t) < rteps
                viw_λ_Δ[idx] = 0
            else
                viw_λ_Δ[idx] = t
            end
            idx += 1
        end
        viw_λ_Δ[idx] = 0
        idx += 1
    end

    # setup Δh
    idx = 1
    @inbounds for j in 1:cone.d
        ∇h_j = ∇h[j]
        Δh[j, j] = ∇2h_j = ∇2h[j]
        for i in 1:(j - 1)
            denom = viw_λ_Δ[idx]
            if iszero(denom)
                Δh[i, j] = (∇2h[i] + ∇2h_j) / 2
            else
                Δh[i, j] = (∇h[i] - ∇h_j) / denom
            end
            idx += 1
        end
        idx += 1
    end

    ζivi = cache.ζi / cone.point[2]
    @. cache.θ = ζivi * Δh + w_λi * w_λi'
    copytri!(cache.θ, 'U')

    cone.hess_aux_updated = true
end

function update_hess(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    d = cone.d
    v = cone.point[2]
    H = cone.hess.data
    cache = cone.cache
    rt2 = cache.rt2
    ζi = cache.ζi
    σ = cache.σ
    viw_X = cache.viw_X
    viw_λ = cache.viw_λ
    ∇h = cache.∇h
    ∇2h = cache.∇2h
    wd = cache.wd
    w1 = cache.w1
    w2 = cache.w2
    ζi2 = abs2(ζi)
    ζivi = ζi / v

    # Huu
    H[1, 1] = ζi2

    # Huv
    H[1, 2] = -ζi2 * σ

    # Hvv
    @inbounds sum1 = sum(abs2(viw_λ[j]) * ∇2h[j] for j in 1:d)
    H[2, 2] = v^-2 + abs2(ζi * σ) + ζivi * sum1

    # Huw
    @. wd = -ζi * ∇h
    mul!(w1, Diagonal(wd), viw_X')
    mul!(w2, viw_X, w1)
    @views Hwu = H[3:end, 1] # use later for Hww
    @views smat_to_svec!(Hwu, w2, rt2)
    @. H[1, 3:end] = ζi * Hwu

    # Hvw
    wd .*= -ζi * σ
    @. wd -= ζivi * ∇2h * viw_λ
    mul!(w1, Diagonal(wd), viw_X')
    mul!(w2, viw_X, w1)
    @views smat_to_svec!(H[2, 3:end], w2, rt2)

    # Hww
    @views Hww = H[3:end, 3:end]
    eig_dot_kron!(Hww, cache.θ, viw_X, w1, w2, cache.w3, cache.w4, rt2)
    mul!(Hww, Hwu, Hwu', true, true)

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiPerSepSpectral{<:MatrixCSqr{T}},
    ) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    v = cone.point[2]
    cache = cone.cache
    ζi = cache.ζi
    σ = cache.σ
    viw_X = cache.viw_X
    viw_λ = cache.viw_λ
    w_λi = cache.w_λi
    ∇h = cache.∇h
    Δh = cache.Δh
    r_X = cache.w1
    w_aux = cache.w2
    w3 = cache.w3
    D_λi = Diagonal(w_λi)
    D_viw_λ = Diagonal(viw_λ)
    D_∇h = Diagonal(∇h)
    ζivi = ζi / v

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @views svec_to_smat!(r_X, arr[3:end, j], cache.rt2)
        mul!(w_aux, Hermitian(r_X, :U), viw_X)
        mul!(r_X, viw_X', w_aux)

        sum1 = sum(∇h[i] * real(r_X[i, i]) for i in 1:cone.d)
        c1 = -ζi * (p - σ * q - sum1) * ζi
        @. w_aux = ζivi * Δh * (r_X - q * D_viw_λ)
        c2 = sum(viw_λ[i] * real(w_aux[i, i]) for i in 1:cone.d)

        rmul!(r_X, D_λi)
        @. w_aux += w_λi * r_X + c1 * D_∇h
        mul!(w3, Hermitian(w_aux, :U), viw_X')
        mul!(w_aux, viw_X, w3)

        prod[1, j] = -c1
        prod[2, j] = c1 * σ - c2 + q / v / v
        @views smat_to_svec!(prod[3:end, j], w_aux, cache.rt2)
    end

    return prod
end

function update_inv_hess_aux(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    @assert !cone.inv_hess_aux_updated
    cone.hess_aux_updated || update_hess_aux(cone)
    v = cone.point[2]
    cache = cone.cache
    σ = cache.σ
    viw_X = cache.viw_X
    viw_λ = cache.viw_λ
    ∇h = cache.∇h
    wd = cache.wd
    α = cache.α
    γ = cache.γ
    ζivi = cache.ζi / v

    @views diag_θ = cache.θ[1:(1 + cone.d):end]
    @. wd = ζivi * cache.∇2h
    @. α = ∇h / diag_θ
    wd .*= viw_λ
    @. γ = wd / diag_θ

    ζ2β = abs2(cache.ζ) + dot(∇h, α)
    c0 = σ + dot(∇h, γ)
    c1 = c0 / ζ2β
    @inbounds sum1 = sum((viw_λ[i] + c1 * α[i] - γ[i]) * wd[i] for i in 1:cone.d)
    c3 = v^-2 + σ * c1 + sum1
    c4 = inv(c3 - c0 * c1)
    c5 = ζ2β * c3
    cache.c0 = c0
    cache.c4 = c4
    cache.c5 = c5

    cone.inv_hess_aux_updated = true
end

function update_inv_hess(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    Hi = cone.inv_hess.data
    cache = cone.cache
    rt2 = cache.rt2
    viw_X = cache.viw_X
    c4 = cache.c4
    wT = cache.wT
    w1 = cache.w1
    w2 = cache.w2

    # Hiuu, Hiuv, Hivv
    Hi[1, 1] = c4 * cache.c5
    Hiuv = Hi[1, 2] = c4 * cache.c0
    Hi[2, 2] = c4

    # Hiuw, Hivw
    @views HiuW = Hi[1, 3:end]
    @views γ_vec = Hi[3:end, 2]
    mul!(w2, Diagonal(cache.γ), viw_X')
    mul!(w1, viw_X, w2)
    smat_to_svec!(γ_vec, w1, rt2)
    @. Hi[2, 3:end] = c4 * γ_vec
    mul!(w2, Diagonal(cache.α), viw_X')
    mul!(w1, viw_X, w2)
    smat_to_svec!(HiuW, w1, rt2)
    @. HiuW += Hiuv * γ_vec

    # Hiww
    @views Hiww = Hi[3:end, 3:end]
    @. wT = inv(cache.θ)
    eig_dot_kron!(Hiww, wT, viw_X, w1, w2, cache.w3, cache.w4, rt2)
    mul!(Hiww, γ_vec, γ_vec', c4, true)

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiPerSepSpectral{<:MatrixCSqr{T}},
    ) where T
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    d = cone.d
    cache = cone.cache
    viw_X = cache.viw_X
    α = cache.α
    γ = cache.γ
    c0 = cache.c0
    c4 = cache.c4
    c5 = cache.c5
    r_X = Hermitian(cache.w1, :U)
    w2 = cache.w2

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @views svec_to_smat!(r_X.data, arr[3:end, j], cache.rt2)
        mul!(w2, r_X, viw_X)
        mul!(r_X.data, viw_X', w2)

        qγr = q + sum(γ[i] * r_X[i, i] for i in 1:d)
        cu = c4 * (c5 * p + c0 * qγr)
        cv = c4 * (c0 * p + qγr)

        prod[1, j] = cu + sum(α[i] * r_X[i, i] for i in 1:d)
        prod[2, j] = cv

        w_prod = r_X
        w_prod.data ./= cache.θ
        for i in 1:d
            w_prod.data[i, i] += p * α[i] + cv * γ[i]
        end
        mul!(w2, w_prod, viw_X')
        mul!(w_prod.data, viw_X, w2)
        @views smat_to_svec!(prod[3:end, j], w_prod.data, cache.rt2)
    end

    return prod
end

function update_dder3_aux(cone::EpiPerSepSpectral{<:MatrixCSqr{T}}) where T
    @assert !cone.dder3_aux_updated
    cone.hess_aux_updated || update_hess_aux(cone)
    d = cone.d
    cache = cone.cache
    viw_λ = cache.viw_λ
    ∇3h = cache.∇3h
    viw_λ_Δ = cache.viw_λ_Δ
    Δh = cache.Δh
    Δ2h = cache.Δ2h

    h_der3(∇3h, viw_λ, cone.h)

    # setup Δ2h
    i6 = inv(T(6))
    idx_jk = 0
    @inbounds for k in 1:d
        ∇3h_k = ∇3h[k]
        idx_ij = 1
        kd = d * (k - 1)
        for j in 1:k
            ∇3h_j = ∇3h[j]
            Δh_jk = Δh[j, k]
            ijk = d * (idx_jk - 1 + j)
            jik = d * idx_jk + j
            kij = d * (idx_ij - 1) + k
            jd = d * (j - 1)

            @inbounds for i in 1:j
                denom_ij = viw_λ_Δ[idx_ij]
                if iszero(denom_ij)
                    denom_ik = viw_λ_Δ[idx_jk + i]
                    if iszero(denom_ik)
                        t = (∇3h[i] + ∇3h_j + ∇3h_k) * i6
                    else
                        t = (Δh[jd + i] - Δh_jk) / denom_ik
                    end
                else
                    t = (Δh[kd + i] - Δh_jk) / denom_ij
                end

                Δ2h[ijk + i] = Δ2h[jik] = Δ2h[kij] = t

                idx_ij += 1
                jik += d
                kij += d
            end
        end
        idx_jk += k
    end

    cone.dder3_aux_updated = true
    return
end

function dder3(
    cone::EpiPerSepSpectral{<:MatrixCSqr{T}},
    dir::AbstractVector{T},
    ) where T
    cone.dder3_aux_updated || update_dder3_aux(cone)
    d = cone.d
    v = cone.point[2]
    dder3 = cone.dder3
    cache = cone.cache
    ζi = cache.ζi
    viw_X = cache.viw_X
    viw_λ = cache.viw_λ
    w_λi = cache.w_λi
    σ = cache.σ
    ∇h = cache.∇h
    ∇2h = cache.∇2h
    ∇3h = cache.∇3h
    Δh = cache.Δh
    Δ2h = cache.Δ2h
    vi = inv(v)

    r_X = cache.w1
    ξ_X = cache.w2
    ξb = cache.w3
    wd = cache.wd

    p = dir[1]
    q = dir[2]
    @views svec_to_smat!(r_X, dir[3:end], cache.rt2)
    mul!(ξ_X, Hermitian(r_X, :U), viw_X)
    mul!(r_X, viw_X', ξ_X)
    LinearAlgebra.copytri!(r_X, 'U', true)

    viq = vi * q
    D = Diagonal(viw_λ)
    @. ξ_X = vi * r_X - viq * D
    @. ξb = ζi * Δh * ξ_X
    @inbounds sum1 = sum(∇h[i] * real(r_X[i, i]) for i in 1:d)
    ζiχ = ζi * (p - σ * q - sum1)
    ξbξ = dot(Hermitian(ξb, :U), Hermitian(ξ_X, :U)) / 2
    c1 = -ζi * (ζiχ^2 + v * ξbξ)

    w_aux = ξb
    w_aux .*= ζiχ + viq
    ξ_X_views = [view(ξ_X, :, j) for j in 1:d]
    col = 1
    @inbounds for j in 1:d
        ξ_Xj = ξ_X_views[j]
        for i in 1:j
            ξ_Xi = ξ_X_views[i]
            @views Δ2h_ij = Δ2h[:, col]
            w_aux[i, j] -= ζi * dot(ξ_Xi, Diagonal(Δ2h_ij), ξ_Xj)
            col += 1
        end
    end
    @inbounds c2 = sum(viw_λ[i] * real(w_aux[i, i]) for i in 1:d)

    @. wd = sqrt(w_λi)
    lmul!(Diagonal(w_λi), r_X)
    rmul!(r_X, Diagonal(wd))
    mul!(w_aux, r_X, r_X', true, true)
    D_∇h = Diagonal(∇h)
    @. w_aux += c1 * D_∇h
    mul!(ξ_X, Hermitian(w_aux, :U), viw_X')
    mul!(w_aux, viw_X, ξ_X)

    dder3[1] = -c1
    @inbounds dder3[2] = c1 * σ - c2 + ξbξ + viq^2 / v
    @views smat_to_svec!(dder3[3:end], w_aux, cache.rt2)

    return dder3
end


text:
# Start version of HelloWorld for Learning Julia course

# Define a function that prints out a message
# This creates a function with a single "method"
function helloWorld()
    println("Hello world!")
end

# By defining another version of the function with different
# arguments, we create another "method" - Julia will choose
# the right one based on the arguments. This is called
# "multiple dispatch"
function helloWorld(count::Int)
    for i in 1:count
        println("Hello world ", i, "!")
    end
end

# Call the version with no arguments
helloWorld()

# Call the version that takes a number
helloWorld(3)

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
using RunMetrics

cd(@__DIR__)

path = joinpath("..", "sample_data", "sample.tcx")

rs = parse_tcx(path)
urs = RunMetrics.unit_run_sum(rs)
# X = RunMetrics.make_design(urs, 100.0)

lags = 1.0:1.0:600.0
n = length(lags)
B = zeros(3, n)
ssrs = zeros(n)
for i = 1:n
    X = RunMetrics.make_design(urs, lags[i])

    β, ssr = RunMetrics.regress(urs.hr, X)
    B[:, i] .= β
    ssrs[i] = ssr
end
min_ind = findmin(ssrs)[2]
min_lag = lags[min_ind]
β = B[:, min_ind]


text:
module MultiQuad

using QuadGK, Cuba, Unitful

export quad, dblquad, tplquad

@doc raw"""
    quad(arg::Function, x1, x2; method = :quadgk, kwargs...)

Performs the integral ``\int_{x1}^{x2}f(x)dx``
See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.

# Examples
```jldoctest
func(x) = x^2*exp(-x)
integral, error = quad(func, 0, 4)
```

`quad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:

```jldoctest
using Unitful

func(x) = x^2
integral, error = quad(func, 1u"m", 5u"m")
```
"""
function quad(arg::Function, x1, x2; method = :quadgk, kwargs...)

    if method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    elseif method == :quadgk
        integrate = quadgk
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    if method == :quadgk
        return quadgk(arg, x1, x2; kwargs...)
    end

    units = unit(arg(x1)) * unit(x1)

    arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64

    function integrand(x, f)
        f[1] = arg2(x[1])
    end

    result, err = integrate(integrand, 1, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

@doc raw"""
    dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)

Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx``
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.

# Examples
```jldoctest
func(y,x) = sin(x)*y^2
integral, error = dblquad(func, 1, 2, x->0, x->x^2, rtol=1e-9)
```

`dblquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:

```jldoctest
using Unitful

func(y,x) = x*y^2
integral, error = dblquad(func, 1u"m", 2u"m", x->0u"m^2", x->x^2, rtol=1e-9)
```
"""
function dblquad(
    arg::Function,
    x1,
    x2,
    y1::Function,
    y2::Function;
    method = :cuhre,
    kwargs...,
)

    if method == :cuhre
        integrate = cuhre
    elseif method == :divonne
        integrate = divonne
    elseif method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))

    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)

    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64

    function integrand(x, f)
        f[1] = arg2(x[1], x[2])
    end

    result, err = integrate(integrand, 2, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

@doc raw"""
    tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)

Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx``
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.

# Examples
```jldoctest
func(z,y,x) = sin(z)*y*x
integral, error = tplquad(func, 0, 4, x->x, x->x^2, (x,y)->2, (x,y)->3*x)
```

`tplquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:

```jldoctest
using Unitful

func(z,y,x) = sin(z)*y*x
integral, error = tplquad(func, 0u"m", 4u"m", x->0u"m^2", x->x^2, (x,y)->0, (x,y)->3)
```
"""
function tplquad(
    arg::Function,
    x1,
    x2,
    y1::Function,
    y2::Function,
    z1::Function,
    z2::Function;
    method = :cuhre,
    kwargs...,
)

    if method == :cuhre
        integrate = cuhre
    elseif method == :divonne
        integrate = divonne
    elseif method == :suave
        integrate = suave
    elseif method == :vegas
        integrate = vegas
    else
        ex = ErrorException("Integration method $method is not supported!")
        throw(ex)
    end

    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
            unit(z1(y1(x1), x1))

    arg0(a, y, x) =
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)

    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)

    arg2(a, b, c) =
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64

    function integrand(x, f)
        f[1] = arg2(x[1], x[2], x[3])
    end

    result, err = integrate(integrand, 3, 1; kwargs...)

    if units == Unitful.NoUnits
        return result[1], err[1]
    else
        return (result[1], err[1]) .* units
    end
end

end # module


text:
using ExtendedExtremes
using Extremes, Random, Test

# Set the seed for reproductible test results
Random.seed!(12)

@testset "ExtendedExtremes.jl" begin
    include("reproducingNaveauResults.jl")
    include("parameterestimation.jl")
    include("EGPpower.jl")
    include("EGPbeta.jl")
    include("EGPbetapower.jl")
end;


text:
using FactorGraph
using OMEinsum: bpcheck, einsum
using Zygote
using Test

Base.conj(arr::AbstractArray{<:Any, 0}) = conj!(copy(arr))

@testset "einsum bp" begin
    a = randn(ComplexF64, 3,3)
    f2(a) = ein"ij,ik,il->jkl"(a, conj(a), a) |> sum |> real

    @test gradient(f2, a)[1] !== nothing
    @test bpcheck(f2, a)
    @test bpcheck(a->treecontract(((1,2),3), ((1,2), (2,3), (3,1)), (a, a, a), ())[] |> real, a)
end

@testset "optcontract" begin
    a = randn(ComplexF64, 3,3)
    @test optcontract(((1,2), (2,3)), (a, a), (1,3)) ≈ a*a
    @test bpcheck(a->ein"ij,ji->"(a, a) |> sum |> real, a)
    @test bpcheck(a->optcontract(((1,2), (2,3)), (a, a), (1,3)) |> sum |> real, a)
    @test bpcheck(a->optcontract(((1,2), (2,3), (3,4)), (a,a,a), (4,1)) |> sum |> real, a)
end


text:
module ClpCInterface

using SparseArrays

export
    # Types
    CoinBigIndex,
    CoinBigDouble,
    ClpModel,
    ClpSolve,

    # Methods
    load_problem,
    load_quadratic_objective,
    read_mps,
    copy_in_integer_information,
    delete_integer_information,
    resize,
    delete_rows,
    add_rows,
    add_row,
    delete_columns,
    add_column,
    add_columns,
    chg_row_lower,
    chg_row_upper,
    chg_column_lower,
    chg_column_upper,
    chg_obj_coefficients,
    drop_names,
    copy_names,
    get_num_rows,
    get_num_cols,
    number_rows,
    number_cols,
    primal_tolerance,
    set_primal_tolerance,
    dual_tolerance,
    set_dual_tolerance,
    dual_objective_limit,
    set_dual_objective_limit,
    objective_offset,
    set_objective_offset,
    problem_name,
    set_problem_name,
    number_iterations,
    get_iteration_count,
    set_number_iterations,
    maximum_iterations,
    set_maximum_iterations,
    maximum_seconds,
    set_maximum_seconds,
    hit_maximum_iterations,
    status,
    set_problem_status,
    secondary_status,
    set_secondary_status,
    optimization_direction,
    set_optimization_direction,
    get_row_activity,
    get_col_solution,
    primal_column_solution,
    set_col_solution,
    primal_row_solution,
    get_row_price,
    dual_row_solution,
    get_reduced_cost,
    dual_column_solution,
    get_row_lower,
    row_lower,
    get_row_upper,
    row_upper,
    get_obj_coefficients,
    objective,
    get_col_lower,
    column_lower,
    get_col_upper,
    column_upper,
    get_num_elements,
    get_vector_starts,
    get_indices,
    get_constraint_matrix,
    get_vector_lengths,
    get_elements,
    get_obj_value,
    objective_value,
    integer_information,
    infeasibility_ray,
    unbounded_ray,
    statusExists,
    status_array,
    copyin_status,
    get_column_status,
    get_row_status,
    set_column_status,
    set_row_status,
    register_call_back,
    clear_call_back,
    log_level,
    set_log_level,
    length_names,
    row_name,
    column_name,
    initial_solve,
    initial_solve_with_options,
    initial_dual_solve,
    initial_primal_solve,
    initial_barrier_solve,
    initial_barrier_no_cross_solve,
    dual,
    primal,
    idiot,
    scaling,
    scaling_flag,
    crash,
    primal_feasible,
    dual_feasible,
    dual_bound,
    set_dual_bound,
    infeasibility_cost,
    set_infeasibility_cost,
    perturbation,
    set_perturbation,
    algorithm,
    set_algorithm,
    sum_dual_infeasibilities,
    number_dual_infeasibilities,
    sum_primal_infeasibilities,
    number_primal_infeasibilities,
    save_model,
    restore_model,
    check_solution,
    is_abandoned,
    is_proven_optimal,
    is_proven_primal_infeasible,
    is_proven_dual_infeasible,
    is_primal_objective_limit_reached,
    is_dual_objective_limit_reached,
    is_iteration_limit_reached,
    get_obj_sense,
    set_obj_sense,
    print_model,
    get_small_element_value,
    set_small_element_value,

    # ClpSolve methods
    set_special_option,
    get_special_option,
    set_presolve_type,
    get_presolve_type,
    set_solve_type,
    get_solve_type,
    get_presolve_passes,
    get_extra_info,
    set_infeasible_return,
    infeasible_return,
    do_dual,
    set_do_dual,
    do_singleton,
    set_do_singleton,
    do_doubleton,
    set_do_doubleton,
    do_tripleton,
    set_do_tripleton,
    do_tighten,
    set_do_tighten,
    do_forcing,
    set_do_forcing,
    do_implied_free,
    set_do_implied_free,
    do_dupcol,
    set_do_dupcol,
    do_duprow,
    set_do_duprow,
    do_singleton_column,
    set_do_singleton_column,
    presolve_actions,
    set_presolve_actions,
    substitution,
    set_substitution

import Base.pointer

if haskey(ENV,"JULIA_CLP_LIBRARY_PATH") || VERSION < v"1.3"
    if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
        include("../deps/deps.jl")
    else
        error("Clp not properly installed. Please run import Pkg; Pkg.build(\"Clp\")")
    end
else
    import Clp_jll: libClp
end
## Shared library interface setup
#{{{

macro clp_ccall(func, args...)
    args = [esc(ex) for ex in args]
    f = "Clp_$(func)"
    quote
        ccall(($f,libClp), $(args...))
    end
end

macro clpsolve_ccall(func, args...)
    args = [esc(ex) for ex in args]
    f = "ClpSolve_$(func)"
    quote
        ccall(($f,libClp), $(args...))
    end
end

# Note: we assume COIN_BIG_INDEX and COIN_BIG_DOUBLE
# were not defined when compiling Clp (which is the
# default)
const CoinBigIndex  = Int32
const CoinBigDouble = Float64

#}}}

## Main types definitions
#{{{
mutable struct ClpModel
    p::Ptr{Cvoid}
    function ClpModel()
        p = @clp_ccall newModel Ptr{Cvoid} ()
        prob = new(p)
        finalizer(delete_model, prob)
        return prob
    end
end

function delete_model(model::ClpModel)
    if model.p == C_NULL
        return
    end
    @clp_ccall deleteModel Cvoid (Ptr{Cvoid},) model.p
    model.p = C_NULL
    return
end

mutable struct ClpSolve
    p::Ptr{Cvoid}
    function ClpSolve()
        p = @clpsolve_ccall new Ptr{Cvoid} ()
        prob = new(p)
        finalizer(delete_solve, prob)
        return prob
    end
end

function delete_solve(solve::ClpSolve)
    if solve.p == C_NULL
        return
    end
    @clpsolve_ccall delete Cvoid (Ptr{Cvoid},) solve.p
    solve.p = C_NULL
    return
end

pointer(model::ClpModel) = model.p
#}}}

## Check functions for internal use
#{{{
# Functions which perform all sorts of
# sanity checks on input parameters and
# throw exceptions in case of errors.
# Ideally, it should never be possible
# to pass an invalid parameter to the
# underlying Clp API.

function _jl__check_model(model::ClpModel)
    if model.p == C_NULL
        error("Invalid ClpModel")
    end
    return true
end

function _jl__check_solve(solve::ClpSolve)
    if solve.p == C_NULL
        error("Invalid ClpSolve")
    end
    return true
end

function _jl__check_row_is_valid(model::ClpModel, row::Integer)
    num_rows = @clp_ccall getNumRows Int32 (Ptr{Cvoid},) model.p
    if !(1 <= row <= num_rows)
        error("Invalid row $row (must be 1 <= row <= $num_rows)")
    end
    return true
end

function _jl__check_col_is_valid(model::ClpModel, col::Integer)
    num_cols = @clp_ccall getNumCols Int32 (Ptr{Cvoid},) model.p
    if !(1 <= col <= num_cols)
        error("Invalid col $col (must be 1 <= col <= $num_cols)")
    end
    return true
end

function _jl__check_file_is_readable(filename::String)
    try
        f = open(filename, "r")
        close(f)
    catch err
        error("file $filename not readable")
    end
    return true
end
#}}}

## CLP functions
#{{{

# inspired by GLPK interface
const VecOrNothing = Union{Vector,Cvoid}
function vec_or_null(::Type{T}, a::VecOrNothing, len::Integer) where T
    if isequal(a, nothing) || isa(a, Array{Cvoid}) # on 0.3, [] is Array{Cvoid}
        return C_NULL
    else # todo: helpful message if convert fails
        if length(a) != len
            error("Expected vector to have length $len")
        end
        return convert(Vector{T},a)
    end
end

# Load model - loads some stuff and initializes others
# Loads a problem (the constraints on the
# rows are given by lower and upper bounds). If a pointer is NULL then the
# following values are the default:
#
# col_ub: all columns have upper bound infinity
# col_lb: all columns have lower bound 0
# row_ub: all rows have upper bound infinity
# row_lb: all rows have lower bound -infinity
# obj: all variables have 0 objective coefficient

# Just like the other load_problem() method except that the matrix is
# given in a standard column major ordered format (without gaps).
function load_problem(model::ClpModel,  num_cols::Integer, num_rows::Integer,
        start::Vector{CoinBigIndex}, index::Vector{Int32},
        value::Vector{Float64},
        col_lb::VecOrNothing, col_ub::VecOrNothing,
        obj::VecOrNothing,
        row_lb::VecOrNothing, row_ub::VecOrNothing)
    _jl__check_model(model)
    @clp_ccall loadProblem Cvoid (Ptr{Cvoid},Int32,Int32,Ptr{CoinBigIndex},Ptr{Int32},
    Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64}) model.p num_cols num_rows start index value vec_or_null(Float64,col_lb,num_cols) vec_or_null(Float64,col_ub,num_cols) vec_or_null(Float64,obj,num_cols) vec_or_null(Float64,row_lb,num_rows) vec_or_null(Float64,row_ub,num_rows)
end

function load_problem(model::ClpModel,  constraint_matrix::AbstractMatrix,
    col_lb::VecOrNothing, col_ub::VecOrNothing,
    obj::VecOrNothing, row_lb::VecOrNothing,
    row_ub::VecOrNothing)
    mat = convert(SparseMatrixCSC{Float64,Int32},constraint_matrix)
    # We need to convert to zero-based, but
    # TODO: don't make extra copies of arrays
    # TODO: check dimensions match
    load_problem(model, mat.n, mat.m, mat.colptr .- convert(Int32,1),
                 mat.rowval .- convert(Int32, 1), mat.nzval,
                 col_lb, col_ub,obj, row_lb, row_ub)
end

# Read quadratic part of the objective (the matrix part).
function load_quadratic_objective(model::ClpModel,
        num_cols::Integer, start::Vector{CoinBigIndex},
        col::Vector{Int32}, element::Vector{Float64})
    _jl__check_model(model)
    @clp_ccall loadQuadraticObjective Cvoid (Ptr{Cvoid},Int32,Ptr{CoinBigIndex},Ptr{Int32},Ptr{Float64}) model.p num_cols start column element
end

function load_quadratic_objective(model::ClpModel,
    hessian_matrix::SparseMatrixCSC{Float64,Int32})
    load_quadratic_objective(model, hessian_matrix.n, hessian_matrix.colptr-convert(Int32,1), hessian_matrix.rowval-convert(Int32,1),hessian_matrix.nzval)
end

# Read an mps file from the given filename.
function read_mps(model::ClpModel, mpsfile::AbstractString, keep_names::Bool, ignore_errors::Bool)
    _jl__check_model(model)
    _jl__check_file_is_readable(mpsfile)

    status = @clp_ccall readMps Int32 (Ptr{Cvoid}, Cstring, Int32, Int32) model.p mpsfile keep_names ignore_errors
    if status != 0
        error("read_mps: error reading file $mpsfile")
    end
    return true
end
read_mps(model::ClpModel, mpsfile::AbstractString) = read_mps(model, mpsfile, true, false)

# Copy in integer information.
function copy_in_integer_information(model::ClpModel, information::Vector{UInt8})
    _jl__check_model(model)
    if length(information) != num_cols(model)
        error("Length of integer information must match number of columns")
    end
    @clp_ccall copyInIntegerInformation Cvoid (Ptr{Cvoid},Ptr{UInt8}) model.p information
end

# Drop integer information.
function delete_integer_information(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall deleteIntegerInformation Cvoid (Ptr{Cvoid},) model.p
end

# Resize rim part of model.
function resize(model::ClpModel, new_num_rows::Integer, new_num_cols::Integer)
    _jl__check_model(model)
    @clp_ccall resize Cvoid (Ptr{Cvoid},Int32,Int32) model.p new_num_rows new_num_cols
end

# Delete rows.
function delete_rows(model::ClpModel, which::Vector{Int32})
    _jl__check_model(model)
    @clp_ccall deleteRows Cvoid (Ptr{Cvoid},Int32,Ptr{Int32}) model.p length(which) which
end

# Add rows.
function add_rows(model::ClpModel, number::Integer, row_lower::Vector{Float64},
        row_upper::Vector{Float64},
        row_starts::Vector{Int32}, columns::Vector{Int32},
        elements::Vector{Float64})
    _jl__check_model(model)
    @clp_ccall addRows Cvoid (Ptr{Cvoid}, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}) model.p number row_lower row_upper row_starts columns elements
end

#This function exists in cpp but not c interface
function add_row(model::ClpModel, nnz::Cint, columns::Vector{Int32},
                 elements::Vector{Float64}, row_lower::Float64, row_upper::Float64)
    add_rows(model, 1, [row_lower], [row_upper], [Cint(0), nnz], columns,
             elements)
end

# Delete columns.
function delete_columns(model::ClpModel, which::Vector{Int32})
    _jl__check_model(model)
    @clp_ccall deleteColumns Cvoid (Ptr{Cvoid},Int32,Ptr{Int32}) model.p length(which) which
end

#This function exists in cpp but not c interface
function add_column(model::ClpModel, nnz::Int32,  rows::Vector{Int32},
                elements::Vector{Float64}, column_lower::Float64,
                column_upper::Float64, objective::Float64)
    add_columns(model, 1, [column_lower], [column_upper], [objective],
                [Int32(0), nnz], rows, elements)
end

# Add columns.
function add_columns(model::ClpModel, number::Integer, column_lower::Vector{Float64},
        column_upper::Vector{Float64},
        objective::Vector{Float64},
        column_starts::Vector{Int32}, rows::Vector{Int32},
        elements::Vector{Float64})
    _jl__check_model(model)

    @clp_ccall addColumns Cvoid (Ptr{Cvoid}, Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}) model.p number column_lower column_upper objective column_starts rows elements

end

function add_columns(model::ClpModel, column_lower::Vector{Float64},
        column_upper::Vector{Float64},
        objective::Vector{Float64},
        new_columns::SparseMatrixCSC{Float64,Int32})
    add_columns(model, new_columns.n, column_lower, column_upper, objective, new_columns.colptr-convert(Int32,1), new_columns.rowval-convert(Int32,1), new_columns.nzval)
end

# Change row lower bounds.
function chg_row_lower(model::ClpModel, row_lower::Vector{Float64})
    _jl__check_model(model)
    if length(row_lower) != get_num_rows(model)
        error("Array length must match number of rows in the model")
    end

    @clp_ccall chgRowLower Cvoid (Ptr{Cvoid}, Ptr{Float64}) model.p row_lower
end

# Change row upper bounds.
function chg_row_upper(model::ClpModel, row_upper::Vector{Float64})
    _jl__check_model(model)
    if length(row_upper) != get_num_rows(model)
        error("Array length must match number of rows in the model")
    end

    @clp_ccall chgRowUpper Cvoid (Ptr{Cvoid}, Ptr{Float64}) model.p row_upper
end

# Change column lower bounds.
function chg_column_lower(model::ClpModel, column_lower::Vector{Float64})
    _jl__check_model(model)
    if length(column_lower) != get_num_cols(model)
        error("Array length must match number of columns in the model")
    end

    @clp_ccall chgColumnLower Cvoid (Ptr{Cvoid}, Ptr{Float64}) model.p column_lower
end

# Change column upper bounds.
function chg_column_upper(model::ClpModel, column_upper::Vector{Float64})
    _jl__check_model(model)
    if length(column_upper) != get_num_cols(model)
        error("Array length must match number of columns in the model")
    end

    @clp_ccall chgColumnUpper Cvoid (Ptr{Cvoid}, Ptr{Float64}) model.p column_upper
end

# Change objective coefficients.
function chg_obj_coefficients(model::ClpModel, obj_in::Vector{Float64})
    _jl__check_model(model)
    if length(obj_in) != get_num_cols(model)
        error("Array length must match number of columns in the model")
    end

    @clp_ccall chgObjCoefficients Cvoid (Ptr{Cvoid},Ptr{Float64}) model.p obj_in
end

# Drops names - makes lengthnames 0 and names empty.
function drop_names(model::ClpModel)
    _jl__check_model(model)

    @clp_ccall dropNames Cvoid (Ptr{Cvoid},) model.p
end

# Copy in names.
# TODO: Need a nice way to accept a vector of Strings
function copy_names(model::ClpModel, row_names::Vector{Vector{UInt8}},
        columnNames::Vector{Vector{UInt8}})
    # TODO
    error("TODO")
    return
end

# Number of rows.
function get_num_rows(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall getNumRows Int32 (Ptr{Cvoid},) model.p
end
const number_rows = get_num_rows

# Number of columns.
function get_num_cols(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall getNumCols Int32 (Ptr{Cvoid},) model.p
end
const number_cols = get_num_cols

# Get primal tolerance.
function primal_tolerance(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall primalTolerance Float64 (Ptr{Cvoid},) model.p
end

# Set primal tolerance to use.
function set_primal_tolerance(model::ClpModel, value::Float64)
    _jl__check_model(model)
    @clp_ccall setPrimalTolerance Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

# Get dual tolerance.
function dual_tolerance(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall dualTolerance Float64 (Ptr{Cvoid},) model.p
end

# Set dual tolerance to use.
function set_dual_tolerance(model::ClpModel, value::Float64)
    _jl__check_model(model)
    @clp_ccall setDualTolerance Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

# Get dual objective limit.
function dual_objective_limit(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall dualObjectiveLimit Float64 (Ptr{Cvoid},) model.p
end

# Set dual objective limit.
function set_dual_objective_limit(model::ClpModel, value::Float64)
    @clp_ccall setDualObjectiveLimit Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

# Get objective offset.
function objective_offset(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall dualObjectiveLimit Float64 (Ptr{Cvoid},) model.p
end

# Set objective offset.
function set_objective_offset(model::ClpModel, value::Float64)
    @clp_ccall setObjectiveOffset Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

#=
These methods are untested. If you need them, open a PR to get them working
and add new unit tests.
# Fill in array with problem name.
function problem_name(model::ClpModel)
    _jl__check_model(model)
    problem_name = Array{UInt8}(undef, 1000)
    @clp_ccall problemName Cvoid (Ptr{Cvoid}, Int32, Ptr{UInt8}) model.p 1000 problem_name
    return bytestring(problem_name)
end

# Set problem name.  Must have \0 at end.
function set_problem_name(model::ClpModel, name::String)
    _jl__check_model(model)
    @assert isascii(name)

    @clp_ccall setProblemName Cvoid (Ptr{Cvoid}, Int32, Ptr{UInt8}) model.p (length(name)+1) bytestring(name)
end
=#

# Get number of iterations
function number_iterations(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall numberIterations Int32 (Ptr{Cvoid},) model.p
end
const get_iteration_count = number_iterations

# Set number of iterations
function set_number_iterations(model::ClpModel, iters::Integer)
    _jl__check_model(model)
    @clp_ccall setNumberIterations Cvoid (Ptr{Cvoid}, Int32) model.p iters
    return
end

# Get maximum number of iterations
function maximum_iterations(model::ClpModel)
    _jl__check_model(model)
    # The C function's name does not follow the `Clp_` convention
    # See https://github.com/JuliaOpt/Clp.jl/issues/67
    ccall((:maximumIterations, libClp), Int32, (Ptr{Cvoid},), model.p)
end

# Set maximum number of iterations
function set_maximum_iterations(model::ClpModel, max_iters::Integer)
    _jl__check_model(model)
    @clp_ccall setMaximumIterations Cvoid (Ptr{Cvoid}, Int32) model.p max_iters
    return
end

# Get maximum time in seconds (from when set is called)
function maximum_seconds(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall maximumSeconds Float64 (Ptr{Cvoid},) model.p
end

# Set maximum time in seconds (from when set is called)
function set_maximum_seconds(model::ClpModel, max_secs::Real)
    _jl__check_model(model)
    @clp_ccall setMaximumSeconds Cvoid (Ptr{Cvoid}, Float64) model.p max_secs
    return
end

# Query whether maximum iterations or time were hit
function hit_maximum_iterations(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall hitMaximumIterations Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Get the status of the problem:
#   0 - optimal
#   1 - primal infeasible
#   2 - dual infeasible
#   3 - stopped on iterations etc
#   4 - stopped due to errors
function status(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall status Int32 (Ptr{Cvoid},) model.p
end

# Set the status of the problem.
function set_problem_status(model::ClpModel, status::Integer)
    _jl__check_model(model)
    @clp_ccall setProblemStatus Cvoid (Ptr{Cvoid}, Int32) model.p status
    return
end

# Get the secondary status of problem - may get extended:
#   0 - none
#   1 - primal infeasible because dual limit reached
#   2 - scaled problem optimal - unscaled has primal infeasibilities
#   3 - scaled problem optimal - unscaled has dual infeasibilities
#   4 - scaled problem optimal - unscaled has both dual and primal infeasibilities
function secondary_status(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall secondaryStatus Int32 (Ptr{Cvoid},) model.p
end

# Set the secondary status of the problem.
function set_secondary_status(model::ClpModel, status::Integer)
    _jl__check_model(model)
    @clp_ccall setSecondaryStatus Cvoid (Ptr{Cvoid}, Int32) model.p status
    return
end

# Get the direction of optimization:
#   1 - minimize
#   -1 - maximize
#   0 - ignore
# XXX: for some reason this returns a floating point (???)
function optimization_direction(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall optimizationDirection Float64 (Ptr{Cvoid},) model.p
end

# Set the direction of optimization.
# XXX: for some reason this takes a floating point argument (???)
function set_optimization_direction(model::ClpModel, value::Real)
    _jl__check_model(model)
    @clp_ccall setOptimizationDirection Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

# TODO: do we actually want to make a copy of the result?
# More efficient to not and let the sufficiently warned user do so if desired.
macro def_get_row_property(fname,clpname)
    return quote
        function $(esc(fname))(model::ClpModel)
            _jl__check_model(model)
            row_property_p = @clp_ccall $clpname Ptr{Float64} (Ptr{Cvoid},) model.p
            num_rows = convert(Int,get_num_rows(model))
            return copy(unsafe_wrap(Array,row_property_p,(num_rows,)))
        end
    end
end

macro def_get_col_property(fname,clpname)
    return quote
        function $(esc(fname))(model::ClpModel)
            _jl__check_model(model)
            col_property_p = @clp_ccall $clpname Ptr{Float64} (Ptr{Cvoid},) model.p
            num_cols = convert(Int,get_num_cols(model))
            return copy(unsafe_wrap(Array,col_property_p,(num_cols,)))
        end
    end
end

# Get primal row solution.
@def_get_row_property get_row_activity getRowActivity

const primal_row_solution = get_row_activity

# Get primal column solution.
@def_get_col_property get_col_solution getColSolution
const primal_column_solution = get_col_solution

function set_col_solution(model::ClpModel, input::Vector{Float64})
    _jl__check_model(model)

    @clp_ccall setColSolution Cvoid (Ptr{Cvoid},Ptr{Float64}) model.p input
end

# Get dual row solution.
@def_get_row_property get_row_price getRowPrice
const dual_row_solution = get_row_price

# Get dual column solution (i.e. the reduced costs).
@def_get_col_property get_reduced_cost getReducedCost
const dual_column_solution = get_reduced_cost

# Get row lower bounds.
@def_get_row_property get_row_lower getRowLower
const row_lower = get_row_lower

# Get row upper bounds.
@def_get_row_property get_row_upper getRowUpper
const row_upper = get_row_upper

@def_get_col_property get_obj_coefficients getObjCoefficients
const objective = get_obj_coefficients

# Get column lower bounds.
@def_get_col_property get_col_lower getColLower
const column_lower = get_col_lower

# Get column upper bounds.
@def_get_col_property get_col_upper getColUpper
const column_upper = get_col_upper

# Get the number of elements in matrix.
function get_num_elements(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall getNumElements Int32 (Ptr{Cvoid},) model.p
end

# Get column starts in matrix.
function get_vector_starts(model::ClpModel)
    _jl__check_model(model)
    vec_starts_p = @clp_ccall getVectorStarts Ptr{CoinBigIndex} (Ptr{Cvoid},) model.p
    num_cols = Int(get_num_cols(model))
    return copy(unsafe_wrap(Array,vec_starts_p, (num_cols+1,)))
end

# Get row indices in matrix.
function get_indices(model::ClpModel)
    _jl__check_model(model)
    # getIndices returns an "int*", how do we know it's Int32??
    row_indices_p = @clp_ccall getIndices Ptr{Int32} (Ptr{Cvoid},) model.p
    num_elts = Int(get_num_elements(model))
    return copy(unsafe_wrap(Array,row_indices_p,(num_elts,)))
end

function get_constraint_matrix(model::ClpModel)
    _jl__check_model(model)
    @assert CoinBigIndex == Int32
    # SparseMatrixCSC requires same integer type for colptr and rowval
    num_cols = convert(Int, get_num_cols(model))
    num_rows = convert(Int, get_num_rows(model))
    colptr = get_vector_starts(model) .+ convert(Int32,1)
    rowval = get_indices(model) .+ convert(CoinBigIndex,1)
    nzval = get_elements(model)

    return SparseMatrixCSC{Float64,CoinBigIndex}(num_rows,num_cols,colptr,rowval,nzval)

end

# Get column vector lengths in matrix.
@def_get_col_property get_vector_lengths getVectorLengths

# Get element values in matrix.
function get_elements(model::ClpModel)
    _jl__check_model(model)
    elements_p = @clp_ccall getElements Ptr{Float64} (Ptr{Cvoid},) model.p
    num_elts = Int(get_num_elements(model))
    return copy(unsafe_wrap(Array,elements_p,(num_elts,)))
end

# Get objective value.
function get_obj_value(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall getObjValue Float64 (Ptr{Cvoid},) model.p
end
const objective_value = get_obj_value

# Get integer information.
function integer_information(model::ClpModel)
    # TODO
    error("TODO")
    return # vector of chars (?)
end

# Get an infeasibility ray (empty vector returned if none/wrong).
function infeasibility_ray(model::ClpModel)
    _jl__check_model(model)
    infeas_ray_p = @clp_ccall(
        infeasibilityRay, Ptr{Float64}, (Ptr{Cvoid},), model.p
    )
    if infeas_ray_p == C_NULL
        return Float64[]
    end
    num_rows = convert(Int, get_num_rows(model))
    return unsafe_wrap(Array, infeas_ray_p, (num_rows,); own = true)
end

# Get an unbounded ray (empty vector returned if none/wrong).
function unbounded_ray(model::ClpModel)
    _jl__check_model(model)
    unbd_ray_p = @clp_ccall(unboundedRay, Ptr{Float64}, (Ptr{Cvoid},), model.p)
    if unbd_ray_p == C_NULL
        return Float64[]
    end
    num_cols = convert(Int, get_num_cols(model))
    return unsafe_wrap(Array, unbd_ray_p, (num_cols,); own = true)
end

# Query whether the status array exists (partly for OsiClp).
function statusExists(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall statusExists Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Return the address of the status array (char[num_rows+num_cols]).
# Status values are:
#   0 - free
#   1 - basic
#   2 - at upper
#   3 - at lower
#   4 - superbasic
#   (5 - fixed)
#   [???]
function status_array(model::ClpModel)
    # TODO
    error("TODO")
    return # pointer to unsigned chars (?)
end

# Copy in the status vector
# XXX copyin rather than copy_in, WTF?
# XXX input vector of unsigned chars (?)
function copyin_status(model::ClpModel, status_array::Vector{UInt8})
    # TODO
    error("TODO")
    return
end

# Get variable basis info.
# Note that we adjust from one-based indices
function get_column_status(model::ClpModel, sequence::Integer)
    _jl__check_model(model)
    _jl__check_col_is_valid(model, sequence)
    return @clp_ccall getColumnStatus Int32 (Ptr{Cvoid},Int32) model.p (sequence-1)
end

# Get row basis info.
function get_row_status(model::ClpModel, sequence::Integer)
    _jl__check_model(model)
    _jl__check_row_is_valid(model, sequence)
    return @clp_ccall getRowStatus Int32 (Ptr{Cvoid},Int32) model.p (sequence-1)
end

# Set variable basis info (and value if at bound).
function set_column_status(model::ClpModel, sequence::Integer, value::Integer)
    _jl__check_model(model)
    _jl__check_col_is_valid(model, sequence)
    @clp_ccall setColumnStatus Cvoid (Ptr{Cvoid},Int32,Int32) model.p (sequence-1) value
end

# Set row basis info (and value if at bound).
function set_row_status(model::ClpModel, sequence::Integer, value::Integer)
    _jl__check_model(model)
    _jl__check_row_is_valid(model, sequence)
    @clp_ccall setRowStatus Cvoid (Ptr{Cvoid},Int32,Int32) model.p (sequence-1) value
end

# No reason to expose these
# Set user pointer (used for generic purposes)
#function get_user_pointer(model::ClpModel)
# Get user pointer (used for generic purposes)
#function set_user_pointer(model::ClpModel, pointer::Ptr{Cvoid})

# Pass in c-pointer to callback function. (from cfunction())
# See Coin_C_defines.h for signature
# TODO: test this
function register_call_back(model::ClpModel, callback::Ptr{Cvoid})
    _jl__check_model(model)
    @clp_ccall registerCallBack Cvoid (Ptr{Cvoid},Ptr{Cvoid}) model.p callback
end

# Unset callback function.
function clear_call_back(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall clearCallBack Cvoid (Ptr{Cvoid},) model.p
    return
end

# Get amount of print out:
#   0 - none
#   1 - just final
#   2 - just factorizations
#   3 - as 2 plus a bit more
#   4 - verbose
#   above that: 8,16,32 etc just for selective debug.
function log_level(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall logLevel Int32 (Ptr{Cvoid},) model.p
end

# Set amount of print out.
function set_log_level(model::ClpModel, value::Integer)
    _jl__check_model(model)
    @clp_ccall setLogLevel Cvoid (Ptr{Cvoid}, Int32) model.p value
    return
end

# Get length of names (0 means no names).
function length_names(model::ClpModel)
    _jl__check_model(model)
    return @clp_ccall lengthNames Int32 (Ptr{Cvoid},) model.p
end

# Return an array with a row name.
# XXX deviation from Clp API syntax, which fills a user-provided
# array of chars:
#   void Clp_rowName(Clp_Simplex * model, int iRow, char * name);
function row_name(model::ClpModel, row::Integer)
    _jl__check_model(model)
    _jl__check_row_is_valid(model, row)
    size = @clp_ccall lengthNames Int32 (Ptr{Cvoid},) model.p
    row_name = Array{UInt8}(undef, size+1)
    @clp_ccall rowName Cvoid (Ptr{Cvoid}, Int32, Ptr{UInt8}) model.p (row-1) row_name
    return unsafe_string(row_name)
end

# Return an array with a column name.
# XXX deviation from Clp API syntax, which fills a user-provided
# array of chars:
#   void Clp_columnName(Clp_Simplex * model, int iColumn, char * name);
function column_name(model::ClpModel, col::Integer)
    _jl__check_model(model)
    _jl__check_col_is_valid(model, col)
    size = @clp_ccall lengthNames Int32 (Ptr{Cvoid},) model.p
    col_name = Array{UInt8}(undef, size+1)
    @clp_ccall columnName Cvoid (Ptr{Cvoid}, Int32, Ptr{UInt8}) model.p (col-1) col_name
    return unsafe_string(col_name)
end

# General solve algorithm which can do presolve.
function initial_solve(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall initialSolve Int32 (Ptr{Cvoid},) model.p
end

function initial_solve_with_options(model::ClpModel,solve::ClpSolve)
    _jl__check_model(model)
    _jl__check_solve(solve)
    @clp_ccall initialSolveWithOptions Int32 (Ptr{Cvoid},Ptr{Cvoid}) model.p solve.p
end

initial_solve(model::ClpModel, solve::ClpSolve) = initial_solve_with_options(model,solve)

# Dual initial solve.
function initial_dual_solve(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall initialDualSolve Int32 (Ptr{Cvoid},) model.p
end

# Primal initial solve.
function initial_primal_solve(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall initialPrimalSolve Int32 (Ptr{Cvoid},) model.p
end

# Barrier initial solve.
function initial_barrier_solve(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall initialBarrierSolve Int32 (Ptr{Cvoid},) model.p
end

# Barrier initial solve, no crossover.
function initial_barrier_no_cross_solve(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall initialBarrierNoCrossSolve Int32 (Ptr{Cvoid},) model.p
end

# Dual algorithm.
function dual(model::ClpModel, ifValuesPass::Integer)
    _jl__check_model(model)
    @clp_ccall dual Int32 (Ptr{Cvoid},Int32) model.p ifValuesPass
end

# Primal algorithm.
function primal(model::ClpModel, ifValuesPass::Integer)
    _jl__check_model(model)
    @clp_ccall primal Int32 (Ptr{Cvoid},Int32) model.p ifValuesPass
end

# Solve the problem with the idiot code.
function idiot(model::ClpModel, tryhard::Integer)
    _jl__check_model(model)
    @clp_ccall idiot Int32 (Ptr{Cvoid},Int32) model.p tryhard
end

# Set or unset scaling:
#  0 - off
#  1 - equilibrium
#  2 - geometric
#  3 - auto
#  4 - dynamic (later)
function scaling(model::ClpModel, mode::Integer)
    _jl__check_model(model)
    if !(0 <= mode <= 4)
        error("Invalid clp scaling mode $mode (must be between 0 and 4)")
    end
    @clp_ccall scaling Cvoid (Ptr{Cvoid}, Int32) model.p mode
    return
end

# Get scaling flag.
function scaling_flag(model::ClpModel)
    _jl__check_model(model)
    return @clp_ccall scalingFlag Int32 (Ptr{Cvoid},) model.p
end

# Crash (at present just aimed at dual); returns:
#   -2 if dual preferred and crash basis created
#   -1 if dual preferred and all slack basis preferred
#   0 if basis going in was not all slack
#   1 if primal preferred and all slack basis preferred
#   2 if primal preferred and crash basis created.
#
#   If gap between bounds <="gap" variables can be flipped
#
#   If "pivot" is:
#     0 - No pivoting (so will just be choice of algorithm)
#     1 - Simple pivoting e.g. gub
#     2 - Mini iterations
function crash(model::ClpModel, gap::Float64, pivot::Int32)
    _jl__check_model(model)
    if !(0 <= pivot <= 2)
        error("Invalid clp crash pivot value $pivot (must be between 0 and 2)")
    end
    return @clp_ccall crash Int32 (Ptr{Cvoid}, Float64, Int32) model.p gap pivot
end

# Query whether the problem is primal feasible.
function primal_feasible(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall primalFeasible Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether the problem is dual feasible.
function dual_feasible(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall dualFeasible Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Get dual bound.
function dual_bound(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall dualBound Float64 (Ptr{Cvoid},) model.p
end

# Set dual bound.
function set_dual_bound(model::ClpModel, value::Real)
    _jl__check_model(model)
    @clp_ccall setDualBound Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

# Get infeasibility cost.
function infeasibility_cost(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall infeasibilityCost Float64 (Ptr{Cvoid},) model.p
end

# Set infeasibility cost.
function set_infeasibility_cost(model::ClpModel, value::Real)
    _jl__check_model(model)
    @clp_ccall setInfeasibilityCost Cvoid (Ptr{Cvoid}, Float64) model.p value
    return
end

# Get Perturbation. Values are:
#   50  - switch on perturbation
#   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
#   101 - we are perturbed
#   102 - don't try perturbing again
# The default is 100; others are for playing.
function perturbation(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall perturbation Int32 (Ptr{Cvoid},) model.p
end

# Set Perturbation.
function set_perturbation(model::ClpModel, value::Integer)
    _jl__check_model(model)
    if !(value == 50 || value == 100 || value == 101 || value == 102)
        error("Invalid clp perturbation value: $value (must be one of 50,100,101,102)")
    end
    @clp_ccall setPerturbation Cvoid (Ptr{Cvoid}, Int32) model.p value
    return
end

# Get current (or last) algorithm.
function algorithm(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall algorithm Int32 (Ptr{Cvoid},) model.p
end

# Set algorithm.
function set_algorithm(model::ClpModel, value::Integer)
    _jl__check_model(model)
    # XXX which values of the algorithm are valid ???
    @clp_ccall setAlgorithm Cvoid (Ptr{Cvoid}, Int32) model.p value
    return
end

# Get the sum of dual infeasibilities.
function sum_dual_infeasibilities(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall sumDualInfeasibilities Float64 (Ptr{Cvoid},) model.p
end

# Get the number of dual infeasibilities.
function number_dual_infeasibilities(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall numberDualInfeasibilities Int32 (Ptr{Cvoid},) model.p
end

# Get the sum of primal infeasibilities.
function sum_primal_infeasibilities(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall sumPrimalInfeasibilities Float64 (Ptr{Cvoid},) model.p
end

# Get the number of primal infeasibilities.
function number_primal_infeasibilities(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall numberPrimalInfeasibilities Int32 (Ptr{Cvoid},) model.p
end

# Save model to file, returns 0 if success.  This is designed for
# use outside algorithms so does not save iterating arrays etc.
# It does not save any messaging information.
# Does not save scaling values.
# It does not know about all types of virtual functions.
function save_model(model::ClpModel, file_name::String)
    _jl__check_model(model)
    @assert isascii(file_name)

    @clp_ccall saveModel Int32 (Ptr{Cvoid}, Cstring) model.p file_name
end

# Restore model from file, returns 0 if success,
# deletes current model.
function restore_model(model::ClpModel, file_name::String)
    _jl__check_model(model)
    @assert isascii(file_name)

    @clp_ccall restoreModel Int32 (Ptr{Cvoid}, Cstring) model.p file_name
end

# Just check solution (for external use) - sets sum of
# infeasibilities etc.
function check_solution(model::ClpModel);
    _jl__check_model(model)
    @clp_ccall checkSolution Cvoid (Ptr{Cvoid},) model.p
    return
end

# Query whether there are numerical difficulties.
function is_abandoned(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isAbandoned Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether optimality is proven.
function is_proven_optimal(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isProvenOptimal Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether primal infeasiblity is proven.
function is_proven_primal_infeasible(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isProvenPrimalInfeasible Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether dual infeasiblity is proven.
function is_proven_dual_infeasible(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isProvenDualInfeasible Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether the given primal objective limit is reached.
function is_primal_objective_limit_reached(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isPrimalObjectiveLimitReached Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether the given dual objective limit is reached.
function is_dual_objective_limit_reached(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isDualObjectiveLimitReached Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Query whether the iteration limit is reached
function is_iteration_limit_reached(model::ClpModel)
    _jl__check_model(model)
    ret = @clp_ccall isIterationLimitReached Int32 (Ptr{Cvoid},) model.p
    return ret != 0
end

# Get the direction of optimization:
#   1  - minimize
#   -1 - maximize
#   0  - ignore
function get_obj_sense(model::ClpModel)
    _jl__check_model(model)
    @clp_ccall getObjSense Float64 (Ptr{Cvoid},) model.p
end

# Set the direction of optimization.
function set_obj_sense(model::ClpModel, objsen::Real)
    _jl__check_model(model)
    if !(objsen == -1 || objsen == 0 || objsen == 1)
        error("Invalid clp objsense $objsense (should be -1,0 or 1)")
    end
    @clp_ccall setObjSense Cvoid (Ptr{Cvoid}, Float64) model.p objsen
    return
end

# Print model for debugging purposes
function print_model(model::ClpModel, prefix::String)
    _jl__check_model(model)
    @assert isascii(prefix)

    @clp_ccall printModel Cvoid (Ptr{Cvoid}, Cstring) model.p prefix
end

function get_small_element_value(model::ClpModel)
    _jl__check_model(model)

    @clp_ccall getSmallElementValue Float64 (Ptr{Cvoid},) model.p
end

function set_small_element_value(model::ClpModel,value::Float64)
    _jl__check_model(model)

    @clp_ccall setSmallElementValue Cvoid (Ptr{Cvoid},Float64) model.p value
end

# ClpSolve functions

function set_special_option(solve::ClpSolve, which::Integer, value::Integer, extraInfo::Integer)
    _jl__check_solve(solve)

    @clpsolve_ccall setSpecialOption Cvoid (Ptr{Cvoid},Int32,Int32,Int32) solve.p which value extraInfo
end

function get_special_option(solve::ClpSolve, which::Integer)
    _jl__check_solve(solve)

    @clpsolve_ccall getSpecialOption Int32 (Ptr{Cvoid},Int32) solve.p which
end

macro def_get_int_property(fname,clpname)
    quote
        function $(esc(fname))(solve::ClpSolve)
            _jl__check_solve(solve)
            @clpsolve_ccall $clpname Int32 (Ptr{Cvoid},) solve.p
        end
    end
end

macro def_set_int_property(fname,clpname)
    quote
        function $(esc(fname))(solve::ClpSolve, val::Integer)
            _jl__check_solve(solve)
            @clpsolve_ccall $clpname Cvoid (Ptr{Cvoid},Int32) solve.p val
        end
    end
end

# 0 - presolve on
# 1 - presolve off
# 2 - presolve number
# 3 - presolve number cost
function set_presolve_type(solve::ClpSolve, value::Integer)
    _jl__check_solve(solve)

    @clpsolve_ccall setPresolveType Cvoid (Ptr{Cvoid},Int32,Int32) solve.p value -1
end

@def_get_int_property get_presolve_type getPresolveType

# 0 - dual simplex
# 1 - primal simplex
# 2 - primal or sprint
# 3 - barrier
# 4 - barrier no crossover
# 5 - automatic
function set_solve_type(solve::ClpSolve, value::Integer)
    _jl__check_solve(solve)

    @clpsolve_ccall setSolveType Cvoid (Ptr{Cvoid},Int32,Int32) solve.p value -1
end

@def_get_int_property get_solve_type getSolveType

@def_get_int_property get_presolve_passes getPresolvePasses

function get_extra_info(solve::ClpSolve, which::Integer)
    _jl__check_solve(solve)

    @clpsolve_ccall getExtraInfo Int32 (Ptr{Cvoid},Int32) solve.p which
end

@def_set_int_property set_infeasible_return setInfeasibleReturn
@def_get_int_property infeasible_return infeasibleReturn

@def_get_int_property do_dual doDual
@def_set_int_property set_do_dual setDoDual

@def_get_int_property do_singleton doSingleton
@def_set_int_property set_do_singleton setDoSingleton

@def_get_int_property do_doubleton doDoubleton
@def_set_int_property set_do_doubleton setDoDoubleton

@def_get_int_property do_tripleton doTripleton
@def_set_int_property set_do_tripleton setDoTripleton

@def_get_int_property do_tighten doTighten
@def_set_int_property set_do_tighten setDoTighten

@def_get_int_property do_forcing doForcing
@def_set_int_property set_do_forcing setDoForcing

@def_get_int_property do_implied_free doImpliedFree
@def_set_int_property set_do_implied_free setDoImpliedFree

@def_get_int_property do_dupcol doDupcol
@def_set_int_property set_do_dupcol setDoDupcol

@def_get_int_property do_duprow doDuprow
@def_set_int_property set_do_duprow setDoDuprow

@def_get_int_property do_singleton_column doSingletonColumn
@def_set_int_property set_do_singleton_column setDoSingletonColumn

@def_get_int_property presolve_actions presolveActions
@def_set_int_property set_presolve_actions setPresolveActions

@def_get_int_property substitution substitution
@def_set_int_property set_substitution setSubstitution

end


text:
#=== compiler apis ===#
@ocl_func(clUnloadCompiler, CL_int, ())

#=== memory object apis ===#
@ocl_func(clCreateSubBuffer, CL_mem,
              (CL_mem, CL_mem_flags, CL_buffer_create_type, Ptr{Nothing}, Ptr{CL_int}))

@ocl_func(clSetMemObjectDestructorCallback, CL_int,
              (CL_mem, CL_callback, Ptr{Nothing}))

@ocl_func(clCreateImage2D, CL_mem,
               (CL_context, CL_mem_flags, Ptr{CL_image_format}, Csize_t, Csize_t, Csize_t,
                Ptr{Nothing}, Ptr{CL_int}))

@ocl_func(clCreateImage3D, CL_mem,
               (CL_context, CL_mem_flags, Ptr{CL_image_format}, Csize_t, Csize_t, Csize_t,
                Csize_t, Ptr{Nothing}, Ptr{CL_int}))

#=== program object apis ===#
@ocl_func(clGetProgramInfo, CL_int,
              (CL_program, CL_program_info, Csize_t, Ptr{Nothing}, Ptr{Csize_t}))

#=== event object apis ===#
@ocl_func(clCreateUserEvent, CL_event,
              (CL_context, Ptr{CL_int}))

@ocl_func(clSetUserEventStatus, CL_int, (CL_event, CL_int))

@ocl_func(clSetEventCallback, CL_int,
              (CL_event, CL_int, CL_callback, CL_user_data))

#=== enqueued commands apis ===#
@ocl_func(clEnqueueReadBufferRect, CL_int,
              (CL_command_queue, CL_mem, CL_bool,
               Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t},
               Csize_t, Csize_t, Csize_t, Csize_t,
               Ptr{Nothing}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueWriteBufferRect, CL_int,
              (CL_command_queue, CL_mem, CL_bool,
               Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t},
               Csize_t, Csize_t, Csize_t, Csize_t,
               Ptr{Nothing}, CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueCopyBufferRect, CL_int,
              (CL_command_queue, CL_mem, CL_mem,
               Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t},
               Csize_t, Csize_t, Csize_t, Csize_t,
               CL_uint, Ptr{CL_event}, Ptr{CL_event}))

@ocl_func(clEnqueueMarker, CL_int,
               (CL_command_queue, Ptr{CL_event}))

@ocl_func(clEnqueueWaitForEvents, CL_int,
               (CL_command_queue, CL_uint, Ptr{CL_event}))

@ocl_func(clEnqueueBarrier, CL_int,
               (CL_command_queue,))

#=== extension function access ===#
@ocl_func(clGetExtensionFunctionAddress, Ptr{Nothing}, (Ptr{Cchar},))

#=== opengl interop functions ===#

@ocl_func(clGetGLContextInfoKHR, CL_int,
              (Ptr{CL_context_properties}, CL_gl_context_info, Csize_t, Ptr{Nothing}, Ptr{Csize_t}))

@ocl_func(clCreateEventFromGLsyncKHR, CL_event,
              (CL_context, GL_sync, Ptr{CL_int}))


text:
# ---
# title: 1310. XOR Queries of a Subarray
# id: problem1310
# author: Indigo
# date: 2021-05-14
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/xor-queries-of-a-subarray/description/>
# hidden: true
# ---
#
# Given the array `arr` of positive integers and the array `queries` where
# `queries[i] = [Li, Ri]`, for each query `i` compute the **XOR** of elements
# from `Li` to `Ri` (that is, `arr[Li] **xor** arr[Li+1] **xor** ... **xor**
# arr[Ri]` ). Return an array containing the result for the given `queries`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,3,4,8], queries = [[0,1],[1,2],[0,3],[3,3]]
#     Output: [2,7,14,8]
#     Explanation:
#     The binary representation of the elements in the array are:
#     1 = 0001
#     3 = 0011
#     4 = 0100
#     8 = 1000
#     The XOR values for queries are:
#     [0,1] = 1 xor 3 = 2
#     [1,2] = 3 xor 4 = 7
#     [0,3] = 1 xor 3 xor 4 xor 8 = 14
#     [3,3] = 8
#
#
# **Example 2:**
#
#
#
#     Input: arr = [4,8,2,10], queries = [[2,3],[1,3],[0,0],[0,3]]
#     Output: [8,0,4,4]
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 3 * 10^4`
#   * `1 <= arr[i] <= 10^9`
#   * `1 <= queries.length <= 3 * 10^4`
#   * `queries[i].length == 2`
#   * `0 <= queries[i][0] <= queries[i][1] < arr.length`
#
#
## @lc code=start
using LeetCode

function xor_queries(arr::Vector{Int}, queries::Vector{Vector{Int}})
    pushfirst!(arr, 0)
    off_arr = OffsetArray(arr, -1)
    for i in 1:(length(arr) - 1)
        @inbounds off_arr[i] = off_arr[i] ⊻ off_arr[i - 1]
    end
    res = fill(0, length(queries))
    for (idx, q) in enumerate(queries)
        @inbounds res[idx] = off_arr[q[2] + 1] ⊻ off_arr[q[1]]
    end
    return res
end
## @lc code=end


text:
try Pkg.installed("Revise")
    using Revise
end

__precompile__(true)

module Flora

include("mlp.jl")
include("rnn.jl")

end # module Flora


text:
module DyadicKDE

using LinearAlgebra
using Distributions
using Statistics
using Convex
using Mosek
using MosekTools
using COSMO

include("kernels.jl")

include("estimation.jl")
export DyadicKernelDensityEstimator
export fit
export estimate_ROT_bandwidth

include("helpers.jl")
export make_data
export mean
export get_RIMSE
export get_ucb_coverage
export get_pci_coverage
export get_bci_coverage
export get_ucb_average_width
export get_pci_average_width
export get_bci_average_width

end


text:
RecipesBase.@recipe function f(y::Vector, forec::Forecast)
    n = length(y)
    steps_ahead = length(forec.observation_forecast)
    forec_idx = collect(n + 1: n + steps_ahead)
    fillrange_color = :steelblue
    # Plot the series
    @series begin
        seriestype := :path
        seriescolor := "black"
        return y
    end
    # Plot the forecast
    @series begin
        seriescolor := "blue"
        return forec_idx, forec.observation_forecast
    end
    # Plot the prediction interval
    q975 = Vector{Float64}(undef, steps_ahead)
    q025 = Vector{Float64}(undef, steps_ahead)
    q900 = Vector{Float64}(undef, steps_ahead)
    q100 = Vector{Float64}(undef, steps_ahead)
    for i in 1:steps_ahead
        q975[i] = quantile(forec.observation_scenarios[i, :], 0.975)
        q025[i] = quantile(forec.observation_scenarios[i, :], 0.025)
        q900[i] = quantile(forec.observation_scenarios[i, :], 0.9)
        q100[i] = quantile(forec.observation_scenarios[i, :], 0.1)
    end
    @series begin
        seriestype := :path
        seriescolor := fillrange_color
        fillcolor := fillrange_color
        fillalpha := 0.5
        fillrange := q975
        label := ""
        return forec_idx, q025
    end
    @series begin
        seriestype := :path
        seriescolor := fillrange_color
        fillcolor := fillrange_color
        fillalpha := 0.5
        fillrange := q900
        label := ""
        return forec_idx, q100
    end
end


text:
using Keys

import Documenter
Documenter.makedocs(
    modules = [Keys],
    format = :html,
    sitename = "Keys.jl",
    root = joinpath(dirname(dirname(@__FILE__)), "docs"),
    pages = Any["Home" => "index.md"],
    strict = true,
    checkdocs = :exports,
    authors = "Brandon Taylor"
)

using Test

x = (1, 2, 3)
t0 = ()
switch = (True(), True(), False())
new = (4, 5)

@test_throws ArgumentError reduce_unrolled(+, t0)

@test getindex_unrolled(x, t0) == t0
@test getindex_unrolled(t0, switch) == t0

@test setindex_unrolled(t0, t0, switch) == t0
@test setindex_unrolled(t0, new, t0) == t0
@test setindex_unrolled(t0, new, switch) == t0
@test setindex_unrolled(x, t0, t0) == x
@test setindex_unrolled(x, t0, switch) == x
@test setindex_unrolled(x, new, t0) == t0
@test setindex_unrolled(t0, t0, t0) == t0


text:
"
 Web socket client sample

 author: Atsushi Sakai
"

using HTTP

function main()
    println(PROGRAM_FILE, " start!!")

    HTTP.WebSockets.open("ws://127.0.0.1:8081") do ws
        write(ws, "Hello")

        while true
            x = readavailable(ws)
            @show x
            println(String(x))
        end
    end

    println(PROGRAM_FILE, " Done!!")
end

if abspath(PROGRAM_FILE) == @__FILE__
    @time main()
end

max_stars_repo_name: AtsushiSakai/JSONServerAndClient

max_issues_repo_name: AtsushiSakai/JSONServerAndClient

max_forks_repo_name: AtsushiSakai/JSONServerAndClient


text:
using NeuroCore.NeuroMetadata: ContrastIngredient

using NeuroCore.AnatomicalAPI

include("contrast_ingredient_tests.jl")
include("encoding_directions.jl")

@testset "Semantic positions" begin
    for (f,vars) in ((AnatomicalAPI.is_left, (:left, :L)),
                     (AnatomicalAPI.is_right, (:right, :R)),
                     (AnatomicalAPI.is_anterior, (:anterior, :A)),
                     (AnatomicalAPI.is_posterior, (:posterior, :P)),
                     (AnatomicalAPI.is_superior, (:superior, :S)),
                     (AnatomicalAPI.is_inferior, (:inferior, :I)),
                     (AnatomicalAPI.is_sagittal, (:sagittal, :left, :L, :right, :R)),
                     (AnatomicalAPI.is_coronal, (:coronal, :anterior, :A, :posterior, :P)),
                     (AnatomicalAPI.is_axial, (:axial, :superior, :S, :inferior, :I)),
                     (AnatomicalAPI.is_white_matter, (:white_matter,)),
                     (AnatomicalAPI.is_gyrus, (:gyrus,)),
                     (AnatomicalAPI.is_cortical, (:cortical,)),
                     (AnatomicalAPI.is_sulcus, (:sulcus,)),
                     #(AnatomicalAPI.is_csp, (:csp,)),
                    )
        for v in vars
            @test @inferred(f(v))
        end
    end
end


text:
using ArgParse
using BioSequences
using CodecZlib
using FASTX

cli = ArgParseSettings()
@add_arg_table cli begin
	"--window"
		help = "window size"
		arg_type = Int
		default = 100
	"fasta"
		help = "fasta file, required"
		required = true
end
arg = parse_args(ARGS, cli)

reader = nothing
if splitext(arg["fasta"])[2] == ".gz"
	reader = FASTA.Reader(GzipDecompressorStream(open(arg["fasta"])))
else
	reader = FASTA.Reader(open(arg["fasta"]))
end

## Main loop ##
w = arg["window"]
for record in reader
	chrom = FASTA.identifier(record)
	dna_seq = FASTA.sequence(record)
	println(length(dna_seq))
	for i = 1:length(dna_seq)-w
		g = 0
		c = 0
		for j = i:i+w
			if 	   dna_seq[j] == DNA_C c+=1
			elseif dna_seq[j] == DNA_G g+=1
			end
		end
		println("$chrom\t$i\t", (g-c)/(g+c))
	end
end

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
module OIFit

export +, ==, getindex, length, Vis2
export annulus, disk, gauss, pseudolorentzian

import Base: +, ==, getindex, length, merge

using CMPFit
using Measurements
using OIFITS
using Random
using SpecialFunctions
using Statistics

include("data.jl")
include("model.jl")
include("fit.jl")

end # module