text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

function tryParseIdxSet(arg::Expr)
    # :in appears as arg.head only prior to 0.5
    if arg.head === :(=) || arg.head === :in
        @assert length(arg.args) == 2
        return true, arg.args[1], arg.args[2]
    elseif isexpr(arg, :call) && arg.args[1] === :in
        return true, arg.args[2], arg.args[3]
    else
        return false, nothing, nothing
    end
end

function parseIdxSet(arg::Expr)
    parse_done, idxvar, idxset = tryParseIdxSet(arg)
    if parse_done
        return idxvar, idxset
    end
    error("Invalid syntax: $arg")
end

# addtoexpr(ex::Number, c::Number) = ex + c

addtoexpr(ex::Number, c::Number, x::Number) = ex + c*x

addtoexpr(ex::Number, c::Number, x::Variable) = AffExpr([x],[c],ex)

function addtoexpr{T<:GenericAffExpr}(ex::Number, c::Number, x::T)
    # It's only safe to mutate the first argument.
    if c == 0
        T(ex)
    else
        x = copy(x)
        scale!(x.coeffs, c)
        x.constant *= c
        x.constant += ex
        x
    end
end

function addtoexpr{T<:GenericQuadExpr}(ex::Number, c::Number, x::T)
    # It's only safe to mutate the first argument.
    if c == 0
        T(ex)
    else
        x = copy(x)
        scale!(x.qcoeffs, c)
        scale!(x.aff.coeffs, c)
        x.aff.constant *= c
        x.aff.constant += ex
        x
    end
end

addtoexpr(ex::Number, c::Variable, x::Variable) = QuadExpr([c],[x],[1.0],zero(AffExpr))

function addtoexpr{T<:GenericAffExpr}(ex::Number, c::T, x::T)
    q = c*x
    q.aff.constant += ex
    q
end

function addtoexpr{C,V}(ex::Number, c::GenericAffExpr{C,V}, x::V)
    q = c*x
    q.aff.constant += ex
    q
end

function addtoexpr{T<:GenericQuadExpr}(ex::Number, c::T, x::Number)
    if x == 0
        T(ex)
    else
        q = c*x
        q.aff.constant += ex
        q
    end
end

function addtoexpr(aff::GenericAffExpr, c::Number, x::Number)
    aff.constant += c*x
    aff
end

function addtoexpr{C,V}(aff::GenericAffExpr{C,V}, c::Number, x::V)
    if c != 0
        push!(aff.vars,   x)
        push!(aff.coeffs, c)
    end
    aff
end

function addtoexpr{C,V}(aff::GenericAffExpr{C,V},c::Number,x::GenericAffExpr{C,V})
    if c != 0
        append!(aff.vars, x.vars)
        sizehint!(aff.coeffs, length(aff.coeffs)+length(x.coeffs))
        for i in 1:length(x.coeffs)
            push!(aff.coeffs, c*x.coeffs[i])
        end
        aff.constant += c*x.constant
    end
    aff
end

# help w/ ambiguity
addtoexpr{C,V<:Number}(aff::GenericAffExpr{C,V}, c::Number, x::Number) = aff + c*x

function addtoexpr{C,V}(aff::GenericAffExpr{C,V}, c::V, x::Number)
    if x != 0
        push!(aff.vars,   c)
        push!(aff.coeffs, x)
    end
    aff
end

addtoexpr{C,V}(aff::GenericAffExpr{C,V},c::V,x::V) =
    GenericQuadExpr{C,V}([c],[x],[one(C)],aff)

# TODO: add generic versions of following two methods
addtoexpr(aff::AffExpr,c::AffExpr,x::Variable) =
    QuadExpr(c.vars,
             fill(x,length(c.vars)),
             c.coeffs,
             addtoexpr(aff,c.constant,x))

addtoexpr(aff::AffExpr,c::Variable,x::AffExpr) =
    QuadExpr(fill(c,length(x.vars)),
             x.vars,
             x.coeffs,
             addtoexpr(aff,c,x.constant))

function addtoexpr{C,V}(aff::GenericAffExpr{C,V},c::GenericAffExpr{C,V},x::Number)
    if x != 0
        append!(aff.vars, c.vars)
        append!(aff.coeffs, c.coeffs * x)
        aff.constant += c.constant * x
    end
    aff
end

function addtoexpr{C,V}(aff::GenericAffExpr{C,V}, c::GenericQuadExpr{C,V}, x::Number)
    if x == 0
        GenericQuadExpr{C,V}(aff)
    else
        GenericQuadExpr{C,V}(copy(c.qvars1),
                            copy(c.qvars2),
                            c.qcoeffs*x,
                            addtoexpr(aff,c.aff,x))
    end
end

function addtoexpr{C,V}(aff::GenericAffExpr{C,V}, c::Number, x::GenericQuadExpr{C,V})
    if c == 0
        GenericQuadExpr{C,V}(aff)
    else
        GenericQuadExpr{C,V}(copy(x.qvars1),
                            copy(x.qvars2),
                            c*x.qcoeffs,
                            addtoexpr(aff,c,x.aff))
    end
end

function addtoexpr{C,V}(ex::GenericAffExpr{C,V}, c::GenericAffExpr{C,V}, x::GenericAffExpr{C,V})
    q = convert(GenericQuadExpr{C,V}, ex)
    addtoexpr(q, one(C), c*x)
    q
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::Number,x::V)
    if c != 0
        push!(quad.aff, convert(C,c), x)
    end
    quad
end

function addtoexpr(quad::GenericQuadExpr,c::Number,x::Number)
    quad.aff.constant += c*x
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},x::V,y::V)
    push!(quad.qvars1, x)
    push!(quad.qvars2, y)
    push!(quad.qcoeffs, one(C))
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::Number,x::GenericAffExpr{C,V})
    if c != 0
        append!(quad.aff.vars, x.vars)
        sizehint!(quad.aff.coeffs, length(quad.aff.coeffs)+length(x.coeffs))
        for i in 1:length(x.coeffs)
            push!(quad.aff.coeffs, c*x.coeffs[i])
        end
        quad.aff.constant += c*x.constant
    end
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::GenericAffExpr{C,V},x::Number)
    if x != 0
        addtoexpr(quad.aff,c,x)
    end
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::GenericAffExpr{C,V},x::V)
    append!(quad.qvars1, c.vars)
    append!(quad.qvars2, fill(x,length(c.vars)))
    append!(quad.qcoeffs, c.coeffs)
    addtoexpr(quad.aff,c.constant,x)
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::V,x::GenericAffExpr{C,V})
    append!(quad.qvars1, fill(c,length(x.vars)))
    append!(quad.qvars2, x.vars)
    append!(quad.qcoeffs, x.coeffs)
    addtoexpr(quad.aff,c,x.constant)
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::GenericQuadExpr{C,V},x::Number)
    if x != 0
        append!(quad.qvars1,c.qvars1)
        append!(quad.qvars2,c.qvars2)
        sizehint!(quad.qcoeffs, length(quad.qcoeffs)+length(c.qcoeffs))
        for i in 1:length(c.qcoeffs)
            push!(quad.qcoeffs, c.qcoeffs[i]*x)
        end
        addtoexpr(quad,c.aff,x)
    end
    quad
end

function addtoexpr{C,V}(quad::GenericQuadExpr{C,V},c::Number,x::GenericQuadExpr{C,V})
    if c != 0
        append!(quad.qvars1,x.qvars1)
        append!(quad.qvars2,x.qvars2)
        sizehint!(quad.qcoeffs, length(quad.qcoeffs)+length(x.qcoeffs))
        for i in 1:length(x.qcoeffs)
            push!(quad.qcoeffs, c*x.qcoeffs[i])
        end
        addtoexpr(quad,c,x.aff)
    end
    quad
end

function addtoexpr{C,V}(ex::GenericQuadExpr{C,V}, c::GenericAffExpr{C,V}, x::GenericAffExpr{C,V})
    q = c*x
    addtoexpr(ex, 1.0, q)
    ex
end

# Catch nonlinear expressions and parameters being used in addconstraint, etc.

typealias _NLExpr Union{NonlinearExpression,NonlinearParameter}
_nlexprerr() = error("""Cannot use nonlinear expression or parameter in @constraint or @objective.
                        Use @NLconstraint or @NLobjective instead.""")
# Following three definitions avoid ambiguity warnings
addtoexpr{C,V<:_NLExpr}(expr::GenericQuadExpr{C,V}, c::GenericAffExpr{C,V}, x::V) = _nlexprerr()
addtoexpr{C,V<:_NLExpr}(expr::GenericQuadExpr{C,V}, c::Number, x::V) = _nlexprerr()
addtoexpr{C,V<:_NLExpr}(expr::GenericQuadExpr{C,V}, c::V, x::GenericAffExpr{C,V}) = _nlexprerr()
for T1 in (GenericAffExpr,GenericQuadExpr), T2 in (Number,Variable,GenericAffExpr,GenericQuadExpr)
    @eval addtoexpr(::$T1, ::$T2, ::_NLExpr) = _nlexprerr()
    @eval addtoexpr(::$T1, ::_NLExpr, ::$T2) = _nlexprerr()
end

addtoexpr(ex, c, x) = ex + c*x

@generated addtoexpr_reorder(ex, arg) = :(addtoexpr(ex, 1.0, arg))

@generated function addtoexpr_reorder(ex, x, y)
    if x <: Union{Variable,AffExpr} && y <: Number
        :(addtoexpr(ex, y, x))
    else
        :(addtoexpr(ex, x, y))
    end
end

@generated function addtoexpr_reorder(ex, args...)
    n = length(args)
    @assert n ≥ 3
    varidx = find(t -> (t == Variable || t == AffExpr), collect(args))
    allscalar = all(t -> (t <: Number), args[setdiff(1:n, varidx)])
    idx = (allscalar && length(varidx) == 1) ? varidx[1] : n
    coef = Expr(:call, :*, [:(args[$i]) for i in setdiff(1:n,idx)]...)
    :(addtoexpr(ex, $coef, args[$idx]))
end

function parseCurly(x::Expr, aff::Symbol, lcoeffs, rcoeffs, newaff=gensym())
    header = x.args[1]
    if length(x.args) < 3
        error("Need at least two arguments for $header")
    end
    if issum(header)
        parseSum(x, aff, lcoeffs, rcoeffs, newaff)
    elseif header ∈ [:norm1, :norm2, :norminf, :norm∞]
        parseNorm(header, x, aff, lcoeffs, rcoeffs, newaff)
    else
        error("Expected sum or norm2 outside curly braces; got $header")
    end
end

function parseSum(x::Expr, aff::Symbol, lcoeffs, rcoeffs, newaff)
    # we have a filter condition
    if isexpr(x.args[2],:parameters)
        cond = x.args[2]
        if length(cond.args) != 1
            error("No commas after semicolon allowed in sum expression, use && for multiple conditions")
        end
        # generate inner loop code first and then wrap in for loops
        inneraff, innercode = parseExpr(x.args[3], aff, lcoeffs, rcoeffs, aff)
        code = quote
            if $(esc(cond.args[1]))
                $innercode
            end
        end
        for level in length(x.args):-1:4
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
        end
    else # no condition
        inneraff, code = parseExpr(x.args[2], aff, lcoeffs, rcoeffs, aff)
        for level in length(x.args):-1:3
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
        end
        len = :len
        # precompute the number of elements to add
        # this is unncessary if we're just summing constants
        _, lastidxset = parseIdxSet(x.args[length(x.args)]::Expr)
        preblock = :($len += length($(esc(lastidxset))))
        for level in (length(x.args)-1):-1:3
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            preblock = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $preblock
                end
            end)
        end
        preblock = quote
            $len = 0
            $preblock
            _sizehint_expr!($aff,$len)
        end
        code = :($preblock;$code)
    end
    :($code; $newaff=$aff)
end

function parseNorm(normp::Symbol, x::Expr, aff::Symbol, lcoeffs, rcoeffs, newaff)
    @assert string(x.args[1])[1:4] == "norm"
    # we have a filter condition
    finalaff = gensym()
    gennorm  = gensym()
    len      = gensym()
    normexpr = gensym()
    if isexpr(x.args[2],:parameters)
        cond = x.args[2]
        if length(cond.args) != 1
            error("No commas after semicolon allowed in sum expression, use && for multiple conditions")
        end
        # generate inner loop code first and then wrap in for loops
        inneraff, innercode = parseExprToplevel(x.args[3], :normaff)
        code = quote
            if $(esc(cond.args[1]))
                normaff = 0.0
                $innercode
                push!($normexpr, $inneraff)
            end
        end
        for level in length(x.args):-1:4
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
        end
        preblock = :($normexpr = GenericAffExpr[])
    else # no condition
        inneraff, code = parseExprToplevel(x.args[2], :normaff)
        code = :(normaff = 0.0; $code; push!($normexpr, $inneraff))
        _, lastidxset = parseIdxSet(x.args[length(x.args)]::Expr)
        preblock = :($len += length($(esc(lastidxset))))
        for level in length(x.args):-1:3
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
            preblock = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $preblock
                end
            end)
        end
        preblock = quote
            $len = 0
            $preblock
            $normexpr = GenericAffExpr[]
            _sizehint_expr!($normexpr, $len)
        end
    end
    if normp == :norm2
        param = 2
    elseif normp == :norm1
        param = 1
    elseif normp ∈ [:norminf, :norm∞]
        param = Inf
    else
        error("Unrecognized norm: $normp")
    end
    quote
        $preblock
        $code
        $gennorm = _build_norm($param,$normexpr)
        $newaff = $(Expr(:call, :addtoexpr_reorder, aff,lcoeffs...,gennorm,rcoeffs...))
    end
end

function parseGenerator(x::Expr, aff::Symbol, lcoeffs, rcoeffs, newaff=gensym())
    @assert isexpr(x,:call)
    @assert length(x.args) > 1
    @assert isexpr(x.args[2],:generator)
    header = x.args[1]
    if issum(header)
        parseGeneratorSum(x.args[2], aff, lcoeffs, rcoeffs, newaff)
    elseif header ∈ [:norm1, :norm2, :norminf, :norm∞]
        parseGeneratorNorm(header, x.args[2], aff, lcoeffs, rcoeffs, newaff)
    else
        error("Expected sum or norm2 outside generator expression; got $header")
    end
end

function parseGeneratorSum(x::Expr, aff::Symbol, lcoeffs, rcoeffs, newaff)
    # we have a filter condition
    if isexpr(x.args[2],:filter)
        cond = x.args[2].args[1]
        # generate inner loop code first and then wrap in for loops
        inneraff, innercode = parseExpr(x.args[1], aff, lcoeffs, rcoeffs, aff)
        code = quote
            if $(esc(cond))
                $innercode
            end
        end
        for level in length(x.args[2].args):-1:2
            _idxvar, idxset = parseIdxSet(x.args[2].args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
        end
    else # no condition
        inneraff, code = parseExpr(x.args[1], aff, lcoeffs, rcoeffs, aff)
        for level in length(x.args):-1:2
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
        end
        len = :len
        # precompute the number of elements to add
        # this is unncessary if we're just summing constants
        _, lastidxset = parseIdxSet(x.args[length(x.args)]::Expr)
        preblock = :($len += length($(esc(lastidxset))))
        for level in (length(x.args)-1):-1:2
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            preblock = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $preblock
                end
            end)
        end
        preblock = quote
            $len = 0
            $preblock
            _sizehint_expr!($aff,$len)
        end
        code = :($preblock;$code)
    end
    :($code; $newaff=$aff)
end

function parseGeneratorNorm(normp::Symbol, x::Expr, aff::Symbol, lcoeffs, rcoeffs, newaff)
    @assert string(normp)[1:4] == "norm"
    finalaff = gensym()
    gennorm  = gensym()
    len      = gensym()
    normexpr = gensym()
    # we have a filter condition
    if isexpr(x.args[2],:filter)
        cond = x.args[2].args[1]
        # generate inner loop code first and then wrap in for loops
        inneraff, innercode = parseExprToplevel(x.args[1], :normaff)
        code = quote
            if $(esc(cond))
                normaff = 0.0
                $innercode
                push!($normexpr, $inneraff)
            end
        end
        for level in length(x.args[2].args):-1:2
            _idxvar, idxset = parseIdxSet(x.args[2].args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
        end
        preblock = :($normexpr = GenericAffExpr[])
    else # no condition
        inneraff, code = parseExprToplevel(x.args[1], :normaff)
        code = :(normaff = 0.0; $code; push!($normexpr, $inneraff))
        _, lastidxset = parseIdxSet(x.args[length(x.args)]::Expr)
        preblock = :($len += length($(esc(lastidxset))))
        for level in length(x.args):-1:2
            _idxvar, idxset = parseIdxSet(x.args[level]::Expr)
            idxvar = esc(_idxvar)
            code = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $code
                end
            end)
            preblock = :(let
                $(localvar(idxvar))
                for $idxvar in $(esc(idxset))
                    $preblock
                end
            end)
        end
        preblock = quote
            $len = 0
            $preblock
            $normexpr = GenericAffExpr[]
            _sizehint_expr!($normexpr, $len)
        end
    end
    if normp == :norm2
        param = 2
    elseif normp == :norm1
        param = 1
    elseif normp ∈ [:norminf, :norm∞]
        param = Inf
    else
        error("Unrecognized norm: $normp")
    end
    quote
        $preblock
        $code
        $gennorm = _build_norm($param,$normexpr)
        $newaff = $(Expr(:call, :addtoexpr_reorder, aff,lcoeffs...,gennorm,rcoeffs...))
    end
end

is_complex_expr(ex) = isa(ex,Expr) && !isexpr(ex,:ref)

parseExprToplevel(x, aff::Symbol) = parseExpr(x, aff, [], [])

# output is assigned to newaff
function parseExpr(x, aff::Symbol, lcoeffs::Vector, rcoeffs::Vector, newaff::Symbol=gensym())
    if !isa(x,Expr)
        # at the lowest level
        callexpr = Expr(:call,:addtoexpr_reorder,aff,lcoeffs...,esc(x),rcoeffs...)
        return newaff, :($newaff = $callexpr)
    else
        if x.head == :call && x.args[1] == :+
            b = Expr(:block)
            aff_ = aff
            for arg in x.args[2:(end-1)]
                aff_, code = parseExpr(arg, aff_, lcoeffs, rcoeffs)
                push!(b.args, code)
            end
            newaff, code = parseExpr(x.args[end], aff_, lcoeffs, rcoeffs, newaff)
            push!(b.args, code)
            return newaff, b
        elseif x.head == :call && x.args[1] == :-
            if length(x.args) == 2 # unary subtraction
                return parseExpr(x.args[2], aff, vcat(-1.0, lcoeffs), rcoeffs, newaff)
            else # a - b - c ...
                b = Expr(:block)
                aff_, code = parseExpr(x.args[2], aff, lcoeffs, rcoeffs)
                push!(b.args, code)
                for arg in x.args[3:(end-1)]
                    aff_,code = parseExpr(arg, aff_, vcat(-1.0, lcoeffs), rcoeffs)
                    push!(b.args, code)
                end
                newaff,code = parseExpr(x.args[end], aff_, vcat(-1.0, lcoeffs), rcoeffs, newaff)
                push!(b.args, code)
                return newaff, b
            end
        elseif x.head == :call && x.args[1] == :*
            # we might need to recurse on multiple arguments, e.g.,
            # (x+y)*(x+y)
            n_expr = mapreduce(is_complex_expr, +, x.args)
            if n_expr == 1 # special case, only recurse on one argument and don't create temporary objects
                which_idx = 0
                for i in 2:length(x.args)
                    if is_complex_expr(x.args[i])
                        which_idx = i
                    end
                end
                return parseExpr(x.args[which_idx], aff, vcat(lcoeffs, [esc(x.args[i]) for i in 2:(which_idx-1)]),
                                                         vcat(rcoeffs, [esc(x.args[i]) for i in (which_idx+1):length(x.args)]),
                                                         newaff)
            else
                blk = Expr(:block)
                for i in 2:length(x.args)
                    if is_complex_expr(x.args[i])
                        s = gensym()
                        newaff_, parsed = parseExprToplevel(x.args[i], s)
                        push!(blk.args, :($s = 0.0; $parsed))
                        x.args[i] = newaff_
                    else
                        x.args[i] = esc(x.args[i])
                    end
                end
                callexpr = Expr(:call,:addtoexpr_reorder,aff,lcoeffs...,x.args[2:end]...,rcoeffs...)
                push!(blk.args, :($newaff = $callexpr))
                return newaff, blk
            end
        elseif x.head == :call && x.args[1] == :^ && is_complex_expr(x.args[2])
            if x.args[3] == 2
                blk = Expr(:block)
                s = gensym()
                newaff_, parsed = parseExprToplevel(x.args[2], s)
                push!(blk.args, :($s = 0.0; $parsed))
                push!(blk.args, :($newaff = $aff + $(Expr(:call,:*,lcoeffs...,newaff_,newaff_,rcoeffs...))))
                return newaff, blk
            elseif x.args[3] == 1
                return parseExpr(x.args[2], aff, lcoeffs, rcoeffs)
            elseif x.args[3] == 0
                return parseExpr(1, aff, lcoeffs, rcoeffs)
            else
                blk = Expr(:block)
                s = gensym()
                newaff_, parsed = parseExprToplevel(x.args[2], s)
                push!(blk.args, :($s = 0.0; $parsed))
                push!(blk.args, :($newaff = $aff + $(Expr(:call,:*,lcoeffs...,Expr(:call,:^,newaff_,esc(x.args[3])),rcoeffs...))))
                return newaff, blk
            end
        elseif x.head == :call && x.args[1] == :/
            @assert length(x.args) == 3
            numerator = x.args[2]
            denom = x.args[3]
            return parseExpr(numerator, aff, lcoeffs,vcat(esc(:(1/$denom)),rcoeffs),newaff)
        elseif isexpr(x,:call) && length(x.args) >= 2 && isexpr(x.args[2],:generator)
            return newaff, parseGenerator(x,aff,lcoeffs,rcoeffs,newaff)
        elseif isexpr(x,:call) && length(x.args) >= 2 && isexpr(x.args[2],:flatten)
            flatten_error(x.args[2])
        elseif x.head == :curly
            return newaff, parseCurly(x,aff,lcoeffs,rcoeffs,newaff)
        else # at lowest level?
            !isexpr(x,:comparison) || error("Unexpected comparison in expression $x")
            callexpr = Expr(:call,:addtoexpr_reorder,aff,lcoeffs...,esc(x),rcoeffs...)
            return newaff, :($newaff = $callexpr)
        end
    end
end

# Semi-automatically generated precompilation hints
include("precompile.jl")


text:
function buildTree(img, capacity)

  root = initQuadTree(img, capacity, 1)

  q=Queue{QuadTree}()
  enqueue!(q, root)
  leaves = []

  while ~isempty(q)
    quadTree = dequeue!(q)

    if quadTree.leaf == true
      continue

    elseif var(Float64.(
        img[quadTree.boundary.x -
          quadTree.boundary.width + 1: quadTree.boundary.x +
          quadTree.boundary.width,
          quadTree.boundary.y -
          quadTree.boundary.height + 1: quadTree.boundary.y +
          quadTree.boundary.height])) < capacity
      quadTree.leaf = true

      xs = quadTree.boundary.x - quadTree.boundary.width +
         1:quadTree.boundary.x + quadTree.boundary.width
      ys = quadTree.boundary.y - quadTree.boundary.height +
         1:quadTree.boundary.y + quadTree.boundary.height
      pxls = collect(zip(repeat(ys, inner=length(xs)),
          repeat(xs, outer=length(ys))))

      for i in pxls
        push!(quadTree.pixels, Pixel(i[1], i[2]))
      end

      push!(leaves, quadTree)

    else
      children = branchTree(quadTree)

      for child in children
        enqueue!(q, child)
      end

      continue
    end
  end

  return leaves
end

ext: jl

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
using PyCall

# Load preprocessed data from a single PyTorch Pickle file.
# As Julia does not support parsing PyTorch files itself,
# delegate to the Python package instead.
function load_torch(path::String)
	torch = pyimport_conda("torch", "pytorch")
    @info "Loading PyTorch file '$path'."
    dataset = torch.load(abspath(path))
    return dataset
end

# Create an iterator for rows of the preprocessed data.
function load_summaries(paths::Array{String})::Channel{SummaryPair}
	return Channel{SummaryPair}() do channel
		for path ∈ paths
			data = load_torch(path)
			@info "Loaded $(length(data)) rows."
			for dict ∈ data
				push!(channel, SummaryPair(dict))
			end
		end
	end
end

# Check if the file is of the given corpus type.
function is_corpus_type(path::String, type::String)
	@assert type ∈ ["train", "valid", "test"]
	isfile(path) && occursin(type, basename(path))
end

function filter_corpus_type!(paths::Array{String}, type::String)
    filter!(path -> is_corpus_type(path, type), paths)
end

# Create an iterator for rows of the preprocessed training/test/validation data.
function summaries(path::String, type::String)::Channel{SummaryPair}
	paths = readdir(path, join=true, sort=true)
	filter_corpus_type!(paths, type)

	return Channel{SummaryPair}() do channel
		loader = load_summaries(paths)
		for row ∈ loader
			push!(channel, row)
		end
	end
end


text:
function formulation(
    i::BipartiteMatchingInstance{Int},
    ::DefaultLinearFormulation;
    solver=nothing,
)
    indices = collect((src(e), dst(e)) for e in edges(i.graph))

    model = Model(solver)
    x = @variable(solver.model, [indices], binary = true)

    for i in 1:n_arms # Left nodes.
        @constraint(solver.model, sum(solver.x[(i, j)] for j in 1:n_arms) == 1)
    end
    for j in 1:n_arms # Right nodes.
        @constraint(solver.model, sum(solver.x[(i, j)] for i in 1:n_arms) == 1)
    end

    @objective(
        model,
        Max,
        sum(rewards[(i, j)] * solver.x[(i, j)] for (i, j) in keys(rewards))
    )
    @constraint(model, sum(x) <= i.m)

    set_silent(model)

    return m, x
end


text:
using ObjectiveC
using Test

@test startswith(read(`hostname`, String), hostname())


text:
using DataInterpolations, Test
using Random

# Linear Interpolation
u = 2.0collect(1:10)
t = 1.0collect(1:10)
A = LinearInterpolation(u,t)

@test A(1) == u[1]
@test A(5) == u[5]
@test A(5.5) == 11.0

u = vcat(2.0collect(1:10)', 3.0collect(1:10)')
A = LinearInterpolation(u,t)

@test A(1) == u[:,1]
@test A(5) == u[:,5]
@test A(5.5) == [11.0,16.5]

# Quadratic Interpolation
u = [1.0, 4.0, 9.0, 16.0]
t = [1.0, 2.0, 3.0, 4.0]
A = QuadraticInterpolation(u,t)

@test A(2.0) == 4.0
@test A(1.5) == 2.25
@test A(3.5) == 12.25
@test A(2.5) == 6.25

u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
A = QuadraticInterpolation(u,t)

@test A(2.0) == [4.0,4.0]
@test A(1.5) == [2.25,2.25]
@test A(3.5) == [12.25,12.25]
@test A(2.5) == [6.25,6.25]

# Lagrange Interpolation
u = [1.0, 4.0, 9.0]
t = [1.0, 2.0, 3.0]
A = LagrangeInterpolation(u,t)

@test A(2.0) == 4.0
@test A(1.5) == 2.25

u = [1.0, 8.0, 27.0, 64.0]
t = [1.0, 2.0, 3.0, 4.0]
A = LagrangeInterpolation(u,t)

@test A(2.0) == 8.0
@test A(1.5) ≈ 3.375
@test A(3.5) ≈ 42.875

u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
A = LagrangeInterpolation(u,t)

@test A(2.0) == [4.0,4.0]
@test A(1.5) ≈ [2.25,2.25]
@test A(3.5) ≈ [12.25,12.25]

# Akima Interpolation
u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
t = collect(0.0:10.0)
A = AkimaInterpolation(u, t)

@test A(0.0) ≈ 0.0
@test A(0.5) ≈ 1.375
@test A(1.0) ≈ 2.0
@test A(1.5) ≈ 1.5
@test A(2.5) ≈ 1.953125
@test A(3.5) ≈ 2.484375
@test A(4.5) ≈ 4.1363636363636366866103344
@test A(5.1) ≈ 5.9803623910336236590978842
@test A(6.5) ≈ 5.5067291516462386624652936
@test A(7.2) ≈ 5.2031367459745245795943447
@test A(8.6) ≈ 4.1796554159017080820603951
@test A(9.9) ≈ 3.4110386597938129327189927
@test A(10.0) ≈ 3.0

@testset "ConstantInterpolation" begin

    t = [1.0, 2.0, 3.0, 4.0]

    @testset "Vector case" for u in
        [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[2]
        @test A(2.0) == u[2]
        @test A(2.5) == u[3]
        @test A(3.0) == u[3]
        @test A(3.5) == u[1]
        @test A(4.0) == u[1]
        @test A(4.5) == u[1]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[1]
        @test A(2.0) == u[2]
        @test A(2.5) == u[2]
        @test A(3.0) == u[3]
        @test A(3.5) == u[3]
        @test A(4.0) == u[1]
        @test A(4.5) == u[1]
    end

    @testset "Matrix case" for u in
        [[1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0], ["B" "C" "A" "B"; "B" "C" "A" "B"]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[:,1]
        @test A(1.0) == u[:,1]
        @test A(1.5) == u[:,2]
        @test A(2.0) == u[:,2]
        @test A(2.5) == u[:,3]
        @test A(3.0) == u[:,3]
        @test A(3.5) == u[:,1]
        @test A(4.0) == u[:,1]
        @test A(4.5) == u[:,1]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[:,1]
        @test A(1.0) == u[:,1]
        @test A(1.5) == u[:,1]
        @test A(2.0) == u[:,2]
        @test A(2.5) == u[:,2]
        @test A(3.0) == u[:,3]
        @test A(3.5) == u[:,3]
        @test A(4.0) == u[:,1]
        @test A(4.5) == u[:,1]
    end
end

# QuadraticSpline Interpolation
u = [0.0, 1.0, 3.0]
t = [-1.0, 0.0, 1.0]

A = QuadraticSpline(u,t)

#             Solution ->
#             f(x) = (x+1)^2 for x -> [-1.0, 0.0]
#             f(x) = 1+2x    for x -> [0.0, 1.0]

@test A(-0.5) == 0.25
@test A(0.7) == 2.4
@test A(-1.0) == 0.0
@test A(0.0) == 1.0
@test A(1.0) == 3.0

# CubicSpline Interpolation
u = [0.0, 1.0, 3.0]
t = [-1.0, 0.0, 1.0]

A = CubicSpline(u,t)

#             Solution ->
#             f(x) = 1 + 1.5x + x^2 + 0.5x^3 for x -> [-1.0, 0.0]
#             f(x) = 1 + 1.5x + x^2 - 0.5x^3   for x -> [0.0, 1.0]

@test A(-0.5) == 0.4375
@test A(0.5) == 1.9375
@test A(-0.7) == 0.2685
@test A(0.3) == 1.5265
@test A(-1.0) == 0.0
@test A(0.0) == 1.0
@test A(1.0) == 3.0

# BSpline Interpolation and Approximation
t = [0,62.25,109.66,162.66,205.8,252.3]
u = [14.7,11.51,10.41,14.95,12.24,11.22]

A = BSplineInterpolation(u,t,2,:Uniform,:Uniform)

@test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
@test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251]

A = BSplineInterpolation(u,t,2,:ArcLen,:Average)

@test [A(25.0), A(80.0)] == [13.363814458968486, 10.685201117692609]
@test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463]

A = BSplineApprox(u,t,2,4,:Uniform,:Uniform)

@test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178]
@test [A(190.0), A(225.0)] ≈ [13.851245975109263, 12.963685868886575]

# Curvefit Interpolation
Random.seed!(12345)
model(x, p) = @. p[1]/(1+exp(x-p[2]))
t = range(-10, stop=10, length=40)
u = model(t, [1.0, 2.0]) + 0.01*randn(length(t))
p0 = [0.5, 0.5]

A = Curvefit(u,t,model,p0,LBFGS())

ts = [-7.0,-2.0,0.0,2.5,5.0]
vs = [1.0039795744162028,0.9854877725868618, 0.881099402277441, 0.3717861289542075, 0.04623053035113763]
us = A.(ts)

@test vs ≈ us

# missing values handling tests
u = [1.0, 4.0, 9.0, 16.0, 25.0, missing, missing]
t = [1.0, 2.0, 3.0, 4.0, missing, 6.0, missing]
A = QuadraticInterpolation(u,t)

@test A(2.0) == 4.0
@test A(1.5) == 2.25
@test A(3.5) == 12.25
@test A(2.5) == 6.25

u = copy(hcat(u, u)')
A = QuadraticInterpolation(u,t)

@test A(2.0) == [4.0, 4.0]
@test A(1.5) == [2.25, 2.25]
@test A(3.5) == [12.25,12.25]
@test A(2.5) == [6.25, 6.25]


text:
compress_then_write(b::StringVector{T}, io) where {T} = begin

    buffer_meta  = (type=eltype(b.buffer) ,  len=write(io, Blosc.compress(b.buffer)))
    offsets_meta = (type=eltype(b.offsets),  len=write(io, Blosc.compress(b.offsets)))
    lengths_meta = (type=eltype(b.lengths),  len=write(io, Blosc.compress(b.lengths)))

    (
        metadata = (buffer = buffer_meta, offsets = offsets_meta, lengths = lengths_meta),
        type = typeof(b),
    )
end

column_loader(::Type{StringVector{T}}, io, metadata) where {T} = begin
    buffer  = column_loader(metadata.metadata.buffer.type , io, metadata.metadata.buffer)
    offsets = column_loader(metadata.metadata.offsets.type, io, metadata.metadata.offsets)
    lengths = column_loader(metadata.metadata.lengths.type, io, metadata.metadata.lengths)

    metadata.type(buffer, offsets, lengths)
end

# tests at test/test-stringarray.jl


text:
@testset "Display" begin
    let
        foo() = foo()
        @test Atom.render(Atom.Inline(), Atom.@errs foo()) != nothing
    end

    # issue: https://github.com/JunoLab/Juno.jl/issues/376
    @test Atom.render(Atom.Inline(), :abc)[:contents] == [":abc"]
    @test Atom.render(Atom.Inline(), Symbol("A B"))[:contents] == ["Symbol(\"A B\")"]

    is_content_in(contents, s) = any(c->occursin(s, c), contents)

    # function display
    # basic
    let d = Atom.render(Atom.Inline(), sin)
        @test d[:type] === :lazy
        @test is_content_in(d[:head][:contents], "sin")
    end

    # in a specific module
    m = eval(:(module $(gensym(:atomjltestmodule)) end))

    let d = Atom.render(Atom.Inline(), Core.eval(m, :(mock_func() = nothing)))
        @test d[:type] === :lazy
        @test is_content_in(d[:head][:contents], "$(m).mock_func")
    end

    # macro
    let d = Atom.render(Atom.Inline(), Core.eval(m, :(macro mock_macro() end)))
        @test d[:type] === :lazy
        @test (is_content_in(d[:head][:contents], "$(m).@mock_macro") ||
               is_content_in(d[:head][:contents], "$(m).var\"@mock_macro\"")) # version 1.6
    end

    @test length(Atom.trim(rand(50), 20)) == 20

    @test Atom.isanon(sin) == false
    @test Atom.isanon(x -> x) == true # the `== true` might fix a weird travis error

    @test Atom.pluralize([1,2,3], "entry", "entries") == "3 entries"
    @test Atom.pluralize([1], "entry", "entries") == "1 entry"
end


text:
function getdocs(ctx::Context, cursor::CXCursor, nested = [])
	getblock(ctx).flags.nodocs && return nothing

	comment = clang_Cursor_getParsedComment(cursor)
	kind = clang_Comment_getKind(comment)
	if kind == CXComment_FullComment
		return Markdown.MD(ctx, cursor, comment)
	elseif kind == CXComment_Null
		return Markdown.MD(ctx, cursor)
	else
		@warn "Expected a full comment, but got $(kind)"
		return nothing
	end
end

function Markdown.MD(ctx::Context, cursor::CXCursor)
	tokens = tokenize(ctx.tu[], cursor)
	tokens = join(map(token -> clang_getTokenKind(token) == CXToken_Comment ? "" : string(ctx.tu[], token), tokens), ' ')
	tokens = replace(tokens, "{" => "{\n")
	tokens = replace(tokens, " ;" => ";\n")
	tokens = replace(tokens, " }" => "\n}")
	tokens = replace(tokens, " ," => ",")
	tokens = replace(tokens, ": " => ":")
	tokens = replace(tokens, " :" => ":")
	tokens = replace(tokens, "( " => "(")
	tokens = replace(tokens, " )" => ")")
	tokens = replace(tokens, "[ " => "[")
	tokens = replace(tokens, " [" => "[")
	tokens = replace(tokens, " ]" => "]")
	tokens = replace(tokens, " (" => "(")
	tokens = replace(tokens, r"(\w)(\(\*)\s+" => s"\1 \2")
	tokens = replace(tokens, "  " => " ")
	tokens = replace(tokens, r"\n\s+" => "\n")

	indents = 0
	result = ""
	isParens = false
	wasComma = false
	wasEOL = false
	for c in tokens
		if c == '}'
			indents -= 1
		elseif c == '{'
			indents += 1
		elseif c == '('
			isParens = true
		elseif c == ')'
			isParens = false
		end

		if wasEOL
			result *= repeat("  ", indents)
		end
		if !(!isParens && wasComma)
			result *= c
		end

		wasEOL = c == '\n'
		wasComma = c == ','
		if !isParens && wasComma
			result *= '\n'
			wasEOL = true
		end
	end

	return Markdown.MD(Markdown.Code(String(language(ctx)), strip(result)))
end

function Markdown.MD(ctx::Context, cursor::CXCursor, comment::CXComment)
	hasDetails = false
	hasParams  = false
	hasReturns = false

	contents = Markdown.MD(ctx, cursor).content
	for ind in 1:clang_Comment_getNumChildren(comment)
		child = clang_Comment_getChild(comment, ind-1)
		kind = clang_Comment_getKind(child)

		if !hasDetails && !isempty(contents)
			hasDetails = true
			push!(contents, Markdown.Header("Details", 2))
		end

		if kind == CXComment_Paragraph
			Bool(clang_Comment_isWhitespace(child)) || push!(contents, Markdown.Paragraph(child))
		elseif kind == CXComment_BlockCommand
			para = clang_BlockCommandComment_getParagraph(child)
			para = Markdown.Paragraph(para)

			num = clang_BlockCommandComment_getNumArgs(child)
			cmd = _string(clang_BlockCommandComment_getCommandName, child)
			if cmd == "brief" || cmd == "par" || cmd == "paragraph"
				push!(contents, para)
			elseif cmd == "note" || cmd == "warning" || cmd == "deprecated"
				para = Markdown.Paragraph(["$(uppercase(cmd)):", para.content...])
				push!(contents, para)
			elseif cmd == "sa" || cmd == "see"
				para = isempty(para.content) ? para : Markdown.Paragraph("See also: [`$(strip(first(para.content)))`](@ref)")
				push!(contents, para)
			elseif cmd == "returns" || cmd == "return"
				if !hasReturns
					hasReturns = true

					hdr = Markdown.Header("Returns", 2)
					if !isempty(contents) && contents[end] isa Markdown.Header
						contents[end] = hdr
					else
						push!(contents, hdr)
					end
				end

				push!(contents, para)
			elseif cmd == "li"
				if !isempty(contents) && contents[end] isa Markdown.List
					push!(contents[end].items, para)
				else
					push!(contents, Markdown.List(para))
				end
			else
				@warn "Unhandled comment block-command comment: $(cmd)"
			end
		elseif kind == CXComment_ParamCommand
			if !hasParams
				hasParams = true

				hdr = Markdown.Header("Parameters", 2)
				if !isempty(contents) && contents[end] isa Markdown.Header
					contents[end] = hdr
				else
					push!(contents, hdr)
				end
			end

			addParameter(contents, child)
		elseif kind == CXComment_VerbatimBlockCommand
			push!(contents, Markdown.Code(child))
		elseif kind == CXComment_VerbatimLine
			text = _string(clang_VerbatimLineComment_getText, child)
			push!(contents, Markdown.Paragraph(text))
		else
			@warn "Unhandled comment document child: $(kind)"
		end
	end

	if length(contents) > 1 && contents[end] isa Markdown.Header
		contents = contents[1:end-1]
	end

	return Markdown.MD(contents)
end

function Markdown.Paragraph(cxcomment::CXComment)
	contents = []

	for ind in 1:clang_Comment_getNumChildren(cxcomment)
		child = clang_Comment_getChild(cxcomment, ind-1)
		kind = clang_Comment_getKind(child)

		if kind == CXComment_Text
			text = _string(clang_TextComment_getText, child)
			push!(contents, text)
		elseif kind == CXComment_InlineCommand
			addInline(contents, child)
		elseif kind == CXComment_HTMLStartTag || kind == CXComment_HTMLEndTag
			# TODO: handle HTML stuff...
		else
			@warn "Unhandled comment paragraph child: $(kind)"
		end
	end

	return Markdown.Paragraph(contents)
end

function Markdown.Code(cxcomment::CXComment)
	lines = []

	for ind in 1:clang_Comment_getNumChildren(cxcomment)
		child = clang_Comment_getChild(cxcomment, ind-1)
		kind = clang_Comment_getKind(child)

		if kind == CXComment_VerbatimBlockLine
			line = _string(clang_VerbatimBlockLineComment_getText, child)

			# a line with newlines is probably incorrectly parsed, so only keep before a newline
			push!(lines, first(split(line, '\n', limit=2)))
		else
			@warn "Unhandled comment code-block child: $(kind)"
		end
	end

	return Markdown.Code(join(lines, '\n'))
end

function addInline(contents, cxcomment)
	num = clang_InlineCommandComment_getNumArgs(cxcomment)
	cmd = _string(clang_InlineCommandComment_getCommandName, cxcomment)

	text = _string(clang_InlineCommandComment_getArgText, cxcomment, 0)
	kind = clang_InlineCommandComment_getRenderKind(cxcomment)
	if kind == CXCommentInlineCommandRenderKind_Bold
		push!(contents, Markdown.Bold(text))
	elseif kind == CXCommentInlineCommandRenderKind_Emphasized
		push!(contents, Markdown.Italic(text))
	elseif kind == CXCommentInlineCommandRenderKind_Monospaced
		push!(contents, Markdown.Code(text))  # TODO: maybe need Code(language, text)
	else
		push!(contents, text)
	end
end

function addParameter(contents, cxcomment)
	num = clang_Comment_getNumChildren(cxcomment)
	num == 1 || error("Incorrect number of parameter-command comment children")

	para = clang_Comment_getChild(cxcomment, 0)
	clang_Comment_getKind(para) == CXComment_Paragraph || error("Expected a parameter-command comment paragraph")

	param = _string(clang_ParamCommandComment_getParamName, cxcomment)
	param = Markdown.Paragraph([Markdown.Code(param), ":", Markdown.Paragraph(para).content...,])

	if !isempty(contents) && contents[end] isa Markdown.List
		push!(contents[end].items, param)
	else
		push!(contents, Markdown.List(param))
	end
end


text:
module c

include("../deps/deps.jl")

export
  # Types
  KCDBPtr, KCCURPtr, CString, KCREC, KCSTR,

  # Error codes
  KCESUCCESS, KCENOIMPL, KCEINVALID, KCENOREPOS, KCENOPERM, KCEBROKEN, KCEDUPREC, KCENOREC, KCELOGIC, KCESYSTEM, KCEMISC,

  # Open modes
  KCOREADER, KCOWRITER, KCOCREATE, KCOTRUNCATE, KCOAUTOTRAN, KCOAUTOSYNC, KCONOLOCK, KCOTRYLOCK, KCONOREPAIR,

  # General functions
  kcfree,

  # DB functions
  kcdbnew, kcdbdel, kcdbopen, kcdbclose, kcdbecode, kcdbemsg, kcdbset, kcdbget,
  kcdbcursor, kcdbcount, kcdbcheck, kcdbclear, kcdbremove, kcdbseize,
  kcdbpath, kcdbcas,
  kcdbsetbulk, kcdbremovebulk,

  # Cursor functions
  kccurdel, kccurjump, kccurstep, kccurget, kccurecode, kccuremsg, kccurdb

# C API types
typealias KCDBPtr Ptr{Void}
typealias KCCURPtr Ptr{Void}
typealias CString Ptr{Uint8}

immutable KCSTR
  buf::CString
  size::Csize_t

  KCSTR(buf::CString, size::Csize_t) = new(buf, size)
  KCSTR(buf, size) = new(convert(CString, buf), convert(Csize_t, size))
  KCSTR(buf::ASCIIString) = new(convert(CString, buf), convert(Csize_t, length(buf)))
end

immutable KCREC
  key::KCSTR
  value::KCSTR
end

# Error codes
const KCESUCCESS = convert(Cint,  0) # success
const KCENOIMPL  = convert(Cint,  1) # not implemented
const KCEINVALID = convert(Cint,  2) # invalid operation
const KCENOREPOS = convert(Cint,  3) # no repository
const KCENOPERM  = convert(Cint,  4) # no permission
const KCEBROKEN  = convert(Cint,  5) # broken file
const KCEDUPREC  = convert(Cint,  6) # record duplication
const KCENOREC   = convert(Cint,  7) # no record
const KCELOGIC   = convert(Cint,  8) # logical inconsistency
const KCESYSTEM  = convert(Cint,  9) # system error
const KCEMISC    = convert(Cint, 15) # miscellaneous error

# Open modes
const KCOREADER   = convert(Uint, 1 << 0) # open as a reader
const KCOWRITER   = convert(Uint, 1 << 1) # open as a writer
const KCOCREATE   = convert(Uint, 1 << 2) # writer creating
const KCOTRUNCATE = convert(Uint, 1 << 3) # writer truncating
const KCOAUTOTRAN = convert(Uint, 1 << 4) # auto transaction
const KCOAUTOSYNC = convert(Uint, 1 << 5) # auto synchronization
const KCONOLOCK   = convert(Uint, 1 << 6) # open without locking
const KCOTRYLOCK  = convert(Uint, 1 << 7) # lock without blocking
const KCONOREPAIR = convert(Uint, 1 << 8) # open without auto repair

# Release a region allocated in the library.
kcfree(ptr) = ccall((:kcfree, libkyotocabinet), Void, (Ptr{Void},), ptr)

# Create a polymorphic database object.
kcdbnew() = ccall((:kcdbnew, libkyotocabinet), KCDBPtr, ())

# Destroy a database object.
function kcdbdel(db::KCDBPtr)
  ccall((:kcdbdel, libkyotocabinet), Void, (KCDBPtr,), db)
end

# Open a database file.
function kcdbopen(db::KCDBPtr, path, mode)
  ccall((:kcdbopen, libkyotocabinet), Cint, (KCDBPtr, CString, Cuint), db, path, mode)
end

# Close the database file.
function kcdbclose(db::KCDBPtr)
  ccall((:kcdbclose, libkyotocabinet), Cint, (KCDBPtr,), db)
end

# Get the code of the last happened error.
function kcdbecode(db::KCDBPtr)
  ccall((:kcdbecode, libkyotocabinet), Cint, (KCDBPtr,), db)
end

# Get the supplement message of the last happened error.
function kcdbemsg(db::KCDBPtr)
  ccall((:kcdbemsg, libkyotocabinet), CString, (KCDBPtr,), db)
end

# Set the value of a record.
function kcdbset(db::KCDBPtr, kbuf, ksize, vbuf, vsize)
  ccall((:kcdbset, libkyotocabinet), Cint, (KCDBPtr, CString, Cuint, CString, Cuint),
    db, kbuf, ksize, vbuf, vsize)
end

# Retrieve the value of a record.
function kcdbget(db::KCDBPtr, kbuf, ksize, vsize)
  ccall((:kcdbget, libkyotocabinet), CString, (KCDBPtr, CString, Cuint, Ptr{Cuint}),
    db, kbuf, ksize, vsize)
end

# Check the existence of a record.
function kcdbcheck(db::KCDBPtr, kbuf, ksize)
  ccall((:kcdbcheck, libkyotocabinet), Cint, (KCDBPtr, CString, Csize_t),
    db, kbuf, ksize)
end

# Get the number of records.
function kcdbcount(db::KCDBPtr)
  ccall((:kcdbcount, libkyotocabinet), Clonglong, (KCDBPtr,), db)
end

# Remove all records.
function kcdbclear(db::KCDBPtr)
  ccall((:kcdbclear, libkyotocabinet), Cint, (KCDBPtr,), db)
end

# Remove a record.
function kcdbremove(db::KCDBPtr, kbuf, ksize)
  ccall((:kcdbremove, libkyotocabinet), Cint, (KCDBPtr, CString, Csize_t),
    db, kbuf, ksize)
end

# Retrieve the value of a record and remove it atomically.
function kcdbseize(db::KCDBPtr, kbuf, ksize, vsize)
  ccall((:kcdbseize, libkyotocabinet), CString, (KCDBPtr, CString, Cuint, Ptr{Cuint}),
    db, kbuf, ksize, vsize)
end

# Get the path of the database file.
function kcdbpath(db::KCDBPtr)
  ccall((:kcdbpath, libkyotocabinet), CString, (KCDBPtr,), db)
end

# Perform compare-and-swap.
function kcdbcas(db::KCDBPtr, kbuf, ksize, ovbuf, ovsize, nvbuf, nvsize)
  ccall((:kcdbcas, libkyotocabinet), Cint,
    (KCDBPtr, CString, Csize_t, CString, Csize_t, CString, Csize_t),
    db, kbuf, ksize, ovbuf, ovsize, nvbuf, nvsize)
end

# Store records at once.
function kcdbsetbulk(db::KCDBPtr, recs, rnum, atomic)
  ccall((:kcdbsetbulk, libkyotocabinet), Int64,
    (KCDBPtr, Ptr{KCREC}, Csize_t, Cint),
    db, convert(Ptr{KCREC}, pointer(recs)), rnum, atomic)
end

# Remove records at once.
function kcdbremovebulk(db::KCDBPtr, keys, knum, atomic)
  ccall((:kcdbremovebulk, libkyotocabinet), Int64,
    (KCDBPtr, Ptr{KCSTR}, Csize_t, Cint),
    db, convert(Ptr{KCSTR}, pointer(keys)), knum, atomic)
end

# Create a polymorphic cursor object.
function kcdbcursor(db::KCDBPtr)
  ccall((:kcdbcursor, libkyotocabinet), KCCURPtr, (KCDBPtr,), db)
end

# Destroy a cursor object.
function kccurdel(cursor::KCCURPtr)
  ccall((:kccurdel, libkyotocabinet), Void, (KCCURPtr,), cursor)
end

# Get a pair of the key and the value of the current record.
function kccurget(cursor::KCCURPtr, ksize, v, vsize, step)
  ccall((:kccurget, libkyotocabinet), CString,
    (KCCURPtr, Ptr{Cuint}, Ptr{CString}, Ptr{Cuint}, Cint),
    cursor, ksize, v, vsize, step)
end

# Jump the cursor to the first record for forward scan.
function kccurjump(cursor::KCCURPtr)
  ccall((:kccurjump, libkyotocabinet), Cint, (KCCURPtr,), cursor)
end

# Step the cursor to the next record.
function kccurstep(cursor::KCCURPtr)
  ccall((:kccurstep, libkyotocabinet), Cint, (KCCURPtr,), cursor)
end

# Get the database object.
function kccurdb(cursor::KCCURPtr)
  ccall((:kccurdb, libkyotocabinet), KCDBPtr, (KCCURPtr,), cursor)
end

# Get the code of the last happened error.
function kccurecode(cursor::KCCURPtr)
  ccall((:kccurecode, libkyotocabinet), Cint, (KCCURPtr,), cursor)
end

# Get the supplement message of the last happened error.
function kccuremsg(cursor::KCCURPtr)
  ccall((:kccuremsg, libkyotocabinet), CString, (KCCURPtr,), cursor)
end

end # module c


text:
function add(obj::PolynomialFunctionMatrix, arg0::PolynomialFunctionMatrix)
    return jcall(obj, "add", PolynomialFunctionMatrix, (PolynomialFunctionMatrix,), arg0)
end

function get_elem(obj::PolynomialFunctionMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getElem", PolynomialFunction, (jint, jint), arg0, arg1)
end

function get_matrix_line(obj::PolynomialFunctionMatrix, arg0::jint)
    return jcall(obj, "getMatrixLine", Vector{PolynomialFunction}, (jint,), arg0)
end

function multiply(obj::PolynomialFunctionMatrix, arg0::PolynomialFunctionMatrix)
    return jcall(obj, "multiply", PolynomialFunctionMatrix, (PolynomialFunctionMatrix,), arg0)
end

function set_elem(obj::PolynomialFunctionMatrix, arg0::jint, arg1::jint, arg2::PolynomialFunction)
    return jcall(obj, "setElem", void, (jint, jint, PolynomialFunction), arg0, arg1, arg2)
end

function set_matrix(obj::PolynomialFunctionMatrix, arg0::Vector{Vector{PolynomialFunction}})
    return jcall(obj, "setMatrix", void, (Vector{Vector{PolynomialFunction}},), arg0)
end

function set_matrix_line(obj::PolynomialFunctionMatrix, arg0::jint, arg1::Vector{PolynomialFunction})
    return jcall(obj, "setMatrixLine", void, (jint, Vector{PolynomialFunction}), arg0, arg1)
end


text:
# This script runs all the PolyhedralRelaxation examples.
# Usage: julia --project=examples run_examples.jl

using Test

const EXAMPLES = filter(
    ex -> endswith(ex, ".jl") && ex != "run_examples.jl",
    readdir(@__DIR__),
)

@testset "run_examples.jl" begin
    @testset "$(example)" for example in EXAMPLES
        include(example)
    end
end


text:
function convertshape(box::ConstrainedDynamics.Box)
    x,y,z = Tuple(box.xyz)
    return GeometryBasics.HyperRectangle(Vec(-x/2,-y/2,-z/2),Vec(x,y,z))
end

function convertshape(cylinder::ConstrainedDynamics.Cylinder)
    r,h = Tuple(cylinder.rh)
    return GeometryBasics.Cylinder(Point(0.0,0.0,-h/2),Point(0.0,0.0,h/2), r)
end

function convertshape(sphere::ConstrainedDynamics.Sphere)
    r = sphere.r
    return GeometryBasics.Sphere(Point(0.0,0.0,0.0), r)
end

function convertshape(pyramid::ConstrainedDynamics.Pyramid)
    w, h = Tuple(pyramid.wh)
    return GeometryBasics.Pyramid(Point(0.0,0.0,-h/4), h, w)
end

function convertshape(mesh::ConstrainedDynamics.Mesh)
    obj = MeshFileObject(mesh.path)
    if obj.mtl_library == "" # no material file (mtl) available -> render as geometry with uniform color
        return MeshFileGeometry(obj.contents, obj.format)
    end
    return obj
end

function convertshape(::ConstrainedDynamics.EmptyShape)
    return nothing
end


text:
"""
The parent type of all graphs.
"""
abstract type Graph end

function get_source_calculator_wrapper(graph::Graph)::CalculatorWrapper
    # FIXME: Change name to wrapper
    graph.source_calculator
end

function get_calculator_wrappers(graph::Graph)
    graph.calculator_wrappers
end

function get_cw_tasks(graph::Graph)
    graph.cw_tasks
end

function set_cw_tasks(graph::Graph, cw_tasks::Vector{Task})
    graph.cw_tasks = cw_tasks
end

function get_next_timestamp(graph::Graph)
    graph.next_timestamp
end

function set_next_timestamp(graph::Graph, next_ts)
    graph.next_timestamp = next_ts
end

function get_sinks_lock(graph::Graph)
    graph.sinks_lock
end

function is_last_sink_not_init(graph::Graph)
    is_last = nothing
    lk = get_sinks_lock(graph)

    lock(lk)
    # TODO: Use getter
    graph.num_sinks_not_init -= 1
    if graph.num_sinks_not_init == 0
        is_last = true
    else
        is_last = false
    end
    unlock(lk)

    return is_last
end

function get_num_sinks_not_init(graph::Graph)
    lk = get_sinks_lock(graph)
    lock(lk)
    ret = graph.num_sinks_not_init
    unlock(lk)

    ret
end

function get_flow_limiter_bootstrap_lock(graph::Graph)
    graph.bootstrap_lock
end

function get_flow_limiter_bootstrap_cond(graph::Graph)
    graph.bootstrap_cond
end

"""
    stop(graph)

Stop the graph gracefully.

The graph will stop pulling new data packets out of the *source* calculator, and all the tasks running calculators will exit after they are done processing the last generated packet.

[`Tissue.wait_until_done(graph)`](@ref) can be called in the main thread to block it until every calculator is fully terminated and closed.

# Examples
```julia
using Tissue
using MyGraphs: CoolestGraph

graph = CoolestGraph()

start(graph)
sleep(5)
stop(graph)
wait_until_done(graph)
```
"""
function stop(graph::Graph)
    graph.done[] = true
end

function is_done(graph::Graph)
    graph.done[]
end

"""
Gets the period (in the sense of inverse frequency) at which the source
calculator should be called to generate new graph input stream packets.
"""
function get_source_period(graph::Graph)::Float64
    graph.gen_period[]
end

function set_source_period(graph::Graph, period::Float64)
    MIN_SLEEPABLE_TIME = 0.001
    graph.gen_period[] = max(MIN_SLEEPABLE_TIME, period)
end

"""
    wait_until_done(graph)

Block the main thread until the graph is done.

Can only be called after [`start(graph)`](@ref) was called. This waits for all the calculators to be finished processing the last packet, and calls [`close(calculator)`](@ref) for each calculator. This can occur either because [`stop(graph)`](@ref) was called, or the source calculator indicated that it is done generating data by returning `nothing`.

"""
function wait_until_done(graph::Graph)
    # Wait on all tasks
    for calculator_wrapper_task in get_cw_tasks(graph)
        wait(calculator_wrapper_task)
    end

    set_cw_tasks(graph, Vector{Task}())

    # cleanup
    # TODO: do this concurrently
    for calculator in map(get_calculator, get_calculator_wrappers(graph))
        close(calculator)
    end
end


text:
#### General Assumptions of these OTS Models ####
#
# - if the branch status is 0 in the input, it is out of service and forced to 0 in OTS
# - the network will be maintained as one connected component (i.e. at least n-1 edges)
#

export run_ots

""
function run_ots(file, model_constructor, solver; kwargs...)
    return run_generic_model(file, model_constructor, solver, post_ots; solution_builder = get_ots_solution, kwargs...)
end

""
function post_ots(pm::GenericPowerModel)
    variable_branch_indicator(pm)
    variable_voltage_on_off(pm)
    variable_generation(pm)
    variable_branch_flow(pm)
    variable_dcline_flow(pm)

    objective_min_fuel_cost(pm)

    constraint_voltage_on_off(pm)

    for i in ids(pm, :ref_buses)
        constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_kcl_shunt(pm, i)
    end

    for i in ids(pm, :branch)
        constraint_ohms_yt_from_on_off(pm, i)
        constraint_ohms_yt_to_on_off(pm, i)

        constraint_voltage_angle_difference_on_off(pm, i)

        constraint_thermal_limit_from_on_off(pm, i)
        constraint_thermal_limit_to_on_off(pm, i)
    end

    for i in ids(pm, :dcline)
        constraint_dcline(pm, i)
    end
end

""
function get_ots_solution(pm::GenericPowerModel, sol::Dict{String,Any})
    add_bus_voltage_setpoint(sol, pm)
    add_generator_power_setpoint(sol, pm)
    add_branch_flow_setpoint(sol, pm)
    add_branch_status_setpoint(sol, pm)
end


text:
using Polyhedra
using CDDLib
using GLVisualize
using GeometryTypes

ext1 = SimpleVRepresentation([0 1 2;0 2 1; 1 0 2; 1 2 0; 2 0 1; 2 1 0])
ext2 = SimpleVRepresentation([0 0 0;0 0 1; 0 1 0; 1 0 0])
ext3 = SimpleVRepresentation([0 0 0;0 0 1; 0 1 0; 1 0 0], IntSet([1])) # BUG in CDD : 0x + 0y + 0z <= 1 -> useless
poly1 = CDDPolyhedron{3,Rational{BigInt}}(ext1)
poly2 = CDDPolyhedron{3,Rational{BigInt}}(ext2)
poly3 = CDDPolyhedron{3,Rational{BigInt}}(ext3)

window = glscreen()

#view(visualize(GeometryTypes.GLNormalMesh(poly1, RGBA(1f0,0f0,0f0,1f0))), window)
view(visualize(GLNormalMesh(poly1)), window)
#view(visualize(poly2), window)
#view(visualize(GeometryTypes.GLNormalMesh(poly3, RGBA(0f0,1f0,0f0,1f0))), window)
view(visualize(GLNormalMesh(poly3)), window)

include("axes.jl")

renderloop(window)


text:
function PiGPIOC.gpioWrite(vector::AbstractVector{Int}, data::Int)
    for x in vector
        x >= 0 && PiGPIOC.gpioWrite(x, data)
    end
end

function PiGPIOC.gpioWrite(vector::NTuple{8,Int}, data::Int)
    for x in vector
        x >= 0 && PiGPIOC.gpioWrite(x, data)
    end
end

function PiGPIOC.gpioWrite(vector::NTuple{5,Int}, data::Int)
    for x in vector
        x >= 0 && PiGPIOC.gpioWrite(x, data)
    end
end

function PiGPIOC.gpioSetMode(vector::AbstractVector{Int}, data::Int)
    for x in vector
        x >= 0 && PiGPIOC.gpioSetMode(x, data)
    end
end

function PiGPIOC.gpioSetMode(vector::NTuple{8,Int}, data::Int)
    for x in vector
        x >= 0 && PiGPIOC.gpioSetMode(x, data)
    end
end

function PiGPIOC.gpioSetMode(vector::NTuple{5,Int}, data::Int)
    for x in vector
        x >= 0 && PiGPIOC.gpioSetMode(x, data)
    end
end

# stop previous (if active) and run
function run_wave(wave_id::Int)
    # get old
    old_wave_id = PiGPIOC.gpioWaveTxAt()

    # run new wave when possible
    ret_code = PiGPIOC.gpioWaveTxSend(wave_id, PiGPIOC.PI_WAVE_MODE_REPEAT_SYNC)
    if ret_code < 0
        # return the number of DMA control blocks in the waveform if OK,
        # otherwise PI_BAD_WAVE_ID, or PI_BAD_WAVE_MODE
        throw("Error in 'PiGPIOC.gpioWaveCreate()' with code: $(ret_code)")
    end

    # safe delete old wave
    no_old_wave = old_wave_id == PiGPIOC.PI_NO_TX_WAVE || # 9999
        old_wave_id == PiGPIOC.PI_WAVE_NOT_FOUND # 9998
    if !no_old_wave
        while PiGPIOC.gpioWaveTxAt() != wave_id
            # wait for starting new wave
        end
        PiGPIOC.gpioWaveDelete(old_wave_id)
    end

    return nothing
end


text:
@testset "Spin-1/2 Fock" begin
    # test six states with N=2
    # (UP and DOWN)
    n = [([1, 1], [0, 0]), ([0, 1], [0, 1]), ([1, 0], [0, 1]), ([0, 1], [1, 0]), ([1, 0], [1, 0]), ([0, 0], [1, 1])]
    idx = [3, 5, 6, 9, 10, 12]

    for i in 1:length(n)
    @test Hilbert.state2idx(n[i]...) == idx[i]
    @test Hilbert.idx2state(idx[i], 2) == n[i]
    @test Hilbert.occupation(idx[i], 2, 1, UP) == n[i][1][1]
    @test Hilbert.occupation(idx[i], 2, 1, DOWN) == n[i][2][1]
    @test Hilbert.occupation(idx[i], 2, 2, UP) == n[i][1][2]
    @test Hilbert.occupation(idx[i], 2, 2, DOWN) == n[i][2][2]
    @test Hilbert.occupationTotal(idx[i], 2) == (sum(n[i][1]), sum(n[i][2]))

end
    basis = Hilbert.BinaryFock(2, 2, :all)
    @test basis.idx == idx
    # basis = BinaryFock(1, :all, :all)
    # println(creation(basis, 1, UP))
    # println(creation(basis, 1, UP))
end


text:
function save_basic_features_instance!(db::SQLite.DB, name, scenario, source_type, source_path)
    println("Saving basic instance features: ($name, $scenario)")
    _nbus = nbus(source_path, source_type)
    _nbranch_unique = nbranch(source_path, source_type; distinct_pair=true)
    _nbranch = nbranch(source_path, source_type)
    _ngen = ngen(source_path, source_type)
    query = """
    UPDATE instances
    SET nbus = $_nbus, nbranch_unique = $_nbranch_unique, nbranch = $_nbranch, ngen = $_ngen
    WHERE name = '$name' AND scenario = $scenario;
    """
    execute_query(db, query)
end

function save_basic_features_instance_dfrow!(db::SQLite.DB, row)
    save_basic_features_instance!(db, row[:name], row[:scenario], row[:source_type], row[:source_path])
end

function save_basic_features_instances!(db::SQLite.DB; recompute=false, subset=nothing)
    query = "SELECT name, scenario, source_type, source_path FROM instances"
    if !recompute
        query *= " WHERE (nbus IS NULL OR nbranch_unique IS NULL OR nbranch IS NULL OR ngen IS NULL)"
    end
    if !isnothing(subset)
        if !recompute
            query *= " AND id IN ($(join(subset, ',')))"
        else
            query *= " WHERE id IN ($(join(subset, ',')))"
        end
    end
    results = DBInterface.execute(db, query) |> DataFrame
    println("Saving basic instance features config: recompute=$recompute")
    println("subset\n$subset\nend subset")
    save_function!(row) = save_basic_features_instance_dfrow!(db, row)
    save_function!.(eachrow(results[!, [:name, :scenario, :source_type, :source_path]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
@testset "Test 3D" begin

    # function get_3D_model(p::Vector{Float64})
    #     r = p[1] + p[2] + p[3]
    #     return r
    # end
    #
    # function get_3D_model_grad()
    #     return [1.0; 1.0; 1.0]
    # end

    p = [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0] # function nodes
    u = [1.0; 1.0; 1.0; 0.0]                                # function values in nodes
    n_1 = size(p, 2)
    dp = Matrix{Float64}(undef, 3, 3*n_1)
    es = Matrix{Float64}(undef, 3, 3*n_1)
    du = Vector{Float64}(undef, 3*n_1)
    grad = [1.0; 1.0; 1.0]

    k = 0
    for i = 1:n_1
        k += 1
        dp[1,k] = p[1,i]
        dp[2,k] = p[2,i]
        dp[3,k] = p[3,i]
        du[k] = grad[1]
        es[1,k] = 1.0
        es[2,k] = 0.0
        es[3,k] = 0.0
        k += 1
        dp[1,k] = p[1,i]
        dp[2,k] = p[2,i]
        dp[3,k] = p[3,i]
        du[k] = grad[2]
        es[1,k] = 0.0
        es[2,k] = 1.0
        es[3,k] = 0.0
        k += 1
        dp[1,k] = p[1,i]
        dp[2,k] = p[2,i]
        dp[3,k] = p[3,i]
        du[k] = grad[3]
        es[1,k] = 0.0
        es[2,k] = 0.0
        es[3,k] = 1.0
    end
    ####

    @testset "Test 3D-RK_H0 kernel" begin
        s = interpolate(p, u, RK_H0(0.00001))
        σ1 = evaluate_one(s, [1.0; 0.0; 0.01])
        @test isapprox(σ1, u[1], atol = 1e-3)
        σ = evaluate(s, [1.0 0.5; 0.0 0.5; 0.01 0.5])
        @test isapprox(σ[1], u[1], atol = 1e-3)

        s = prepare(p, RK_H0(0.00001)) # prepare spline
        s = construct(s, u) # construct spline
        σ1 = evaluate_one(s, [1.0; 0.0; 0.01])
        @test isapprox(σ1, u[1], atol = 1e-3)
        σ = evaluate(s, [1.0 0.5; 0.0 0.5; 0.01 0.5])
        @test isapprox(σ[1], u[1], atol = 1e-3)

        cond = estimate_cond(s)
        @test cond ≈ 1.0e6

        iq = estimate_accuracy(s)
        @test iq ≈ 15.0

        eps = get_epsilon(s)
        @test eps ≈ 1.0e-5

        est_eps = estimate_epsilon(p)
        @test isapprox(est_eps, 1.036, atol = 1e-3)
    end

    @testset "Test 3D-RK_H1 kernel" begin
        s = interpolate(p, u, dp, es, du, RK_H1(0.1))
        σ1 = evaluate_one(s, [1.0; 0.0; 0.01])
        @test isapprox(σ1, u[1], atol = 1e-2)
        σ = evaluate(s, [1.0 0.5; 0.0 0.5; 0.01 0.5])
        @test isapprox(σ[1], u[1], atol = 1e-2)

        s = prepare(p, dp, es, RK_H1(0.1)) # prepare spline
        s = construct(s, u, du) # construct spline
        σ1 = evaluate_one(s, [1.0; 0.0; 0.01])
        @test isapprox(σ1, u[1], atol = 1e-2)
        σ = evaluate(s, [1.0 0.5; 0.0 0.5; 0.01 0.5])
        @test isapprox(σ[1], u[1], atol = 1e-2)

        g = evaluate_gradient(s, [0.01; 0.2; 0.3])
        @test all(isapprox.(g, grad, atol = 0.01))

        cond = estimate_cond(s)
        @test cond ≈ 1.0e5

        iq = estimate_accuracy(s) + 1.0
        @test iq ≈ 14.0 atol = 1.0

        eps = get_epsilon(s)
        @test eps ≈ 0.1

        est_eps = estimate_epsilon(p)
        @test isapprox(est_eps, 1.036, atol = 1e-3)

        est_eps = estimate_epsilon(p, dp)
        @test isapprox(est_eps, 1.415, atol = 1e-3)
    end

    @testset "Test 3D-RK_H2 kernel" begin
        s = interpolate(p, u, dp, es, du, RK_H2(0.1))
        σ1 = evaluate_one(s, [1.0; 0.0; 0.01])
        @test isapprox(σ1, u[1], atol = 2e-2)
        σ = evaluate(s, [1.0 0.5; 0.0 0.5; 0.01 0.5])
        @test isapprox(σ[1], u[1], atol = 2e-2)

        s = prepare(p, dp, es, RK_H2(0.1)) # prepare spline
        s = construct(s, u, du) # construct spline
        σ1 = evaluate_one(s, [1.0; 0.0; 0.01])
        @test isapprox(σ1, u[1], atol = 2e-2)
        σ = evaluate(s, [1.0 0.5; 0.0 0.5; 0.01 0.5])
        @test isapprox(σ[1], u[1], atol =2e-2)

        g = evaluate_gradient(s, [0.01; 0.2; 0.3])
        @test all(isapprox.(g, grad, atol = 0.01))

        cond = estimate_cond(s)
        @test cond ≈ 1.0e8

        iq = estimate_accuracy(s) + 1.0
        @test iq ≈ 13.0 atol = 1.0

        eps = get_epsilon(s)
        @test eps ≈ 0.1

        est_eps = estimate_epsilon(p)
        @test isapprox(est_eps, 1.036, atol = 1e-3)

        est_eps = estimate_epsilon(p, dp)
        @test isapprox(est_eps, 1.415, atol = 1e-3)
    end

end


text:
using Test

include("example.jl")
# add here the file with your unit tests


text:
using Base: min, max

export
    Polygon,
    add_vertex!, build!, is_point_inside

import .Geometry.add_vertex!
import .Geometry.build!

mutable struct Polygon{T <: AbstractFloat}
    mesh::Mesh

    function Polygon{T}() where {T <: AbstractFloat}
        new(Mesh())
    end
end

function add_vertex!(poly::Polygon{T}, x::T, y::T) where {T <: AbstractFloat}
    add_vertex!(poly.mesh, x, y)
end

function build!(poly::Polygon)
    build!(poly.mesh)
end

function is_point_inside(poly::Polygon{T}, p::Point{T}) where {T <: AbstractFloat}
    i = 1
    c = false
    vertices = poly.mesh.vertices
    nvert = length(vertices)
    j = 2

    while i < nvert
        if ((vertices[i].y > p.y) ≠ (vertices[j].y > p.y)) &&
            (p.x < (vertices[j].x - vertices[i].x) * (p.y - vertices[i].y) /
              (vertices[j].y - vertices[i].y) + vertices[i].x)
            c = !c
        end
        i += 1
        j += 1
    end

    # Last edge to close loop
    i = j - 1
    j = 1
    if ((vertices[i].y > p.y) ≠ (vertices[j].y > p.y)) &&
        (p.x < (vertices[j].x - vertices[i].x) * (p.y - vertices[i].y) /
          (vertices[j].y - vertices[i].y) + vertices[i].x)
        c = !c
    end

    c
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Test
using Ceres
@testset "Fit curve test" begin
  data = [
    0.000000e+00, 1.133898e+00,
    7.500000e-02, 1.334902e+00,
    1.500000e-01, 1.213546e+00,
    2.250000e-01, 1.252016e+00,
    3.000000e-01, 1.392265e+00,
    3.750000e-01, 1.314458e+00,
    4.500000e-01, 1.472541e+00,
    5.250000e-01, 1.536218e+00,
    6.000000e-01, 1.355679e+00,
    6.750000e-01, 1.463566e+00,
    7.500000e-01, 1.490201e+00,
    8.250000e-01, 1.658699e+00,
    9.000000e-01, 1.067574e+00,
    9.750000e-01, 1.464629e+00,
    1.050000e+00, 1.402653e+00,
    1.125000e+00, 1.713141e+00,
    1.200000e+00, 1.527021e+00,
    1.275000e+00, 1.702632e+00,
    1.350000e+00, 1.423899e+00,
    1.425000e+00, 1.543078e+00,
    1.500000e+00, 1.664015e+00,
    1.575000e+00, 1.732484e+00,
    1.650000e+00, 1.543296e+00,
    1.725000e+00, 1.959523e+00,
    1.800000e+00, 1.685132e+00,
    1.875000e+00, 1.951791e+00,
    1.950000e+00, 2.095346e+00,
    2.025000e+00, 2.361460e+00,
    2.100000e+00, 2.169119e+00,
    2.175000e+00, 2.061745e+00,
    2.250000e+00, 2.178641e+00,
    2.325000e+00, 2.104346e+00,
    2.400000e+00, 2.584470e+00,
    2.475000e+00, 1.914158e+00,
    2.550000e+00, 2.368375e+00,
    2.625000e+00, 2.686125e+00,
    2.700000e+00, 2.712395e+00,
    2.775000e+00, 2.499511e+00,
    2.850000e+00, 2.558897e+00,
    2.925000e+00, 2.309154e+00,
    3.000000e+00, 2.869503e+00,
    3.075000e+00, 3.116645e+00,
    3.150000e+00, 3.094907e+00,
    3.225000e+00, 2.471759e+00,
    3.300000e+00, 3.017131e+00,
    3.375000e+00, 3.232381e+00,
    3.450000e+00, 2.944596e+00,
    3.525000e+00, 3.385343e+00,
    3.600000e+00, 3.199826e+00,
    3.675000e+00, 3.423039e+00,
    3.750000e+00, 3.621552e+00,
    3.825000e+00, 3.559255e+00,
    3.900000e+00, 3.530713e+00,
    3.975000e+00, 3.561766e+00,
    4.050000e+00, 3.544574e+00,
    4.125000e+00, 3.867945e+00,
    4.200000e+00, 4.049776e+00,
    4.275000e+00, 3.885601e+00,
    4.350000e+00, 4.110505e+00,
    4.425000e+00, 4.345320e+00,
    4.500000e+00, 4.161241e+00,
    4.575000e+00, 4.363407e+00,
    4.650000e+00, 4.161576e+00,
    4.725000e+00, 4.619728e+00,
    4.800000e+00, 4.737410e+00,
    4.875000e+00, 4.727863e+00,
    4.950000e+00, 4.669206e+00];

  function mycostfunction(user_data,
                          parameters)
      x = user_data[1];
      y = user_data[2];
      m = parameters[1][1];
      c = parameters[1][2];
      # show(user_data)
      #res = Tracker.collect([(y - exp(m * x + c))])
      res = (y .- exp(m .* x .+ c))

      return [res]
  end

  kNumObservations = 67
  params = [[0.0, 0.0]]

  problem = Problem()

  # function AddResidualBlock!(problem::Problem,
  #                            cost_function_jl,
  #                            cost_function_data::Vector{Float64},
  #                            params::Vector{Vector{Float64}},
  #                            residuals_num::Int64)

  for idx = 1:kNumObservations
    # costf = CostFunction(mycostfunction, copy(data[idx*2-1:idx*2]), 1, [2])
    AddResidualBlock!(problem, mycostfunction, copy(data[idx*2-1:idx*2]), params, 1)
  end

  solve(problem)

  @test params ≈ [[0.291861, 0.131439]] atol=1e-4

  function mycostfunction2(user_data,
                          parameters)
      x = user_data[1];
      y = user_data[2];
      m = parameters[1][1];
      c = parameters[2][1];
      # show(user_data)
      #res = Tracker.collect([(y - exp(m * x + c))])
      res = (y .- exp(m .* x .+ c))

      return [res]
  end

  kNumObservations = 67
  params = [[0.0], [0.0]]

  problem = Problem()

  for idx = 1:kNumObservations
    # costf = CostFunction(mycostfunction2, copy(data[idx*2-1:idx*2]), 1, [1, 1])
    AddResidualBlock!(problem, mycostfunction2, copy(data[idx*2-1:idx*2]), params, 1)
  end

  solve(problem)

  @test params ≈ [[0.291861], [0.131439]] atol=1e-4

  function mycostfunction3(user_data,
                          parameters)
      x = user_data[1];
      y = user_data[2];
      m = parameters[1];
      c = parameters[2];
      # show(user_data)
      #res = Tracker.collect([(y - exp(m * x + c))])
      res1 = (y .- exp(m[1] .* x .+ c[1]))
      res2 = (y .- exp((m[2] * 10.0) .* x .+ (c[2]-0.5)))

      return [res1, res2]
  end

  kNumObservations = 67
  params = [[0.0, 0.0], [0.0, 0.0]]

  loss_func = HuberLoss(5.0)

  problem = Problem()

  for idx = 1:kNumObservations
    # costf = CostFunction(mycostfunction3, copy(data[idx*2-1:idx*2]), 2, [2, 2])
    AddResidualBlock!(problem, mycostfunction3, copy(data[idx*2-1:idx*2]), params, 2,
        loss_function=loss_func
    )
  end

  solve(problem)

  @test params ≈ [[0.291861, 0.291861 / 10.0], [0.131439, 0.131439 + 0.5]] atol=1e-4
end


text:
using Test
using IncrementalInference

@testset "Testing getCliqVarInitOrderUp" begin

fg = generateCanonicalFG_lineStep(3;
                                  poseEvery=1,
                                  landmarkEvery=5,
                                  posePriorsAt=[0],
                                  landmarkPriorsAt=[0,2],
                                  sightDistance=3)

fg.solverParams.useMsgLikelihoods = true
# addVariable!(subfg, :x0, Con)

@test getCliqVarInitOrderUp(fg) == [:x0, :lm0, :x3, :x2, :x1]

solveTree!(fg)

# construct a message on :x3
msg = LikelihoodMessage(status=IIF.UPSOLVED)
seps = [:x3]
for vid in seps
  var = DFG.getVariable(fg, vid)
  if isInitialized(var)
    msg.belief[var.label] = TreeBelief(var)
  end
end
addMsgFactors!(fg, msg, IIF.UpwardPass)

# less sensitive to order given msgjoint message updates
@test intersect(getCliqVarInitOrderUp(fg), [:x3, :x0, :lm0, :x2, :x1]) |> length == 5

# construct a message on :lm0,:x2
msg = LikelihoodMessage(status=IIF.UPSOLVED, hasPriors=false)
seps = [:lm0, :x2]
for vid in seps
    var = DFG.getVariable(fg, vid)
    if isInitialized(var)
        msg.belief[var.label] = TreeBelief(var)
    end
end
addMsgFactors!(fg, msg, IIF.UpwardPass)

@warn("disabled getCliqVarInitOrderUp test until #1010 is closed")
# @test getCliqVarInitOrderUp(fg) == [:x3, :x0, :lm0, :x1, :x2]

# test order with mixture prior #998
cv = 3.0
doorPrior = Mixture(Prior,
                    [Normal(-100,cv);Normal(0,cv);Normal(100,cv);Normal(300,cv)],
                    [1/4;1/4;1/4;1/4] )

fg = initfg()
fg.solverParams.graphinit = false
addVariable!(fg, :x0, ContinuousScalar)

addFactor!(fg, [:x0], doorPrior)

addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg,[:x0; :x1], LinearRelative( Normal(0.0,1.0)))

@test lsfPriors(fg) == [:x0f1]

@test getCliqVarInitOrderUp(fg) == [:x0, :x1]

end


text:
using Base
using LinearAlgebra
using SparseArrays

function _expmat2(mat)
    res = one(mat)
    resprec = zero(mat)
    mult = one(mat)
    i = 1
    cpt = 0
    while cpt < 4
        resprec .= res
        mult *= mat
        mult /= i
        res += mult
        i += 1
        cpt = resprec == res ? cpt + 1 : 0
    end
    return res
end
function _expmat1(mat)
    valnorm = norm(mat)

    if valnorm > 1
        b = mat / valnorm
        coef_int = Integer(floor(valnorm))
        coef_mantissa = valnorm - coef_int
        return _expmat2(b)^coef_int * _expmat2(coef_mantissa * b)
    else
        return _expmat2(mat)
    end
end
function _expmat0(mat)
    res = setprecision(precision(BigFloat) + 32) do
        _expmat1(mat)
    end
    return BigFloat.(res)
end
Base.exp(mat::Array{Complex{BigFloat},2}) = _expmat0(mat)
Base.exp(mat::Array{BigFloat,2}) = _expmat0(mat)
Base.exp(mat::Array{Integer,2}) = _expmat1(mat)
Base.exp(mat::SparseMatrixCSC{Complex{BigFloat},Int64}) = _expmat0(mat)
Base.exp(mat::SparseMatrixCSC{BigFloat,Int64}) = _expmat0(mat)
Base.exp(mat::SparseMatrixCSC{Float64,Int64}) = _expmat1(mat)
Base.exp(mat::Array{Rational,2}) = Base.exp(float(mat))
Base.exp(mat::SparseMatrixCSC{Rational,Integer}) = Base.exp(float(mat))


text:
using Test
using UnPack
using MOONS.Grids
const G = Grids

@testset "Coordinates init" begin
  for FT in (Float32, Float64)
    a = FT(0.0)
    b = FT(1.0)
    n = 10
    c = Coordinates(a,b,n)
    @test c.c.Δh[1] ≈ FT(1/10)
    @test c.n.Δh[1] ≈ FT(1/10)
    @test h_extreme(c.n, Min()) ≈ FT(0)
    @test h_extreme(c.n, Max()) ≈ FT(1)
    sprint(show, c)
    @test c.c == coords(c, Center1D())
    @test c.n == coords(c, Vertex1D())
    # FT==Float64 && @show c
  end
end

@testset "Helper funcs" begin
    @test n_hat(Min()) == -1
    @test binary(Min()) == 0
    @test ghost_vec(Min()) == (1, 0, 0)
    @test ghost_dual(Min()) == (1, 0)

    @test n_hat(Max()) == 1
    @test binary(Max()) == 1
    @test ghost_vec(Max()) == (0, 0, 1)
    @test ghost_dual(Max()) == (0, 1)
    @test orthogs(1) == (2, 3)
    @test orthogs(2) == (3, 1)
    @test orthogs(3) == (1, 2)
end

@testset "Stretch funcs" begin
    FT = Float64
    Re = FT(100)
    Ha = FT(200)
    hmin = FT(0)
    hmax = FT(1)
    δ = hmax*99/100
    @test all(G.ReynoldsBL(Re, FT[0, 0, 0], FT[1, 1, 1]) .≈ [1.0540925533894598, 1.0540925533894598, 1.0540925533894598])
    @test all(G.Re_Ha_BL(Re, Ha, FT[0, 0, 0], FT[1, 1, 1]) .≈ [1.002509414234171, 1.002509414234171, 1.002509414234171])
    @test G.Re_Ha_BL_1D(Re, Ha, hmin, hmax) ≈ 1.002509414234171
    @test G.Reynolds_BL_1D(Re, hmin, hmax) ≈ 1.0540925533894598
    @test G.Hartmann_BL_1D(Ha, hmin, hmax) ≈ 1.002509414234171
    @test G.RobertsBL1D(δ, hmax) ≈ FT(1000)
    @test G.β_Δh_big(hmin, hmax, 100, FT(0.02)) ≈ 1.0443565558255983
    @test G.β_Δh_both(hmin, hmax, 100, FT(0.001)) ≈ 1.0211975072167816
    @test G.β_Δh_small(hmin, hmax, 100, FT(0.001)) ≈ 1.0218301816802464

    Δh = FT(0.1)
    @test all(uniform_left(hmin, Δh, 3) .≈ [-0.3, -0.2, -0.1, 0.0])
    @test all(uniform_right(hmin, Δh, 3) .≈ [0.0, 0.1, 0.2, 0.3])
end

@testset "Grid init" begin
  for FT in (Float32, Float64)

    a = FT(0.0)
    b = FT(1.0)
    n = 10
    c = Coordinates(a,b,n)
    @test length(c.c.h) == c.c.s
    @test length(c.n.h) == c.n.s
    @test c.n.h[1] < a
    @test c.n.h[end] > b
    @test c.n.h[2] ≈ a
    @test c.n.h[end-1] ≈ b
    @test c.c.h[1] < a
    @test c.c.h[end] > b

    a = FT(0.0)
    b = FT(1.0)
    n = 10
    β = FT(1.1)
    c1 = Coordinates(a,b,n; warpfun=Roberts_both, args=(β,))
    c2 = Coordinates(a,b,n; warpfun=Roberts_both, args=(β,))
    c3 = Coordinates(a,b,n; warpfun=Roberts_both, args=(β,))
    grid = Grid((c1,c2,c3))

    yc = FT(0.5)
    τ = FT(10)
    c1 = Coordinates(a,b,n; warpfun=cluster, args=(yc,τ))
    c2 = Coordinates(a,b,n; warpfun=cluster, args=(yc,τ))
    c3 = Coordinates(a,b,n; warpfun=cluster, args=(yc,τ))
    grid = Grid((c1,c2,c3))

    c1 = Coordinates(a,b,n; warpfun=Roberts_right, args=(β,))
    c2 = Coordinates(a,b,n; warpfun=Roberts_left, args=(β,))
    c3 = Coordinates(a,b,n; warpfun=Roberts_both, args=(β,))
    grid = Grid((c1,c2,c3))
    sprint(show, grid)
    # FT==Float64 && @show grid
    @unpack xn, yn, zn, xc, yc, zc = all_coords(grid)
    @test xn == grid.c[1].n.h
    @test yn == grid.c[2].n.h
    @test zn == grid.c[3].n.h
    @test xc == grid.c[1].c.h
    @test yc == grid.c[2].c.h
    @test zc == grid.c[3].c.h

    # Test failures
    if FT === Float32
        @test_throws ArgumentError("Coordinates are not finite.") Coordinates(a,b,n; warpfun=Roberts_left, args=(FT(1.000000000000001),))
    end

    @test over_points(c1, All(), Vertex1D()) == 1:13
    @test over_points(c1, Interior(), Vertex1D()) == 3:11
    @test over_points(c1, Boundary(), Vertex1D()) == 2:12
    @test over_points(c1, Ghost(), Vertex1D()) == (1:1, 13:13)

    @test over_points(c1, All(), Center1D()) == 1:12
    @test over_points(c1, Interior(), Center1D()) == 2:11
    @test over_points(c1, Ghost(), Center1D()) == (1:1, 12:12)
  end
end

@testset "Stretched parameter matching" begin
  for FT in (Float32, Float64)
    a = FT(0.0)
    b = FT(1.0)
    n = 10
    β = FT(1.1)
    c = Coordinates(a,b,n; warpfun=Roberts_both, args=(β,))
    β = G.β_Δh_small(a, b, n, c.n.Δh[1])
    @test β ≈ 1.1284008533596548
    @test length(c.c.h) == c.c.s
    @test length(c.n.h) == c.n.s
    @test c.n.h[1] < a
    @test c.n.h[end] > b
    @test c.n.h[2] ≈ a
    @test c.n.h[end-1] ≈ b
    @test c.c.h[1] < a
    @test c.c.h[end] > b
  end
end

@testset "Field-like indexing" begin
    FT = Float64
    a = FT(0.0)
    b = FT(1.0)
    n = 10
    β = FT(1.1)
    grid = Cube(a,b,n)
    @test all(grid.c[1].c.h .≈ [grid.cent.x[i,1,1] for i in 1:grid.c[1].c.s])
    @test all(grid.c[2].c.h .≈ [grid.cent.y[1,j,1] for j in 1:grid.c[2].c.s])
    @test all(grid.c[3].c.h .≈ [grid.cent.z[1,1,k] for k in 1:grid.c[3].c.s])

    @test all(grid.c[1].n.h .≈ [grid.corn.x[i,1,1] for i in 1:grid.c[1].n.s])
    @test all(grid.c[2].n.h .≈ [grid.corn.y[1,j,1] for j in 1:grid.c[2].n.s])
    @test all(grid.c[3].n.h .≈ [grid.corn.z[1,1,k] for k in 1:grid.c[3].n.s])
end


text:
module Space

export sample, contains,
    Box, Discrete, TupleSpace, DictSpace, MultiBinary, MultiDiscrete

abstract type AbstractSpace end

include("box.jl")
include("discrete.jl")
include("tuple-space.jl")
include("dict-space.jl")
include("multi-binary.jl")
include("multi-discrete.jl")

Base.in(x, space_obj::AbstractSpace) = contains(space_obj, x)
end #module


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

const GSW_INVALID_VALUE = 9.0e15
const GSW_ERROR_LIMIT = 1.0e10
const INTERP_METHOD_LINEAR = 1
const INTERP_METHOD_PCHIP = 2


text:
my_f(x,y) = 2x+3y
derivative_of_my_f(x,y) = ForwardDiff.derivative(x->my_f(x,y),x)


text:
export statistics, log_bayes_score_component, log_bayes_score_component_BDeu, net_score, net_score_BDeu

function statistics(
    targetind::Int,
    parents::AbstractVector{Int},
    bincounts::AbstractVector{Int},
    discrete_data::AbstractMatrix{Int}
    )

    q = 1
    if !isempty(parents)
        Np = length(parents)
        q  = prod(bincounts[parents])
        stridevec = fill(1, Np)
        for k = 2:Np
            stridevec[k] = stridevec[k-1] * bincounts[parents[k-1]]
        end
        js = (discrete_data[parents,:] - 1)' * stridevec + 1
    else
        js = fill(1, size(discrete_data,2))
    end
    full(sparse(vec(discrete_data[targetind,:]), vec(js), 1, bincounts[targetind], q))
end
function statistics(
    parents :: Vector{Vector{Int}},
    bincounts :: AbstractVector{Int},
    discrete_data :: AbstractMatrix{Int},
    )

    #=
    Computes sufficient statistics from a discrete dataset
    for a Bayesian Net structure

    INPUT:
        parents:
            list of lists of parent indices
            A variable with index i has bincounts[i]
            and row in discrete_data[i,:]
            No acyclicity checking is done
        bincounts:
            list of variable bin counts, or number of
            discrete values the variable can take on, v ∈ {1 : bincounts[i]}
        discrete_data:
            table of discrete values [n×m]
            where n is the number of nodes
            and m is the number of samples

    OUTPUT:
        N :: Vector{Matrix{Int}}
            a sufficient statistics table for each variable
            Variable with index i has statistics table N[i],
            which is r × q where
            r = bincounts[i] is the number of variable instantiations and
            q is the number of parental instantiations of variable i

            The r-values are ordered from 1 → bincounts[i]
            The q-values are ordered in the same ordering as ind2sub() in Julia Base
                Thus the instantiation of the first parent (by order given in parents[i])
                is varied the fastest.

            ex:
                Variable 1 has parents 2 and 3, with r₁ = 2, r₂ = 2, r₃ = 3
                q for variable 1 is q = r₂×r₃ = 6
                N[1] will be a 6×2 matrix where:
                    N[1][1,1] is the number of time v₁ = 1, v₂ = 1, v₃ = 1
                    N[1][2,1] is the number of time v₁ = 1, v₂ = 2, v₃ = 1
                    N[1][3,1] is the number of time v₁ = 1, v₂ = 1, v₃ = 2
                    N[1][4,1] is the number of time v₁ = 1, v₂ = 2, v₃ = 2
                    N[1][5,1] is the number of time v₁ = 1, v₂ = 1, v₃ = 3
                    N[1][6,1] is the number of time v₁ = 1, v₂ = 2, v₃ = 3
                    N[1][6,2] is the number of time v₁ = 2, v₂ = 1, v₃ = 1
                    ...

    This function uses sparse matrix black magic and was
    mercilessly stolen from Ed Schmerling.
    =#

    n, m = size(discrete_data)
    N = Array(Matrix{Int}, n)
    for i = 1 : n
        # p = parents[i]
        # q = 1
        # if !isempty(p)
        #     Np = length(p)
        #     q  = prod(bincounts[parents[i]])
        #     stridevec = fill(1, Np)
        #     for k = 2:Np
        #         stridevec[k] = stridevec[k-1] * bincounts[p[k-1]]
        #     end
        #     js = (discrete_data[p,:] - 1)' * stridevec + 1
        # else
        #     js = fill(1, m)
        # end
        # N[i] = full(sparse(vec(discrete_data[i,:]), vec(js), 1, bincounts[i], q))
        N[i] = statistics(i, parents[i], bincounts, discrete_data)
    end
    N
end
function statistics(net::Network, bincounts::AbstractVector{Int}, discrete_data::AbstractMatrix{Int})

    n = length(bincounts)
	parentList = [convert(Vector{Int}, collect(get_parents(net, i-1)+1)) for i in 1:n]
	statistics(parentList, bincounts, discrete_data)
end
function statistics(adj::BitMatrix, bincounts::AbstractVector{Int}, discrete_data::AbstractMatrix{Int})

    n = length(bincounts)
    @assert(size(adj,1) == size(adj,2) == n)
    parentList = [find(adj[:,i]) for i in 1:n]
    statistics(parentList, bincounts, discrete_data)
end

function log_bayes_score_component{I<:Integer}(i::Int, parents::AbstractVector{I}, r::AbstractVector{Int}, d::AbstractMatrix{Int})

    #=
    Computes the bayesian score component for the given target variable index
        This assumes a unit dirichlet prior (alpha)

    INPUT:
        i       - index of the target variable
        parents - list of indeces of parent variables (should not contain self)
        r       - list of instantiation counts accessed by variable index
                  r[1] gives number of discrete states variable 1 can take on
        d       - matrix of sufficient statistics / counts
                  d[j,k] gives the number of times the target variable took on its kth instantiation
                  given the jth parental instantiation

    OUTPUT:
        the log bayesian score, Float64
    =#

    alpha = fill(1, r[i], isempty(parents) ? 1 : prod(r[parents]))
    (n, m) = size(d)
    if !isempty(parents)
        Np = length(parents)
        stridevec = fill(1, Np)
        for k = 2:Np
            stridevec[k] = stridevec[k-1] * r[parents[k-1]]
        end
        js = (d[parents,:] - 1)' * stridevec + 1
    else
        js = fill(1, m)
    end
    N = sparse(vec(d[i,:]), vec(js), 1, size(alpha)...) # note: duplicates are added together
    return sum(lgamma(alpha + N)) - sum(lgamma(alpha)) + sum(lgamma(sum(alpha,1))) - sum(lgamma(sum(alpha,1) + sum(N,1)))::Float64
end
function log_bayes_score_component_BDeu{I<:Integer}(i::Int, parents::AbstractVector{I}, r::AbstractVector{Int}, d::AbstractMatrix{Int})

	# obtain BDeu alpha
	if isempty(parents)
		α = 1.0 / r[i]
	else
		α = 1.0 / (r[i] * prod(r[parents]))
	end

    alpha = fill(α, r[i], isempty(parents) ? 1 : prod(r[parents]))
    (n, m) = size(d)
    if !isempty(parents)
        Np = length(parents)
        stridevec = fill(1, Np)
        for k = 2:Np
            stridevec[k] = stridevec[k-1] * r[parents[k-1]]
        end
        js = (d[parents,:] - 1)' * stridevec + 1
    else
        js = fill(1, m)
    end
    N = sparse(vec(d[i,:]), vec(js), 1, size(alpha)...) # note: duplicates are added together
    return sum(lgamma(alpha + N)) - sum(lgamma(alpha)) + sum(lgamma(sum(alpha,1))) - sum(lgamma(sum(alpha,1) + sum(N,1)))
end
function log_bayes_score_component{I<:Integer}(
    i       :: Int,
    parents :: AbstractVector{I},
    r       :: AbstractVector{Int},
    d       :: AbstractMatrix{Int},
    alpha   :: AbstractMatrix{Float64}, # r[i]×prod(r[parents])
    )

    (n, m) = size(d)
    if !isempty(parents)
        Np = length(parents)
        stridevec = fill(1, Np)
        for k = 2:Np
            stridevec[k] = stridevec[k-1] * r[parents[k-1]]
        end
        js = (d[parents,:] - 1)' * stridevec + 1
    else
        js = fill(1, m)
    end
    N = sparse(vec(d[i,:]), vec(js), 1, size(alpha)...) # note: duplicates are added together
    return sum(lgamma(alpha + N)) - sum(lgamma(alpha)) + sum(lgamma(sum(alpha,1))) - sum(lgamma(sum(alpha,1) + sum(N,1)))
end
function log_bayes_score_component{I<:Integer}(i::Int, parents::AbstractVector{I}, r::SharedArray{Int, 1}, d::SharedArray{Int, 2}, cache::Dict{Vector{Int}, Float64})
	if haskey(cache, parents)
        return cache[parents]
    end
    return (cache[parents] = log_bayes_score_component(i, parents, r, d))
end
function log_bayes_score_component{I<:Integer}(i::Int, parents::AbstractVector{I}, r::AbstractVector{Int}, d::Matrix{Int}, cache::Dict{Vector{Int}, Float64})
	if haskey(cache, parents)
        return cache[parents]
    end
    return (cache[parents] = log_bayes_score_component(i, parents, r, d))
end
function log_bayes_score_component{I<:Integer}(
    i       :: Int,
    parents :: AbstractVector{I},
    r       :: AbstractVector{Int},
    d       :: AbstractMatrix{Int},
    alpha   :: AbstractMatrix{Float64}, # r[i]×prod(r[parents])
    cache   :: Dict{Vector{Int}, Float64}
    )

    if haskey(cache, parents)
        return cache[parents]
    end
    return (cache[parents] = log_bayes_score_component(i, parents, r, d, alpha))
end

function net_score( net::Network, r::SharedArray{Int, 1}, d::SharedArray{Int, 2} )
	# compute the log bayes score for the given network on the dataset
	s = @parallel (+) for n in get_all_nodes(net)
		log_bayes_score_component(n + 1, get_parents(net, n)+1, r, d)
	end
	s
end
function net_score( net::Network, r::AbstractVector{Int}, d::Matrix{Int} )
	# compute the log bayes score for the given network on the dataset
	s = @parallel (+) for n in get_all_nodes(net)
		log_bayes_score_component(n + 1, get_parents(net, n)+1, r, d)
	end
	s
end
function net_score_BDeu( net::Network, r::SharedArray{Int, 1}, d::SharedArray{Int, 2} )
	# compute the log bayes score for the given network on the dataset using the BDeu prior
	s = @parallel (+) for n in get_all_nodes(net)
		log_bayes_score_component_BDeu(n + 1, get_parents(net, n)+1, r, d)
	end
	s
end
function net_score_BDeu( net::Network, r::AbstractVector{Int}, d::Matrix{Int} )
	# compute the log bayes score for the given network on the dataset using the BDeu prior
	s = @parallel (+) for n in get_all_nodes(net)
		log_bayes_score_component_BDeu(n + 1, get_parents(net, n)+1, r, d)
	end
	s
end


text:
function test_FixedPolicy()
    pol = FixedPolicy(2)
    a = choose(rand(2), pol)
    update!(pol, rand(5), a, 1.0)
    @test a == 2
end

test_FixedPolicy()


text:
module LinearTimeVaryingModelsBase

using Parameters, Statistics, LinearAlgebra
# Interface exports
export Trajectory, AbstractModel, AbstractCost, ModelAndCost,f,
dc,calculate_cost,calculate_final_cost, predict, simulate, df,costfun, LTVStateSpaceModel,
SimpleLTVModel, covariance, whiten!

export rms, sse, nrmse, modelfit, aic

rms(x::AbstractVector) = sqrt(mean(abs2,x))
sse(x::AbstractVector) = x⋅x

rms(x::AbstractMatrix) = sqrt.(mean(abs2.(x),dims=2))[:]
sse(x::AbstractMatrix) = sum(abs2,x,dims=2)[:]
modelfit(y,yh) = 100 * (1 .-rms(y.-yh)./rms(y.-mean(y)))
aic(x::AbstractVector,d) = log(sse(x)) .+ 2d/size(x,2)
const nrmse = modelfit

# Trajectory =========================================
@with_kw mutable struct Trajectory
    x::Matrix{Float64}
    u::Matrix{Float64}
    y::Matrix{Float64}
    xu::Matrix{Float64}
    yu::Matrix{Float64}
    nx::Int
    nu::Int
    ny::Int
end
function Trajectory(x,u)
    T = size(x,2)
    @assert T ∈ (size(u,2),size(u,2)+1) ||  "The second dimension (time) of x must be equal to or one greater than the second dimension of u "
    x,u,y = x[:,1:T-1],u[:,1:T-1],x[:,2:T]
    Trajectory(x,u,y,[x;u],[y;u],size(x,1), size(u,1), size(y,1))
end
function Trajectory(x,u,y)
    @assert size(x,2) == size(u,2) == size(y,2) "The second dimension of x,u and y (time) must be the same"
    Trajectory(x,u,y,[x;u],[y;u],size(x,1), size(u,1), size(y,1))
end
Base.length(t::Trajectory) = size(t.u,2)

function Base.iterate(t::Trajectory, state=1)
    state == length(t) && return nothing
    (t.x[:,state], t.u[:,state], t.y[:,state]), state+1
end

function whiten!(t::Trajectory)
    whiten!(t.x, 2)
    whiten!(t.u, 2)
    whiten!(t.y, 2)
end

function whiten!(m::Matrix, d=2)
    m .-= mean(m,d)
    m ./= std(m,d)
end

# Model interface ====================================
"""
Model interface, implement the following functions\n
see also `AbstractCost`, `ModelAndCost`
```
fit_model(::Type{AbstractModel}, batch::Batch)::AbstractModel

predict(model::AbstractModel, x, u)

function df(model::AbstractModel, x, u, I::UnitRange)
    return fx,fu,fxx,fxu,fuu
end
```
"""
abstract type AbstractModel end

abstract type LTVModel <: AbstractModel end

abstract type LTVStateSpaceModel <: LTVModel end

mutable struct SimpleLTVModel{T} <: LTVStateSpaceModel
    At::Array{T,3}
    Bt::Array{T,3}
    extended::Bool
    function SimpleLTVModel(At::Array{T,3},Bt::Array{T,3},extend::Bool) where T
        if extend
            At = cat(At,At[:,:,end], dims=3)
            Bt = cat(Bt,Bt[:,:,end], dims=3)
        end
        return new{T}(At,Bt,extend)
    end
end

function df(model::SimpleLTVModel, x, u)
    return model.At,model.Bt,[],[],[]
end

function covariance(model::SimpleLTVModel, x, u)
    n = size(model.At,1)
    zeros(n,n)
end

# SimpleLTVModel(At,Bt,extend::Bool) = SimpleLTVModel{eltype(At)}(At,Bt,extend)

"""
    model = fit_model(::Type{AbstractModel}, x,u)::AbstractModel

Fits a model to data
"""
function fit_model(::Type{AbstractModel}, x,u)::AbstractModel end

"""
    fit_model!(model::AbstractModel, x,u)::AbstractModel

Refits a model to new data
"""
function fit_model!(model::AbstractModel, x,u)::AbstractModel end

"""
    xnew = predict(model::AbstractModel, t::Trajectory [, i])
    xnew = predict(model::AbstractModel, x, u [, i])

Predict the next state given the current state and action
"""
function predict(model::AbstractModel, x, u, i) end

predict(model::AbstractModel, t::Trajectory, args...) = predict(model, t.x, t.u, args...)

"""
    xnew = simulate(model::AbstractModel, t::Trajectory)
    xnew = simulate(model::AbstractModel, x0, u)

Simulate system forward in time given the initial state and actions
"""
function simulate(model::AbstractModel, x0, u) end

simulate(model::AbstractModel, t::Trajectory) = simulate(model, t.x[:,1], t.u)

"""
    fx,fu,fxx,fxu,fuu = df(model::AbstractModel, x, u)

Get the linearized dynamics at `x`,`u`
"""
function df(model::AbstractModel, x, u) end

function covariance(model::AbstractModel, x, u)
    cov(x[:,2:end]-predict(model, x, u)[:,1:end-1], dims=2)
end
# Model interface ====================================

# Cost interface ====================================
"""
Cost interface, implement the following functions\n
see also `AbstractModel`, `ModelAndCost`
```
function calculate_cost(::Type{AbstractCost}, x::AbstractVector, u)::Number

function calculate_cost(::Type{AbstractCost}, x::AbstractMatrix, u)::AbstractVector

function calculate_final_cost(::Type{AbstractCost}, x::AbstractVector)::Number

function dc(::Type{AbstractCost}, x, u)
    return cx,cu,cxx,cuu,cxu
end
```
"""
abstract type AbstractCost end

function calculate_cost(c::AbstractCost, x::AbstractVector, u)::Number end

function calculate_cost(c::AbstractCost, x::AbstractMatrix, u)::AbstractVector end

function calculate_final_cost(c::AbstractCost, x::AbstractVector)::Number end

function dc(c::AbstractCost, x, u) end
# Cost interface ====================================

"""
1. Define types that implement the interfaces `AbstractModel` and `AbstractCost`.
2. Create object modelcost = ModelAndCost(model, cost)
3. Run macro @define_modelcost_functions(modelcost). This macro defines the following functions
```
f(x, u, i)  = f(modelcost, x, u, i)
fT(x)       = fT(modelcost, x)
df(x, u, I) = df(modelcost, x, u, I)
```
see also `AbstractModel`, `AbstractCost`
"""
mutable struct ModelAndCost
    model::AbstractModel
    cost::AbstractCost
end

function f(modelcost::ModelAndCost, x, u, i)
    predict(modelcost.model, x, u, i)
end

function costfun(modelcost::ModelAndCost, x, u)
    calculate_cost(modelcost.cost, x, u)
end

"""
    fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu = df(modelcost::ModelAndCost, x, u)

Get the linearized dynamics and cost at `x`,`u`
"""
function df(modelcost::ModelAndCost, x, u)
    fx,fu,fxx,fxu,fuu = df(modelcost.model, x, u)
    cx,cu,cxx,cuu,cxu = dc(modelcost.cost, x, u)
    return fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu
end

end # module