text:
module WhenCases
using MLStyle
using MLStyle.Render

export @when, @otherwise, gen_when

@active WhenSpliter{s::String}(x) begin
    @match x begin
        Expr(:macrocall,
            &(Symbol("@", s)),
            ln::LineNumberNode,
            Expr(:block, elts...) || a && Do(elts = [a])
        ) => (ln, elts)
        Expr(:macrocall,
            &(Symbol("@", s)),
            ln::LineNumberNode
            # no args
        ) => (ln, [:(_)])
        _ => nothing
    end
end

function split_case_and_block(stmts, first_bindings, first_source)
    blocks :: Vector{Any} = []
    binding_seqs :: Vector{Any} = [first_bindings]
    sources :: Vector{LineNumberNode} = [first_source]
    current_block = []

    function make_block!()

        # avoid setting LineNumberNode in the end of block.
        block_size = length(current_block)
        take_size = block_size
        for i in block_size:-1:(block_size - 1)
            take_size = i
            if !(current_block[i] isa LineNumberNode)
                break
            end
        end
        # push current_block(cache) to current_block, then clear cache
        push!(blocks, Expr(:block, view(current_block, take_size)...))
        empty!(current_block)
        nothing
    end

    function append_block!(a)
        push!(current_block, a)
        nothing
    end

    for stmt in stmts
        @match stmt begin
            WhenSpliter{"when"}(source, bindings) =>
                begin
                    push!(binding_seqs, bindings)
                    push!(sources, source)
                    make_block!()
                end
            WhenSpliter{"otherwise"}(source, _) =>
                begin
                    push!(binding_seqs, [])
                    push!(sources, source)
                    make_block!()
                end
            a => append_block!(a)
        end
    end
    make_block!()
    # thus, the length of `sources`, `binding_seqs` and `blocks`
    # are guaranteed to be the same.
    collect(zip(sources, binding_seqs, blocks))
end

"""
Used in the `@when` block:

```
@when (a, b) = x begin
    a + b
@otherwise
    0
end
```
"""
macro otherwise()
    throw(SyntaxError("@otherwise is only used inside @when block, as a token to indicate default case."))
end

"""
Code generation for `@when`.
You should pass an `Expr(:let, ...)` as the first argument.
"""
function gen_when(let_expr, source :: LineNumberNode, mod :: Module)
    @match let_expr begin
        Expr(:let,
            Expr(:block, bindings...) || a && Do(bindings = [a]),
            Expr(:block, stmts...) || b && Do(stmts = [b])) =>

            begin
                sources_cases_blocks = split_case_and_block(stmts, bindings, source)
                # pprint(sources_cases_blocks)
                foldr(sources_cases_blocks, init=:nothing) do (source, bindings, block), last_block
                    foldr(bindings, init=block) do each, last_ret
                        @match each begin
                            :($a = $b) =>
                                let cbl = @format [source, a, last_ret, last_block] quote
                                            source
                                            a => last_ret
                                            _ => last_block
                                        end
                                    gen_match(b, cbl, source, mod)
                                    # :(@match $b $cbl)
                                end

                            # like:
                            # let a; a = 1; a end
                            new_source::LineNumberNode => begin
                                    source = new_source
                                    last_ret
                                end
                            :(if $a; $(_...) end) ||
                            :($a.?) => @format [source, a, last_ret, last_block] quote
                                    source
                                    a ? last_ret : last_block
                                end
                            a => @format [source, a, last_ret] quote
                                    source
                                    let a
                                        last_ret
                                    end
                                end
                        end
                    end
                end
            end

        a => let short_msg = SubString(string(a), 1, 20)
                throw(SyntaxError("Expected a let expression, got a `$short_msg` at $(string(source))."))
            end
    end
end

"""
1. Allow destructuring in binding sequences of `let` syntax.

In binding sequences,
- For bindings in form of `a = b` and `f(x) = y`, it's regarded as pattern matching here.
- For others like `@inline f(x) = 1`, it's the same as the original let binding(not pattern matching).

```julia
    @when let (a, 1) = x,
          [b, c, 5] = y
        (a, b, c)
    end
```

2. For a regular assignment, like

```julia
    @when (a, 2) = x begin
        # dosomething
    end
```

It's nothing different from

```julia
    @match x begin
        (a, 2) => # dosomething
        _ => nothing
    end
```

3. For multiple branches, we can have

```julia
x = 1
@when let (_, _) = x
    :tuple
@when begin ::Float = x end
    :float
@when ::Int = x
    :int
@otherwise
    :unknown
end
```

4. Also, you can use predicates when matching, in the form
of `cond.?` or `if cond end`:

```julia
x = 1
y = (1, 2)
cond1 = true
cond2 = true
@when let cond1.?,
          (a, b) = x
    a + b
@when begin if cond2 end
            (a, b) = y
      end
    a + b
end
```
"""
macro when(let_expr)
    gen_when(let_expr, __source__, __module__) |> esc
end

macro when(assignment, ret)
    @match assignment begin
        :($_ = $_) =>
            let let_expr = Expr(:let, Expr(:block, assignment), ret)
                gen_when(let_expr, __source__, __module__) |> esc
            end
        _ => @syntax_err "Not match the form of `@when a = b expr`"
    end
end

end

ext: jl


text:
# root(x::EzXML.Document) exported from EzXML
xhtmlroot(x::EzXML.Document) = root(x)
xhtmlhead(x::EzXML.Document) = elements(root(x))[1]
xhtmlbody(x::EzXML.Document) = elements(root(x))[2]

function xhtml2html(x::EzXML.Document)
    str = string(x)
    idx = findfirst("<!DOCTYPE", str)[1]
    return Gumbo.parsehtml(SubString(str,idx,length(str)))
end

function nelements(x::EzXML.Node)
    if hasparentelement(x)
        x = parentelement(x)
    end
    return EzXML.countelements(x)
end

haschildnode(node) = EzXML.hasnode(node)
haschildelement(node) = EzXML.haselement(node)

#=
    root = ElementNode("root")
    x = ElementNode("x")
    first_child!(root, x)
    y = ElementNode("y")
    first_child!(root, y)
    z = ElementNode("z")
    first_child!(root, z)
=#

"""
    first_child!(parent::Node, child::Node)

    Link `child` as the first of the children of `parent`.
"""
function first_child!(parent::EzXML.Node, child::EzXML.Node)
    EzXML.check_topmost(child)
    if EzXML.hasnode(parent)
        EzXML.linkprev!(EzXML.firstnode(parent), child)
    else
        EzXML.link!(parent, child)
    end
    return child
end

"""
     last_child!(parent::Node, child::Node)

     Link `child` as the last of the children of `parent`.
"""
function last_child!(parent::EzXML.Node, child::EzXML.Node)
    EzXML.check_topmost(child)
    EzXML.link!(parent, child)
    return child
end

"""
    next_child!(curr_node, next_child)

    Link `next_child` as the child that follows the parent-level `curr_node`.
"""
function next_child!(parent::EzXML.Node, child::EzXML.Node)
    EzXML.check_topmost(child)
    if EzXML.hasnode(parent)
        EzXML.linknext!(EzXML.firstnode(parent), child)
    else
        EzXML.link!(parent, child)
    end
    return child
end

"""
    prev_child!(curr_node, prev_child)

    Link `prev_child` as the child that preceeds the parent-level `curr_node`.
"""
#=
function prev_child!(curr_node::EzXML.Node, prev_child::EzXML.Node)
    if EzXML.hasprevnode(curr_node)
        curr_node = EzXML.prevnode(curr_node)
        if EzXML.hasnode(curr_node)
            last_sibling!(EzXML.lastnode(curr_node), prev_child)
        else
            throw(ErrorException("no previous parent node"))
        end
    else
        throw(ErrorException("no previous parent node"))
    end
    return prev_child
end
=#

"""
    next_sibling!(curr_sibling, next_sibling)

    Link `next_sibling` as the next sibling of `curr_sibling`.
"""
function next_sibling!(curr_sibling::EzXML.Node, next_sibling::EzXML.Node)
    EzXML.check_topmost(next_sibling)
    linknext!(curr_sibling, next_sibling)
    return next_sibling
end

"""
    prev_sibling!(curr_sibling, prev_sibling)

    Link `prev_sibling` as the prev sibling of `curr_sibling`.
"""
function prev_sibling!(curr_sibling::EzXML.Node, prev_sibling::EzXML.Node)
    EzXML.check_topmost(prev_sibling)
    linkprev!(curr_sibling, prev_sibling)
    return prev_sibling
end

"""
    first_sibling!(curr_sibling::EzXML.Node, first_sibling::EzXML.Node)

    Link `first_sibling` as the first sibling of the siblings to which `curr_sibling` belongs.
"""
function first_sibling!(curr_sibling::EzXML.Node, first_sibling::EzXML.Node)
    EzXML.check_topmost(first_sibling)
    initial_sibling_node = EzXML.firstnode(parent(curr_sibling))
    EzXML.linkprev!(initial_sibling_node, first_sibling)    # link first_sibling as prev sibling of initial_sibling_node
    return first_sibling
end

"""
    last_sibling!(curr_sibling::EzXML.Node, last_sibling::EzXML.Node)

    Link `last_sibling` as the last sibling of the siblings to which `curr_sibling` belongs.
"""
function last_sibling!(curr_sibling::EzXML.Node, last_sibling::EzXML.Node)
    EzXML.check_topmost(last_sibling)
    final_sibling_node = EzXML.lastnode(parent(curr_sibling))
    EzXML.linknext!(final_sibling_node, last_sibling)    # link last_sibling as next sibling of final_sibling_node
    return last_sibling
end

#=
    manipulation

        for Nodes, for Elements
          first, next, prev, last
          first_ptr, next_ptr, prev_ptr, last_ptr

        for Nodes
          link!: introduce to a parent a new, final child
          parent_ptr(Node)

        to provide
            linkfirstchild!, linklastchild!
            linkprevsibling!, linknextsibling!
            linkprevchild!, linknextchild!
            linkfirstsibling!, linknextsibling!

    link!(parent::Node, child::Node)
    Link `child` at the end of children of `parent`.
    linklast!(parent::Node, child::Node)

    linkfirst!(parent::Node, child::Node)
    Link `child` at the start of children of `parent`.

    linknext!(target::Node, node::Node)
    Link `node` as the next sibling of `target`.

    linkprev!(target::Node, node::Node)
    Link `node` as the prev sibling of `target`.

    unlink!(node::Node)
    Unlink `node` from its context.

    addelement!(parent::Node, name::AbstractString)
    Add a new child element of `name` with no content to `parent`
      and return the new child element.

    addelement!(parent::Node, name::AbstractString, content::AbstractString)
    Add a new child element of `name` with `content` to `parent`
      and return the new child element.

    nodecontent(node::Node)
    Return the node content of `node`.

    setnodecontent!(node::Node, content::AbstractString)
    Replace the content of `node`.

    # create a vector of _
    nodes(node::Node, [backward=false])
    Create a vector of child nodes.

    elements(node::Node, [backward=false])
    Create a vector of child elements.

    attributes(node::Node)
    Create a vector of attributes.

    # attributes

    haskey(node::Node, attr::AbstractString)
    getindex(node::Node, attr::AbstractString)
    setindex!(node::Node, val, attr::AbstractString)
    delete!(node::Node, attr::AbstractString)

=#

#=
prettyxml(x::EzXML.Document) = prettyxml(root(x))

function prettyxml(x::EzXML.Node)
    io = IOBuffer()
    EzXML.print(io, x)
    str = String(take!(io))
    retuprevrettyxml(str)
e

    Unlink `n` from its context.nd

function prettyxml(x::AbstractString)
    s = replace(x, "><" => ">\n<")
    strs = strip.(String.(split(s, "\n")))
    prettystrs = Vector{String}(undef, length(strs))
    indent = ""
    idx = 1
    for str in strs
        open_tag = openxmltag(str)
        bothtags = !isnothing(open_tag) && (str[end-1:end] == "/>" || occursin(string("</",open_tag,">"), str))
        if bothtags
           s = string(indent, str)
        else
           close_tag = closexmltag(str)
           if !isnothing(open_tag)
               s = string(indent, str)
               indent = string(indent, "  ")
           elseif !isnothing(close_tag)
               indent = indent[1:end-2]
               s = string(indent, str)
           else
               s  = string(indent, str)
           end
       end
       prettystrs[idx] = s
       idx += 1
    end
    return join(prettystrs, "\n")
end

function openxmltag(str::AbstractString)
   (length(str) < 2 || str[1] !== '<' || str[2] === '/') && return nothing
   firstspace = findfirst(" ", str)
   if isnothing(firstspace)
       firstclose = findfirst(">", str)
       if isnothing(firstclose)
           tag = str[2:end]
       else
           tag = str[2:first(firstclose)-1]
       end
   else
       tag = str[2:first(firstspace)-1]
   end
   return tag
end

function closexmltag(str::AbstractString)
   (length(str) < 2 || str[1] !== '<' || str[2] !== '/') && return nothing
   firstspace = findfirst(" ", str)
   if isnothing(firstspace)
       if str[end] === '>'
           tag = str[3:end-1]
       else
           tag = str[3:end]
       end
   else
       tag = str[3:first(firstspace)-1]
   end
   return tag
end
=#

#=
xhtmlstr = """
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">

  <meta charset="utf-8" />
  <meta name="viewport" content="width=device-width, initial-scale=1" />
  <title>Title</title>
  <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/bulma/0.7.5/css/bulma.min.css" />

      these words

      © 2019 Jeffrey Sarnoff. All Rights Reserved.

"""
=#


text:
@testset "MFA" begin
    Random.seed!(28374544);

    # generate dataset
    d = 2
    theta = [0,3,5,7] * (pi/4)
    k = length(theta)
    X = hcat([MultivariateMixtures.clusters(t, 250, d) for t in theta]...)

    FA1 = fit(FactorAnalysis, X[:,1:250])
    MFA = fit_mm(FactorAnalysis, X, k, tol=1e-5)
    @test MultivariateMixtures.component_type(MFA) <: FactorAnalysis
    @test ncomponents(MFA) == k
    @test probs(MFA) ≈ fill(1/4, k) atol=1e-2
    FA2 = component(MFA,3)
    @test indim(FA2) == d
    @test outdim(FA2) == 1
    @test mean(FA1) ≈ mean(FA2) atol=1e-2
    @test cov(FA1) ≈ cov(FA2) atol=1e-2
    @test_broken length(MFA) == d
    @test_broken size(MFA) == (d,)
    @test_broken mean(MFA) ≈ mean(X, dims=2) atol=1e-2
    @test_broken cov(MFA) ≈ cov(X, dims=2) atol=1e-1

    X = hcat([MultivariateMixtures.clusters(t, 250, 10) for t in theta]...)
    FA1 = fit(FactorAnalysis, X[:,1:250], method=:em)
    MFA = fit_mm(FactorAnalysis, X, k, tol=1e-5)
    FA2 = component(MFA,1)
    @test mean(FA1) ≈ mean(FA2) atol=0.01
    @test cov(FA1) ≈ cov(FA2) atol=0.1

    MFA = fit_mm(FactorAnalysis, X, k, factors=3, tol=1e-5, maxiter=10)
    @test outdim(component(MFA,3)) == 3
end


text:
using SageMath
using Test

@testset "SageMath.jl" begin
    # Write your tests here.
end


text:
struct CrystalsCrystal{T}
    data::Set{T}
end
function Base.show(io::IO, cr::CrystalsCrystal)
    print(io, "(crystals.crystal "*join([aetherlang_repr(x) for x in cr.data], ' ')*")")
end

# overload indexin to make `contains?` from base lib work
Base.indexin(v, cr::CrystalsCrystal)::Vector{Union{Nothing, Int64}} = v in cr.data ? [1] : [nothing]
# overload length to make `len` from base lib work
Base.length(cr::CrystalsCrystal)::Int = length(cr.data)

crystaltype(::Type{CrystalsCrystal{T}}) where {T} = T

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


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

name = "glmnet"
version = v"5.0.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JuliaStats/GLMNet.jl.git", "d8cab55556a8c184a879446476efcc35e68a3eee"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/GLMNet.jl/deps/

flags="-fdefault-real-8 -ffixed-form -shared -O3"
if [[ ${target} != *mingw* ]]; then
    flags="${flags} -fPIC"
fi
if [[ ${target} != aarch64* ]] && [[ ${target} != arm* ]]; then
    flags="${flags} -m${nbits}"
fi

${FC} ${LDFLAGS} ${flags} glmnet5.f90 -o libglmnet.${dlext}

mkdir -p $libdir
mv libglmnet.${dlext} $libdir/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libglmnet", :libglmnet)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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


text:
using Base.Test
using Lora

println("    Testing BasicUnvVariableNState constructors and methods...")

nstatev = Float64[1.25, -4.4, 7.5]
nstate = BasicUnvVariableNState(nstatev)

@test eltype(nstate) == Float64
@test nstate.value == nstatev
@test nstate.n == length(nstatev)

statev = Float64(5.2)
savei = 2

copy!(nstate, BasicUnvVariableState(statev), savei)
@test nstate.value == [nstatev[1:savei-1]; statev; nstatev[savei+1:end]]

nstatev = Float32[-2.17, 1.92, -0.15, -0.65]
nstate = BasicUnvVariableNState(nstatev)

@test eltype(nstate) == Float32
@test nstate.value == nstatev
@test nstate.n == length(nstatev)

statev = Float32(-2.12)
savei = 4

copy!(nstate, BasicUnvVariableState(statev), savei)
@test nstate.value == [nstatev[1:savei-1]; statev; nstatev[savei+1:end]]

nstaten = 10
nstate = BasicUnvVariableNState(nstaten, BigFloat)

@test eltype(nstate) == BigFloat
@test nstate.n == nstaten

println("    Testing BasicMuvVariableNState constructors and methods...")

nstatev = Float64[1.35 3.7 4.5; 5.6 8.81 9.2]
nstate = BasicMuvVariableNState(nstatev)

@test eltype(nstate) == Float64
@test nstate.value == nstatev
@test nstate.size == size(nstatev, 1)
@test nstate.n == size(nstatev, 2)

statev = Float64[5.2, 3.31]
savei = 2

copy!(nstate, BasicMuvVariableState(statev), savei)
@test nstate.value == [nstatev[:, 1:savei-1] statev nstatev[:, savei+1:end]]

nstatev= BigFloat[
  0.41257   0.106756   0.817916   0.569789  0.54802;
  0.630804  0.0212354  0.0729593  0.483741  0.596365;
  0.82968   0.4872     0.185226   0.354095  0.944551;
  0.634923  0.448942   0.300905   0.243899  0.126606
]
nstate = BasicMuvVariableNState(nstatev)

@test eltype(nstate) == BigFloat
@test nstate.value == nstatev
@test nstate.size == size(nstatev, 1)
@test nstate.n == size(nstatev, 2)

statev = BigFloat[0.0646775, 0.379354, 0.0101067, 0.821756]
savei = 1

copy!(nstate, BasicMuvVariableState(BigFloat[0.0646775, 0.379354, 0.0101067, 0.821756]), 1)
@test nstate.value == [nstatev[:, 1:savei-1] statev nstatev[:, savei+1:end]]

nstatesize = 3
nstaten = 10
nstate = BasicMuvVariableNState(nstatesize, nstaten, Float16)

@test eltype(nstate) == Float16
@test nstate.size == nstatesize
@test nstate.n == nstaten

println("    Testing BasicMavVariableNState constructors and methods...")

nstatesize = (2, 4)
nstaten = 5
nstatev = Array(Float64, nstatesize..., nstaten)
nstatev[:, :, 1] = Float64[
   0.680789  -0.194683   1.86498    0.490497;
  -0.730417   0.305873  -0.0434663  0.879241
]
nstatev[:, :, 2] = Float64[
  -1.63194  -0.257043  -0.981173  0.524005;
  -1.61526   0.173245  -0.677052  1.88443
]
nstatev[:, :, 3] = Float64[
  -0.131442  -1.42516   -0.267594  -1.45806;
   0.544324  -0.109355  -0.219908   1.07273
]
nstatev[:, :, 4] = Float64[
  0.134486   2.0883    -0.455804  -2.20976;
  0.468841  -0.552023  -0.837046  -0.183255
]
nstatev[:, :, 5] = Float64[
  -1.3907   -1.18854  0.985564   0.373107;
   1.21794  -1.04891  0.611239  -0.526945
]
nstate = BasicMavVariableNState(nstatev)

@test eltype(nstate) == Float64
@test nstate.value == nstatev
@test nstate.size == nstatesize
@test nstate.n == nstaten

statev = Float64[
   1.36032    0.200834  -1.83856  -1.3039;
  -0.641921  -1.31766    1.35137   0.938878
]
savei = 3

copy!(nstate, BasicMavVariableState(statev), savei)
for i in [1:savei-1; savei; savei+1:nstaten]
  @test nstate.value[:, :, i] == nstatev[:, :, i]
end
@test nstate.value[:, :, savei] == statev

nstatesize = (3, 4)
nstaten = 2
nstatev = Array(Float16, nstatesize..., nstaten)
nstatev[:, :, 1] = Float16[
  -0.372612   0.00410037  -1.11811   -0.278956;
  -0.744771  -1.63419     -2.12376    1.161;
   1.0475    -1.9372       0.198637  -1.45228
]
nstatev[:, :, 2] = Float16[
  -0.124304  -1.60967   -0.0316585  -1.06222;
  -1.18789    0.849805   0.442735    0.454251;
   0.356338   0.60818    0.0970741   0.433813
]
nstate = BasicMavVariableNState(nstatev)

@test eltype(nstate) == Float16
@test nstate.value == nstatev
@test nstate.size == nstatesize
@test nstate.n == nstaten

statev = Float16[
   0.538734  -2.18586     1.4828     0.0479034;
   0.225186   0.28757    -0.580649   0.535829;
  -1.92379    0.0638309  -0.205978  -0.604045
]
savei = 2

copy!(nstate, BasicMavVariableState(statev), savei)
for i in [1:savei-1; savei; savei+1:nstaten]
  @test nstate.value[:, :, i] == nstatev[:, :, i]
end
@test nstate.value[:, :, savei] == statev

nstatesize = (3, 5)
nstaten = 12
nstate = BasicMavVariableNState(nstatesize, nstaten, BigFloat)

@test eltype(nstate) == BigFloat
@test nstate.size == nstatesize
@test nstate.n == nstaten


text:
function main()
    A, B, C = parse.(Float64, split(readline()))
    g(t; A=A, B=B, C=C) = A * t + B * sin(C * t * π) - 100
    t̲, t̄ = 0, (B + 100) / A
    ε = 1e-6

    t = (t̄ + t̲) / 2
    while abs(g(t)) > ε
        if g(t) > 0
            t̄ = t
        else
            t̲ = t
        end
        t = (t̄ + t̲) / 2
    end

    println(t)
    println()
end

main()


text:
include("node.jl")
include("arg.jl")
include("dash.jl")
include("optional.jl")
include("options.jl")
include("vararg.jl")


text:
export UA

mutable struct UA

    ntau :: Int64
    ε    :: Float64
    tau  :: Vector{Float64}
    ltau :: Vector{Float64}
    ftau :: Array{ComplexF64,1}
    ptau :: FFTW.cFFTWPlan{ComplexF64,-1,false,1}
    pl   :: Array{ComplexF64,2}
    ql   :: Array{ComplexF64,2}
    r    :: Array{ComplexF64,2}
    r̃    :: Array{ComplexF64,2}
    rtau :: FFTW.cFFTWPlan{ComplexF64,-1,false,2}

    function UA( ntau, ε, nbpart )

        dtau = 2π / ntau

        ltau  = zeros(Float64, ntau)
        ltau .= vcat(0:ntau÷2-1, -ntau÷2:-1)

        tau   = zeros(Float64, ntau)
        tau  .= [ i*dtau for i=0:ntau-1 ]

        ftau  = zeros(ComplexF64, ntau)

        ptau  = FFTW.plan_fft(ftau)

        pl = zeros(ComplexF64, (ntau, nbpart))
        ql = zeros(ComplexF64, (ntau, nbpart))

        r  = zeros(ComplexF64, (ntau,2))
        r̃  = zeros(ComplexF64, (ntau,2))

        rtau = plan_fft(r,1)

        new( ntau, ε, tau, ltau, ftau, ptau, pl, ql, r, r̃, rtau )

    end

end


text:
# Implement GreedyMotifSearch
using DataStructures

const NUCLEOTIDES = ('A', 'C', 'G', 'T')

struct ProfileColumn
    A::Float64
    C::Float64
    G::Float64
    T::Float64
end

function Base.argmax(column::ProfileColumn)
    ((getfield(column, symbol), String(symbol)) for symbol in fieldnames(ProfileColumn)) |> maximum |> last
end

convert_to_intlist(line) = map(x -> parse(Int, x), split(line))

calc_hamming_distance(s1, s2) = count(((c1, c2),) -> c1 != c2, zip(s1, s2))

function create_profile_matrix_from_motifs(motifs)
    nr_motifs = length(motifs)
    profile_matrix = ProfileColumn[]
    for column in zip(motifs...)
        counts = counter(column)
        vals = (get(counts, nucleotide, 0) / nr_motifs for nucleotide in NUCLEOTIDES)
        push!(profile_matrix, ProfileColumn(vals...))
    end
    profile_matrix
end

function calc_profile_matrix_score(motifs)
    profile_matrix = create_profile_matrix_from_motifs(motifs)
    consensus = join(map(argmax, profile_matrix))
    foldl((score, motif) -> score + calc_hamming_distance(motif, consensus), motifs; init=0)
end

calc_kmer_probability(kmer, profile_matrix) =
    foldl((p, (nucleotide, column)) -> p * getfield(column, Symbol(nucleotide)), zip(kmer, profile_matrix); init=1.0)

function profile_most_probable_kmer(text, profile_matrix, k)
    max_probability = 0.0
    most_probable_kmer = text[1:k]
    for ix in 1:(length(text)-k+1)
        kmer = text[ix:ix+k-1]
        p = calc_kmer_probability(kmer, profile_matrix)
        if p > max_probability
            max_probability = p
            most_probable_kmer = kmer
        end
    end
    most_probable_kmer
end

function greedy_motif_search(texts, k)
    best_motifs = map(text -> text[1:k], texts)
    best_score = calc_profile_matrix_score(best_motifs)
    for ix in 1:(length(texts[1])-k+1)
        kmer = texts[1][ix:ix+k-1]
        motifs = [kmer]
        for text in texts[2:end]
            profile = create_profile_matrix_from_motifs(motifs)
            motif = profile_most_probable_kmer(text, profile, k)
            push!(motifs, motif)
        end
        score = calc_profile_matrix_score(motifs)
        if score < best_score
            best_score = score
            best_motifs = motifs
        end
    end
    best_motifs
end

function main()
    k, t = convert_to_intlist(readline())
    texts = map(_ -> readline(), 1:t)
    result = greedy_motif_search(texts, k)
    foreach(println, result)
end

main()


text:
import Base.isempty

export ConvexHullArray, CHArray,
       array

# ================================
#  Convex hull of an array of sets
# ================================

"""
    ConvexHullArray{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the symbolic convex hull of a finite number of convex sets.

### Fields

- `array` -- array of sets

### Notes

The `EmptySet` is the neutral element for `ConvexHullArray`.

Constructors:

- `ConvexHullArray(array::Vector{<:LazySet})` -- default constructor

- `ConvexHullArray([n]::Int=0, [N]::Type=Float64)`
  -- constructor for an empty hull with optional size hint and numeric type

### Examples

Convex hull of 100 two-dimensional balls whose centers follows a sinusoidal:

```jldoctest
julia> b = [Ball2([2*pi*i/100, sin(2*pi*i/100)], 0.05) for i in 1:100];

julia> c = ConvexHullArray(b);
```
"""
struct ConvexHullArray{N<:Real, S<:LazySet{N}} <: LazySet{N}
    array::Vector{S}
end

isoperationtype(::Type{<:ConvexHullArray}) = true
isconvextype(::Type{<:ConvexHullArray}) = true

# constructor for an empty hull with optional size hint and numeric type
function ConvexHullArray(n::Int=0, N::Type=Float64)
    a = Vector{LazySet{N}}()
    sizehint!(a, n)
    return ConvexHullArray(a)
end

# EmptySet is the neutral element for ConvexHullArray
@neutral(ConvexHullArray, EmptySet)

# Universe is the absorbing element for ConvexHullArray
@absorbing(ConvexHullArray, Universe)

# add functions connecting ConvexHull and ConvexHullArray
@declare_array_version(ConvexHull, ConvexHullArray)

"""
    CHArray

Alias for `ConvexHullArray`.
"""
const CHArray = ConvexHullArray

"""
    array(cha::ConvexHullArray{N, S}) where {N<:Real, S<:LazySet{N}}

Return the array of a convex hull of a finite number of convex sets.

### Input

- `cha` -- convex hull array

### Output

The array of a convex hull of a finite number of convex sets.
"""
function array(cha::ConvexHullArray{N, S}) where {N<:Real, S<:LazySet{N}}
    return cha.array
end

"""
    dim(cha::ConvexHullArray)

Return the dimension of the convex hull of a finite number of convex sets.

### Input

- `cha` -- convex hull array

### Output

The ambient dimension of the convex hull of a finite number of convex sets.
"""
function dim(cha::ConvexHullArray)
    @assert !isempty(cha.array)
    return dim(cha.array[1])
end

"""
    σ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}

Return the support vector of a convex hull array in a given direction.

### Input

- `d`   -- direction
- `cha` -- convex hull array
"""
function σ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}
    s = σ(d, cha.array[1])
    ri = dot(d, s)
    rmax = ri
    for (i, chi) in enumerate(cha.array[2:end])
        si = σ(d, chi)
        ri = dot(d, si)
        if ri > rmax
            rmax = ri
            s = si
        end
    end
    return s
end

"""
    ρ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}

Return the support function of a convex hull array in a given direction.

### Input

- `d`   -- direction
- `cha` -- convex hull array

### Output

The support function of the convex hull array in the given direction.

### Algorithm

This algorihm calculates the maximum over all ``ρ(d, X_i)`` where the
``X_1, …, X_k`` are the sets in the array `cha`.
"""
function ρ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}
    return maximum([ρ(d, Xi) for Xi in array(cha)])
end

"""
    isbounded(cha::ConvexHullArray)

Determine whether a convex hull of a finite number of convex sets is
bounded.

### Input

- `cha` -- convex hull of a finite number of convex sets

### Output

`true` iff all wrapped sets are bounded.
"""
function isbounded(cha::ConvexHullArray)
    return all(isbounded, cha.array)
end

"""
    isempty(cha::ConvexHullArray)

Return if a convex hull array is empty or not.

### Input

- `cha` -- convex hull array

### Output

`true` iff all wrapped sets are empty.
"""
function isempty(cha::ConvexHullArray)
    return all(isempty, array(cha))
end

"""
    vertices_list(X::ConvexHullArray{N, Singleton{N, VT}}) where {N, VT}

Return the list of vertices of the convex hull array of singletons.

### Input

- `X` -- convex hull array of singletons

### Output

The list of elements in the array that defines `X`.
"""
function vertices_list(X::ConvexHullArray{N, Singleton{N, VT}}) where {N, VT}
    m = length(X.array)
    vertices = Vector{VT}(undef, m)
    @inbounds for i in 1:m
        vertices[i] = X.array[i].element
    end
    return vertices
end

# list of constraints of the convex hull array of singletons
function constraints_list(X::ConvexHullArray{N, Singleton{N, VT}}) where {N, VT}
    n = dim(X)
    ST = n == 2 ? VPolygon : VPolytope
    V = convert(ST, X)
    return constraints_list(V)
end

# membership in convex hull array of singletons
function ∈(x::AbstractVector{N}, X::ConvexHullArray{N, Singleton{N, VT}}) where {N, VT}
    n = length(x)
    ST = n == 2 ? VPolygon : VPolytope
    V = convert(ST, X)
    return x ∈ V
end


text:
# This Map has non-trivial domain, thus we need the define testargs
"""
    PosNegReindex(values_pos,values_neg)
"""
struct PosNegReindex{A,B} <: Map
  values_pos::A
  values_neg::B
end

function testargs(k::PosNegReindex,i::Integer)
  @check length(k.values_pos) !=0 || length(k.values_neg) != 0 "This map has empty domain"
  @check eltype(k.values_pos) == eltype(k.values_neg) "This map is type-instable"
  length(k.values_pos) !=0 ? (one(i),) : (-one(i))
end

function return_value(k::PosNegReindex,i::Integer)
  if length(k.values_pos)==0 && length(k.values_neg)==0
    @check eltype(k.values_pos) == eltype(k.values_neg) "This map is type-instable"
    testitem(k.values_pos)
  else
    evaluate(k,testargs(k,i)...)
  end
end

function return_cache(k::PosNegReindex,i::Integer)
  c_p = array_cache(k.values_pos)
  c_n = array_cache(k.values_neg)
  c_p, c_n
end

function evaluate!(cache,k::PosNegReindex,i::Integer)
  c_p, c_n = cache
  i>0 ? getindex!(c_p,k.values_pos,i) : getindex!(c_n,k.values_neg,-i)
end

function evaluate(k::PosNegReindex,i::Integer)
  i>0 ? k.values_pos[i] : k.values_neg[-i]
end

function lazy_map(k::PosNegReindex,::Type{T},i_to_iposneg::AbstractArray) where T
  if all_pos(i_to_iposneg)
    lazy_map(Reindex(k.values_pos),get_array(i_to_iposneg))
  elseif all_neg(i_to_iposneg)
    lazy_map(Reindex(k.values_neg),lazy_map(iposneg->-iposneg,get_array(i_to_iposneg)))
  else
    s = size(i_to_iposneg)
    lazy_map(evaluate,T,Fill(k, s), i_to_iposneg)
  end
end

function lazy_map(::typeof(evaluate),::Type{T},a::LazyArray{<:Fill{<:PosNegReindex}}...) where T
  i_to_iposneg = a[1].args[1]
  if all(map( ai-> is_exhaustive(a[1].args[1]),a)) && all( map( ai-> i_to_iposneg==a[1].args[1],a) )
    bpos = map(ai->ai.maps.value.values_pos,a)
    bneg = map(ai->ai.maps.value.values_neg,a)
    cpos = lazy_map(evaluate,T,bpos...)
    cneg = lazy_map(evaluate,T,bneg...)
    lazy_map(PosNegReindex(cpos,cneg),T,i_to_iposneg)
  else
    LazyArray(T,a...)
  end
end

function lazy_map(::typeof(evaluate),::Type{T},a::LazyArray{<:Fill{<:PosNegReindex}},x::AbstractArray) where T
  i_to_iposneg = a.args[1]
  _lazy_map_evaluate_posneg(i_to_iposneg,T,a,x)
end

function _lazy_map_evaluate_posneg(i_to_iposneg,::Type{T},a,x) where T
  apos = a.maps.value.values_pos
  aneg = a.maps.value.values_neg
  if is_exhaustive(i_to_iposneg)
    ipos_to_i, ineg_to_i = pos_and_neg_indices(i_to_iposneg)
    xpos = lazy_map(Reindex(x),ipos_to_i)
    xneg = lazy_map(Reindex(x),ineg_to_i)
    cpos = lazy_map(evaluate,apos,xpos)
    cneg = lazy_map(evaluate,aneg,xneg)
    lazy_map(PosNegReindex(cpos,cneg),T,i_to_iposneg)
  else
    LazyArray(T,a,x)
  end
end

function lazy_map(::typeof(evaluate),::Type{T},b::Fill,a::LazyArray{<:Fill{<:PosNegReindex}}...) where T
  i_to_iposneg = a[1].args[1]
  if all(map( ai-> is_exhaustive(ai.args[1]),a)) && all( map( ai-> i_to_iposneg==ai.args[1],a) )
    k = b.value
    bpos = map(ai->ai.maps.value.values_pos,a)
    bneg = map(ai->ai.maps.value.values_neg,a)
    cpos = lazy_map(k,T,bpos...)
    cneg = lazy_map(k,T,bneg...)
    lazy_map(PosNegReindex(cpos,cneg),T,i_to_iposneg)
  else
    LazyArray(T,b,a...)
  end
end

# Helper functions to work with arrays representing a binary partition

function pos_and_neg_length(i_to_iposneg)
  Npos = maximum(i_to_iposneg)
  Nneg = minimum(i_to_iposneg)
  Nneg = Nneg ==  1  ? 0 : -Nneg
  Npos == -1  ? 0 : Npos , Nneg
end

function is_exhaustive(i_to_iposneg)
  Npos, Nneg = pos_and_neg_length(i_to_iposneg)
  if length(i_to_iposneg) != Npos+Nneg
    return false
  end
  ipos_to_touched = fill(false,Npos)
  ineg_to_touched = fill(false,Nneg)
  for iposneg in i_to_iposneg
    iposneg>0 ? ipos_to_touched[iposneg] = true : ineg_to_touched[-iposneg] = true
  end
  all(ipos_to_touched) && all(ineg_to_touched)
end

function pos_and_neg_indices(i_to_iposneg)
  @check is_exhaustive(i_to_iposneg)
  Npos, Nneg = pos_and_neg_length(i_to_iposneg)
  ipos_to_i = zeros(Int,Npos)
  ineg_to_i = zeros(Int,Nneg)
  @inbounds for (i,iposneg) in enumerate(i_to_iposneg)
    iposneg>0 ? ipos_to_i[iposneg] = i : ineg_to_i[-iposneg] = i
  end
  ipos_to_i, ineg_to_i
end

function aligned_with_pos(i_to_iposneg,j_to_i,npos)
  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
  j_to_iposneg == 1:npos
end

function aligned_with_neg(i_to_iposneg,j_to_i,nneg)
  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
  j_to_iposneg == -(1:nneg)
end

function all_in_pos(i_to_iposneg,j_to_i)
  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
  all_pos(j_to_iposneg)
end

function all_in_neg(i_to_iposneg,j_to_i)
  j_to_iposneg = lazy_map(Reindex(i_to_iposneg),j_to_i)
  all_neg(j_to_iposneg)
end

all_pos(i_to_iposneg) = all( lazy_map(iposneg->iposneg>0, i_to_iposneg))
all_neg(i_to_iposneg) = all( lazy_map(iposneg->iposneg<0, i_to_iposneg))

# This is important to do optimizations associated with ExtendedFESpace
"""
struct representing a binary partition of a range of indices

Using this allows one to do a number of important optimizations when working with `PosNegReindex`
"""
struct PosNegPartition{T,V,Vp,Vn} <: AbstractVector{T}
  i_to_iposneg::V
  ipos_to_i::Vp
  ineg_to_i::Vn
  function PosNegPartition(
    i_to_iposneg::AbstractVector,
    ipos_to_i::AbstractVector,
    ineg_to_i::AbstractVector)

    @check eltype(ipos_to_i) == eltype(ineg_to_i)
    @check all((ipos_to_i, ineg_to_i) .== pos_and_neg_indices(i_to_iposneg))
    T = eltype(i_to_iposneg)
    @check T <: Integer
    V = typeof(i_to_iposneg)
    Vp = typeof(ipos_to_i)
    Vn = typeof(ineg_to_i)
    new{T,V,Vp,Vn}(i_to_iposneg,ipos_to_i,ineg_to_i)
  end
end

function PosNegPartition(ipos_to_i::AbstractArray,Ni::Integer)
  i_to_iposneg = zeros(typeof(Ni),Ni)
  i_to_iposneg[ipos_to_i] .= 1:length(ipos_to_i)
  Nneg = count(k->k==0,i_to_iposneg)
  ineg_to_i = similar(ipos_to_i,eltype(ipos_to_i),(Nneg,))
  ineg = 1
  for (i,iposneg) in enumerate(i_to_iposneg)
    if iposneg==0
      i_to_iposneg[i] = -ineg
      ineg_to_i[ineg] = i
      ineg += 1
    end
  end
  PosNegPartition(i_to_iposneg,ipos_to_i,ineg_to_i)
end

Base.IndexStyle(::Type{<:PosNegPartition}) = IndexLinear()
@propagate_inbounds Base.getindex(a::PosNegPartition,i::Integer) = a.i_to_iposneg[i]
Base.size(a::PosNegPartition) = size(a.i_to_iposneg)
get_array(a::PosNegPartition) = a.i_to_iposneg

pos_and_neg_length(a::PosNegPartition) = (length(a.ipos_to_i),length(a.ineg_to_i))
is_exhaustive(a::PosNegPartition) = true
pos_and_neg_indices(a::PosNegPartition) = (a.ipos_to_i,a.ineg_to_i)
function aligned_with_pos(a::PosNegPartition,j_to_i,npos)
  @check length(a.ipos_to_i) == npos
  a.ipos_to_i === j_to_i || a.ipos_to_i == j_to_i
end
function aligned_with_neg(a::PosNegPartition,j_to_i,nneg)
  @check length(a.ineg_to_i) == nneg
  a.ineg_to_i === j_to_i || a.ineg_to_i == j_to_i
end
all_pos(a::PosNegPartition) = length(a.ineg_to_i)==0
all_neg(a::PosNegPartition) = length(a.ipos_to_i)==0

# @propagate_inbounds function Base.setindex!(a::LazyArray{<:Fill{<:PosNegReindex}},v,j::Integer)
#   k = a.map.value
#   i_to_v = k.values
#   j_to_i, = a.f
#   i = j_to_i[j]
#   i_to_v[i]=v
# end


text:
using Xpress
using MathOptInterface
using Random
using Test

const MOI = MathOptInterface

function callback_simple_model()
    model = Xpress.Optimizer(
        PRESOLVE = 0,
        CUTSTRATEGY = 0,
        HEURSTRATEGY = 0,
        SYMMETRY = 0,
        OUTPUTLOG = 0
    )

    MOI.Utilities.loadfromstring!(model, """
        variables: x, y
        maxobjective: y
        c1: x in Integer()
        c2: y in Integer()
        c3: x in Interval(0.0, 2.5)
        c4: y in Interval(0.0, 2.5)
    """)
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    return model, x, y
end

function callback_knapsack_model()
    model = Xpress.Optimizer(
        PRESOLVE = 0,
        CUTSTRATEGY = 0,
        HEURSTRATEGY = 0,
        SYMMETRY = 0,
        OUTPUTLOG = 0
    )
    MOI.set(model, MOI.NumberOfThreads(), 2)

    N = 30
    x = MOI.add_variables(model, N)
    MOI.add_constraints(model, x, MOI.ZeroOne())
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    Random.seed!(1)
    item_weights, item_values = rand(N), rand(N)
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0)
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_values, x), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    return model, x, item_weights
end

@testset "LazyConstraintCallback" begin
    @testset "LazyConstraint" begin
        model, x, y = callback_simple_model()
        global lazy_called = false
        MOI.set(model, MOI.LazyConstraintCallback(), cb_data -> begin
            global lazy_called = true
            x_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
            y_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), y)
            status = MOI.get(model, MOI.CallbackNodeStatus(cb_data))::MOI.CallbackNodeStatusCode
            if round.(Int, [x_val, y_val]) ≈ [x_val, y_val] atol=1e-6
                @test status == MOI.CALLBACK_NODE_STATUS_INTEGER
            else
                @test status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
            end
            @test MOI.supports(model, MOI.LazyConstraint(cb_data))
            if y_val - x_val > 1 + 1e-6
                MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                        0.0
                    ),
                    MOI.LessThan{Float64}(1.0)
                )
            elseif y_val + x_val > 3 + 1e-6
                MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                        0.0
                    ), MOI.LessThan{Float64}(3.0)
                )
            end
        end)
        @test MOI.supports(model, MOI.LazyConstraintCallback())
        MOI.optimize!(model)
        @test lazy_called
        @test MOI.get(model, MOI.VariablePrimal(), x) == 1
        @test MOI.get(model, MOI.VariablePrimal(), y) == 2
    end
    @testset "OptimizeInProgress" begin
        model, x, y = callback_simple_model()
        MOI.set(model, MOI.LazyConstraintCallback(), cb_data -> begin
            @test_throws(
                MOI.OptimizeInProgress(MOI.VariablePrimal()),
                MOI.get(model, MOI.VariablePrimal(), x)
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveValue()),
                MOI.get(model, MOI.ObjectiveValue())
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveBound()),
                MOI.get(model, MOI.ObjectiveBound())
            )
        end)
        MOI.optimize!(model)
    end

    @testset "HeuristicSolution" begin
        model, x, y = callback_simple_model()
        cb = nothing
        MOI.set(model, MOI.LazyConstraintCallback(), cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                [x],
                [2.0]
            )
        end)
        @test_throws(
            MOI.InvalidCallbackUsage(
                MOI.LazyConstraintCallback(),
                MOI.HeuristicSolution(cb)
            ),
            MOI.optimize!(model)
        )
   end
end

@testset "UserCutCallback" begin
    @testset "UserCut" begin
        model, x, item_weights = callback_knapsack_model()
        global user_cut_submitted = false
        MOI.set(model, MOI.UserCutCallback(), cb_data -> begin
            terms = MOI.ScalarAffineTerm{Float64}[]
            accumulated = 0.0
            for (i, xi) in enumerate(x)
                if MOI.get(model, MOI.CallbackVariablePrimal(cb_data), xi) > 0.0
                    push!(terms, MOI.ScalarAffineTerm(1.0, xi))
                    accumulated += item_weights[i]
                end
            end
            @test MOI.supports(model, MOI.UserCut(cb_data))
            if accumulated > 10.0
                MOI.submit(
                    model,
                    MOI.UserCut(cb_data),
                    MOI.ScalarAffineFunction{Float64}(terms, 0.0),
                    MOI.LessThan{Float64}(length(terms) - 1)
                )
                global user_cut_submitted = true
            end
        end)
        @test MOI.supports(model, MOI.UserCutCallback())
        MOI.optimize!(model)
        @test user_cut_submitted
    end
    @testset "HeuristicSolution" begin
        model, x, item_weights = callback_knapsack_model()
        cb = nothing
        MOI.set(model, MOI.UserCutCallback(), cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                [x[1]],
                [0.0]
            )
        end)
        @test_throws(
            MOI.InvalidCallbackUsage(
                MOI.UserCutCallback(),
                MOI.HeuristicSolution(cb)
            ),
            MOI.optimize!(model)
        )
    end
end

@testset "HeuristicCallback" begin
    @testset "HeuristicSolution" begin
        model, x, item_weights = callback_knapsack_model()
        global callback_called = false
        MOI.set(model, MOI.HeuristicCallback(), cb_data -> begin
            x_vals = MOI.get.(model, MOI.CallbackVariablePrimal(cb_data), x)
            status = MOI.get(model, MOI.CallbackNodeStatus(cb_data))::MOI.CallbackNodeStatusCode
            if round.(Int, x_vals) ≈ x_vals atol=1e-6
                @test status == MOI.CALLBACK_NODE_STATUS_INTEGER
            else
                @test status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
            end
            @test MOI.supports(model, MOI.HeuristicSolution(cb_data))
            @test MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                floor.(x_vals)
            ) == MOI.HEURISTIC_SOLUTION_UNKNOWN
            global callback_called = true
        end)
        @test MOI.supports(model, MOI.HeuristicCallback())
        MOI.optimize!(model)
        @test callback_called
    end
    @testset "LazyConstraint" begin
        model, x, item_weights = callback_knapsack_model()
        cb = nothing
        MOI.set(model, MOI.HeuristicCallback(), cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0)
            )
        end)
        @test_throws(
            MOI.InvalidCallbackUsage(
                MOI.HeuristicCallback(),
                MOI.LazyConstraint(cb)
            ),
            MOI.optimize!(model)
        )
    end
    @testset "UserCut" begin
        model, x, item_weights = callback_knapsack_model()
        cb = nothing
        MOI.set(model, MOI.HeuristicCallback(), cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.UserCut(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0)
            )
        end)
        @test_throws(
            MOI.InvalidCallbackUsage(
                MOI.HeuristicCallback(),
                MOI.UserCut(cb)
            ),
            MOI.optimize!(model)
        )
    end
end

@testset "Xpress.CallbackFunction" begin
    @testset "OptimizeInProgress" begin
        model, x, y = callback_simple_model()
        MOI.set(model, Xpress.CallbackFunction(), (cb_data) -> begin
            @test_throws(
                MOI.OptimizeInProgress(MOI.VariablePrimal()),
                MOI.get(model, MOI.VariablePrimal(), x)
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveValue()),
                MOI.get(model, MOI.ObjectiveValue())
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveBound()),
                MOI.get(model, MOI.ObjectiveBound())
            )
        end)
        @test MOI.supports(model, Xpress.CallbackFunction())
        MOI.optimize!(model)
    end
    @testset "LazyConstraint" begin
        model, x, y = callback_simple_model()
        cb_calls = Int32[]
        global generic_lazy_called = false
        function callback_function(cb_data)
            push!(cb_calls, 1)
            Xpress.get_cb_solution(model, cb_data.model)
            x_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
            y_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), y)
            if y_val - x_val > 1 + 1e-6
                MOI.submit(model, MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                        0.0
                    ),
                    MOI.LessThan{Float64}(1.0)
                )
            elseif y_val + x_val > 3 + 1e-6
                MOI.submit(model, MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                        0.0
                    ),
                    MOI.LessThan{Float64}(3.0)
                )
            end
        end
        MOI.set(model, Xpress.CallbackFunction(), callback_function)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), x) == 1
        @test MOI.get(model, MOI.VariablePrimal(), y) == 2
        @test length(cb_calls) > 0
    end
    @testset "UserCut" begin
        model, x, item_weights = callback_knapsack_model()
        user_cut_submitted = false
        cb_calls = Int32[]
        MOI.set(model, Xpress.CallbackFunction(), (cb_data) -> begin
            push!(cb_calls)
            if Xpress.getintattrib(cb_data.model, Xpress.Lib.XPRS_CALLBACKCOUNT_OPTNODE) > 1
                return
            end
            Xpress.get_cb_solution(model, cb_data.model)
            terms = MOI.ScalarAffineTerm{Float64}[]
            accumulated = 0.0
            for (i, xi) in enumerate(x)
                if MOI.get(model, MOI.CallbackVariablePrimal(cb_data), xi) > 0.0
                    push!(terms, MOI.ScalarAffineTerm(1.0, xi))
                    accumulated += item_weights[i]
                end
            end
            if accumulated > 10.0
                MOI.submit(
                    model,
                    MOI.UserCut(cb_data),
                    MOI.ScalarAffineFunction{Float64}(terms, 0.0),
                    MOI.LessThan{Float64}(length(terms) - 1)
                )
                user_cut_submitted = true
            end
            return
        end)
        MOI.optimize!(model)
        @test user_cut_submitted
    end
    @testset "HeuristicSolution" begin
        model, x, item_weights = callback_knapsack_model()
        callback_called = false
        cb_calls = Int32[]
        MOI.set(model, Xpress.CallbackFunction(), (cb_data) -> begin
            if Xpress.getintattrib(cb_data.model, Xpress.Lib.XPRS_CALLBACKCOUNT_OPTNODE) > 1
                return
            end
            Xpress.get_cb_solution(model, cb_data.model)
            x_vals = MOI.get.(model, MOI.CallbackVariablePrimal(cb_data), x)
            @test MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                floor.(x_vals)
            ) == MOI.HEURISTIC_SOLUTION_UNKNOWN
            callback_called = true
            return
        end)
        MOI.optimize!(model)
        @test callback_called
    end
end

@testset "Xpress.CallbackFunction.CallbackNodeStatus" begin
    model, x, item_weights = callback_knapsack_model()
    global unknown_reached = false
    MOI.set(model, Xpress.CallbackFunction(), (cb_data) -> begin
        if MOI.get(model, MOI.CallbackNodeStatus(cb_data)) == MOI.CALLBACK_NODE_STATUS_UNKNOWN
            global unknown_reached = true
        end
    end)
    MOI.optimize!(model)
    @test unknown_reached
end


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

if !isdefined(@__MODULE__, :DefsAppId)
    const DefsAppId = String
else
    @warn("Skipping redefinition of DefsAppId to String")
end


text:
include("orbit/orbit.jl")
include("external_model/external_model.jl")
include("external_model/disturbance_torque.jl")
include("internal_model/magnetic_torque.jl")
#include("dynamic_model/dynamic_model.jl")
include("dynamics/dynamics.jl")
#include("satellite/satellite.jl")
include("satellite/attitude_control.jl")
include("satellite/attitude_determination.jl")
include("satellite/target_decision.jl")
include("satellite/control_algorithm.jl")
include("plot/plot_plots.jl")
#include("plot/plot_makie.jl")
include("coordinates.jl")

using SatelliteToolbox
using Dates
using LinearAlgebra

using HDF5

function main()
	#=
	設定パラメータ
	=#
	DataNum = 5000 #シミュレータ反復回数
	dt = 1 ##シミュレータの計算間隔 [s]
	start_time = DateTime(2021, 05, 23, 18, 00, 00)	#シミュレート開始時刻
	TLEFileName = "./orbit/ISS_TLE.txt"

	# 初期姿勢，角速度の設定
	sat_attqua_Initial = [cosd(15), sind(15)/sqrt(3), sind(15)/sqrt(3), sind(15)/sqrt(3)]
	sat_ω_Initial = [0., 0., 0.]

	# 撮影用パラメータの設定
	limit_time = DataNum                 # 計算開始時刻から "limit_time × dt" sec の間に撮影を行う
 	targetpos_geod = [-40., -110., 25.7]    # 撮影対象の位置@Geodetic
	cam_viewangle = 40                   # カメラの視野角
	sat_axisval = 80                     # 撮影を許可する角度範囲（直下向きを0°とし，それと撮影時カメラ方向との角度差に対する制限）
	cam_origindir = [0., 0., 1.]         # カメラの方向ベクトル@SCSF
	picture_aspectratio = [4, 3]         # 撮影画像の縦横比
	num = 1                              # ターゲット軌跡記憶用配列のどこまでデータが入ったかを記憶する

	# 制御用パラメータの設定
	kp = 0.00000030                       # クロスプロダクト則比例ゲイン
	kr = 0.0000030                         # クロスプロダクト則微分ゲイン
	mtq_maxcurrent = 0.002              # 磁気トルカの最大駆動電流
	mtq_scutter = 255                    # 磁気トルカの駆動電流分割数（" ± mtq_scutter" 段階で行う）
	Tmax = 1.0*10^(-7)                   # 出力トルクの最大値
	t_scatternum = 255                   # 出力トルクの分割数 (" ± t_scatternum" 段階で行う)
	I = [(0.1^2)/6 0.        0.;
		 0.        (0.1^2)/6 0.;
		 0.        0.        (0.1^2)/6]  # 衛星の慣性テンソル
	sat_size = [0.1 , 0.1 , 0.1]         # 衛星各辺長さ [x,y,z][m]
	cm = [0.005, 0.005, 0.005]           # 衛星体心から見た重心の位置ベクトル [x,y,z][m]
	Cd = 2.0                             # 抗力係数（宇宙空間では通常2~3）
	target_updatefreq = 12               # 目標姿勢の更新頻度 [step/回]
	target_updaterange = 120             # 目標姿勢の更新を行う時間範囲（"撮影時刻 ± target_updaterange" の間は目標姿勢の更新を行う）
	CP2Bdot_delay = -1                    # CP制御の後 nステップはBdot制御を行わない
	CP2Bdot_count = 0                    # CP制御からのカウント数

	# 誤差検証用パラメータ
	x_ecef_error = [0., 0., 0. ]
	targetqua_error = SatelliteToolbox.Quaternion(cos(deg2rad(0)), 0, sin(deg2rad(0)), 0)

	#=
	データ保存用変数
	=#
	# 衛星位置，進行方向
	dotvs = zeros(DataNum)                    # 進行方向とSCOFx軸とのずれ（内積）
	direct_on_SCOFs = zeros(DataNum,3)        # SCOF上での衛星の進行方向
	# 外環境
	mag_vecs = zeros(DataNum,3)               # 地磁場@SEOF
	sun_vecs = zeros(DataNum,3)               # 太陽方向@ECEF
	atoms_denses = zeros(DataNum)             # 大気密度
	# 姿勢
	sat_attqua_elements = zeros(DataNum+1, 4) # 衛星姿勢クォータニオンの4要素
	sat_ω = zeros(DataNum+1, 3)               # 衛星の各軸回り角速度@SCSF
	sat_eularangle = zeros(DataNum, 3)      # 衛星姿勢（オイラー角）
	# トルク
	disturbance = zeros(DataNum, 3)           # 擾乱トルクの合計値@SCSF
	airtorques = zeros(DataNum, 3)            # 空力トルク

	suntorques = zeros(DataNum, 3)            # 太陽輻射圧トルク
	magtorques = zeros(DataNum, 3)            # 磁気トルク
	# 制御値
	mtq_currentlog = zeros(DataNum, 3)        # 磁気トルカ入力電流
	M_reqs = zeros(DataNum, 3)                # 所望磁気モーメント
	T_reqs = zeros(DataNum, 3)                # 所望トルク
	Terrors = zeros(DataNum, 3)               # 所望トルクの内，出せないトルク
	controlmethod = zeros(DataNum)            # どの制御を行っているかの判定
	# 撮影
	cam_dir = zeros(DataNum, 3)               # カメラ方向ベクトル@SEOF
	tar_dir = zeros(DataNum, 3)               # 目標方向ベクトル@SEOF
	target_diffs = zeros(DataNum)             # 目標方向ベクトルとカメラ方向ベクトルとの角度誤差
	tarpos_log = zeros(target_updaterange*2, 2)  # 撮影画像上でのターゲット位置の軌跡

	# 初期値の代入
	sat_attqua_elements[1,:] = sat_attqua_Initial
	sat_ω[1, :] = sat_ω_Initial

	# 軌道・太陽方向計算については先に行い、全時間分を配列に保存する
	JD_log, x_ecef_log, v_ecef_log, x_geod_log, tles, eop_IAU2000A = orbit_cal(DataNum,dt,start_time,TLEFileName)
	sunvector_index = sunvector_model(JD_log[1], tles)

	h5open("orbit.h5", "w") do file
		write(file, "JD_log", JD_log)  # alternatively, say "@write file A"
		write(file, "x_ecef_log", x_ecef_log)
		write(file, "v_ecef_log", v_ecef_log)
		write(file, "x_geod_log", x_geod_log)
		write(file, "sunvector_index", sunvector_index)
	end

	# 初期状態の衛星位置@ECI
	x_eci_1st = rECEFtoECI(ITRF(), GCRF(), JD_log[1], eop_IAU2000A)*x_ecef_log[1,:]

	# 目標姿勢の計算
	targetpos_ecef = GeodetictoECEF(deg2rad(targetpos_geod[1]), deg2rad(targetpos_geod[2]), targetpos_geod[3])
	satqua = SatelliteToolbox.Quaternion(sat_attqua_elements[1,1], sat_attqua_elements[1,2], sat_attqua_elements[1,3], sat_attqua_elements[1,4])
	shoot_time = shootingtime_decision2(x_ecef_log[1:DataNum, :], targetpos_ecef, 1, limit_time, cam_viewangle, sat_axisval)
	println("shoottime", shoot_time)
	targetqua = targetqua_dicision(targetpos_ecef, x_ecef_log[shoot_time-target_updaterange, :], v_ecef_log[shoot_time-target_updaterange, :], sat_axisval)
	println("targetqua:", targetqua)
	# targetqua = SatelliteToolbox.Quaternion(cosd(0), sind(0)/sqrt(3), sind(0)/sqrt(3), sind(0)/sqrt(3))

	for i=1:DataNum

		# 時刻表示
		current_time = start_time + Second(dt) * i
		println("")
		println("Current DateTime:",current_time)
		println("       loop time:", i)
		println("              JD:",JD_log[i])
		println("      x_ecef_log:",x_ecef_log[i,:])
		println("      v_ecef_log:",v_ecef_log[i,:])
		println("      x_geod_log:",x_geod_log[i,:])
		dotv = dot(x_ecef_log[i,:] , v_ecef_log[i,:])
		dotvs[i] = dotv
		println("            dotv:",dotv)

		# 衛星外環境モデル計算
		mag_vec, sun_vec, atoms_dens = external_model(current_time,x_ecef_log[i,:],x_geod_log[i,:], x_eci_1st, eop_IAU2000A, sunvector_index)
		mag_vec_seof = ecef_to_DCM(x_ecef_log[i,:],v_ecef_log[i,:],true)*mag_vec
		mag_vecs[i,:] = mag_vec_seof
		sun_vecs[i,:] = sun_vec
		atoms_denses[i] = atoms_dens
		println("         mag_vec:",mag_vec)
		println("         sun_vec:",sun_vec)
		println("      atoms_dens:",atoms_dens)

		direct_on_SCOFs[i,:] = ecef_to_DCM(x_ecef_log[i,:],v_ecef_log[i,:],true) * v_ecef_log[i,:]

		# 衛星姿勢
		qua = SatelliteToolbox.Quaternion(sat_attqua_elements[i,1], sat_attqua_elements[i,2], sat_attqua_elements[i,3], sat_attqua_elements[i,4])

		#擾乱の計算
		v_scof = [norm(v_ecef_log[i,:]), 0., 0.]
		Ts = sun_pressure(sun_vecs[i,:], qua, sat_size, cm)
		Ta = air_pressure(atoms_denses[i], v_scof, qua, sat_size, cm, Cd)
		suntorques[i,:] = Ts
		airtorques[i,:] = Ta
		disturbance[i,:] = Ts + Ta

		# 地磁場ベクトル @SCSF の計算
		magvec_scsfqua = qua \ mag_vec_seof * qua
		magvec_scsf = [magvec_scsfqua.q1*10^(-9), magvec_scsfqua.q2*10^(-9), magvec_scsfqua.q3*10^(-9)]
		println("  scsf_mag_vec:", magvec_scsf)

		# 姿勢決定による衛星姿勢情報の取得
		qua_ad, vel_ad, omega_ad = attitude_empty(qua, v_ecef_log[i, :], sat_ω[i, :])

		#=
		# B-dot法による回転抑制制御
		if i==1
			M = B_dot(magvec_scsf, sat_ω[i, :], sat_ω[i, :])
		else
			M = B_dot(magvec_scsf, sat_ω[i, :], sat_ω[i-1, :])
		end
		=#

		# クロスプロダクト則による指向制御

		tarqua = targetqua_error * targetqua
		Treq, M = cross_product(tarqua,qua_ad, kp, kr, omega_ad, magvec_scsf)
		# Treq, M = pseudo_inverse(tarqua, qua_ad,kp, kr, omega_ad, magvec_scsf)
		# Treq, M, n = crossproduct_adj(targetqua, qua_ad, kp, kr, omega_ad, magvec_scsf)
		# controlmethod[i] = n

		#=
		# 所望トルクと地磁場のなす角が80°~100°の時はCrossProduct，それ以外はB-dot制御を行う
		tarqua = targetqua_error * targetqua
		Treq, M, n = cp_and_bdot(tarqua, qua_ad, omega_ad, kp, kr, magvec_scsf, boundangle)
		Terrors[i, :] = dot(Treq, magvec_scsf)/(norm(magvec_scsf)^2) * magvec_scsf
		=#

		M_reqs[i, :] = M
		println("request_Moment:", M)

		i_m = mm2current_theory(M)
		# i_m = digitizing_3elesvec(i_m, mtq_maxcurrent, mtq_scutter, true)
		Tm = magnetic_torque(i_m, magvec_scsf)
		mtq_currentlog[i, :] = i_m
		magtorques[i,:] = Tm
		println("  current_mag:",i_m)

		# Treqの離散化
		# Treq = digitizing_3elesvec(Treq, Tmax, t_scatternum, true)
		T_reqs[i, :] = Treq

		# Tm = collect(inv(ecef_to_DCM(x_ecef_log[i,:],v_ecef_log[i,:],true)) * Tm)

		# ダイナミクス
		# next_qua, ω = dynamics(qua, sat_ω[i,:], Ta+Ts+Tm, I, dt)
		# next_qua, ω = dynamics(qua, sat_ω[i,:], Tm, I, dt)
		# next_qua, ω = dynamics(qua, sat_ω[i,:], Treq, I, dt)
		next_qua, ω = dynamics(qua, sat_ω[i,:], Treq+Ta+Ts, I, dt)

		sat_attqua_elements[i+1, :] = [next_qua.q0, next_qua.q1, next_qua.q2, next_qua.q3]
		sat_ω[i+1, :] = ω
		println("quaternion:", qua)
		println("onega : ", ω)

		# カメラ方向
		cam_dir[i, :], target_diffs[i] = anglediffs(qua, targetqua, cam_origindir)

		tar_dir_q = targetqua * cam_origindir / targetqua
		tar_dir[i, :] = [tar_dir_q.q1, tar_dir_q.q2, tar_dir_q.q3]

		# オイラー角
		roll, pitch, yaw = Euler_Angles(qua \ targetqua)
		sat_eularangle[i, 1] = roll
		sat_eularangle[i, 2] = pitch
		sat_eularangle[i, 3] = yaw

		# 撮影地点の画像上の軌跡計算, 目標姿勢の更新
		if norm(i - shoot_time) < target_updaterange

			tarpos_log[num, :] = cal_targetlocus(targetpos_ecef, x_ecef_log, x_geod_log, v_ecef_log, x_ecef_error, qua, i)
			num = num + 1

			if norm(i-shoot_time) % target_updatefreq == 0
				# targetqua = targetqua_dicision(targetpos_ecef, x_ecef_log[i+target_updatefreq, :], v_ecef_log[i+target_updatefreq, :], sat_axisval)
			end

		end

	end

	plot = true

	if plot == true
		plot_2scalar([1:DataNum],dotvs,"dotvs")
		plot_vec(direct_on_SCOFs,"direct_on_SCOFs")

		plot_2scalar([1:DataNum],atoms_denses,"atoms_dens")
		plot_2scalar([1:DataNum],x_geod_log[:,3],"height")
		plot_vec_range(mag_vecs,[-50000, 50000], [-50000, 50000], [-50000, 50000],"mag_vec")
		plot_3scalar([1:DataNum], mag_vecs[1:DataNum, 1], mag_vecs[1:DataNum, 2], mag_vecs[1:DataNum, 3], ["x", "y", "z"], "magvec_elements")
		plot_vec(sun_vecs,"sun_vec")
		plot_3scalar([1:DataNum], sun_vecs[1:DataNum, 1], sun_vecs[1:DataNum, 2], sun_vecs[1:DataNum, 3], ["x", "y", "z"], "sunvec_elements")
		plot_vec(x_ecef_log,"x_ecef_log")
		plot_vec(v_ecef_log,"v_ecef_log")
		plot_vec(x_geod_log,"x_geod_log")
		plot_2scalar(rad2deg.(x_geod_log[:,2]),rad2deg.(x_geod_log[:,1]),"x_geod_log_2d")
		plot_2scalar(rad2deg.(x_geod_log[shoot_time-target_updaterange:shoot_time+target_updaterange,2]),rad2deg.(x_geod_log[shoot_time-target_updaterange:shoot_time+target_updaterange,1]),"x_geod_log_2d_2")
		plot_vec(disturbance, "disturbance")
		plot_3scalar([1:DataNum], sat_ω[1:DataNum, 1], sat_ω[1:DataNum, 2], sat_ω[1:DataNum, 3], ["x", "y", "z"], "omega")
		plot_3scalar([1:DataNum], airtorques[1:DataNum, 1], suntorques[1:DataNum, 1], magtorques[1:DataNum, 1], ["air", "sun", "mag"], "x_torques")
		plot_3scalar([1:DataNum], airtorques[1:DataNum, 2], suntorques[1:DataNum, 2], magtorques[1:DataNum, 2], ["air", "sun", "mag"], "y_torques")
		plot_3scalar([1:DataNum], airtorques[1:DataNum, 3], suntorques[1:DataNum, 3], magtorques[1:DataNum, 3], ["air", "sun", "mag"], "z_torques")
		plot_3scalar([1:DataNum], T_reqs[1:DataNum, 1], T_reqs[1:DataNum, 2], T_reqs[1:DataNum, 3], ["x", "y", "z"], "request_torques")
		plot_vec(T_reqs, "request_torques_vec")
		plot_3scalar([1:DataNum], mtq_currentlog[:, 1], mtq_currentlog[:, 2], mtq_currentlog[:, 3], ["x", "y", "z"], "MTQcurrent")
		plot_3scalar([1:DataNum], sat_attqua_elements[1:DataNum, 2], sat_attqua_elements[1:DataNum, 3], sat_attqua_elements[1:DataNum, 4], ["q1", "q2", "q3"], "attqua_1")
		plot_2scalar([1:DataNum], sat_attqua_elements[1:DataNum, 1], "attqua_2")
		plot_vec(cam_dir, "cam_dir")
		plot_vec(cam_dir[ Int(floor(DataNum*0.9)) : DataNum, : ], "cam_dir2")
		plot_vecs(cam_dir, tar_dir, ["cam", "target"], "cam&tar_dir")
		plot_2scalar([1:DataNum], target_diffs, "target_diffs")
		plot_2scalar([shoot_time-target_updaterange:shoot_time+target_updaterange], target_diffs[shoot_time-target_updaterange:shoot_time+target_updaterange], "target_diffs2")
		plot_2scalar(tarpos_log[:, 1], tarpos_log[:, 2], "tarpos")
		plot_2scalar([1:DataNum], controlmethod[:], "controlmethod")
		plot_3scalar([1:DataNum], Terrors[:, 1], T_reqs[1:DataNum, 1], magtorques[1:DataNum, 1], ["-", "reqs", "out"], "x_torques_check")
		plot_3scalar([1:DataNum], Terrors[:, 2], T_reqs[1:DataNum, 2], magtorques[1:DataNum, 2], ["-", "reqs", "out"], "y_torques_check")
		plot_3scalar([1:DataNum], Terrors[:, 3], T_reqs[1:DataNum, 3], magtorques[1:DataNum, 3], ["-", "reqs", "out"], "z_torques_check")

		plot_3scalar([1:DataNum], target_diffs, controlmethod*180, zeros(DataNum), ["targetdiffs", "controlmethod", "-"], "diffs_method_check")

		plot_3scalar([1:DataNum], sat_eularangle[:, 1], sat_eularangle[:, 2], sat_eularangle[:, 3], ["roll", "pitch", "yaw"], "EularAngleDiffs")

		# 撮影限界	画像サイズ4:3として範囲換算
		picture_radius = x_geod_log[shoot_time, 3] * tand(cam_viewangle/2)
		picture_xratio = picture_aspectratio[1] / norm(picture_aspectratio)
		picture_yratio = picture_aspectratio[2] / norm(picture_aspectratio)
		picture_xmax = picture_xratio * picture_radius
		picture_ymax = picture_yratio * picture_radius
		plot_2scalar_range(tarpos_log[:, 1], tarpos_log[:, 2], [-1*picture_xmax, picture_xmax], [-1*picture_ymax, picture_ymax], "targetlocus")

		# plot_2scalar(targetlocus_picture[:, 1], targetlocus_picture[:, 2], "targetlocus")

	end

end

main()

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
# MIP01 from CPLEX.jl
function int1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # an example on mixed integer programming
    #
    #   maximize 1.1x + 2 y + 5 z
    #
    #   s.t.  x + y + z <= 10
    #         x + 2 y + z <= 15
    #
    #         x is continuous: 0 <= x <= 5
    #         y is integer: 0 <= y <= 10
    #         z is binary

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}),(MOI.SingleVariable,MOI.GreaterThan{Float64}), (MOI.SingleVariable, MOI.ZeroOne), (MOI.SingleVariable, MOI.Integer)])

    MOI.empty!(model)
    @test MOI.isempty(model)

    MOI.canaddvariable(model)
    v = MOI.addvariables!(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3

    cf = MOI.ScalarAffineFunction(v, [1.0,1.0,1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf), MOI.LessThan{Float64})
    c = MOI.addconstraint!(model, cf, MOI.LessThan(10.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    cf2 = MOI.ScalarAffineFunction(v, [1.0,2.0,1.0], 0.0)
    @test MOI.canaddconstraint(model, typeof(cf2), MOI.LessThan{Float64})
    c2 = MOI.addconstraint!(model, cf2, MOI.LessThan(15.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 2

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.Interval{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(v[1]), MOI.Interval(0.0, 5.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Interval{Float64}}()) == 1

    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.Interval{Float64})
    MOI.addconstraint!(model, MOI.SingleVariable(v[2]), MOI.Interval(0.0, 10.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Interval{Float64}}()) == 2
    @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(v[2])), MOI.Integer)
    MOI.addconstraint!(model, MOI.SingleVariable(v[2]), MOI.Integer())
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Integer}()) == 1

    @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(v[3])), MOI.ZeroOne)
    MOI.addconstraint!(model, MOI.SingleVariable(v[3]), MOI.ZeroOne())
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.ZeroOne}()) == 1

    objf = MOI.ScalarAffineFunction(v, [1.1, 2.0, 5.0], 0.0)
    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.ResultCount())
        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FeasiblePoint, MOI.NearlyFeasiblePoint ]

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 19.4 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4,5,1] atol=atol rtol=rtol

        @test MOI.canget(model, MOI.ConstraintPrimal(), typeof(c))
        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 10 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.ConstraintPrimal(), typeof(c2))
        @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 15 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.DualStatus()) == false

        if MOI.canget(model, MOI.ObjectiveBound())
            @test MOI.get(model, MOI.ObjectiveBound()) >= 19.4
        end
        if MOI.canget(model, MOI.RelativeGap())
            @test MOI.get(model, MOI.RelativeGap()) >= 0.0
        end
        if MOI.canget(model, MOI.SolveTime())
            @test MOI.get(model, MOI.SolveTime()) >= 0.0
        end
        if MOI.canget(model, MOI.SimplexIterations())
            @test MOI.get(model, MOI.SimplexIterations()) >= 0
        end
        if MOI.canget(model, MOI.BarrierIterations())
            @test MOI.get(model, MOI.BarrierIterations()) >= 0
        end
        if MOI.canget(model, MOI.NodeCount())
            @test MOI.get(model, MOI.NodeCount()) >= 0
        end
    end
end

Base.isapprox(a::T, b::T; kwargs...) where T <: Union{MOI.SOS1, MOI.SOS2} = isapprox(a.weights, b.weights; kwargs...)

# sos from CPLEX.jl" begin
function int2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [ (MOI.VectorOfVariables, MOI.SOS1),
    #                                                                    (MOI.VectorOfVariables, MOI.SOS2) ])
    @testset "SOSI" begin
        #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.VectorOfVariables,MOI.SOS1), (MOI.SingleVariable,MOI.LessThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        v = MOI.addvariables!(model, 3)
        @test MOI.get(model, MOI.NumberOfVariables()) == 3
        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(v[1])), MOI.LessThan{Float64})
        MOI.addconstraint!(model, MOI.SingleVariable(v[1]), MOI.LessThan(1.0))
        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(v[2])), MOI.LessThan{Float64})
        MOI.addconstraint!(model, MOI.SingleVariable(v[2]), MOI.LessThan(1.0))
        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(v[3])), MOI.LessThan{Float64})
        MOI.addconstraint!(model, MOI.SingleVariable(v[3]), MOI.LessThan(2.0))

        @test MOI.canaddconstraint(model, MOI.VectorOfVariables, MOI.SOS1{Float64})
        c1 = MOI.addconstraint!(model, MOI.VectorOfVariables([v[1], v[2]]), MOI.SOS1([1.0, 2.0]))
        @test MOI.canaddconstraint(model, MOI.VectorOfVariables, MOI.SOS1{Float64})
        c2 = MOI.addconstraint!(model, MOI.VectorOfVariables([v[1], v[3]]), MOI.SOS1([1.0, 2.0]))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SOS1}()) == 2

        @test MOI.canget(model, MOI.ConstraintSet(), typeof(c2))
        @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c2))
        #=
            To allow for permutations in the sets and variable vectors
            we're going to sort according to the weights
        =#
        cs_sos = MOI.get(model, MOI.ConstraintSet(), c2)
        cf_sos = MOI.get(model, MOI.ConstraintFunction(), c2)
        p = sortperm(cs_sos.weights)
        @test cs_sos.weights[p] ≈ [1.0, 2.0] atol=atol rtol=rtol
        @test cf_sos.variables[p] == v[[1,3]]

        objf = MOI.ScalarAffineFunction(v, [2.0, 1.0, 1.0], 0.0)
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.ResultCount())
            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0,1,2] atol=atol rtol=rtol

            @test MOI.canget(model, MOI.DualStatus()) == false
        end

        @test MOI.candelete(model, c1)
        MOI.delete!(model, c1)
        @test MOI.candelete(model, c2)
        MOI.delete!(model, c2)

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.ResultCount())
            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 5 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1,1,2] atol=atol rtol=rtol
        end
    end
    @testset "SOSII" begin
        #@test MOI.supportsproblem(model,
        #    MOI.ScalarAffineFunction{Float64},
        #    [
        #        (MOI.VectorOfVariables,MOI.SOS1),
        #        (MOI.VectorOfVariables,MOI.SOS2),
        #        (MOI.SingleVariable, MOI.ZeroOne),
        #        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})
        #        ]
        #)

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        v = MOI.addvariables!(model, 10)
        @test MOI.get(model, MOI.NumberOfVariables()) == 10

        bin_constraints = []
        for i in 1:8
            @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.Interval{Float64})
            MOI.addconstraint!(model, MOI.SingleVariable(v[i]), MOI.Interval(0.0, 2.0))
            @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(v[i])), MOI.ZeroOne)
            push!(bin_constraints, MOI.addconstraint!(model, MOI.SingleVariable(v[i]), MOI.ZeroOne()))
        end

        MOI.addconstraint!(model,
            MOI.ScalarAffineFunction(v[[1,2,3,9]], [1.0,2.0,3.0,-1.0], 0.0),
            MOI.EqualTo(0.0)
        )

        MOI.addconstraint!(model,
            MOI.ScalarAffineFunction(v[[4,5,6,7,8,10]], [5.0,4.0,7.0,2.0,1.0,-1.0], 0.0),
            MOI.EqualTo(0.0)
        )

        MOI.addconstraint!(model,
            MOI.VectorOfVariables(v[[1, 2, 3]]),
            MOI.SOS1([1.0, 2.0, 3.0])
        )

        vv   = MOI.VectorOfVariables(v[[4,5,6,7,8]])
        sos2 = MOI.SOS2([5.0, 4.0, 7.0, 2.0, 1.0])
        @test MOI.canaddconstraint(model, typeof(vv), MOI.SOS2{Float64})
        c = MOI.addconstraint!(model, vv, sos2)

        @test MOI.canget(model, MOI.ConstraintSet(), typeof(c))
        @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c))
        #=
            To allow for permutations in the sets and variable vectors
            we're going to sort according to the weights
        =#
        cs_sos = MOI.get(model, MOI.ConstraintSet(), c)
        cf_sos = MOI.get(model, MOI.ConstraintFunction(), c)
        p = sortperm(cs_sos.weights)
        @test cs_sos.weights[p] ≈ [1.0, 2.0, 4.0, 5.0, 7.0] atol=atol rtol=rtol
        @test cf_sos.variables[p] == v[[8,7,5,4,6]]

        objf = MOI.ScalarAffineFunction([v[9], v[10]], [1.0, 1.0], 0.0)
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.ResultCount())
            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 15.0 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 3.0, 12.0] atol=atol rtol=rtol

            @test MOI.canget(model, MOI.DualStatus()) == false
        end

        for cref in bin_constraints
            @test MOI.candelete(model, cref)
            MOI.delete!(model, cref)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.ResultCount())
            @test MOI.get(model, MOI.ResultCount()) >= 1

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 30.0 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.0, 0.0, 2.0, 2.0, 0.0, 2.0, 0.0, 0.0, 6.0, 24.0] atol=atol rtol=rtol

            @test MOI.canget(model, MOI.DualStatus()) == false
        end
    end
end

# CPLEX #76
function int3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # integer knapsack problem
    # max   z - 0.5 ( b1 + b2 + b3) / 40
    # s.t.  0 <= z - 0.5 eᵀ b / 40 <= 0.999
    #       b1, b2, ... b10 ∈ {0, 1}
    #       z in {0, 1, 2, ..., 100}

    MOI.empty!(model)
    @test MOI.isempty(model)

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64},
    #    [
    #        (MOI.SingleVariable,MOI.ZeroOne),
    #        (MOI.SingleVariable,MOI.Integer),
    #        (MOI.SingleVariable,MOI.Interval{Float64}),
    #        (MOI.ScalarAffineFunction{Float64},MOI.Interval{Float64})
    #    ]
    #)

    MOI.canaddvariable(model)
    z = MOI.addvariable!(model)
    @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(z)), MOI.Integer)
    MOI.addconstraint!(model, MOI.SingleVariable(z), MOI.Integer())
    @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.Integer)
    MOI.addconstraint!(model, MOI.SingleVariable(z), MOI.Interval(0.0, 100.0))

    MOI.canaddvariable(model)
    b = MOI.addvariables!(model, 10)

    for bi in b
        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(bi)), MOI.ZeroOne)
        MOI.addconstraint!(model, MOI.SingleVariable(bi), MOI.ZeroOne())
    end

    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
    c = MOI.addconstraint!(model, MOI.ScalarAffineFunction(vcat(z, b), vcat(1.0, fill(-0.5 / 40, 10)), 0.0), MOI.Interval(0.0, 0.999))

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(vcat(z, b[1:3]), vcat(1.0, fill(-0.5 / 40, 3)), 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 atol=atol rtol=rtol

        # test for CPLEX.jl #76
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1 atol=atol rtol=rtol
    end
end

# Mixed-integer linear problems

function knapsacktest(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # integer knapsack problem
    # max 5a + 3b + 2c + 7d + 4e
    # st  2a + 8b + 4c + 2d + 5e <= 10
    #                  a,b,c,d,e ∈ binary

    MOI.empty!(model)
    @test MOI.isempty(model)

    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.SingleVariable,MOI.ZeroOne),(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64})])

    MOI.canaddvariable(model)
    v = MOI.addvariables!(model, 5)
    @test MOI.get(model, MOI.NumberOfVariables()) == 5

    for vi in v
        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(vi)), MOI.ZeroOne)
        MOI.addconstraint!(model, MOI.SingleVariable(vi), MOI.ZeroOne())
    end
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.ZeroOne}()) == 5
    @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
    c = MOI.addconstraint!(model, MOI.ScalarAffineFunction(v, [2.0, 8.0, 4.0, 2.0, 5.0], 0.0), MOI.LessThan(10.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}()) == 1

    @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(v, [5.0, 3.0, 2.0, 7.0, 4.0], 0.0))
    @test MOI.canset(model, MOI.ObjectiveSense())
    MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

    if MOI.canset(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
        MOI.set!(model, MOI.VariablePrimalStart(), v, [0.0, 0.0, 0.0, 0.0, 0.0])
    end

    if config.solve
        MOI.optimize!(model)

        @test MOI.canget(model, MOI.TerminationStatus())
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

        @test MOI.canget(model, MOI.PrimalStatus())
        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FeasiblePoint, MOI.NearlyFeasiblePoint ]

        @test MOI.canget(model, MOI.ObjectiveValue())
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 16 atol=atol rtol=rtol

        @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [1, 0, 0, 1, 1] atol=atol rtol=rtol
    end
end

const intlineartests = Dict("knapsack" => knapsacktest,
                            "int1"     => int1test,
                            "int2"     => int2test,
                            "int3"     => int3test)

@moitestset intlinear


text:
module Flatland

using LinearAlgebra

abstract type Shape end

abstract type AbstractRectangle <: Shape end

mutable struct Rectangle <: AbstractRectangle
    x::Float64
    y::Float64
    l::Float64
    w::Float64
    function Rectangle((x, y); l=1.0, w=1.0)
        return new(x, y, l, w)
    end
end

function Rectangle((xmin, xmax), (ymin, ymax))
    @assert xmin < xmax && ymin < ymax "Corners have to be ordered: `xmin < xmax && ymin < ymax `"
    x = (xmin + xmax) / 2
    y = (ymin + ymax) / 2
    l = xmax - xmin
    w = ymax - ymin
    return Rectangle((x, y), l=l, w=w)
end

mutable struct Square <: AbstractRectangle
    x::Float64
    y::Float64
    l::Float64
    function Square((x, y); l=1.0)
        return new(x, y, l)
    end
end

lw(shape::Rectangle) = shape.l, shape.w
lw(shape::Square) = shape.l, shape.l

mutable struct RegularPolygon{N} <: Shape
    x::Float64
    y::Float64
    R::Float64
    θ::Float64  # angle
    function RegularPolygon((x, y), n::Int; R=1.0, θ=0.0)
        @assert n ≥ 3 "polygons need a minimum of three corners"
        return new{n}(x, y, R, θ)
    end
end

#=
`Circle`s are pretty straightforward, having a center and a radius.
=#

mutable struct Circle <: Shape
    x::Float64
    y::Float64
    R::Float64
    function Circle((x, y); R=1.0)
        return new(x, y, R)
    end
end

abstract type AbstractTriangle <: Shape end

mutable struct Triangle <: AbstractTriangle
    x1::Float64
    x2::Float64
    x3::Float64
    y1::Float64
    y2::Float64
    y3::Float64
    Triangle((x1, y1), (x2, y2), (x3, y3)) = new(x1, x2, x3, y1, y2, y3)
end

ncorners(::AbstractRectangle) = 4
ncorners(::AbstractTriangle) = 3
ncorners(::Circle) = 0
ncorners(::RegularPolygon{N}) where {N} = N

function corners(shape::AbstractRectangle)
    x, y = center(shape)
    l, w = lw(shape)
    return [(x+l/2, y+w/2), (x-l/2, y+w/2), (x-l/2, y-w/2), (x+l/2, y-w/2)]
end

function corners(shape::RegularPolygon)
    x, y = center(shape)
    θ = shape.θ
    R = shape.R
    n = ncorners(shape)
    return [(x+R*cos(t+θ), y+R*sin(t+θ)) for t in range(0, step=2π/n, length=n)]
end

function corners(shape::Triangle)
    return [(shape.x1, shape.y1), (shape.x2, shape.y2), (shape.x3, shape.y3)]
end

center(shape::Shape) = shape.x, shape.y

function center(shape::Triangle)
    (x1, y1), (x2, y2), (x3, y3) = corners(shape)
    return ((x1 + x2 + x3) / 3, (y1 + y2 + y3) / 3)
end

# xycoords returns two of the outline vectors: ycoords and ycoords, for `Circle`, you can specify the number of points to take
# (50 by default).

xycoords(s::Shape) = [first(p) for p in corners(s)], [last(p) for p in corners(s)]

function xycoords(shape::Circle; n=50)
    ts = range(0, 2π, length=n)
    return shape.x .+ shape.R * cos.(ts), shape.y .+ shape.R * sin.(ts)
end

# x,y-bounding

#=
The fuctions below yield the outer limits of the x and y axes of your shape. Can you complete the method as a oneliner?

Hint: The function `extrema` could be useful here...
=#

xlim(shape::Shape) = extrema(xycoords(shape)[1])
xlim(shape::Circle) = (shape.x - shape.R, shape.x + shape.R)
xlim(shape::AbstractRectangle) = (shape.x - 0.5shape.l, shape.x + 0.5shape.l)

ylim(shape::Shape) = extrema(xycoords(shape)[2])
ylim(shape::Circle) = (shape.y - shape.R, shape.y + shape.R)
ylim(shape::AbstractRectangle) = (shape.y - 0.5lw(shape)[2], shape.y + 0.5lw(shape)[2])

boundingbox(shape::Shape) = Rectangle(xlim(shape), ylim(shape))

# # Area

#=
Next, we compute the area of our shapes.
=#

area(shape::AbstractRectangle) = prod(lw(shape))

function area(shape::RegularPolygon)
    R = shape.R
    n = ncorners(shape)
    dθ = 2π / n  # angle of a triangle
    h = R * cos(dθ/2)
    b = R * sin(dθ/2)
    return b * h * n
end

function area(shape::Triangle)
    (x1, y1), (x2, y2), (x3, y3) = corners(shape)
    return abs(x1*y2 + x2*y3 + x3*y1 - y1*x2 - y2*x3 - y3*x1) / 2
end

area(shape::Circle) = shape.R^2 * π

# move, rotate and scale

#=
Moving, rotating and scaling should also work.
Important, the functions work in-place, meaning that the modify your structure (that is why use use `mutable` structures).

For `Circle` and `AbstractRectangle` types, `rotate!` leaves them unchanged.

Hint: Rotations are in radials, so between $0$ and $2\pi$.
=#

function move!(shape::Shape, (dx, dy))
    shape.x += dx
    shape.y += dy
end

function move!(shape::Triangle, (dx, dy))
    shape.x1 += dx
    shape.x2 += dx
    shape.x3 += dx
    shape.y1 += dy
    shape.y2 += dy
    shape.y3 += dy
    shape
end

rotate!(shape::Union{Circle,AbstractRectangle}, dθ) = shape

function rotate!(shape::RegularPolygon, dθ)
    shape.θ += dθ
    shape
end

function rotate!(shape::Triangle, dθ)
    xc, yc = center(shape)
    # center triangle
    move!(shape, (-xc, -yc))
    (x1, y1), (x2, y2), (x3, y3) = corners(shape)
    cosdθ = cos(dθ)
    sindθ = sin(dθ)
    shape.x1 = x1 * cos(dθ) - y1 * sin(dθ)
    shape.x2 = x2 * cos(dθ) - y2 * sin(dθ)
    shape.x3 = x3 * cos(dθ) - y3 * sin(dθ)
    shape.y1 = y1 * cos(dθ) + x1 * sin(dθ)
    shape.y2 = y2 * cos(dθ) + x2 * sin(dθ)
    shape.y3 = y3 * cos(dθ) + x3 * sin(dθ)
    # set to original position
    move!(shape, (xc, yc))
    shape
end

function scale!(shape::Union{Circle,RegularPolygon}, a)
    @assert a > 0 "scaling has to be a positive number"
    shape.R *= a
end

function scale!(shape::Rectangle, a)
    @assert a > 0 "scaling has to be a positive number"
    shape.w *= a
    shape.l *= a
end

function scale!(shape::Square, a)
    @assert a > 0 "scaling has to be a positive number"
    shape.w *= a
end

function scale!(shape::Triangle, a)
    @assert a > 0 "scaling has to be a positive number"
    xc, yc = center(shape)
    move!(shape, (-xc, -yc))
    shape.x1 *= a
    shape.x2 *= a
    shape.x3 *= a
    shape.y1 *= a
    shape.y2 *= a
    shape.y3 *= a
    move!(shape, (xc, yc))
    return shape
end

# plotting utilities

using RecipesBase

#=
OK, let's take a look at our shapes! We use `RecipesBase` to allow plotting.
This falls back on `xycoords` (can you see how it works?), so make sure this method is operational.
=#

@recipe function f(s::Shape)
    xguide --> "x"
    yguide --> "y"
    label --> ""
    aspect_ratio := :equal
    seriestype := :shape
    x, y = xycoords(s)
    return x, y
end

Base.in((x, y), s::Circle) = (s.x - x)^2 + (s.y - y)^2 ≤ s.R^2

function Base.in((x, y), shape::AbstractRectangle)
    xc, yc = center(shape)
    l, w = lw(shape)
    return (xc - 0.5l ≤ x ≤ xc + 0.5l) && (yc - 0.5w ≤ y ≤ yc + 0.5w)
end

crossprod((x1, y1), (x2, y2)) = x1 * y2 - x2 * y1

"""
    same_side((a, b), p, q)

Given a line described by two points, `a` and `b`, check whether two points
`p` and `q` are on the same side.
"""
function same_side((a, b), p, q)
    # normal vector on the line
    n = (a[2] - b[2], b[1] - a[1])
    # check if they are on both sides by projection
	return sign(n ⋅ (p .- a)) == sign(n ⋅ (q .-a ))
end

"""
    linecross((p1, p2), (q1, q2))

Check whether line segments `(p1, p2)` and `(q1, q2)` intersect.
"""
function linecross((p1, p2), (q1, q2))
    v = p2 .- p1
    w = q2 .- q1
    vw = crossprod(v, w)
    t = crossprod(q1 .- p1, w) / vw
    s = crossprod(q1 .- p1, v) / vw
    return 0.0 ≤ t ≤ 1.0 &&  0.0 ≤ s ≤ 1.0
end

function Base.in(q, s::Triangle)
    p1, p2, p3 = corners(s)
    return same_side((p1, p2), p3, q) &&
            same_side((p2, p3), p1, q) &&
            same_side((p3, p1), p2, q)
end

function Base.in(q, shape::Shape)
    corns = corners(shape)
    n = ncorners(shape)
    c = center(shape)
    # check if q is always on the same side as the center
    for i in 1:n-1
        !same_side((corns[i], corns[i+1]), c, q) && return false
    end
    !same_side((corns[end], corns[1]), c, q) && return false
    return true
end

@inline function boundboxes_overlap(shape1::Shape, shape2::Shape)
    (xmin1, xmax1), (xmin2, xmax2) = xlim(shape1), xlim(shape2)
    (ymin1, ymax1), (ymin2, ymax2)  = ylim(shape1), ylim(shape2)
    # check for x and y overlap
    return (xmin1 ≤ xmin2 ≤ xmax1 || xmin1 ≤ xmax2 ≤ xmax1 || xmin2 ≤ xmin1 ≤ xmax2) &&
            (ymin1 ≤ ymin2 ≤ ymax1 || ymin1 ≤ ymax2 ≤ ymax1 ||ymin2 ≤ ymin1 ≤ ymax2)
end

Base.intersect(shape1::AbstractRectangle, shape2::AbstractRectangle) = boundboxes_overlap(shape1, shape2)

function Base.intersect(shape1::T, shape2::T) where {T<:Shape}
    (shape1.x - shape2.x)^2 + (shape1.y - shape2.y)^2 > (shape1.R + shape2.R)^2 && return false
    return center(shape1) ∈ shape2 ||
            center(shape2) ∈ shape1 ||
            any(c->c ∈ shape2, corners(shape1)) ||
            any(c->c ∈ shape1, corners(shape2))
end

function Base.intersect(shape1::Triangle, shape2::Triangle)
    # first check if the bounding boxes overlap
    boundboxes_overlap(shape1, shape2) || return false
    # yes?
    # now check if one shape is within the other one
    p1, p2, p3 = corners(shape1)
    q1, q2, q3 = corners(shape2)
    (p1 ∈ shape2 || q1 ∈ shape1) && return true
    # if not, we have to check whether two lines intersect
    return linecross((p1, p2), (q1, q2)) ||
            linecross((p1, p2), (q2, q3)) ||
            linecross((p1, p2), (q1, q3)) ||
            linecross((p3, p2), (q1, q2)) ||
            linecross((p3, p2), (q2, q3)) ||
            linecross((p3, p2), (q1, q3))
end

function Base.intersect(shape1::Circle, shape2::Circle)
    c1 = center(shape1)
    c2 = center(shape2)
    d = sum(abs2, c1 .- c2) |> sqrt
    return d < shape1.R + shape2.R
end

# random generation

#=
Finally, `randplace!` takes a shape, rotates it randomly and moves it randomly within the bounds of the limits
`(xmin, xmax)` and `(ymin, ymax)`.
=#

function randplace!(shape::Shape, (xmin, xmax), (ymin, ymax); rotate=true)
    # random rotation
    rotate && rotate!(shape, 2π * rand())
    # random tranlation within bound
    dxmin, dxmax = (xmin, xmax) .- xlim(shape)
    dymin, dymax = (ymin, ymax) .- ylim(shape)
    dx = (dxmax - dxmin) * rand() + dxmin
    dy = (dymax - dymin) * rand() + dymin
    move!(shape, (dx, dy))
    return shape
end

function rejection_sampling!(shapes::Vector{<:Shape}, xlims, ylims)
    trials = 0
    n = length(shapes)
    while true
        trials += 1
        for (i, shape) in enumerate(shapes)
            randplace!(shape, xlims, ylims)
            # any intersection with previous shapes: start again
            overlap = false
            for j in 1:i-1
                if intersect(shape, shapes[j])
                    overlap = true
                    break
                end
            end
            overlap && break
            i==n && return shapes, trials
        end
    end
end

# function to place `n` copies of a given `shape`.

function rejection_sampling(shape, n, xlims, ylims)
    shapes = [deepcopy(shape) for i in 1:n]
    trials = rejection_sampling!(shapes, xlims, ylims)
    return shapes, trials
end

end


text:
# GeometryTypes doesn't define an Ellipsoid type yet, so we'll make one ourselves!
mutable struct HyperEllipsoid{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    radii::Vec{N, T}
end

origin(geometry::HyperEllipsoid{N, T}) where {N, T} = geometry.center
radii(geometry::HyperEllipsoid{N, T}) where {N, T} = geometry.radii
center(geometry::HyperEllipsoid) = origin(geometry)

center(geometry::HyperRectangle) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperCube) = minimum(geometry) + 0.5 * widths(geometry)
center(geometry::HyperSphere) = origin(geometry)
center(geometry::Cylinder) = origin(geometry) + geometry.extremity / 2

struct PointCloud{T, Point <: StaticArray{Tuple{3}, T}} <: AbstractGeometry{3, T}
    points::Vector{Point}
    channels::Dict{Symbol, Any}
end

function PointCloud(points::AbstractVector{Point},
                    channels::Dict=Dict{Symbol, Any}()) where {T, Point <: StaticArray{Tuple{3}, T}}
    PointCloud{T, Point}(points, channels)
end

function PointCloud(points::AbstractVector{V},
                    channels::Dict=Dict{Symbol, Any}()) where {T, V <: AbstractVector{T}}
    PointCloud{T, Point{3, T}}(points, channels)
end

struct Triad <: AbstractGeometry{3, Float64}
    scale::Float64
    tube::Bool

    Triad(scale=1.0, tube=true) = new(scale, tube)
end

struct ArrowHead
    radius::Float64
    length::Float64

    ArrowHead(radius=0.05, length=radius) = new(radius, length)
end

struct PolyLine{T, Point <: StaticArray{Tuple{3}, T}} <: AbstractGeometry{3, T}
    points::Vector{Point}
    radius::Float64
    closed::Bool
    start_head::Union{ArrowHead, Nothing}
    end_head::Union{ArrowHead, Nothing}
end

PolyLine(points::AbstractVector{Point}, radius, closed, start_head, end_head) where {T, Point <: StaticArray{Tuple{3}, T}} =
    PolyLine{T, Point}(points, radius, closed, start_head, end_head)

PolyLine(points::AbstractVector{V}, args...) where {T, V <: AbstractVector{T}} =
    PolyLine(convert.(Point{3, Float64}, points), args...)

function PolyLine(points::AbstractVector{V};
    radius=0.0,
    closed=false,
    start_head=nothing,
    end_head=nothing) where {T, V <: AbstractVector{T}}
    PolyLine(points, radius, closed,
             start_head,
             end_head)
end

"""
A MeshFile is a wrapper around a mesh that indicates that the mesh should be
transmitted to the viewer by saving its data to a file and sending that file's
path to the viewer. This can avoid issues with very large meshes (which cannot
otherwise be sent over LCM), but will not work if the viewer is not running
locally.
"""
struct MeshFile <: AbstractGeometry{3, Float64}
    filename::String
end

function MeshFile(mesh::AbstractMesh)
    path = string(tempname(), ".ply")
    save(path, mesh)
    MeshFile(path)
end