text:
"""
    function symbols_(abm,exp,l=Symbol[])

Finds all symbols that have been used in a piece of code and show how they have been employed in the model.
"""
function symbols_(abm::Agent,exp::Expr,
                symb = DataFrame(Symbol = Symbol[], updated=Bool[], assigned=Bool[], referenced=Bool[],
                                called=Bool[], placeDeclaration = Symbol[], type = Symbol[]),
                called = false)

    for (pos,a) in enumerate(exp.args)
        if called
            aux = whereDeclared_(abm,a)
            for i in ENDSYMBOLS a = Meta.parse(split(string(a),i)[1]) end
            push!(symb,(a,false,false,false,true,aux[1],aux[2]))
            called = false
        elseif exp.head in UPDATINGOPERATORS && pos == 1
            if typeof(a) == Symbol
                aux = whereDeclared_(abm,a)
                for i in ENDSYMBOLS a = Meta.parse(split(string(a),i)[1]) end
                push!(symb,(a,true,false,false,false,aux[1],aux[2]))
            elseif typeof(a) == Expr && a.head == :ref
                aux = whereDeclared_(abm,a.args[1])
                push!(symb,(a.args[1],true,false,true,false,aux[1],aux[2]))
            end
        elseif exp.head in [:(=)] && pos == 1
            if typeof(a) == Symbol
                aux = whereDeclared_(abm,a)
                for i in ENDSYMBOLS a = Meta.parse(split(string(a),i)[1]) end
                push!(symb,(a,false,true,false,false,aux[1],aux[2]))
            elseif typeof(a) == Expr && a.head == :ref
                aux = whereDeclared_(abm,a.args[1])
                push!(symb,(a.args[1],false,true,true,false,aux[1],aux[2]))
            end
        elseif typeof(a) == Expr && a.head == :ref
            aux = whereDeclared_(abm,a.args[1])
            push!(symb,(a.args[1],false,false,true,false,aux[1],aux[2]))
        elseif typeof(a) == Expr && a.head == :call
            symbols_(abm,a,symb,true)
        elseif typeof(a) == Symbol
            aux = whereDeclared_(abm,a)
            for i in ENDSYMBOLS a = Meta.parse(split(string(a),i)[1]) end
            push!(symb,(a,false,false,false,false,aux[1],aux[2]))
        elseif typeof(a) == Expr
            symbols_(abm,a,symb)
        end
    end

    return symb
end

# """
#     function vectorize_(abm::Agent,code::Expr)

# Function to subtitute all the declared symbols of the agents in the expression into vector form.
# optional arguments base and update add some intermediate names to the vectorized variables and to updated ones.
# """
# function validSymbolsAgent_(abm::Agent)

#     use =  Symbol[]
#     declaration = Symbol[]

#     #Vectorisation changes
#     for (i,v) in enumerate(abm.declaredSymbols["Variable"])

#         push!(use["Variable"],v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))
#         push!(declaration["Variable"],v,Meta.parse(string(v,"₁")))
#         push!(allSymbols,v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))

#     end

#     for (i,v) in enumerate(abm.declaredSymbols["Local"])

#         push!(use["Local"],v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))
#         push!(declaration["Local"],v,Meta.parse(string(v,"₁")))
#         push!(allSymbols,v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))

#     end

#     for (i,v) in enumerate(abm.declaredSymbols["Interaction"])

#         push!(use["Interaction"],v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))
#         push!(declaration["Interaction"],v,Meta.parse(string(v,"₁")))
#         push!(allSymbols,v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))

#     end

#     for (i,v) in enumerate(abm.declaredSymbols["Global"])

#         push!(use["Global"],v)
#         push!(declaration["Global"],v)
#         push!(allSymbols,v)

#     end

#     for (i,v) in enumerate(abm.declaredSymbols["GlobalArray"])

#         push!(use["GlobalArray"],v)
#         push!(declaration["GlobalArray"],v)
#         push!(allSymbols,v)

#     end

#     for (i,v) in enumerate(abm.declaredSymbols["Identity"])

#         push!(use["Identity"],v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))
#         push!(declaration["Identity"],v,Meta.parse(string(v,"₁")))
#         push!(allSymbols,v,Meta.parse(string(v,"₁")),Meta.parse(string(v,"₂")))

#     end

#     return code
# end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
macro scalloptest()
  println("")
  reload(Pkg.dir("Scallop", "run_tests.jl"))
end

ext: jl


text:
# read the output plimit
export psp_readplimit

const READ_PLIMIT = [0x50, 0x0d]

"""
    psp_readplimit(io)

Read the power limit.
"""
function psp_readplimit(io_psp::IO)
  write(io_psp, READ_PLIMIT)
  plimitstring = bytestring(read(io_psp,UInt8, 7))
  if  uppercase(plimitstring[1]) != 'P' ||
      plimitstring[5:7] != "\r\r\n"
        error("PSP Invalid Response: ",plimitstring)
  end
  parse(Float64,plimitstring[2:4])
end


text:
# Intersect
# =========

type BAMIntersectionIterator
    reader::BAMReader
    chunks::Vector{Chunk}
    refid::Int
    interval::UnitRange{Int}
end

function Base.start(iter::BAMIntersectionIterator)
    if !isempty(iter.chunks)
        seek(iter.reader, first(iter.chunks).start)
    end
    rec = BAMRecord()
    return advance!(iter, rec, 1)
end

function Base.done(iter::BAMIntersectionIterator, i_rec)
    i, _ = i_rec
    return i > endof(iter.chunks)
end

function Base.next(iter::BAMIntersectionIterator, i_rec)
    i, rec = i_rec
    retrec = copy(rec)
    return retrec, advance!(iter, rec, i)
end

function advance!(iter, rec, i)
    while i ≤ endof(iter.chunks)
        chunk = iter.chunks[i]
        while virtualoffset(iter.reader.stream) < chunk.stop
            read!(iter.reader, rec)
            if isoverlapping(rec, iter.refid, iter.interval)
                return i, rec
            end
        end
        i += 1
        if i ≤ endof(iter.chunks)
            seek(iter.reader, iter.chunks[i].start)
        end
    end
    return i, rec
end

function isoverlapping(rec, refid_, interval)
    return ismapped(rec) &&
        refid(rec) == refid_ &&
        position(rec) ≤ last(interval) &&
        rightmost_position(rec) ≥ first(interval)
end

function Base.intersect(reader::BAMReader, refid::Integer, interval::UnitRange)
    if isnull(reader.index)
        error("no index")
    end
    chunks = overlapchunks(get(reader.index), refid, interval)
    return BAMIntersectionIterator(reader, chunks, refid, interval)
end

function Base.intersect(reader::BAMReader, refname::AbstractString, interval::UnitRange)
    refid = findfirst(reader.refseqnames, refname)
    if refid == 0
        error("sequence name $refname is not in the header")
    end
    return intersect(reader, refid, interval)
end


text:
using Sctp
import Sctp: SockaddrIn
import Sockets: IPv4
import Sockets: @ip_str
sctp_init(true)
# s1addr = SockaddrIn(ip"0.0.0.0")
@show s = SCTPSocket()
@show b = bind(s, ip"0.0.0.0")
if b < 0
  error("bind failed")
end
# listen(s1, 10)

# @show s2 = SCTPSocket()

payload = Cint[1, 2, 3, 4, 5]

addr = SockaddrIn(ip"127.0.0.1", 7)
#function test()
#  global t1 = @async begin
#    @show sendv(s2, s1addr, payload)
#  end

#  global t2 = @async begin
#    @show recvv(s1)
#  end
#  #
#end

function stop()
  # wait(t1)
  # wait(t2)

  close(s)
  # close(s2)
  sctp_finish()
end


text:
function deflection_model(loc, turbine::Turbine)

    dx = loc[1]-turbine.coord.x
    dy = loc[2]-(turbine.coord.y+deflection_y)
    dz = loc[3]-(turbine.coord.z+turbine.hub_height+deflection_z)

    del = sqrt(dy^2+dz^2)

end


text:
using ONNX, Flux, ProtoBuf
include("ops_tests.jl")

# test and 3v1d
main_test("$ONNX_TEST_PATH/test_and_bcast3v1d",
    read_output("$ONNX_TEST_PATH/test_and_bcast3v1d"),
        read_input("$ONNX_TEST_PATH/test_and_bcast3v1d")[1],
            read_input("$ONNX_TEST_PATH/test_and_bcast3v1d")[2])

# test and 3v2d
main_test("$ONNX_TEST_PATH/test_and_bcast3v2d",
    read_output("$ONNX_TEST_PATH/test_and_bcast3v2d"),
        read_input("$ONNX_TEST_PATH/test_and_bcast3v2d")[1],
            read_input("$ONNX_TEST_PATH/test_and_bcast3v2d")[2])

# test and 4v2d
main_test("$ONNX_TEST_PATH/test_and_bcast4v2d",
    read_output("$ONNX_TEST_PATH/test_and_bcast4v2d"),
        read_input("$ONNX_TEST_PATH/test_and_bcast4v2d")[1],
            read_input("$ONNX_TEST_PATH/test_and_bcast4v2d")[2])

# test and 4v3d
main_test("$ONNX_TEST_PATH/test_and_bcast4v3d",
    read_output("$ONNX_TEST_PATH/test_and_bcast4v3d"),
        read_input("$ONNX_TEST_PATH/test_and_bcast4v3d")[1],
            read_input("$ONNX_TEST_PATH/test_and_bcast4v3d")[2])

# test and 4v4d
main_test("$ONNX_TEST_PATH/test_and_bcast4v4d",
    read_output("$ONNX_TEST_PATH/test_and_bcast4v4d"),
        read_input("$ONNX_TEST_PATH/test_and_bcast4v4d")[1],
            read_input("$ONNX_TEST_PATH/test_and_bcast4v4d")[2])

#Test and2d
main_test("$ONNX_TEST_PATH/test_and2d",
          read_output("$ONNX_TEST_PATH/test_and2d"),
          read_input("$ONNX_TEST_PATH/test_and2d")[1],
          read_input("$ONNX_TEST_PATH/test_and2d")[2])

#Test and3d
main_test("$ONNX_TEST_PATH/test_and3d",
          read_output("$ONNX_TEST_PATH/test_and3d"),
          read_input("$ONNX_TEST_PATH/test_and3d")[1],
          read_input("$ONNX_TEST_PATH/test_and3d")[2])

#Test and4d
main_test("$ONNX_TEST_PATH/test_and4d",
          read_output("$ONNX_TEST_PATH/test_and4d"),
          read_input("$ONNX_TEST_PATH/test_and4d")[1],
          read_input("$ONNX_TEST_PATH/test_and4d")[2])

# Test identity
main_test("$ONNX_TEST_PATH/test_identity",
    read_output("$ONNX_TEST_PATH/test_identity"),
        read_input("$ONNX_TEST_PATH/test_identity")[1])

# Test equal
main_test("$ONNX_TEST_PATH/test_equal",
    read_output("$ONNX_TEST_PATH/test_equal"),
        read_input("$ONNX_TEST_PATH/test_equal")[1],
            read_input("$ONNX_TEST_PATH/test_equal")[2])

# Test equal bcast
main_test("$ONNX_TEST_PATH/test_equal_bcast",
    read_output("$ONNX_TEST_PATH/test_equal_bcast"),
        read_input("$ONNX_TEST_PATH/test_equal_bcast")[1],
            read_input("$ONNX_TEST_PATH/test_equal_bcast")[2])

# Test greater
main_test("$ONNX_TEST_PATH/test_greater",
    read_output("$ONNX_TEST_PATH/test_greater"),
        read_input("$ONNX_TEST_PATH/test_greater")[1],
            read_input("$ONNX_TEST_PATH/test_greater")[2])

# Test greater bcast
main_test("$ONNX_TEST_PATH/test_greater_bcast",
    read_output("$ONNX_TEST_PATH/test_greater_bcast"),
        read_input("$ONNX_TEST_PATH/test_greater_bcast")[1],
            read_input("$ONNX_TEST_PATH/test_greater_bcast")[2])

# test xor2d
main_test("$ONNX_TEST_PATH/test_xor2d",
    read_output("$ONNX_TEST_PATH/test_xor2d"),
        read_input("$ONNX_TEST_PATH/test_xor2d")[1],
            read_input("$ONNX_TEST_PATH/test_xor2d")[2])

# test xor3d
main_test("$ONNX_TEST_PATH/test_xor3d",
    read_output("$ONNX_TEST_PATH/test_xor3d"),
        read_input("$ONNX_TEST_PATH/test_xor3d")[1],
            read_input("$ONNX_TEST_PATH/test_xor3d")[2])

# test xor4d
main_test("$ONNX_TEST_PATH/test_xor4d",
    read_output("$ONNX_TEST_PATH/test_xor4d"),
        read_input("$ONNX_TEST_PATH/test_xor4d")[1],
            read_input("$ONNX_TEST_PATH/test_xor4d")[2])

# test xor bcast 3v1d
main_test("$ONNX_TEST_PATH/test_xor_bcast3v1d",
    read_output("$ONNX_TEST_PATH/test_xor_bcast3v1d"),
        read_input("$ONNX_TEST_PATH/test_xor_bcast3v1d")[1],
            read_input("$ONNX_TEST_PATH/test_xor_bcast3v1d")[2])

# test xor bcast 3v2d
main_test("$ONNX_TEST_PATH/test_xor_bcast3v2d",
    read_output("$ONNX_TEST_PATH/test_xor_bcast3v2d"),
        read_input("$ONNX_TEST_PATH/test_xor_bcast3v2d")[1],
            read_input("$ONNX_TEST_PATH/test_xor_bcast3v2d")[2])

# test xor bcast 4v2d
main_test("$ONNX_TEST_PATH/test_xor_bcast4v2d",
    read_output("$ONNX_TEST_PATH/test_xor_bcast4v2d"),
        read_input("$ONNX_TEST_PATH/test_xor_bcast4v2d")[1],
            read_input("$ONNX_TEST_PATH/test_xor_bcast4v2d")[2])

# test xor bcast 4v3d
main_test("$ONNX_TEST_PATH/test_xor_bcast4v3d",
    read_output("$ONNX_TEST_PATH/test_xor_bcast4v3d"),
        read_input("$ONNX_TEST_PATH/test_xor_bcast4v3d")[1],
            read_input("$ONNX_TEST_PATH/test_xor_bcast4v3d")[2])

# test xor bcast 4v4d
main_test("$ONNX_TEST_PATH/test_xor_bcast4v4d",
    read_output("$ONNX_TEST_PATH/test_xor_bcast4v4d"),
        read_input("$ONNX_TEST_PATH/test_xor_bcast4v4d")[1],
            read_input("$ONNX_TEST_PATH/test_xor_bcast4v4d")[2])


text:
using Compat.SparseArrays
using MathOptInterface
using Cbc.CbcCInterface

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const CbcCI = CbcCInterface

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::CbcModel
    Optimizer() = new(CbcModel()) # Initializes with an empty model
end

struct CbcModelFormat
    num_rows::Int
    num_cols::Int
    row_idx::Vector{Int}
    col_idx::Vector{Int}
    values::Vector{Float64}
    col_lb::Vector{Float64}
    col_ub::Vector{Float64}
    obj::Vector{Float64}
    row_lb::Vector{Float64}
    row_ub::Vector{Float64}
    function CbcModelFormat(num_rows::Integer, num_cols::Integer)
        obj = fill(0.0, num_cols)
        row_idx = Int[]
        col_idx = Int[]
        values = Float64[]
        col_lb = fill(-Inf, num_cols)
        col_ub = fill(Inf, num_cols)
        row_lb = fill(-Inf, num_rows)
        row_ub = fill(Inf, num_rows)
        constraint_matrix = Tuple{Int,Int,Float64}[]
        # An `InexactError` might occur if `num_rows` or `num_cols` is too
        # large, e.g., if `num_cols isa Int64` and is larger than 2^31 on a
        # 32-bit hardware
        new(num_rows, num_cols, row_idx, col_idx, values, col_lb, col_ub,
            obj, row_lb, row_ub)
    end
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap,f::MOI.SingleVariable,
                         s::MOI.EqualTo)

    cbc_model_format.col_lb[mapping.varmap[f.variable].value] = s.value
    cbc_model_format.col_ub[mapping.varmap[f.variable].value] = s.value
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.SingleVariable,
                         s::MOI.LessThan)

    cbc_model_format.col_ub[mapping.varmap[f.variable].value] = s.upper
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.SingleVariable,
                         s::MOI.GreaterThan)

    cbc_model_format.col_lb[mapping.varmap[f.variable].value] = s.lower
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.SingleVariable,
                         s::MOI.Interval)

    cbc_model_format.col_lb[mapping.varmap[f.variable].value] = s.lower
    cbc_model_format.col_ub[mapping.varmap[f.variable].value] = s.upper
end

function push_terms(row_idx::Vector{Int}, col_idx::Vector{Int},
                    values::Vector{Float64}, ci::MOI.ConstraintIndex{F,S},
                    terms::Vector{MOI.ScalarAffineTerm{Float64}},
                    mapping::MOIU.IndexMap) where {F,S}

    for term in terms
        push!(row_idx, mapping.conmap[ci].value)
        push!(col_idx, mapping.varmap[term.variable_index].value)
        push!(values, term.coefficient)
    end
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.ScalarAffineFunction,
                         s::MOI.EqualTo)

    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_lb[mapping.conmap[ci].value] = s.value - f.constant
    cbc_model_format.row_ub[mapping.conmap[ci].value] = s.value - f.constant
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.ScalarAffineFunction,
                         s::MOI.GreaterThan)

    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_lb[mapping.conmap[ci].value] = s.lower - f.constant
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.ScalarAffineFunction,
                         s::MOI.LessThan)

    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_ub[mapping.conmap[ci].value] = s.upper - f.constant
end

function load_constraint(ci::MOI.ConstraintIndex,
                         cbc_model_format::CbcModelFormat,
                         mapping::MOIU.IndexMap, f::MOI.ScalarAffineFunction,
                         s::MOI.Interval)

    push_terms(cbc_model_format.row_idx, cbc_model_format.col_idx,
               cbc_model_format.values, ci, f.terms, mapping)
    cbc_model_format.row_ub[mapping.conmap[ci].value] = s.upper - f.constant
    cbc_model_format.row_lb[mapping.conmap[ci].value] = s.lower - f.constant
end

function load_obj(cbc_model_format::CbcModelFormat, mapping::MOIU.IndexMap,
    f::MOI.ScalarAffineFunction)
    # We need to increment values of objective function with += to handle
    # cases like $x_1 + x_2 + x_1$. This is safe because objective function
    # is initialized with zeros in the constructor.
    for term in f.terms
        cbc_model_format.obj[mapping.varmap[term.variable_index].value] +=
            term.coefficient
    end
end

function copy_constraints!(cbc_model_format::CbcModelFormat,
                           user_optimizer::MOI.ModelLike,
                           mapping::MOIU.IndexMap)

    for (F,S) in MOI.get(user_optimizer, MOI.ListOfConstraints())
        if !(S <: Union{MOI.ZeroOne, MOI.Integer})
            for ci in MOI.get(user_optimizer, MOI.ListOfConstraintIndices{F,S}())
                f = MOI.get(user_optimizer, MOI.ConstraintFunction(), ci)
                s = MOI.get(user_optimizer,  MOI.ConstraintSet(), ci)
                load_constraint(ci, cbc_model_format, mapping, f, s)
            end
        end
    end
end

function update_bounds_for_binary_vars!(col_lb::Vector{Float64},
                                        col_ub::Vector{Float64},
                                        zero_one_indices::Vector{Int})

    for idx in zero_one_indices
        if col_lb[idx] < 0.0
            col_lb[idx] = 0.0
        end
        if col_ub[idx] > 1.0
            col_ub[idx] = 1.0
        end
    end
end

function update_indices(user_optimizer::MOI.ModelLike,
                        mapping::MOIU.IndexMap,
                        ci::Vector{MOI.ConstraintIndex{F,S}},
                        list_of_indices::Vector{Int}) where {F, S}

    for i in 1:length(ci)
        f = MOI.get(user_optimizer, MOI.ConstraintFunction(), ci[i])
        push!(list_of_indices, mapping.varmap[f.variable].value)
    end
end

# This function creates MOI.ConstraintIndexes for MOI.SingleVariable constraints
# and updates the mapping with such indices. SingleVariables are seen as bounds by Cbc.
# This function also updates the vectors zero_one_indices and integer_indices.
function create_indices_for_bounds(user_optimizer::MOI.ModelLike, mapping::MOIU.IndexMap,
                                   zero_one_indices::Vector{Int}, integer_indices::Vector{Int})
    indices = MOI.ConstraintIndex[]
    num_bounds = 0
    list_of_constraints = MOI.get(user_optimizer, MOI.ListOfConstraints())
    for (F,S) in list_of_constraints
        if F == MOI.SingleVariable
            list_of_ci = MOI.get(user_optimizer, MOI.ListOfConstraintIndices{F,S}())
            if S == MOI.ZeroOne
                update_indices(user_optimizer, mapping, list_of_ci, zero_one_indices)
            elseif S == MOI.Integer
                update_indices(user_optimizer, mapping, list_of_ci, integer_indices)
            end
            for i in 1:length(list_of_ci)
                mapping.conmap[list_of_ci[i]] = MOI.ConstraintIndex{F,S}(num_bounds + i)
            end
            num_bounds += MOI.get(user_optimizer, MOI.NumberOfConstraints{F,S}())
        end
    end
end

# This function creates MOI.ConstraintIndex's for constraints that are not
# MOI.SingleVariable constraints and updates the mapping with such indices.
function create_indices_for_rows(cbc_optimizer::Optimizer, user_optimizer::MOI.ModelLike,
                                 mapping::MOIU.IndexMap)
    list_of_constraints = MOI.get(user_optimizer, MOI.ListOfConstraints())
    num_rows = 0
    indices = MOI.ConstraintIndex[]
    for (F,S) in list_of_constraints
        if !(MOI.supports_constraint(cbc_optimizer, F, S))
            throw(MOI.UnsupportedConstraint{F,S}("Cbc MOI Interface does not support constraints of type " * (F,S) * "."))
        end
        if F != MOI.SingleVariable
            list_of_ci = MOI.get(user_optimizer, MOI.ListOfConstraintIndices{F,S}())
            for i in 1:length(list_of_ci)
                mapping.conmap[list_of_ci[i]] = MOI.ConstraintIndex{F,S}(num_rows + i)
            end
            num_rows += MOI.get(user_optimizer, MOI.NumberOfConstraints{F,S}())
        end
    end
end

function MOI.copy_to(cbc_optimizer::Optimizer,
    user_optimizer::MOI.ModelLike; copy_names=false)

    mapping = MOIU.IndexMap()

    num_cols = MOI.get(user_optimizer, MOI.NumberOfVariables())
    var_index = MOI.get(user_optimizer, MOI.ListOfVariableIndices())
    for i in 1:num_cols
        mapping.varmap[var_index[i]] = MOI.VariableIndex(i)
    end

    zero_one_indices = Int[]
    integer_indices = Int[]
    create_indices_for_rows(cbc_optimizer, user_optimizer, mapping) # Updates mapping
    num_rows = length(mapping.conmap)
    # Rows are not created for bounds, but the mapping is still updated.
    create_indices_for_bounds(user_optimizer, mapping, zero_one_indices, integer_indices)

    cbc_model_format = CbcModelFormat(num_rows, num_cols)
    copy_constraints!(cbc_model_format, user_optimizer, mapping)
    update_bounds_for_binary_vars!(cbc_model_format.col_lb,
                                   cbc_model_format.col_ub, zero_one_indices)

    # Copy the objective function.
    objF = MOI.get(user_optimizer,
                   MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    load_obj(cbc_model_format, mapping, objF)
    sense = MOI.get(user_optimizer, MOI.ObjectiveSense())
    MOI.set(cbc_optimizer, MOI.ObjectiveSense(), sense)

    # Load the problem into Cbc.
    CbcCI.loadProblem(cbc_optimizer.inner,
                      sparse(cbc_model_format.row_idx, cbc_model_format.col_idx,
                             cbc_model_format.values, cbc_model_format.num_rows,
                             cbc_model_format.num_cols),
                      cbc_model_format.col_lb, cbc_model_format.col_ub,
                      cbc_model_format.obj, cbc_model_format.row_lb,
                      cbc_model_format.row_ub)

    empty!(cbc_model_format.row_idx)
    empty!(cbc_model_format.col_idx)
    empty!(cbc_model_format.values)

    ## Set the integer variables.
    for idx in vcat(integer_indices, zero_one_indices)
        CbcCI.setInteger(cbc_optimizer.inner, idx-1)
    end

    return MOIU.IndexMap(mapping.varmap, mapping.conmap)
end

function MOI.optimize!(cbc_optimizer::Optimizer)
    # Call solve function
    CbcCI.solve(cbc_optimizer.inner)
end

function MOI.add_variable(cbc_optimizer::Optimizer)
    throw(MOI.AddVariableNotAllowed())
end

## supports constraints

MOI.supports_constraint(::Optimizer, ::Type{<:Union{MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable}},
::Type{<:Union{MOI.EqualTo{Float64}, MOI.Interval{Float64}, MOI.LessThan{Float64},
MOI.GreaterThan{Float64}, MOI.ZeroOne, MOI.Integer}}) = true

MOI.supports(cbc_optimizer::Optimizer, object::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true

MOI.supports(cbc_optimizer::Optimizer, object::MOI.ObjectiveSense) = true

## Set functions

function MOI.write_to_file(cbc_optimizer::Optimizer, filename::String)
    if !endswith("filename", "mps")
        error("Optimizer only supports writing .mps files")
    else
        writeMps(cbc_optimizer.inner, filename)
    end
end

# empty!
function MOI.empty!(cbc_optimizer::Optimizer)
    cbc_optimizer.inner = CbcModel()
end

function MOI.set(cbc_optimizer::Optimizer, object::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense == MOI.MAX_SENSE
        CbcCI.setObjSense(cbc_optimizer.inner, -1)
    else ## Other senses are set as minimization (cbc default)
        CbcCI.setObjSense(cbc_optimizer.inner, 1)
    end
end

## Get functions

function MOI.is_empty(cbc_optimizer::Optimizer)
    return (CbcCI.getNumCols(cbc_optimizer.inner) == 0 && CbcCI.getNumRows(cbc_optimizer.inner) == 0)
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.NumberOfVariables)
    return getNumCols(cbc_optimizer.inner)
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.ObjectiveBound)
    return CbcCI.getBestPossibleObjValue(cbc_optimizer.inner)
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.NodeCount)
    return CbcCI.getNodeCount(cbc_optimizer.inner)
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.ObjectiveValue)
    return CbcCI.getObjValue(cbc_optimizer.inner)
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.VariablePrimal,
                 ref::MOI.VariableIndex)

    variablePrimals = CbcCI.getColSolution(cbc_optimizer.inner)
    return variablePrimals[ref.value]
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.VariablePrimal,
                 ref::Vector{MOI.VariableIndex})

    variablePrimals = CbcCI.getColSolution(cbc_optimizer.inner)
    return [variablePrimals[vi.value] for vi in ref]
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.ResultCount)
    if (isProvenInfeasible(cbc_optimizer.inner) ||
        isContinuousUnbounded(cbc_optimizer.inner) ||
        isAbandoned(cbc_optimizer.inner) ||
        CbcCI.getObjValue(cbc_optimizer.inner) >= 1e300)

        return 0
    end
    return 1
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.ObjectiveSense)
    CbcCI.getObjSense(cbc_optimizer.inner) == 1 && return MOI.MIN_SENSE
    CbcCI.getObjSense(cbc_optimizer.inner) == -1 && return MOI.MAX_SENSE
end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.TerminationStatus)

    if isProvenInfeasible(cbc_optimizer.inner)
        return MOI.INFEASIBLE
    elseif isContinuousUnbounded(cbc_optimizer.inner)
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif isNodeLimitReached(cbc_optimizer.inner)
        return MOI.NODE_LIMIT
    elseif isSecondsLimitReached(cbc_optimizer.inner)
        return MOI.TIME_LIMIT
    elseif isSolutionLimitReached(cbc_optimizer.inner)
        return MOI.SOLUTION_LIMIT
    elseif (isProvenOptimal(cbc_optimizer.inner) ||
            isInitialSolveProvenOptimal(cbc_optimizer.inner) ||
            MOI.get(cbc_optimizer, MOI.ResultCount()) == 1)
        return MOI.OPTIMAL
    elseif isAbandoned(cbc_optimizer.inner)
        return MOI.INTERRUPTED
    elseif optimizeNotCalled(cbc_optimizer.inner)
        return MOI.OPTIMIZE_NOT_CALLED
    else
        error("Internal error: Unrecognized solution status:",
              " status = $(status(cbc_optimizer.inner)),",
              " secondaryStatus = $(secondaryStatus(cbc_optimizer.inner))")
    end

end

function MOI.get(cbc_optimizer::Optimizer, object::MOI.PrimalStatus)
    if isProvenOptimal(cbc_optimizer.inner) || isInitialSolveProvenOptimal(cbc_optimizer.inner)
        return MOI.FEASIBLE_POINT
    elseif isProvenInfeasible(cbc_optimizer.inner)
        return MOI.INFEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function MOI.get(cbc_optimizer::Optimizer, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end


text:
function init_pot_Hcoul( LF, center )
  Npoints = LF.Nx * LF.Ny * LF.Nz
  Vpot = zeros( Float64,Npoints )
  for ip = 1:Npoints
    r = norm( LF.lingrid[:,ip] - center[:] )
    Vpot[ip] = -1.0/r
  end
  return Vpot
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using InteractiveUtils, Test

@testset "warntype highlighting" begin
    # Make sure that "expected" unions are highlighted with warning color instead of error color
    io = IOBuffer()
    f(x) = x > 1 ? "foo" : nothing
    c = Base.text_colors[Base.warn_color()]
    InteractiveUtils.highlighting[:warntype] = false
    code_warntype(IOContext(io, :color => true), f, Tuple{Int64})
    @test !occursin(c, String(take!(io)))
    InteractiveUtils.highlighting[:warntype] = true
    code_warntype(IOContext(io, :color => true), f, Tuple{Int64})
    @test occursin(c, String(take!(io)))
end

@testset "no colors" begin
    io = IOBuffer()
    code_llvm(io, sqrt, (Float32,))
    @test !occursin("\e", String(take!(io)))

    InteractiveUtils.print_llvm(io, "; comment")
    @test String(take!(io)) == "; comment\n"

    InteractiveUtils.highlighting[:llvm] = false
    code_llvm(IOContext(io, :color => true), abs, (Float64,))
    @test !occursin("\e", String(take!(io)))
    InteractiveUtils.highlighting[:llvm] = true
    code_llvm(IOContext(io, :color => true), abs, (Float64,))
    @test occursin("\e", String(take!(io)))

    code_native(io, sqrt, (Float32,))
    @test !occursin("\e", String(take!(io)))

    InteractiveUtils.print_native(io, "; comment", :x86)
    @test String(take!(io)) == "; comment\n"

    InteractiveUtils.print_native(io, "; comment", :unknown)
    @test String(take!(io)) == "; comment\n"

    InteractiveUtils.highlighting[:native] = false
    code_native(IOContext(io, :color => true), abs, (Float64,))
    @test !occursin("\e", String(take!(io)))
    InteractiveUtils.highlighting[:native] = true
    code_native(IOContext(io, :color => true), abs, (Float64,))
    @test occursin("\e", String(take!(io)))
end

function hilight_llvm(s)
    io = IOBuffer()
    InteractiveUtils.print_llvm(IOContext(io, :color=>true), s)
    r = String(take!(io))
    get(ENV, "JULIA_DEBUG", "") == "" && return r
    println(stdout, " input: ", s)
    println(stdout, "result: ", r)
    flush(stdout)
    r
end
function hilight_native(s, arch)
    io = IOBuffer()
    InteractiveUtils.print_native(IOContext(io, :color=>true), s, arch)
    r = String(take!(io))
    get(ENV, "JULIA_DEBUG", "") == "" && return r
    println(stdout, " input: ", s)
    println(stdout, "result: ", r)
    flush(stdout)
    r
end
hilight_x86(s) = hilight_native(s, :x86)
hilight_arm(s) = hilight_native(s, :arm)

function esc_code(s)
    io = IOBuffer()
    InteractiveUtils.printstyled_ll(IOContext(io, :color=>true), "!", s)
    split(String(take!(io)), "!")
end

const D, XD = esc_code(:default)
const C, XC = esc_code(:comment)
const L, XL = esc_code(:label)
const I, XI = esc_code(:instruction)
const T, XT = esc_code(:type)
const N, XN = esc_code(:number)
const B, XB = esc_code(:bracket)
const V, XV = esc_code(:variable)
const K, XK = esc_code(:keyword)
const F, XF = esc_code(:funcname)

const COM = D * "," * XD
const COL = D * ":" * XD
const EQU = D * "=" * XD
const P = B * "(" * XB
const XP = B * ")" * XB
const S = B * "[" * XB
const XS = B * "]" * XB
const U = B * "{" * XB
const XU = B * "}" * XB

@testset "LLVM IR" begin
    @testset "comment" begin
        @test hilight_llvm("; comment ; // # ") == "$(C); comment ; // # $(XC)\n"
    end
    @testset "lavel" begin
        @test hilight_llvm("top:") == "$(L)top:$(XL)\n"

        @test hilight_llvm("L7:\t\t; preds = %top") ==
            "$(L)L7:$(XL)\t\t$(C); preds = %top$(XC)\n"
    end
    @testset "define" begin
        @test hilight_llvm("define double @julia_func_1234(float) {") ==
            "$(K)define$(XK) $(T)double$(XT) " *
            "$(F)@julia_func_1234$(XF)$P$(T)float$(XT)$XP $U\n"

        @test hilight_llvm("}") == "$XU\n"
    end

    @testset "declare" begin
        @test hilight_llvm("declare i32 @jl_setjmp(i8*) #2") ==
            "$(K)declare$(XK) $(T)i32$(XT) " *
            "$(F)@jl_setjmp$(XF)$P$(T)i8$(XT)$(D)*$(XD)$XP $(D)#2$(XD)\n"
    end

    @testset "type" begin
        @test hilight_llvm("%jl_value_t = type opaque") ==
            "$(V)%jl_value_t$(XV) $EQU $(K)type$(XK) $(T)opaque$(XT)\n"
    end

    @testset "target" begin
        datalayout = "e-m:e-i64:64-f80:128-n8:16:32:64-S128-ni:10:11:12:13"
        @test hilight_llvm("target datalayout = \"$datalayout\"") ==
            "$(K)target$(XK) $(K)datalayout$(XK) $EQU $(V)\"$datalayout\"$(XV)\n"
    end

    @testset "attributes" begin
        @test hilight_llvm(
            """attributes #1 = { uwtable "frame-pointer"="all" "thunk" }""") ==
            "$(K)attributes$(XK) $(D)#1$(XD) $EQU " *
            "$U $(K)uwtable$(XK) $(V)\"frame-pointer\"$(XV)$EQU" *
            "$(V)\"all\"$(XV) $(V)\"thunk\"$(XV) $XU\n"
    end

    @testset "terminator" begin
        @test hilight_llvm("  ret i8 %12") ==
            "  $(I)ret$(XI) $(T)i8$(XT) $(V)%12$(XV)\n"

        @test hilight_llvm("  br i1 %2, label %L6, label %L4") ==
            "  $(I)br$(XI) $(T)i1$(XT) $(V)%2$(XV)$COM " *
            "$(T)label$(XT) $(L)%L6$(XL)$COM $(T)label$(XT) $(L)%L4$(XL)\n"

        @test hilight_llvm("  br label %L5") ==
            "  $(I)br$(XI) $(T)label$(XT) $(L)%L5$(XL)\n"

        @test hilight_llvm("  unreachable") == "  $(I)unreachable$(XI)\n"
    end

    @testset "arithmetic" begin
        @test hilight_llvm("   %11 = add nuw nsw i64 %value_phi10, 1") ==
            "   $(V)%11$(XV) $EQU $(I)add$(XI) $(K)nuw$(XK) $(K)nsw$(XK) " *
            "$(T)i64$(XT) $(V)%value_phi10$(XV)$COM $(N)1$(XN)\n"

        @test hilight_llvm("   %13 = fadd double %12, -2.000000e+00") ==
            "   $(V)%13$(XV) $EQU $(I)fadd$(XI) " *
            "$(T)double$(XT) $(V)%12$(XV)$COM $(N)-2.000000e+00$(XN)\n"

        @test hilight_llvm("      %21 = fmul contract double %20, 0x0123456789ABCDEF") ==
            "      $(V)%21$(XV) $EQU $(I)fmul$(XI) $(K)contract$(XK) " *
            "$(T)double$(XT) $(V)%20$(XV)$COM $(N)0x0123456789ABCDEF$(XN)\n"
    end

    @testset "bitwise" begin
        @test hilight_llvm("   %31 = shl i64 %value_phi4, 52") ==
            "   $(V)%31$(XV) $EQU " *
            "$(I)shl$(XI) $(T)i64$(XT) $(V)%value_phi4$(XV)$COM $(N)52$(XN)\n"
    end

    @testset "aggregate" begin
        @test hilight_llvm("    %4 = extractvalue { i64, i1 } %1, 0") ==
            "    $(V)%4$(XV) $EQU $(I)extractvalue$(XI) " *
            "$U $(T)i64$(XT)$COM $(T)i1$(XT) $XU $(V)%1$(XV)$COM $(N)0$(XN)\n"
    end

    @testset "memory access" begin
        @test hilight_llvm("  %dims = alloca [1 x i64], align 8") ==
            "  $(V)%dims$(XV) $EQU $(I)alloca$(XI) " *
            "$S$(N)1$(XN) $(D)x$(XD) $(T)i64$(XT)$XS$COM $(K)align$(XK) $(N)8$(XN)\n"

        @test hilight_llvm("    %51 = load i32," *
                           " i32* inttoptr (i64 226995504 to i32*), align 16") ==
            "    $(V)%51$(XV) $EQU $(I)load$(XI) $(T)i32$(XT)$COM " *
            "$(T)i32$(XT)$(D)*$(XD) $(K)inttoptr$(XK) $P$(T)i64$(XT) $(N)226995504$(XN) " *
            "$(K)to$(XK) $(T)i32$(XT)$(D)*$(XD)$XP$COM $(K)align$(XK) $(N)16$(XN)\n"

        @test hilight_llvm("    %53 = load %jl_value_t addrspace(10)*, " *
                           "%jl_value_t addrspace(10)* addrspace(11)* %52, align 8") ==
            "    $(V)%53$(XV) $EQU $(I)load$(XI) $(V)%jl_value_t$(XV) " *
            "$(K)addrspace$(XK)$P$(N)10$(XN)$XP$(D)*$(XD)$COM " *
            "$(V)%jl_value_t$(XV) $(K)addrspace$(XK)$P$(N)10$(XN)$XP$(D)*$(XD) " *
            "$(K)addrspace$(XK)$P$(N)11$(XN)$XP$(D)*$(XD) " *
            "$(V)%52$(XV)$COM $(K)align$(XK) $(N)8$(XN)\n"

        @test hilight_llvm("    store i64 %61, i64 addrspace(11)* %60, align 8") ==
            "    $(I)store$(XI) $(T)i64$(XT) $(V)%61$(XV)$COM " *
            "$(T)i64$(XT) $(K)addrspace$(XK)$P$(N)11$(XN)$XP$(D)*$(XD) " *
            "$(V)%60$(XV)$COM $(K)align$(XK) $(N)8$(XN)\n"

        @test hilight_llvm("  store volatile %jl_value_t addrspace(10)** %62, " *
                           "%jl_value_t addrspace(10)*** %63, align 8") ==
            "  $(I)store$(XI) $(K)volatile$(XK) $(V)%jl_value_t$(XV) " *
            "$(K)addrspace$(XK)$P$(N)10$(XN)$XP$(D)**$(XD) $(V)%62$(XV)$COM " *
            "$(V)%jl_value_t$(XV) $(K)addrspace$(XK)$P$(N)10$(XN)$XP$(D)***$(XD) " *
            "$(V)%63$(XV)$COM $(K)align$(XK) $(N)8$(XN)\n"

        @test hilight_llvm("     %71 = getelementptr i8, i8* %70, i64 8") ==
            "     $(V)%71$(XV) $EQU $(I)getelementptr$(XI) $(T)i8$(XT)$COM " *
            "$(T)i8$(XT)$(D)*$(XD) $(V)%70$(XV)$COM $(T)i64$(XT) $(N)8$(XN)\n"
    end

    @testset "conversion" begin
        @test hilight_llvm("  %22 = zext i1 %21 to i8") ==
            "  $(V)%22$(XV) $EQU $(I)zext$(XI) $(T)i1$(XT) $(V)%21$(XV) " *
            "$(K)to$(XK) $(T)i8$(XT)\n"

        @test hilight_llvm("     %24 = sitofp i64 %23 to double") ==
            "     $(V)%24$(XV) $EQU $(I)sitofp$(XI) $(T)i64$(XT) $(V)%23$(XV) " *
            "$(K)to$(XK) $(T)double$(XT)\n"

        @test hilight_llvm("  %26 = ptrtoint i8* %25 to i64") ==
            "  $(V)%26$(XV) $EQU $(I)ptrtoint$(XI) $(T)i8$(XT)$(D)*$(XD) " *
            "$(V)%25$(XV) $(K)to$(XK) $(T)i64$(XT)\n"

        @test hilight_llvm("  %28 = bitcast %jl_value_t addrspace(10)* %27 " *
                           "to [2 x i16] addrspace(10)*") ==
            "  $(V)%28$(XV) $EQU $(I)bitcast$(XI) $(V)%jl_value_t$(XV) " *
            "$(K)addrspace$(XK)$P$(N)10$(XN)$XP$(D)*$(XD) $(V)%27$(XV) " *
            "$(K)to$(XK) $S$(N)2$(XN) $(D)x$(XD) $(T)i16$(XD)$XS " *
            "$(K)addrspace$(XK)$P$(N)10$(XN)$XP$(D)*$(XD)\n"
    end

    @testset "other" begin
        @test hilight_llvm("  %31 = icmp slt i64 %30, 0") ==
            "  $(V)%31$(XV) $EQU $(I)icmp$(XI) $(I)slt$(XI) " *
            "$(T)i64$(XT) $(V)%30$(XV)$COM $(N)0$(XN)\n"

        @test hilight_llvm("  %value_phi34 = phi double [ %33, %L50 ], [ %32, %L60 ]") ==
            "  $(V)%value_phi34$(XV) $EQU $(I)phi$(XI) $(T)double$(XT) " *
            "$S $(V)%33$(XV)$COM $(L)%L50$(XL) $XS$COM " *
            "$S $(V)%32$(XV)$COM $(L)%L60$(XL) $XS\n"

        @test hilight_llvm("   %.v = select i1 %35, i64 %36, i64 63") ==
            "   $(V)%.v$(XV) $EQU $(I)select$(XI) $(T)i1$(XT) $(V)%35$(XV)$COM " *
            "$(T)i64$(XT) $(V)%36$(XV)$COM $(T)i64$(XT) $(N)63$(XN)\n"

        @test hilight_llvm("   %38 = call i64 @llvm.cttz.i64(i64 %37, i1 false)") ==
            "   $(V)%38$(XV) $EQU $(I)call$(XI) $(T)i64$(XT) " *
            "$(F)@llvm.cttz.i64$(XF)$P$(T)i64$(XT) $(V)%37$(XV)$COM " *
            "$(T)i1$(XT) $(K)false$(XK)$XP\n"
    end
end

@testset "x86 ASM" begin
    @testset "comment" begin
        @test hilight_x86("; comment ; // # ") == "$(C); comment ; // # $(XC)\n"
    end
    @testset "label" begin
        @test hilight_x86("L123:") == "$(L)L123:$(XL)\n"
    end
    @testset "directive" begin
        @test hilight_x86("\t.text") == "\t$(D).text$(XD)\n"
    end

    @testset "0-operand" begin
        # AT&T
        @test hilight_x86("\tretq") == "\t$(I)retq$(XI)\n"

        # Intel
        @test hilight_x86("\tret") == "\t$(I)ret$(XI)\n"
    end
    @testset "1-operand" begin
        # AT&T
        @test hilight_x86("\tpopq\t%rax") == "\t$(I)popq$(XI)\t$(V)%rax$(XV)\n"

        @test hilight_x86("\tpushl\t\$4294967295\t# imm = 0xFFFFFFFF") ==
            "\t$(I)pushl$(XI)\t$(N)\$4294967295$(XN)\t$(C)# imm = 0xFFFFFFFF$(XC)\n"

        @test hilight_x86("\tja\tL234") == "\t$(I)ja$(XI)\t$(L)L234$(XL)\n"

        @test hilight_x86("\tnopw\t%cs:(%rax,%rax)") ==
            "\t$(I)nopw$(XI)\t$(V)%cs$(XV)$COL$P$(V)%rax$(XV)$COM$(V)%rax$(XV)$XP\n"

        # Intel
        @test hilight_x86("\tpop\trax") == "\t$(I)pop$(XI)\t$(V)rax$(XV)\n"

        @test hilight_x86("\tpush\t4294967295") ==
            "\t$(I)push$(XI)\t$(N)4294967295$(XN)\n"

        @test hilight_x86("\tja\tL234") == "\t$(I)ja$(XI)\t$(L)L234$(XL)\n"

        @test hilight_x86("\tnop\tword ptr cs:[rax + rax]") ==
            "\t$(I)nop$(XI)\t$(K)word$(XK) $(K)ptr$(XK) " *
            "$(V)cs$(XV)$COL$S$(V)rax$(XV) $(D)+$(XD) $(V)rax$(XV)$XS\n"
    end
    @testset "2-operand" begin
        # AT&T
        @test hilight_x86("\tshrq\t\$63, %rcx") ==
            "\t$(I)shrq$(XI)\t$(N)\$63$(XN)$COM $(V)%rcx$(XV)\n"

        @test hilight_x86("\tvmovsd\t(%rsi,%rdx,8), %xmm1\t# xmm1 = mem[0],zero") ==
            "\t$(I)vmovsd$(XI)\t$P$(V)%rsi$(XV)$COM$(V)%rdx$(XV)$COM$(N)8$(XN)$XP" *
            "$COM $(V)%xmm1$(XV)\t$(C)# xmm1 = mem[0],zero$(XC)\n"

        @test hilight_x86("\tmovabsq\t\$\"#string#338\", %rax") ==
            "\t$(I)movabsq$(XI)\t$(F)\$\"#string#338\"$(XF)$COM $(V)%rax$(XV)\n"

        # Intel
        @test hilight_x86("\tshr\trcx, 63") ==
            "\t$(I)shr$(XI)\t$(V)rcx$(XV)$COM $(N)63$(XN)\n"

        @test hilight_x86(
            "\tvmovsd\txmm1, dword ptr [rsi + 8*rdx]\t# xmm1 = mem[0],zero") ==
            "\t$(I)vmovsd$(XI)\t$(V)xmm1$(XV)$COM $(K)dword$(XK) $(K)ptr$(XK) " *
            "$S$(V)rsi$(XV) $(D)+$(XD) $(N)8$(XN)$(D)*$(XD)$(V)rdx$(XV)$XS" *
            "\t$(C)# xmm1 = mem[0],zero$(XC)\n"

        @test hilight_x86("\tmovabs\trax, offset \"#string#338\"") ==
            "\t$(I)movabs$(XI)\t$(V)rax$(XV)$COM " *
            "$(K)offset$(XK) $(F)\"#string#338\"$(XF)\n"
    end
    @testset "3-operand" begin
        # AT&T
        @test hilight_x86("\tvaddsd\t(%rax), %xmm0, %xmm0") ==
            "\t$(I)vaddsd$(XI)\t$P$(V)%rax$(XV)$XP$COM " *
            "$(V)%xmm0$(XV)$COM $(V)%xmm0$(XV)\n"

        # Intel
        @test hilight_x86("\tvaddsd\txmm0, xmm0, qword ptr [rax]") ==
            "\t$(I)vaddsd$(XI)\t$(V)xmm0$(XV)$COM $(V)xmm0$(XV)$COM " *
            "$(K)qword$(XK) $(K)ptr$(XK) $S$(V)rax$(XV)$XS\n"
    end
    @testset "4-operand" begin
        # AT&T
        @test hilight_x86("\tvroundsd\t\$4, %xmm1, %xmm1, %xmm1") ==
            "\t$(I)vroundsd$(XI)\t$(N)\$4$(XN)$COM " *
            "$(V)%xmm1$(XV)$COM $(V)%xmm1$(XV)$COM $(V)%xmm1$(XV)\n"

        # Intel
        @test hilight_x86("\tvroundsd\txmm1, xmm1, xmm1, 4") ==
            "\t$(I)vroundsd$(XI)\t" *
            "$(V)xmm1$(XV)$COM $(V)xmm1$(XV)$COM $(V)xmm1$(XV)$COM $(N)4$(XN)\n"
    end
    @testset "AVX-512" begin
        # AT&T
        @test hilight_x86("\tvmovaps\t(%eax), %zmm0") ==
            "\t$(I)vmovaps$(XI)\t$P$(V)%eax$(XV)$XP$COM $(V)%zmm0$(XV)\n"

        @test hilight_x86("\tvpaddd\t%zmm3, %zmm1, %zmm1 {%k1}") ==
            "\t$(I)vpaddd$(XI)\t$(V)%zmm3$(XV)$COM $(V)%zmm1$(XV)$COM " *
            "$(V)%zmm1$(XV) $U$(V)%k1$(XV)$XU\n"

        @test hilight_x86("\tvdivpd\t%zmm3, %zmm1, %zmm0 {%k1} {z}") ==
            "\t$(I)vdivpd$(XI)\t$(V)%zmm3$(XV)$COM $(V)%zmm1$(XV)$COM " *
            "$(V)%zmm0$(XV) $U$(V)%k1$(XV)$XU $U$(K)z$(XK)$XU\n"

        @test hilight_x86("\tvdivps\t(%ebx){1to16}, %zmm5, %zmm4") ==
            "\t$(I)vdivps$(XI)\t$P$(V)%ebx$(XV)$XP$U$(K)1to16$(XK)$XU$COM " *
            "$(V)%zmm5$(XV)$COM $(V)%zmm4$(XV)\n"

        @test hilight_x86("\tvcvtsd2si\t{rn-sae}, %xmm0, %eax") ==
            "\t$(I)vcvtsd2si$(XI)\t$U$(K)rn-sae$(XK)$XU$COM " *
            "$(V)%xmm0$(XV)$COM $(V)%eax$(XV)\n"

        # Intel
        @test hilight_x86("\tvmovaps\tzmm0, zmmword ptr [eax]") ==
            "\t$(I)vmovaps$(XI)\t$(V)zmm0$(XV)$COM " *
            "$(K)zmmword$(XK) $(K)ptr$(XK) $S$(V)eax$(XV)$XS\n"

        @test hilight_x86("\tvpaddd\tzmm1 {k1}, zmm1, zmm3") ==
            "\t$(I)vpaddd$(XI)\t$(V)zmm1$(XV) $U$(V)k1$(XV)$XU$COM " *
            "$(V)zmm1$(XV)$COM $(V)zmm3$(XV)\n"

        @test hilight_x86("\tvdivpd\tzmm0 {k1} {z}, zmm1, zmm3") ==
            "\t$(I)vdivpd$(XI)\t$(V)zmm0$(XV) $U$(V)k1$(XV)$XU $U$(K)z$(XK)$XU$COM " *
            "$(V)zmm1$(XV)$COM $(V)zmm3$(XV)\n"

        @test hilight_x86("\tvdivps\tzmm4, zmm5, dword ptr [ebx]{1to16}") ==
            "\t$(I)vdivps$(XI)\t$(V)zmm4$(XV)$COM $(V)zmm5$(XV)$COM " *
            "$(K)dword$(XK) $(K)ptr$(XK) $S$(V)ebx$(XV)$XS$U$(K)1to16$(XK)$XU\n"

        @test hilight_x86("\tvcvtsd2si\teax, xmm0$(XV), {rn-sae}") ==
            "\t$(I)vcvtsd2si$(XI)\t$(V)eax$(XV)$COM " *
            "$(V)xmm0$(XV)$COM $U$(K)rn-sae$(XK)$XU\n"
    end
end

@testset "ARM ASM" begin
    @testset "comment" begin
        @test hilight_arm("; comment ; // # ") == "$(C); comment ; // # $(XC)\n"
    end
    @testset "label" begin
        @test hilight_arm("L45:") == "$(L)L45:$(XL)\n"
    end
    @testset "directive" begin
        @test hilight_arm("\t.text") == "\t$(D).text$(XD)\n"
    end

    @testset "0-operand" begin
        @test hilight_arm("\tret") == "\t$(I)ret$(XI)\n"
    end
    @testset "1-operand" begin
        @test hilight_arm("\tbl\t0x12") == "\t$(I)bl$(XI)\t$(N)0x12$(XN)\n"

        @test hilight_arm("\tb\tL345") == "\t$(I)b$(XI)\t$(L)L345$(XL)\n"

        @test hilight_arm("\tb.gt\tL67") == "\t$(I)b.gt$(XI)\t$(L)L67$(XL)\n"

        @test hilight_arm("\tpop\t{r11, pc}") ==
            "\t$(I)pop$(XI)\t$U$(V)r11$(XV)$COM $(V)pc$(XV)$XU\n"
    end
    @testset "2-operand" begin
        @test hilight_arm("\tcmp\tx10, #2047\t// =2047") ==
            "\t$(I)cmp$(XI)\t$(V)x10$(XV)$COM $(N)#2047$(XN)\t$(C)// =2047$(XC)\n"

        @test hilight_arm("\tldr\td1, [x10]") ==
            "\t$(I)ldr$(XI)\t$(V)d1$(XV)$COM $S$(V)x10$(XV)$XS\n"

        @test hilight_arm("\tstr\tx30, [sp, #-16]!") ==
            "\t$(I)str$(XI)\t$(V)x30$(XV)$COM " *
            "$S$(V)sp$(XV)$COM $(N)#-16$(XN)$XS$(K)!$(XK)\n"

        @test hilight_arm("\tmov\tv0.16b, v1.16b") ==
            "\t$(I)mov$(XI)\t$(V)v0.16b$(XV)$COM $(V)v1.16b$(XV)\n"
    end
    @testset "3-operand" begin
        @test hilight_arm("\tfmul\td2, d0, d2") ==
            "\t$(I)fmul$(XI)\t$(V)d2$(XV)$COM $(V)d0$(XV)$COM $(V)d2$(XV)\n"

        @test hilight_arm("\tmovk\tx10, #65535, lsl #32") ==
            "\t$(I)movk$(XI)\t$(V)x10$COM $(N)#65535$(XN)$COM $(K)lsl$(XK) $(N)#32$(XN)\n"

        @test hilight_arm("\tcneg\tx8, x8, ge") ==
            "\t$(I)cneg$(XI)\t$(V)x8$(XV)$COM $(V)x8$(XV)$COM $(K)ge$(XK)\n"
    end
    @testset "4-operand" begin
        @test hilight_arm("\tadd\tx8, x9, x8, lsl #52") ==
            "\t$(I)add$(XI)\t$(V)x8$(XV)$COM $(V)x9$(XV)$COM $(V)x8$(XV)$COM " *
            "$(K)lsl$(XK) $(N)#52$(XN)\n"

        @test hilight_arm("\tfcsel\td1, d0, d1, eq") ==
            "\t$(I)fcsel$(XI)\t" *
            "$(V)d1$(XV)$COM $(V)d0$(XV)$COM $(V)d1$(XV)$COM $(K)eq$(XK)\n"
    end
    @testset "NEON" begin
        hilight_arm("\tvmul.f32\tq8, q9, q8") ==
            "\t$(I)vmul.f32$(XI)\t$(V)q8$(XV)$COM $(V)q9$(XV)$COM $(V)q8$(XV)\n"
        hilight_arm("\tvcvt.s32.f64\ts2, d20") ==
            "\t$(I)vcvt.s32.f64$(XI)\t$(V)s2$(XV)$COM $(V)d20$(XV)\n"
        hilight_arm("\tvld1.32\t{d18, d19}, [r1]") ==
            "\t$(I)vld1.32$(XI)\t$U$(V)d18$(XV)$COM $(V)d19$(XV)$XU$COM $S$(V)r1$(XV)$XS\n"
    end
    @testset "SVE" begin
        hilight_arm("\tld1d\tz1.d, p0/z, [x0, x4, lsl #3]") ==
            "\t$(I)ld1d$(XI)\t$(V)z1.d$(XV)$COM " *
            "$(V)p0$(XV)$(K)/z$(XK)$COM " *
            "$S$(V)x0$(XV)$COM $(V)x4$(XV)$COM $(K)lsl$(XK) $(N)#3$(XN)$XS\n"
        hilight_arm("\tb.first\tL123") == "\t$(I)b.first$(XI)\t$(L)L123$(XL)"
    end
end

max_stars_repo_name: tomerarnon/julia

max_issues_repo_name: tomerarnon/julia

max_forks_repo_name: jishnub/julia


text:
# z[1] = xDot
# z[2] = yDot
# z[3] = psiDot
# z[4] = ePsi
# z[5] = eY
# z[6] = s

function simKinModel(z::Array{Float64},u::Array{Float64},dt::Float64,coeff::Array{Float64},modelParams::ModelParams)

    zNext::Array{Float64}
    L_a = modelParams.l_A
    L_b = modelParams.l_B
    c = ([z[1]^8 z[1]^7 z[1]^6 z[1]^5 z[1]^4 z[1]^3 z[1]^2 z[1] 1]*coeff)[1]        # Polynomial

    bta = atan(L_a/(L_a+L_b)*tan(u[2]+abs(u[2])*u[2]))
    dsdt = z[4]*cos(z[3]+bta)/(1-z[2]*c)

    zNext = copy(z)
    zNext[1] = z[1] + dt*dsdt                               # s
    zNext[2] = z[2] + dt*z[4] * sin(z[3] + bta)             # eY
    zNext[3] = z[3] + dt*(z[4]/L_a*sin(bta)-dsdt*c)         # ePsi
    zNext[4] = z[4] + dt*(u[1] - modelParams.c_f*z[4])      # v

    return zNext
end

function simDynModel_exact(z::Array{Float64},u::Array{Float64},dt::Float64,coeff::Array{Float64},modelParams::ModelParams)
    # This function uses smaller steps to achieve higher fidelity than we would achieve using longer timesteps
    z_final = copy(z)
    u[1] = min(u[1],3)
    u[1] = max(u[1],-3)
    u[2] = min(u[2],pi/6)
    u[2] = max(u[2],-pi/6)
    dtn = dt/100
    for i=1:100
        z_final = simDynModel(z_final,u,dtn,coeff,modelParams)
    end
    return z_final
end

function simDynModel(z::Array{Float64},u::Array{Float64},dt::Float64,coeff::Array{Float64},modelParams::ModelParams)

    zNext::Array{Float64}
    L_f = modelParams.l_A
    L_r = modelParams.l_B
    c0  = modelParams.c0
    m   = modelParams.m
    I_z = modelParams.I_z
    c_f = modelParams.c_f

    a_F = 0
    a_R = 0
    if abs(z[1]) >= 0.1
        a_F     = atan((z[2] + L_f*z[3])/abs(z[1])) - z[8]
        a_R     = atan((z[2] - L_r*z[3])/abs(z[1]))
    end
    if max(abs(a_F),abs(a_R))>30/180*pi
        warn("Large tire angles: a_F = $a_F, a_R = $a_R, xDot = $(z[1]), d_F = $(z[8])")
    end

    FyF = -pacejka(a_F)
    FyR = -pacejka(a_R)

    c = ([z[1]^8 z[2]^7 z[3]^6 z[4]^5 z[5]^4 z[6]^3 z[7]^2 z[8] 1]*coeff)[1]                        # Polynomial for curvature

    dsdt = (z[1]*cos(z[4]) - z[2]*sin(z[4]))/(1-z[5]*c)

    zNext = copy(z)
    zNext[1] = z[1] + dt * (z[7] + z[2]*z[3] - c_f*z[1])                    # xDot
    zNext[2] = z[2] + dt * (2/m*(FyF*cos(z[8]) + FyR) - z[3]*z[1])          # yDot
    zNext[3] = z[3] + dt * (2/I_z*(L_f*FyF - L_r*FyR))                      # psiDot
    zNext[4] = z[4] + dt * (z[3]-dsdt*c)                                    # ePsi
    zNext[5] = z[5] + dt * (z[1]*sin(z[4]) + z[2]*cos(z[4]))                # eY
    zNext[6] = z[6] + dt * dsdt                                             # s
    zNext[7] = z[7] + dt * (u[1] - z[7]) * 100                              # a
    zNext[8] = z[8] + dt * (u[2] - z[8]) * 10                               # d_f

    return zNext
end

function pacejka(a)
    B = 1.0             # This value determines the steepness of the curve
    C = 1.25
    mu = 0.8            # Friction coefficient (responsible for maximum lateral tire force)
    m = 1.98
    g = 9.81
    D = mu * m * g/2
    C_alpha_f = D*sin(C*atan(B*a))
    return C_alpha_f
end

function simDynModel_exact_xy(z::Array{Float64},u::Array{Float64},dt::Float64,modelParams::ModelParams)
    dtn = dt/10
    t = 0:dtn:dt
    z_final = copy(z)
    ang = zeros(2)
    for i=1:length(t)-1
        z_final, ang = simDynModel_xy(z_final,u,dtn,modelParams)
    end
    return z_final, ang
end

function simDynModel_xy(z::Array{Float64},u::Array{Float64},dt::Float64,modelParams::ModelParams)

    zNext::Array{Float64}
    L_f = modelParams.l_A
    L_r = modelParams.l_B
    m   = modelParams.m
    I_z = modelParams.I_z

    a_F = 0.0
    a_R = 0.0
    if abs(z[3]) > 0.2
        a_F     = atan((z[4] + L_f*z[6])/abs(z[3])) - z[8]
        a_R     = atan((z[4] - L_r*z[6])/abs(z[3]))
    end

    FyF = -pacejka(a_F)
    FyR = -pacejka(a_R)

    if abs(a_F) > 30/180*pi || abs(a_R) > 30/180*pi
        warn("Large slip angles in simulation: a_F = $a_F, a_R = $a_R")
    end

    zNext = copy(z)
    # compute next state
    zNext[1]        = zNext[1]       + dt * (cos(z[5])*z[3] - sin(z[5])*z[4])               # x
    zNext[2]        = zNext[2]       + dt * (sin(z[5])*z[3] + cos(z[5])*z[4])               # y
    zNext[3]        = zNext[3]       + dt * (z[7] + z[4]*z[6] - 0.5*z[3])                   # v_x
    zNext[4]        = zNext[4]       + dt * (1/m*(FyF*cos(z[8]) + FyR) - z[6]*z[3])         # v_y
    zNext[5]        = zNext[5]       + dt * (z[6])                                          # psi
    zNext[6]        = zNext[6]       + dt * (1/I_z*(L_f*FyF - L_r*FyR))                     # psiDot
    zNext[7]        = zNext[7]       + dt * (u[1]-z[7])*100                                 # a
    zNext[8]        = zNext[8]       + dt * (u[2]-z[8])*100                                 # d_f

    zNext[3] = max(0,zNext[3])              # limit speed to positive values (BARC specific)
    return zNext, [a_F a_R]
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
using ArgParse2
using Artifacts
using DataFrames
using DataStructures
using BioSequences
using CodecZlib
using BioAlignments
using DataFramesMeta
using Statistics
using st16SseqJuliaTools
using XAM
using BioTools.BLAST
using CSV
using ProgressMeter

function main()
    parser = ArgumentParser(description = "Integrate clustering counts and assignment")
    add_argument!(parser, "--assembly-clusters", "-a", help = "Assembly clustering info",type = String)
    add_argument!(parser, "--deuplicated-clusters", "-d", help = "Assembly clustering info",type = String)
    add_argument!(parser, "--kraken", "-k", help = "Kraken output of pseudocontigs ",type = String)
    add_argument!(parser, "--kraken2", "-f", help = "Kraken output of fused contigs",type = String)
    add_argument!(parser, "--label", "-l", help = "Label for outout",type = String)
    args = parse_args(parser)
    taxdb=artifact"taxdump"

    #Parse data of fused clusters assembly
    fused_taxids = Array{String,1}()
    data_on_fused = DataFrame(Cluster_ID = Array{String,1}(),Nuber_of_contigs=Array{Int64,1}())
    for l in readlines(args.kraken2)
        parts = split(l,"\t")
        tid = parts[3]
        tid = split(split(split(tid,"taxid")[2])[1],")")[1]
        cid = split(parts[2],"@")[1]
        ncont = split(split(parts[2],"@")[2],":")[1]
        push!(data_on_fused,[cid,parse(Int64,ncont)])
        push!(fused_taxids,tid)
    end
    data_on_fused.Genus_based_on_fused_contigs = taxid_to_genus(fused_taxids,taxdb)
    data_on_fused.Species_based_on_fused_contigs = taxid_to_species(fused_taxids,taxdb)

    #Integrate deduplicated and clsutered for assembly clusters
    clusters_cf= args.assembly_clusters
    clusters_df= args.deuplicated_clusters
    krakenout= args.kraken
    stem=args.label

    cluster_4_workD,r2cD,c2rD = parse_mapping(clusters_df)
    cluster_4_workC,r2cC,c2rC = parse_mapping(clusters_cf)

    ct = 0
    clusterids = collect(keys(c2rC))
    cl_vs_abund = DataFrame(cluster=Array{String,1}(),cnt=Array{Float64,1}(),size=Array{Int64,1}())
    println(stderr, "Assigning counts")
    @showprogress for c in clusterids
        ct +=length(c2rC[c])
        c_reads=c2rC[c]
        fr = 0
        for cd in keys(c2rD)
            d_reads = c2rD[cd]
            totc = length(d_reads)
            overlapc = length(intersect(d_reads,c_reads))
            fr += overlapc/totc
        end
        push!(cl_vs_abund,[c,fr,length(c_reads)])
    end
    taxids=Array{String,1}()
    clid_spec=Dict{String,String}()
    clids=Array{String,1}()
    println(stderr,"Checking species and genus ssignments...")
    for l in readlines(krakenout)
        id = split(split(l)[2],"_")[1]
        push!(clids,id)
        taxid = split(split(split(l,"taxid")[2])[1],")")[1]
        push!(taxids,taxid)
    end
    species=taxid_to_species(taxids,taxdb)
    genus=taxid_to_genus(taxids,taxdb)
    println(stderr,"...done")
    clid_spec=Dict{String,String}(zip(clids,species))
    clid_genus=Dict{String,String}(zip(clids,genus))
    cl_vs_abund = sort(cl_vs_abund,[:size],rev=true)
    spassignments = Array{String,1}()
    geassignments = Array{String,1}()
    assigned = collect(keys(clid_spec))
    Cluster_ID = Array{String,1}()
    @showprogress for i in  1:nrow(cl_vs_abund)
        is = string(i)
        push!(Cluster_ID,is)
        if is in assigned
            push!(spassignments,clid_spec[is])
            push!(geassignments,clid_genus[is])
        else
            push!(spassignments,"NA")
            push!(geassignments,"NA")
        end
    end
    cl_vs_abund.Species_based_on_pseudocontigs = spassignments
    cl_vs_abund.Genus_based_on_pseudocontigs = geassignments
    cl_vs_abund.Fraction = cl_vs_abund.cnt ./ sum(cl_vs_abund.cnt)
    cl_vs_abund.Sample = fill(split(stem,"/")[end],length(spassignments))
    cl_vs_abund.Cluster_ID = Cluster_ID
    out = join(cl_vs_abund, data_on_fused, on=:Cluster_ID, kind = :outer)
    CSV.write(args.label*"_info_on_pseudocontig_and_contigs.csv",out)
    end

main()


text:
function discrete_phs(N,a,b)
	Nx = N   	#    energy variables: x
	Ne = Nx+1   # co-energy variables: e

	ei,we = lgwt(Ne,a,b)    # discretization of co-energy variables
	#ei = collect(linspace(a,b,Ne,))
	if Nx > 1
		xi,wx = lgwt(Nx,a,b)
		#xi = collect(linspace(a,b,Nx))
	else
		xi,wx = (a+b)/2, (b-a);
	end
	#M = massmatrix(ei,xi,we);
	xquad, wquad = lgwt(Ne,a,b)
	M = massmatrix(ei,xi, xquad, wquad);
	D = dermatrix(ei,xi,1);
	p0 = map(i->leg_pol(a,ei,i), 1:length(ei))
	pL = map(i->leg_pol(b,ei,i), 1:length(ei))

	F = [M 0*M;
	     pL' pL'*0;
		 0*M M;
		 0*p0' p0'];
	E = -[D*0 D;
	      pL'*0 -pL';
		  D D*0;
		  p0' p0'*0];
		  J = E/F
    phJ = -[zeros(N,N) J[1:N,(N+2):end-1];
		          J[N+2:end-1,1:N] zeros(N,N)];
    phB = -[J[[1:N;N+2:2*N+1],N+1] J[[1:N;(N+2):(2*N+1)],end]];
    phD = [J[N+1,N+1] J[N+1,end];
        J[end,N+1] J[end,end]];

	if N > 1
		Q = massmatrix(xi,xi,xquad,wquad)
		#Q = sparse(diagm(wx[:],0));
		Q = blkdiag(Q,Q);
		else
		#Q = diagm([wx; wx][:]);
	end

	pflow = Collocation_points(xi, wx)
	peffort = Collocation_points(ei, we)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

function discrete_phs2(N,a,b)
	Nx = N
	Ne = Nx+2
	ei,wi = lgwt(Ne,a,b)
	if Nx > 2
		xi,wx = lgwt(Nx,a,b)
	else
		xi,wx = (a+b)/2, (b-a);
	end
	M = massmatrix(ei,xi,wi);
	D = dermatrix(ei,xi,2);
	leg_pol_der = i -> (algo_diff(x->leg_pol(x,ei,i), 1))
	p0 = map(i->leg_pol(a,ei,i), 1:length(ei))
	p0d = map(i->leg_pol_der(i)(a), 1:length(ei))
	pL = map(i->leg_pol(b,ei,i), 1:length(ei))
	pLd = map(i->leg_pol_der(i)(b), 1:length(ei))

	F = [M 0*M;
	     pL' pL'*0;
	     pLd' pL'*0
		 0*M M;
	     0*pLd' p0d'
		 0*p0' p0'];
	E = [D*0 D;
	      pL'*0 -pLd';
	      pL'*0  pL';
	       -D D*0;
	       p0' p0'*0;
	      -p0d' p0'*0];
	J = E/F
	neworder = [1:N; (1:N) .+ (N+2); N+1:N+2; 2*N+3:2*N+4];
    Jnew = J[neworder,neworder]
    phJ = - Jnew[1:2*N,1:2*N];
    phB = - Jnew[1:2*N, 2*N .+ (1:4)];
    phD = Jnew[2*N .+ (1:4), 2*N .+ (1:4)];

	if N > 1
		#Q = (diagm(wx[:],0));
		Q = (diagm(0 => wx[:]));
		Q = blkdiag(Q,Q);
		else
		Q = diagm([wx; wx][:]);
	end

	pflow = Collocation_points(xi, wx)
	peffort = Collocation_points(ei, wi)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

function discrete_phs2_distports(N,K,a,b,ad,bd)
	Nx = N
	Ne = Nx+2
	zei,we = lgwt(Ne,a,b)
	zvi,wv = lgwt(K,ad,bd)
	if Nx > 2
		zxi,wx = lgwt(Nx,a,b)
	else
		zxi,wx = (a+b)/2, (b-a);
	end
	M = massmatrix(zei,zxi,we);

	leg_pol_der2 = (z,zi,i) -> algo_diff(x->leg_pol(x,zi,i),2)(z);
	# M_theta_phizz can be exactly obtained by quadrature if K >= N-3:
	if K < N-3
		warn("K < N-3: The distributed ports matrix is not exactly computed by quadrature")
	end
	M_theta_phizz = massmatrix(zvi,zxi,wv, leg_pol_der2);
	M_phi = diagm(0 => wx[:]);
	B_bar = M_phi \ M_theta_phizz;
	D = dermatrix(zei,zxi,2);
	leg_pol_der = i -> (algo_diff(x->leg_pol(x,zei,i), 1))
	p0 = map(i->leg_pol(a,zei,i), 1:length(zei))
	p0d = map(i->leg_pol_der(i)(a), 1:length(zei))
	pL = map(i->leg_pol(b,zei,i), 1:length(zei))
	pLd = map(i->leg_pol_der(i)(b), 1:length(zei))

	F = [M 0*M;
	     p0' pL'*0;
	     p0d' pL'*0
		 0*M M;
	     0*pLd' pLd'
		 0*p0' pL'];
	F = blkdiag(F,Matrix(1I, K, K));
	E = [D*0 D;
	      pL'*0 p0d';
	      pL'*0  -p0';
	       -D D*0;
	       -pL' p0'*0;
	      pLd' p0'*0];
	E = blkdiag(E, zeros(K,K));
	E[1:size(D,1),end-K+1:end] = - B_bar;
    E[end-K+1:end,1:size(D,2)] =  B_bar' * M;
	J = E/F
	neworder = [1:N; (1:N) .+ (N+2); N+1:N+2; 2*N+3:2*N+4; (2*N+4) .+ (1:K)];
    Jnew = J[neworder,neworder]
    phJ = - Jnew[1:2*N,1:2*N];
    phB = - Jnew[1:2*N, 2*N .+ (1:4)];
	phBdistributed = - Jnew[1:2*N, (2*N + 4) .+ (1:K)];
    phD = Jnew[2*N .+ (1:4), 2*N .+ (1:4)];

	if N > 1
		Q = (diagm(0 => wx[:]));
		Q = blkdiag(Q,Q);
		else
		Q = diagm([wx; wx][:]);
	end

	pflow = Collocation_points(zxi, wx)
	peffort = Collocation_points(zei, we)
	pdistributed = Collocation_points(zvi, wv)

	p = Phs(phJ, phB, phD, Q)
	p.Bd = phBdistributed
	#p.Q = Q
	p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
module PaulsRechner

using Dates

# package code goes here
export
        Rechentrainer,
        Emma

include("Rechner.jl")
include("EmmasRechner.jl")

end # module


text:
using Flux, Random, Parameters, Zygote, LinearAlgebra

# Power iteration algorithm for computing the spectral norm
function power_iteration!(W, u, n_iterations, rng = Random.GLOBAL_RNG)
    v = nothing
    for i=1:n_iterations
        Wu = W' * u
        v = Wu ./ (norm(Wu) + eps(Float32))
        Wv = W * v
        u .= Wv ./ (norm(Wv) + eps(Float32))
    end
    u, v
end

# Compute the maximum singular value
msv(u, v, W) = u' * W * v

## Dense Layer with Spectral Normalization
struct DenseSN{F,S<:AbstractArray,T<:AbstractArray, I<:Int, VV<:AbstractArray}
  W::S
  b::T
  σ::F
  n_iterations::I # Number of power iterations for computing max singular value
  u::VV # Left vector for power iteration
end

function DenseSN(in::Integer, out::Integer, σ = identity;
               initW = Flux.glorot_uniform, initb = Flux.zeros, n_iterations = 1)
  return DenseSN(initW(out, in), initb(out), σ, n_iterations, randn(Float32, out, 1))
end

Flux.@functor DenseSN

Flux.trainable(a::DenseSN) = (a.W, a.b)

function (a::DenseSN)(x::AbstractArray)
  W, b, σ = a.W, a.b, a.σ
  u, v = Zygote.ignore(() -> power_iteration!(W, a.u, a.n_iterations))
  σ.((W ./ msv(u, v, W))*x .+ b)
end

function Base.show(io::IO, l::DenseSN)
  print(io, "DenseSN(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

(a::DenseSN{<:Any,W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  invoke(a, Tuple{AbstractArray}, x)

(a::DenseSN{<:Any,W})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  a(T.(x))

function outdims(l::DenseSN, isize)
    first(isize) == size(l.W, 2) || throw(DimensionMismatch("input size should equal to ($(size(l.W, 2)),), got $isize"))
    return (size(l.W, 1),)
end

## Convluational layer with Spectral Normalization
struct ConvSN{N,M,F,A,V, I<:Int, VV<:AbstractArray}
  σ::F
  weight::A
  bias::V
  stride::NTuple{N,Int}
  pad::NTuple{M,Int}
  dilation::NTuple{N,Int}
  n_iterations::I # Number of power iterations for computing max singular value
  u::VV # Left vector for power iteration
end

to2D(W) = reshape(W, :, size(W, ndims(W)))

function ConvSN(w::AbstractArray{T,N}, b::Union{Flux.Zeros, AbstractVector{T}}, σ = identity;
              stride = 1, pad = 0, dilation = 1, n_iterations = 1) where {T,N}
  stride = Flux.expand(Val(N-2), stride)
  dilation = Flux.expand(Val(N-2), dilation)
  pad = Flux.calc_padding(Conv, pad, size(w)[1:N-2], dilation, stride)
  u = randn(Float32, size(to2D(w), 1), 1)
  return ConvSN(σ, w, b, stride, pad, dilation, n_iterations, u)
end

function ConvSN(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
            init = Flux.glorot_uniform,  stride = 1, pad = 0, dilation = 1,
            weight = Flux.convfilter(k, ch, init = init), bias = Flux.zeros(ch[2]), n_iterations = 1) where N
  ConvSN(weight, bias, σ, stride = stride, pad = pad, dilation = dilation, n_iterations = n_iterations)
end

Flux.@functor ConvSN

Flux.trainable(a::ConvSN) = (a.weight, a.bias)

function (c::ConvSN)(x::AbstractArray)
  σ, b = c.σ, reshape(c.bias, ntuple(_->1, length(c.stride))..., :, 1)
  cdims = DenseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  u, v = Zygote.ignore(() -> power_iteration!(to2D(c.weight), c.u, c.n_iterations))
  σ.(conv(x, c.weight ./ msv(u, v, to2D(c.weight)), cdims) .+ b)
end

function Base.show(io::IO, l::ConvSN)
  print(io, "ConvSN(", size(l.weight)[1:ndims(l.weight)-2])
  print(io, ", ", size(l.weight, ndims(l.weight)-1), "=>", size(l.weight, ndims(l.weight)))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

(a::ConvSN{<:Any,<:Any,W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  invoke(a, Tuple{AbstractArray}, x)

(a::ConvSN{<:Any,<:Any,W})(x::AbstractArray{<:Real}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
  a(T.(x))

outdims(l::ConvSN, isize) =
  output_size(DenseConvDims(_paddims(isize, size(l.weight)), size(l.weight); stride = l.stride, padding = l.pad, dilation = l.dilation))

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
using TestImages
using Test

@testset "contrast with kwargs" begin
    # test for the fix in #228
    # Uses 2D images to verify the `axes` and `flipy` kwargs work
    # with a passed CLim argument
    mri = testimage("mri")
    vmin = Float32(minimum(mri))
    vmax = Float32(maximum(mri))
    clim = ImageView.CLim(vmin, vmax)

    ImageView.imshow(mri, clim; axes=(1, 2))
    @test true
    ImageView.imshow(mri, clim; axes=(3, 1), flipy=true)
    @test true
    ImageView.imshow(mri, clim; axes=(3, 2), flipy=true)
    @test true
end


text:
export compute_compton_photon_flux, compute_compton_photons_per_bin

function compute_compton_photons_per_bin(energy_range, parameters)
    photon_number_per_bin = donthcomp(energy_range, parameters)
    return photon_number_per_bin
end


text:
using ImPlot
using Test

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


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    WhiteningAlgorithm

Abstract type for sample whitening algorithms.
"""
abstract type WhiteningAlgorithm end
export WhiteningAlgorithm

"""
NoWhitening

No-op whitening transformation, leaves samples unchanged.
"""
struct NoWhitening <: WhiteningAlgorithm end
export NoWhitening

"""
    CholeskyWhitening

Whitening transformation based on a Cholesky transformation of the empirical
sample covariance matrix.
"""
struct CholeskyWhitening <: WhiteningAlgorithm end
export CholeskyWhitening

"""
    CholeskyPartialWhitening

Whitening transformation based on a Cholesky transformation of the empirical
sample covariance matrix.

Only transforms dimensions (degrees of freedom) for which the marginalized
distribution aymptotically approaches zero in the positive and negative
direction.
"""
struct CholeskyPartialWhitening <: WhiteningAlgorithm end
export CholeskyPartialWhitening

"""
    StatisticalWhitening

Whitening transformation based statistical whitening.
"""
struct StatisticalWhitening <: WhiteningAlgorithm end
export StatisticalWhitening

# ToDo: Move whitening code from AHMI here, add new public API
# function `bat_whiten`. Remove `_amhi_whitening_func` and
# change AHMI to use `bat_whiten` directly.

_amhi_whitening_func(algorithm::NoWhitening) = no_whitening
_amhi_whitening_func(algorithm::CholeskyWhitening) = cholesky_whitening!
_amhi_whitening_func(algorithm::CholeskyPartialWhitening) = cholesky_partial_whitening!
_amhi_whitening_func(algorithm::StatisticalWhitening) = statistical_whitening!


text:
"""
Overview: functions used to manage sparse matrices
in the code.
"""

"""
> init_symCSC(mesh::Grid, ϕ::Field, n2n::NodalConnectivity, naux::Int64, tp::String)
It initializes a symmetric CSC matrix from a field. It eventually add auxiliary equations.
- mesh : Grid
- ϕ : field of interest
- nn : number of queries
- n2n : node-to-node connectivity
- naux : number of auxiliary equations
"""
function init_symCSC(mesh::Grid, ϕ::Field, n2n::NodalConnectivity,
                     naux::Int64, tp::String)
  #
  neq = ϕ.neq + naux
  #
  dofs1 = zeros(Int64,ϕ.cpn)
  dofs2 = zeros(Int64,ϕ.cpn)
  ent = zeros(Int64,neq)
  # count nonzero entries
  nnz = 0
  for n = 1:mesh.nn
    dofs!(ϕ,1,[n],dofs1)
    conn = get_node_connectivity(n2n,n)
    lreg = size(conn)[1]
    for i = 1:ϕ.cpn
      if (dofs1[i]>0)
        for j = 1:ϕ.cpn
          if (dofs1[j]>0 && dofs1[i]>=dofs1[j])
            ent[dofs1[j]] += 1
            nnz += 1
          end
        end
        for j = 1:lreg
          dofs!(ϕ,1,[conn[j]],dofs2)
          for k = 1:ϕ.cpn
            if (dofs2[k]>0 && dofs1[i]>=dofs2[k])
              ent[dofs2[k]] += 1
              nnz += 1
            end
          end
        end
      end
    end
  end
  # add auxiliary
  for i = 1:ϕ.neq
    ent[i] += naux
  end
  for i = 1:naux
    nnz += ϕ.neq+i
    ent[ϕ.neq+i] = naux+1-i
  end
  # allocate arrays
  jcol = zeros(Int64,neq+1)
  irow = zeros(Int64,nnz)
  if (tp=="r")
    val = zeros(Float64,nnz)
  elseif (tp=="c")
    val = zeros(ComplexF64,nnz)
  end
  # fill matrices
  j = 1
  jcol[j] = 1
  for i = 1:neq
    j += ent[i]
    jcol[i+1] = j
  end
  #
  fill!(ent,0.0)
  for n = 1:mesh.nn
    dofs!(ϕ,1,[n],dofs1)
    conn = get_node_connectivity(n2n,n)
    lreg = size(conn)[1]
    for i = 1:ϕ.cpn
      idof = dofs1[i]
      if (idof>0)
        for j = 1:ϕ.cpn
          jdof = dofs1[j]
          if (jdof>0 && dofs1[i]>=dofs1[j])
            iv = ent[jdof] + jcol[jdof]
            ent[jdof] += 1
            irow[iv] = idof
          end
        end
        for j = 1:lreg
          dofs!(ϕ,1,[conn[j]],dofs2)
          for k = 1:ϕ.cpn
            jdof = dofs2[k]
            if (jdof>0 && dofs1[i]>=dofs2[k])
              iv = ent[jdof] + jcol[jdof]
              ent[jdof] += 1
              irow[iv] = idof
            end
          end
        end
      end
    end
  end
  # auxiliary
  for i = 1:ϕ.neq
    iv = ent[i] + jcol[i] - 1
    ent[i] += naux
    for j = 1:naux
      irow[iv+j] = ϕ.neq+j
    end
  end
  for i = 1:naux
    iv = ent[ϕ.neq+i] + jcol[ϕ.neq+i] - 1
    ent[i] += naux+1-i
    for j = i:naux
      irow[iv+j-i+1] = ϕ.neq+j
    end
  end
  #
  for i = 1:neq
    iv1 = jcol[i]
    iv2 = jcol[i+1]-1
    sort!(irow[iv1:iv2])
  end
  #
  return SparseMatrixCSC(neq,neq,jcol,irow,val)
  #
end


text:
pages = [
        "DiffEqBayes.jl: Bayesian Parameter Estimation for Differential Equations" => "index.md",
        "Methods" => "methods.md",
        "Examples" => "examples.md"
    ]


text:
# Extension operator: joConvolve
# Author : Rajiv Kumar
## helper module
module joConvolve_etc
    using JOLI: jo_convert
    using FFTW
    ## forward 1D convolution
# function Convolve1D_fwd{VT<:Number}(x::AbstractVector{VT}, fKernel::AbstractVector,idx::AbstractVector, k::Integer,m::Integer,cflag::Bool,rdt=DataType)
function Convolve1D_fwd(x::AbstractVector, fKernel::AbstractVector,idx::AbstractVector, k::Integer,m::Integer,cflag::Bool,rdt=DataType)
        #fx = fft([full(x);zeros(VT,k-1,1)[:]]);
        #fx = fft([full(x);zeros(k-1,1)[:]]);
         fx = fft([x;zeros(k-1,1)[:]]);
         y = ifft(fKernel.*fx);
         y = y[idx];
         if ~cflag && isreal(y)
             y = real(y);
         end
        return jo_convert(rdt,y,false)
 end
  ## transpose 1D convolution
 #function Convolve1D_tran{VT<:Number}(x::AbstractVector{VT}, fKernel::AbstractVector,idx::AbstractVector, k::Integer,m::Integer,cflag::Bool,rdt=DataType)
 function Convolve1D_tran(x::AbstractVector, fKernel::AbstractVector,idx::AbstractVector, k::Integer,m::Integer,cflag::Bool,rdt=DataType)
        #z      = zeros(VT,m+k-1,1)[:];
        z      = zeros(m+k-1,1)[:];
        #z[idx] = full(x);
        z[idx] = x;
        y      = ifft(conj(fKernel).*fft(z));
        y      = y[1:m];
         if ~cflag && isreal(y)
             y = real(y);
         end
        return jo_convert(rdt,y,false)
 end
end
using .joConvolve_etc

export joConvolve

#function joConvolve{KDT<:Number}(m::Integer, n::Integer, kernel::AbstractVector{KDT}; DDT::DataType=joFloat,RDT::DataType=DDT)
function joConvolve(m::Integer, n::Integer, kernel::AbstractVector; DDT::DataType=joFloat,RDT::DataType=DDT)
                # convolution loop start
                offset = 1
                k      = length(kernel);
                cflag  = (eltype(kernel)<:Complex); #~isreal(kernel);
                # Shift kernel and add internal padding
                # kernel = [kernel[offset:end];zeros(KDT,m-1,1)[:];kernel[1:offset-1]];
                kernel = [kernel[offset:end];zeros(m-1,1)[:];kernel[1:offset-1]];
                   idx = collect(1:m);

                # Precompute kernel in frequency domain
               #fKernel = fft(full(kernel));
               fKernel = fft(kernel);
               # Convolution
               return joLinearFunctionFwd(length(idx),m,
                                        v1->joConvolve_etc.Convolve1D_fwd(v1,fKernel,idx, k,m,cflag,RDT),
                                        v2->joConvolve_etc.Convolve1D_tran(v2,fKernel,idx, k,m,cflag,DDT),
                                        v2->joConvolve_etc.Convolve1D_tran(v2,fKernel,idx, k,m,cflag,DDT),
                                        v1->joConvolve_etc.Convolve1D_fwd(v1,fKernel,idx, k,m,cflag,RDT),
                                        DDT,RDT;fMVok=false,
                                        name="joConvolve1D")
    return nothing
end


text:
__precompile__()

module IntegratedOscillatorModel
# import functions explicitly to override them
import Base.sort
import Base.max
import Base.min

# library imports
using Base.Threads
using DataStructures
using Reexport
using DifferentialEquations
using Parameters
@reexport using LaTeXStrings
@reexport using Plots
@reexport using Plots.Measures

include("helper.jl")
include("meds.jl")
include("params.jl")
include("ode_systems.jl")
include("solver_utils.jl")

# exported functions and variables that are available after using IntegratedOscillatorModel
export simulate
export loopvals
export labels
export y0
export y0_stat
export sys
#export pfk_activity
export params
export sortedvalues
export rowtocolumn
export max
export min
#export vecdiff
export vkf
#export cycle
export Med
export ExpMed

export Activa
export Dz
export Tolb
export Tg
export KCl
export Glucose

export Meds
export trysolve
export sys_const_pfk
export System

end # module


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using BenchmarkTools
using CUDA
using Oceananigans
using Benchmarks

using Oceananigans.Solvers

# Benchmark function

function benchmark_fourier_tridiagonal_poisson_solver(Arch, N, topo)
    grid = VerticallyStretchedRectilinearGrid(architecture=Arch(), topology=topo, size=(N, N, N), x=(0, 1), y=(0, 1), z_faces=collect(0:N))
    solver = FourierTridiagonalPoissonSolver(Arch(), grid)

    solve_poisson_equation!(solver) # warmup

    trial = @benchmark begin
        @sync_gpu solve_poisson_equation!($solver)
    end samples=10

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Ns = [256]
PB = (Periodic, Bounded)
Topologies = collect(Iterators.product(PB, PB, (Bounded,)))[:]

# Run and summarize benchmarks

suite = run_benchmarks(benchmark_fourier_tridiagonal_poisson_solver; Architectures, Ns, Topologies)

df = benchmarks_dataframe(suite)
sort!(df, [:Architectures, :Topologies, :Ns], by=(string, string, identity))
benchmarks_pretty_table(df, title="Fourier-tridiagonal Poisson solver benchmarks")

if GPU in Architectures
    df = gpu_speedups_suite(suite) |> speedups_dataframe
    sort!(df, [:Topologies, :Ns], by=(string, identity))
    benchmarks_pretty_table(df, title="Fourier-tridiagonal Poisson solver CPU to GPU speedup")
end

for Arch in Architectures
    suite_arch = speedups_suite(suite[@tagged Arch], base_case=(Arch, Ns[1], (Periodic, Periodic, Bounded)))
    df_arch = speedups_dataframe(suite_arch, slowdown=true)
    sort!(df_arch, [:Topologies, :Ns], by=string)
    benchmarks_pretty_table(df_arch, title="Fourier-tridiagonal Poisson solver relative performance ($Arch)")
end


text:
using Documenter

makedocs(
format   = :html,
    debug = true,
    doctestfilters = [r"Ptr{0x[0-9]+}"],
    sitename = "Documenter example",
    pages = ["index.md"],
)


text:
using Base.Test
import Ogre

# test resources

module test_eos_resources
    import Ogre

    type SimpleEOS <: Ogre.EOS
        A::Float64
        n::Float64
        C::Float64
    end

    type ComplicatedEOS <: Ogre.EOS
        A::Float64
        B::Float64
    end

    type ConstantEOS <: Ogre.EOS
        C::Float64
    end

    for eos in (SimpleEOS, ComplicatedEOS, ConstantEOS)
        Ogre.@addEOSCall eos
    end

    (eos::SimpleEOS)(P::Number) = eos.A*P^eos.n + eos.C
    (eos::SimpleEOS)(P::Number, T::Number) = eos(P)
    (eos::ComplicatedEOS)(P::Number, T::Number) = eos.A*P + eos.B*T
    (eos::ConstantEOS)(P::Number) = eos.C
    (eos::ConstantEOS)(P::Number, T::Number) = eos(P)

    Ogre.istempdependent(::Union{SimpleEOS, ConstantEOS}) = false
    Ogre.istempdependent(::ComplicatedEOS) = true
end

@testset "Basic equation of state tests" begin
    res = test_eos_resources

    eos1 = res.SimpleEOS(1.,2.,3.)
    eos2 = res.ComplicatedEOS(1.,2.)
    eos3 = res.ConstantEOS(4.)

    vs = Ogre.ValueSet
    vs1 = vs(1,2,3)
    vs2 = vs(1,2,3,4)

    @testset "Calling with physical value sets" begin
        @test eos1(vs1) == 1*3^2 + 3
        @test eos1(vs2) == 1*3^2 + 3
        @test_throws MethodError eos2(vs1)
        @test eos2(vs2) == 1*3 + 2*4
        @test eos3(vs1) == 4
        @test eos3(vs2) == 4
    end

    @testset "Handling mass piecewise EOSes" begin
        mpw = Ogre.MassPiecewiseEOS([eos1, eos2, eos3], [0, 1, 2, 3])
        mpw2 = Ogre.MassPiecewiseEOS([eos1, eos2, eos3], 3, [1/3, 1/3, 1/3])

        @test mpw(vs(0.5, 1, 2)) == eos1(2)
        @test mpw2(vs(0.5, 1, 2)) == eos1(2)
        @test mpw(vs(1.5, 1, 2, 3)) == eos2(2, 3)
        @test mpw2(vs(1.5, 1, 2, 3)) == eos2(2, 3)
        @test mpw(vs(2.5, 1, 2, 3)) == eos3(2, 3)
        @test mpw2(vs(2.5, 1, 2, 3)) == eos3(2, 3)
        @test_throws MethodError mpw(vs(1.5, 2, 3))
        @test_throws MethodError mpw2(vs(1.5, 2, 3))
    end
end


text:
int main() {
        { return true; };
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
"""
    Limbdark

Fast, analytic, and numerically stable transit light curves for stars with
arbitrary order limb darkening.

"""
module Limbdark

include("transit_structure.jl")
include("transit_poly_struct.jl")
include("integrate_lightcurve.jl")

export transit_init, transit_poly, integrate_lightcurve

end


text:
#=

 Test Cases for interval constructors from IEEE Std 1788-2015

 Copyright 2016 Oliver Heimlich

 Copying and distribution of this file, with or without modification,
 are permitted in any medium without royalty provided the copyright
 notice and this notice are preserved.  This file is offered as-is,
 without any warranty.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)

# According to the examples in Section 7.4.2, unbounded intervals can be constructed with non-common inputs.
@testset "IEEE1788.a" begin
    @test Interval(-Inf, Inf) == entireinterval(Float64)
end

# Examples from Sections 9.7.1 and 9.8
@testset "IEEE1788.b" begin
    @test @interval("[1.2345]") == Interval(0x1.3c083126e978dp+0, 0x1.3c083126e978ep+0)
    @test @interval("[1,+infinity]") == Interval(1.0, Inf)
    @test @decorated("[1,1e3]_com") == DecoratedInterval(Interval(1.0, 1000.0), com)
    @test decoration(@decorated("[1,1e3]_com")) == decoration(DecoratedInterval(Interval(1.0, 1000.0), com))
    @test @decorated("[1,1E3]_COM") == DecoratedInterval(Interval(1.0, 1000.0), com)
    @test decoration(@decorated("[1,1E3]_COM")) == decoration(DecoratedInterval(Interval(1.0, 1000.0), com))
end

# Examples from Table 9.4
@testset "IEEE1788.c" begin
    @test @interval("[1.e-3, 1.1e-3]") == Interval(0x4.189374bc6a7ecp-12, 0x4.816f0068db8bcp-12)
    @test @interval("[-0x1.3p-1, 2/3]") == Interval(-0x9.8000000000000p-4, +0xa.aaaaaaaaaaab0p-4)
    @test @interval("[3.56]") == Interval(0x3.8f5c28f5c28f4p+0, 0x3.8f5c28f5c28f6p+0)
    @test @interval("3.56?1") == Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0)
    @test @interval("3.56?1e2") == Interval(355.0, 357.0)
    @test @interval("3.560?2") == Interval(0x3.8ed916872b020p+0, 0x3.8fdf3b645a1ccp+0)
    @test @interval("3.56?") == Interval(0x3.8e147ae147ae0p+0, 0x3.90a3d70a3d70cp+0)
    @test @interval("3.560?2u") == Interval(0x3.8f5c28f5c28f4p+0, 0x3.8fdf3b645a1ccp+0)
    @test @interval("-10?") == Interval(-10.5, -9.5)
    @test @interval("-10?u") == Interval(-10.0, -9.5)
    @test @interval("-10?12") == Interval(-22.0, 2.0)
end

# Examples from Section 10.5.1
@testset "IEEE1788.d" begin
    @test @interval("[1.234e5,Inf]") == Interval(123400.0, Inf)
    @test @interval("3.1416?1") == Interval(0x3.24395810624dcp+0, 0x3.24467381d7dc0p+0)
    @test @interval("[Empty]") == ∅
end

# Example from Section 11.3
@testset "IEEE1788.e" begin
    @test DecoratedInterval(2, 1) == nai()
end

# Examples from Table 12.1
@testset "IEEE1788.e" begin
    @test @decorated("[ ]") == DecoratedInterval(∅, trv)
    @test decoration(@decorated("[ ]")) == decoration(DecoratedInterval(∅, trv))
    @test @decorated("[entire]") == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test decoration(@decorated("[entire]")) == decoration(DecoratedInterval(Interval(-Inf, Inf), dac))
    @test @decorated("[1.e-3, 1.1e-3]") == DecoratedInterval(Interval(0x4.189374bc6a7ecp-12, 0x4.816f0068db8bcp-12), com)
    @test decoration(@decorated("[1.e-3, 1.1e-3]")) == decoration(DecoratedInterval(Interval(0x4.189374bc6a7ecp-12, 0x4.816f0068db8bcp-12), com))
    @test @decorated("[-Inf, 2/3]") == DecoratedInterval(Interval(-Inf, +0xa.aaaaaaaaaaab0p-4), dac)
    @test decoration(@decorated("[-Inf, 2/3]")) == decoration(DecoratedInterval(Interval(-Inf, +0xa.aaaaaaaaaaab0p-4), dac))
    @test @decorated("[0x1.3p-1,]") == DecoratedInterval(Interval(0x1.3p-1, Inf), dac)
    @test decoration(@decorated("[0x1.3p-1,]")) == decoration(DecoratedInterval(Interval(0x1.3p-1, Inf), dac))
    @test @decorated("[,]") == DecoratedInterval(entireinterval(Float64), dac)
    @test decoration(@decorated("[,]")) == decoration(DecoratedInterval(entireinterval(Float64), dac))
    @test @decorated("3.56?1") == DecoratedInterval(Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0), com)
    @test decoration(@decorated("3.56?1")) == decoration(DecoratedInterval(Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0), com))
    @test @decorated("3.56?1e2") == DecoratedInterval(Interval(355.0, 357.0), com)
    @test decoration(@decorated("3.56?1e2")) == decoration(DecoratedInterval(Interval(355.0, 357.0), com))
    @test @decorated("3.560?2") == DecoratedInterval(Interval(0x3.8ed916872b020p+0, 0x3.8fdf3b645a1ccp+0), com)
    @test decoration(@decorated("3.560?2")) == decoration(DecoratedInterval(Interval(0x3.8ed916872b020p+0, 0x3.8fdf3b645a1ccp+0), com))
    @test @decorated("3.56?") == DecoratedInterval(Interval(0x3.8e147ae147ae0p+0, 0x3.90a3d70a3d70cp+0), com)
    @test decoration(@decorated("3.56?")) == decoration(DecoratedInterval(Interval(0x3.8e147ae147ae0p+0, 0x3.90a3d70a3d70cp+0), com))
    @test @decorated("3.560?2u") == DecoratedInterval(Interval(0x3.8f5c28f5c28f4p+0, 0x3.8fdf3b645a1ccp+0), com)
    @test decoration(@decorated("3.560?2u")) == decoration(DecoratedInterval(Interval(0x3.8f5c28f5c28f4p+0, 0x3.8fdf3b645a1ccp+0), com))
    @test @decorated("-10?") == DecoratedInterval(Interval(-10.5, -9.5), com)
    @test decoration(@decorated("-10?")) == decoration(DecoratedInterval(Interval(-10.5, -9.5), com))
    @test @decorated("-10?u") == DecoratedInterval(Interval(-10.0, -9.5), com)
    @test decoration(@decorated("-10?u")) == decoration(DecoratedInterval(Interval(-10.0, -9.5), com))
    @test @decorated("-10?12") == DecoratedInterval(Interval(-22.0, 2.0), com)
    @test decoration(@decorated("-10?12")) == decoration(DecoratedInterval(Interval(-22.0, 2.0), com))
    @test @decorated("-10??u") == DecoratedInterval(Interval(-10.0, Inf), dac)
    @test decoration(@decorated("-10??u")) == decoration(DecoratedInterval(Interval(-10.0, Inf), dac))
    @test @decorated("-10??") == DecoratedInterval(Interval(-Inf, Inf), dac)
    @test decoration(@decorated("-10??")) == decoration(DecoratedInterval(Interval(-Inf, Inf), dac))
    @test @decorated("[nai]") == nai()
    @test @decorated("3.56?1_def") == DecoratedInterval(Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0), def)
    @test decoration(@decorated("3.56?1_def")) == decoration(DecoratedInterval(Interval(0x3.8ccccccccccccp+0, 0x3.91eb851eb8520p+0), def))
end

# Examples from Section 12.11.3
@testset "IEEE1788.f" begin
    @test @interval("[]") == ∅
    @test @interval("[empty]") == ∅
    @test @interval("[ empty ]") == ∅
    @test @interval("[,]") == entireinterval(Float64)
    @test @interval("[ entire ]") == entireinterval(Float64)
end
# FactCheck.exitstatus()


text:
using WAV
using DifferentialEquations, Flux, DiffEqFlux, Plots
using StochasticDiffEq
using DiffEqBase.EnsembleAnalysis
using StatsBase, Statistics

songfn = "/home/sippycups/Music/2019/81 - 9 21 19 2.wav"
wave, sr = wavread(songfn)
start = 5
secs = 5
start_idx = Int(sr*start)
cut_idx = start_idx + Int(sr*secs)
t = float(start_idx):float(cut_idx)
# tspan = (0.0, 1.0)
target_data = wave[start_idx:cut_idx, 2:end]'
u0 = target_data[:, 1]
dim = length(u0)

sde_data_vars = zeros(dim, length(t)) .+ 1e-3

drift_dudt = Chain(
    Dense(dim, 20, tanh),
    # Dense(20, 20, tanh),
    Dense(20, 20, tanh),
    Dense(20, dim)
) #|> gpu

diffusion_dudt = Chain(Dense(dim,dim))

n_sde = NeuralDSDE(drift_dudt,diffusion_dudt,(t[1], t[end]),SOSRI(),saveat=t,reltol=1e-1,abstol=1e-1)
ps = Flux.params(n_sde)
pred = n_sde(u0)

function predict_n_sde()
	Array(n_sde(u0))
end

function loss_n_sde(;n=5)
	samples = [predict_n_sde() for i in 1:n]
	means = reshape(mean.([[samples[i][j] for i in 1:length(samples)] for j in 1:length(samples[1])]),size(samples[1])...)
	vars = reshape(var.([[samples[i][j] for i in 1:length(samples)] for j in 1:length(samples[1])]),size(samples[1])...)
	sum(abs2,target_data - means) + sum(abs2,sde_data_vars - vars)  # try to replace sde_data_vars with flux_err
end

repeated_data = Iterators.repeated((), 1000)
opt = ADAM(0.025)

losses = []
cb = function ()
	cur_pred = predict_n_sde()
	display(string("pred : ", cur_pred[1:2, 2:end]))
	display(string("target : ", target_data[1:2, 2:end]))
	pl = scatter(t[2:end], target_data[1, 2:end], label="flux_target", markersize=5, markercolor=:blue)
	scatter!(pl, t[2:end], target_data[2, 2:end], label="flux_err_target", markersize=5, markercolor=:green)
	scatter!(pl, t[2:end], cur_pred[1, 2:end], label="flux_pred", markersize=5, markercolor=:red) #, markercolor=:red)
	scatter!(pl, t[2:end], cur_pred[2, 2:end], label="flux_err_pred", markersize=5, markercolor=:orange)
	yticks!([-5:5;])
	xticks!(t[2]:t[end])
	plot!(pl, xlabel="normed_mjd", ylabel="normed_param", size=(900, 900))
end
cb()
Flux.train!(loss_n_sde, ps, repeated_data, opt, cb = cb)

end

max_stars_repo_name: anandijain/audio

max_issues_repo_name: anandijain/audio

max_forks_repo_name: anandijain/audio


text:
# packaged version of align.jl. Run align_prep.jl or prep_full.jl first.

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

import NLopt
import PlotlyJS
using Plots; plotly()

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

save_BSON_flag = true

N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )
N_β = sum( NMRCalibrate.getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )

shift_lb = -ones(N_d)
shift_ub = ones(N_d)

Es = collect( NMRSpectraSimulator.καFIDModelType(Bs[i]) for i = 1:length(Bs) )

a_setp, b_setp, minxs,
    rets = NMRCalibrate.setupitpab(0.1, 10, 0.7; optim_algorithm = :LN_BOBYQA)
#

w = ones(Float64, length(As))
#LS_inds = 1:length(U_cost)
loop_range = 1:length(cost_inds_set)

println("Timing:")
@time obj_funcs, minfs, minxs, rets = NMRCalibrate.aligncompoundκ(y,
    U_y,
    P_y,
    As,
    Bs,
    Es,
    fs,
    SW,
    Δsys_cs,
    a_setp, b_setp, κs_β_DOFs, κs_β_orderings,
    shift_lb,
    shift_ub,
    cost_inds_set;
    loop_range = loop_range,
    w = w,
    N_starts = 100,
    #N_starts = 500, # ethanol region 1.
    local_optim_algorithm = NLopt.LN_BOBYQA,
    #xtol_rel = 1e-3,
    xtol_rel = 1e-9,
    maxeval = 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    κ_lb_default = 1e-2,
    κ_ub_default = 1e2,
    β_max_iters = 500,
    #β_max_iters = 1000,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf);

dummy = 1

if save_BSON_flag
    save_path = joinpath(project_folder, "alignment_results.bson")
    BSON.bson(save_path, region_min_dist = region_min_dist,
    minfs = minfs,
    minxs = minxs,
    rets = rets)
end

function plotregion(P, U, q_U, P_y, y, P_cost, y_cost, display_threshold_factor, display_reduction_factor,
    save_folder, title_string, file_name;
    save_plot_flag = true,
    display_plot_flag = true,
    canvas_size = (1000, 400))

    # # reduce the plotting positions for low signal regions. Otherwise the plot store size will be too large, and the time to load the plot will be long.
    # inds, _ = NMRSpectraSimulator.prunelowsignalentries(q_U, display_threshold_factor, display_reduction_factor)
    # P_display = P[inds]
    # U_display = U[inds]
    # q_U_display = q_U[inds]
    P_display = P
    U_display = U
    q_U_display = q_U

    # plot.
    isdir(save_folder) || mkpath(save_folder)

    plots_save_path = joinpath(save_folder, file_name)
    #title_string = "$(project_name) alignment results, region $(r), real part"

    plot_obj = Plots.plot( P_display,
        real.(q_U_display),
        title = title_string,
        label = "model",
        seriestype = :line,
        ticks = :native,
        xlims = (P_display[1],P_display[end]),
        hover = P_display,
        linewidth = 4,
        xlabel = "ppm",
        ylabel = "real part of spectrum",
        xflip = true,
        size = canvas_size)

    Plots.plot!(plot_obj, P_y, real.(y), label = "full data",
        seriestype = :line,
        linestyle = :dot,
        xflip = true,
        linewidth = 4)

    Plots.plot!(plot_obj, P_cost, real.(y_cost), label = "fit data",
        markershape = :circle,
        seriestype = :scatter,
        xflip = true)

    if save_plot_flag
        Plots.savefig(plot_obj, plots_save_path)
    end

    if display_plot_flag
        display(plot_obj)
    end

    return nothing
end

function plotalignmentresults(As, Es, w, save_folder,
    P_y, y,
    region_min_dist,
    obj_funcs, minfs, minxs, rets,
    display_reduction_factor, display_threshold_factor,
    cost_inds_set, loop_range;
    canvas_size = (1000, 400),
    display_flag = false,
    save_plot_flag = true,
    N_viz = 50000)

    U = LinRange(u_min, u_max, N_viz)
    P = hz2ppmfunc.(U)
    U_rad = U .* (2*π)

    for r in loop_range

        y_cost = y[cost_inds_set[r]]
        P_cost = P_y[cost_inds_set[r]]

        q2 = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = w)

        obj_funcs[r](minxs[r])

        # # debug.
        # N_κ, N_κ_singlets = NMRCalibrate.countκ(Es)
        # N_κ_vars = N_κ + N_κ_singlets
        # NMRCalibrate.parseκ!(Es, ones(N_κ_vars))
        # #

        q_U = q2.(U_rad)

        file_name = "results_real_$(r).html"
        title_string = "$(project_name) alignment results, region $(r), real part"

        plotregion(P, U, q_U, P_y, y, P_cost, y_cost, display_threshold_factor, display_reduction_factor,
            save_folder, title_string, file_name;
            save_plot_flag = save_plot_flag,
            display_plot_flag = display_flag,
            canvas_size = canvas_size)
    end

    return nothing
end

#### visualize.
# minxs[1][1] = 0.000
# minxs[2][1] = 0.000
display_reduction_factor = 100
display_threshold_factor = 0.05/10
if "L-Isoleucine" in molecule_names
    display_reduction_factor = 1
    display_threshold_factor = 0.001/10
end
save_plot_flag = true
display_flag = true

plotalignmentresults(As, Es, w, project_folder,
    P_y, y,
    region_min_dist,
    obj_funcs, minfs, minxs, rets,
    display_reduction_factor, display_threshold_factor,
    cost_inds_set, loop_range;
    canvas_size = (1000, 400),
    display_flag = display_flag,
    save_plot_flag = save_plot_flag,
    N_viz = 50000)

for r in loop_range
    println("region $(r):")
    println("objective: $(minfs[r]), return status: $(rets)")
    println("shift variable:")
    display(minxs)
    println()
end


text:
# # Extended formulation

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Extended Formulation.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Extended Formulation.ipynb)

# In this notebook, we show how to work with the extended formulation of a polyhedron.
# The convex hull of the union of polyhedra that are H-represented can be obtained
# as the projection of a H-representation [Theorem 3.3, B85].
# In order to use the resulting polyhedron as a constraint set in an optimization problem,
# there is no need to compute the resulting H-representation of this projection.
# Moreover, other operations such as intersection are also implemented between extended H-representations.
# We illustrate this with a simple example. We start by defining the H-representation of a square with JuMP.

# [B85] Balas, E., 1985.
# *Disjunctive programming and a hierarchy of relaxations for discrete optimization problems*.
# SIAM Journal on Algebraic Discrete Methods, 6(3), pp.466-486.

using JuMP
model = Model()
@variable(model, -1 <= x <= 1)
@variable(model, -1 <= y <= 1)
using Polyhedra
square = hrep(model)

# Note that the names of the JuMP variables are used as the names of the corresponding dimensions for
# the polyhedron `square`.

dimension_names(square) #!jl
@test dimension_names(square) == ["x", "y"] #jl

# In the following, `diagonal` and `antidiag` are projections of extended H-representations of dimension 7
# hence `diamond` is a projection of an extended H-representation of dimensions 12.

diagonal = convexhull(translate(square, [-2, -2]), translate(square, [2, 2]))
@test fulldim(diagonal.set) == 7 #jl
@test dimension_names(diagonal.set) == ["x"; "y"; fill("", 5)] #jl
antidiag = convexhull(translate(square, [-2,  2]), translate(square, [2, -2]))
@test fulldim(antidiag.set) == 7 #jl
@test dimension_names(antidiag.set) == ["x"; "y"; fill("", 5)] #jl
diamond = diagonal ∩ antidiag
@test fulldim(diamond.set) == 12 #jl
@test dimension_names(diamond.set) == ["x"; "y"; fill("", 10)] #jl

# Note that the names the first two dimensions are still identical to the names of the JuMP variables
# and the auxiliary variables have no name.

dimension_names(diamond.set) #!jl
@test dimension_names(diamond.set) == ["x"; "y"; fill("", 10)] #jl

# We don't need to compute the result of the projection to solve an optimization problem
# over `diamond`. For instance, to compute the maximal value that `y` can take over this
# polytope with the GLPK solver, we can do as follows.
# Note that if we use anonymous JuMP variables, the name of the JuMP variables will be
# the names of the corresponding dimensions of the polyhedron.
# Therefore, we can retrieve the JuMP variable according to the corresponding dimension name
# with `variable_by_name`.

import GLPK
model = Model(GLPK.Optimizer)
@variable(model, [1:2] in diamond)
x = variable_by_name(model, "x")
y = variable_by_name(model, "y")
@objective(model, Max, y)
optimize!(model)
value(x), value(y) #!jl
@test value(x) ≈ 0.0 atol=1e-6 #jl
@test value(y) ≈ 2.0           #jl

# In the optimization problem, above, the auxiliary variables of the extended formulation
# are transparently added inside a bridge.
# To manipulate the auxiliary variables, one can use the extended H-representation directly instead of its projection.
# Note that as the auxiliary dimensions have no name, we cannot use `variable_by_name` to retrieve the corresponding
# JuMP variables. We can instead catch the returned value of `@variable` in some variable `v` in order to use anonymous JuMP
# variables while still assigning the created JuMP variables to `v`.

import GLPK
model = Model(GLPK.Optimizer)
v = @variable(model, [1:12] in diamond.set)
y = variable_by_name(model, "y")
@objective(model, Max, y)
optimize!(model)
value.(v) #!jl
@test value.(v) ≈ [0.0, 2.0, -0.75, -0.25, 0.75, 2.25, 0.25, -0.75, 2.25, 0.75, -0.25, 0.75] #jl

max_stars_repo_path: examples/Extended Formulation.jl

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_path: examples/Extended Formulation.jl

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_path: examples/Extended Formulation.jl

max_forks_repo_name: blegat/Polyhedra