text:
# =========================================================================
#    COLUMN GENERATION MECHANISMS
# =========================================================================

function master(n::node;silent::Bool)
    mp = buildMaster(n;silent = silent)
    optimize!(mp)

    return mp
end

function buildMaster(n::node;silent::Bool)
    mp = Model(get_default_optimizer())
    if silent || isnothing(silent)
        set_silent(mp)
    end

    R = Dict(1:length(n.columns) .=> n.columns)

    @variable(mp, θ[keys(R), keys(b().K), b().T] >= 0) #column decision
    @variable(mp, I[keys(b().V), vcat(first(b().T)-1, b().T)]) #inventory level
    @variable(mp, 0 <= slack[i = keys(b().V), t = b().T] <= n.stab.slLim[i,t]) #slack var
    @variable(mp, 0 <= surp[i = keys(b().V), t = b().T] <= n.stab.suLim[i,t]) #surp var

    @objective(
        mp, Min,
        sum(
            θ[r,k,t] * (
                sum(
                    b().dist[i,j] * (
                        b().K[k].vx * R[r].x[i,j,k,t] +
                        b().K[k].vl * R[r].l[i,j,k,t]
                    )
                    for i in b().K[k].cover, j in b().K[k].cover
                ) +
                sum(
                    b().K[k].fd * R[r].u[i,k,t]
                    for i in b().K[k].cover
                ) +
                sum(
                    b().K[k].fp * R[r].z[i,k,t]
                    for i in b().K[k].cover
                )
            )
            for r in keys(R), k in keys(b().K), t in b().T
        ) + #pickup delivery costs
        sum(
            b().V[i].h * I[i,t]
            for i in keys(b().V), t in b().T
        ) + #inventory costs
        sum(
            n.stab.slCoeff * slack[i,t]
            for i in keys(b().V), t in b().T
        ) - #stabilizer
        sum(
            n.stab.suCoeff * surp[i,t]
            for i in keys(b().V), t in b().T
        ) #stabilizer
    )

    @constraint(mp, λ[i = keys(b().V), t = b().T],
        I[i,t-1] + sum(R[r].q[i,k,t] * θ[r,k,t] for r in keys(R), k in keys(b().K)) +
        slack[i,t] - surp[i,t] == b().d[i,t] + I[i,t] #inventory balance
    )

    @constraint(mp, δ[i = keys(b().V), k = keys(b().K), t = b().T],
        sum(R[r].z[i,k,t] * θ[r,k,t] for r in keys(R)) <= b().K[k].freq
    )

    @constraint(mp, [i = keys(b().V), t = b().T],
        b().V[i].MIN <= I[i,t] <= b().V[i].MAX #inventory capacity interval
    )

    @constraint(mp, [i = keys(b().V)],
        I[i,first(b().T)-1] == b().V[i].START #starting inventory level
    )

    # ================================
    #    BOUND GENERATOR
    # ================================

    return mp
end

function getDuals(mp::Model)
    λ = dual.(mp.obj_dict[:λ])
    δ = dual.(mp.obj_dict[:δ])

    return dval(λ,δ)
end

function sub(n::node,duals::dval;silent::Bool)
    sp = buildSub(n,duals;silent=silent)
    optimize!(sp)

    return sp
end

function buildSub(n::node,duals::dval;silent::Bool)
    sp = Model(get_default_optimizer())
    if silent || isnothing(silent)
        set_silent(sp)
    end

    @variable(sp, q[i = keys(b().V), k = keys(b().K), t = b().T])
    @variable(sp, u[i = keys(b().V), k = keys(b().K), t = b().T] >= 0)
    @variable(sp, v[i = keys(b().V), k = keys(b().K), t = b().T] >= 0)
    @variable(sp,
        l[i = collect(keys(b().V)), j = collect(keys(b().V)),
        k = collect(keys(b().K)), t = b().T] >= 0
    ) #quantity variables

    @variable(sp, p[i = keys(b().V), k = keys(b().K), t = b().T] >= 0, Int)
    @variable(sp, y[i = keys(b().V), k = keys(b().K), t = b().T] >= 0, Int)
    @variable(sp, z[i = keys(b().V), k = keys(b().K), t = b().T] >= 0, Int)
    @variable(sp,
        x[i = collect(keys(b().V)), j = collect(keys(b().V)),
        k = collect(keys(b().K)), t = b().T] >= 0, Int
    ) #0-1 variables

    @objective(
        sp, Min,
        sum(
            sum(
                b().dist[i,j] * (
                    b().K[k].vx * x[i,j,k,t] +
                    b().K[k].vl * l[i,j,k,t]
                )
                for i in b().K[k].cover, j in b().K[k].cover
            ) +
            sum(
                b().K[k].fd * u[i,k,t]
                for i in b().K[k].cover
            ) +
            sum(
                b().K[k].fp * z[i,k,t]
                for i in b().K[k].cover
            )
            for k in keys(b().K), t in b().T
        ) - #column cost
        sum(
            sum(
                q[i,k,t] * duals.λ[i,t]
                for i in b().K[k].cover
            )
            for k in keys(b().K), t in b().T
        ) - #dual part 1
        sum(
            sum(
                z[i,k,t] * duals.δ[i,k,t]
                for i in b().K[k].cover
            )
            for k in keys(b().K), t in b().T
        ) #dual part 2
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        q[i,k,t] == u[i,k,t] - v[i,k,t] #q breakdown
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        p[i,k,t] == y[i,k,t] + z[i,k,t] #p breakdown
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        u[i,k,t] <= b().K[k].Q * y[i,k,t] # u - y correlation
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        v[i,k,t] <= b().K[k].Q * z[i,k,t] # v - z correlation
    )

    @constraint(sp,
        [k = keys(b().K), i = b().K[k].cover, j = b().K[k].cover, t = b().T],
        l[i,j,k,t] <= b().K[k].Q * x[i,j,k,t] # l - x correlation
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        sum(x[j,i,k,t] for j in b().K[k].cover) == p[i,k,t] #traverse in to i
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        sum(x[i,j,k,t] for j in b().K[k].cover) == p[i,k,t] #traverse out from i
    )

    @constraint(sp, [k = keys(b().K), i = b().K[k].cover, t = b().T],
        sum(l[j,i,k,t] for j in b().K[k].cover) -
        sum(l[i,j,k,t] for j in b().K[k].cover) == q[i,k,t] #vehicle load balance
    )

    @constraint(sp, [k = keys(b().K), t = b().T],
        sum(z[i,k,t] for i in b().K[k].cover) <= b().K[k].freq #maximum manifest
    )

    @constraint(sp, [k = keys(b().K), t = b().T],
        sum(q[i,k,t] for i in b().K[k].cover) == 0 #all pickup delivered
    )

    # ================================
    #    BOUND GENERATOR
    # ================================

    return sp
end

function getCols(sp::Model)
    q = value.(sp.obj_dict[:q])
    u = value.(sp.obj_dict[:u])
    v = value.(sp.obj_dict[:v])
    l = value.(sp.obj_dict[:l])
    p = value.(sp.obj_dict[:p])
    y = value.(sp.obj_dict[:y])
    z = value.(sp.obj_dict[:z])
    x = value.(sp.obj_dict[:x])

    return col(q,u,v,l,p,y,z,x)
end

function updateStab!(stab::stabilizer,param::Float64)
    for i in first(stab.slLim.axes),t in last(stab.slLim.axes)
        stab.slLim[i,t] = param * stab.slLim[i,t]
        if stab.slLim[i,t] < 1
            stab.slLim[i,t] = 0
        end
    end

    for i in first(stab.suLim.axes),t in last(stab.suLim.axes)
        stab.suLim[i,t] = param * stab.suLim[i,t]
        if stab.suLim[i,t] < 1
            stab.suLim[i,t] = 0
        end
    end

    return stab
end

function checkStab(mp::Model)
    s = sum(value.(mp.obj_dict[:slack])) + sum(value.(mp.obj_dict[:surp]))

    return s
end

function colGen(n::node;silent::Bool,maxCG::Float64,track::Bool)
    terminate = false
    iter = 0

    while !terminate
        if iter < maxCG
            mp = master(n;silent=silent)
            if has_values(mp) && has_duals(mp)
                if track #print master problem obj
                    println("obj: $(objective_value(mp))")
                end

                duals = getDuals(mp)
                sp = sub(n,duals;silent=silent)
                if track #print subproblem price
                    println("price: $(objective_value(sp))")
                end

                if isapprox(objective_value(sp), 0, atol = 1e-8) || objective_value(sp) > 0
                    if isapprox(checkStab(mp), 0, atol = 1e-8)
                        terminate = true #action
                        push!(n.status,"EVALUATED") #report
                        if track
                            println("EVALUATED")
                        end
                    else
                        updateStab!(n.stab,0.5) #action
                        push!(n.status,"STABILIZED") #report
                        if track
                            println("STABILIZED")
                        end
                    end
                else
                    push!(n.columns,getCols(sp)) #action
                    push!(n.status,"ADD_COLUMN") #report
                    if track
                        println("ADD_COLUMN")
                    end
                end

                iter += 1 #iteration update
            else
                terminate = true #action
                push!(n.status,"NO_SOLUTION")
                if track
                    println("NO_SOLUTION")
                end
            end
        else
            terminate = true #action
            push!(n.status,"EVALUATED-TIME OUT") #report
            if track
                println("EVALUATED-TIME OUT")
            end
        end
    end

    if n.status[end] != "NO_SOLUTION"
        println("NODE $(n.self) FINISHED.")
    else
        println("NODE $(n.self) FAILED.")
    end

    return n
end


text:
# Julia wrapper for header: /usr/include/libavfilter/buffersink.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_buffersink_get_frame_flags,
    av_buffersink_params_alloc,
    av_abuffersink_params_alloc,
    av_buffersink_set_frame_size,
    av_buffersink_get_type,
    av_buffersink_get_time_base,
    av_buffersink_get_format,
    av_buffersink_get_frame_rate,
    av_buffersink_get_w,
    av_buffersink_get_h,
    av_buffersink_get_sample_aspect_ratio,
    av_buffersink_get_channels,
    av_buffersink_get_channel_layout,
    av_buffersink_get_sample_rate,
    av_buffersink_get_hw_frames_ctx,
    av_buffersink_get_frame,
    av_buffersink_get_samples

function av_buffersink_get_frame_flags(ctx, frame, flags::Integer)
    ccall((:av_buffersink_get_frame_flags, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVFrame}, Cint), ctx, frame, flags)
end

function av_buffersink_params_alloc()
    ccall((:av_buffersink_params_alloc, libavfilter), Ptr{AVBufferSinkParams}, ())
end

function av_abuffersink_params_alloc()
    ccall((:av_abuffersink_params_alloc, libavfilter), Ptr{AVABufferSinkParams}, ())
end

function av_buffersink_set_frame_size(ctx, frame_size::Integer)
    ccall((:av_buffersink_set_frame_size, libavfilter), Cvoid, (Ptr{AVFilterContext}, UInt32), ctx, frame_size)
end

function av_buffersink_get_type(ctx)
    ccall((:av_buffersink_get_type, libavfilter), AVMediaType, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_time_base(ctx)
    ccall((:av_buffersink_get_time_base, libavfilter), AVRational, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_format(ctx)
    ccall((:av_buffersink_get_format, libavfilter), Cint, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_frame_rate(ctx)
    ccall((:av_buffersink_get_frame_rate, libavfilter), AVRational, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_w(ctx)
    ccall((:av_buffersink_get_w, libavfilter), Cint, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_h(ctx)
    ccall((:av_buffersink_get_h, libavfilter), Cint, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_sample_aspect_ratio(ctx)
    ccall((:av_buffersink_get_sample_aspect_ratio, libavfilter), AVRational, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_channels(ctx)
    ccall((:av_buffersink_get_channels, libavfilter), Cint, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_channel_layout(ctx)
    ccall((:av_buffersink_get_channel_layout, libavfilter), UInt64, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_sample_rate(ctx)
    ccall((:av_buffersink_get_sample_rate, libavfilter), Cint, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_hw_frames_ctx(ctx)
    ccall((:av_buffersink_get_hw_frames_ctx, libavfilter), Ptr{AVBufferRef}, (Ptr{AVFilterContext},), ctx)
end

function av_buffersink_get_frame(ctx, frame)
    ccall((:av_buffersink_get_frame, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVFrame}), ctx, frame)
end

function av_buffersink_get_samples(ctx, frame, nb_samples::Integer)
    ccall((:av_buffersink_get_samples, libavfilter), Cint, (Ptr{AVFilterContext}, Ptr{AVFrame}, Cint), ctx, frame, nb_samples)
end

ext: jl


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const ApiDef_Visibility = (;[
    Symbol("DEFAULT_VISIBILITY") => Int32(0),
    Symbol("VISIBLE") => Int32(1),
    Symbol("SKIP") => Int32(2),
    Symbol("HIDDEN") => Int32(3),
]...)

mutable struct ApiDef_Endpoint <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ApiDef_Endpoint(; kwargs...)
        obj = new(meta(ApiDef_Endpoint), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ApiDef_Endpoint
const __meta_ApiDef_Endpoint = Ref{ProtoMeta}()
function meta(::Type{ApiDef_Endpoint})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ApiDef_Endpoint)
            __meta_ApiDef_Endpoint[] = target = ProtoMeta(ApiDef_Endpoint)
            fnum = Int[1,3,4]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :deprecated => Bool, :deprecation_version => Int32]
            meta(target, ApiDef_Endpoint, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ApiDef_Endpoint[]
    end
end
function Base.getproperty(obj::ApiDef_Endpoint, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :deprecated
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :deprecation_version
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct ApiDef_Arg <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ApiDef_Arg(; kwargs...)
        obj = new(meta(ApiDef_Arg), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ApiDef_Arg
const __meta_ApiDef_Arg = Ref{ProtoMeta}()
function meta(::Type{ApiDef_Arg})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ApiDef_Arg)
            __meta_ApiDef_Arg[] = target = ProtoMeta(ApiDef_Arg)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :rename_to => AbstractString, :description => AbstractString]
            meta(target, ApiDef_Arg, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ApiDef_Arg[]
    end
end
function Base.getproperty(obj::ApiDef_Arg, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :rename_to
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct ApiDef_Attr <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ApiDef_Attr(; kwargs...)
        obj = new(meta(ApiDef_Attr), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ApiDef_Attr
const __meta_ApiDef_Attr = Ref{ProtoMeta}()
function meta(::Type{ApiDef_Attr})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ApiDef_Attr)
            __meta_ApiDef_Attr[] = target = ProtoMeta(ApiDef_Attr)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :rename_to => AbstractString, :default_value => AttrValue, :description => AbstractString]
            meta(target, ApiDef_Attr, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ApiDef_Attr[]
    end
end
function Base.getproperty(obj::ApiDef_Attr, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :rename_to
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :default_value
        return (obj.__protobuf_jl_internal_values[name])::AttrValue
    elseif name === :description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct ApiDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ApiDef(; kwargs...)
        obj = new(meta(ApiDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ApiDef
const __meta_ApiDef = Ref{ProtoMeta}()
function meta(::Type{ApiDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ApiDef)
            __meta_ApiDef[] = target = ProtoMeta(ApiDef)
            fnum = Int[1,12,13,2,3,4,5,11,6,7,8,9,10]
            allflds = Pair{Symbol,Union{Type,String}}[:graph_op_name => AbstractString, :deprecation_message => AbstractString, :deprecation_version => Int32, :visibility => Int32, :endpoint => Base.Vector{ApiDef_Endpoint}, :in_arg => Base.Vector{ApiDef_Arg}, :out_arg => Base.Vector{ApiDef_Arg}, :arg_order => Base.Vector{AbstractString}, :attr => Base.Vector{ApiDef_Attr}, :summary => AbstractString, :description => AbstractString, :description_prefix => AbstractString, :description_suffix => AbstractString]
            meta(target, ApiDef, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ApiDef[]
    end
end
function Base.getproperty(obj::ApiDef, name::Symbol)
    if name === :graph_op_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :deprecation_message
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :deprecation_version
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :visibility
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :endpoint
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ApiDef_Endpoint}
    elseif name === :in_arg
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ApiDef_Arg}
    elseif name === :out_arg
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ApiDef_Arg}
    elseif name === :arg_order
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :attr
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ApiDef_Attr}
    elseif name === :summary
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :description_prefix
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :description_suffix
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct ApiDefs <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ApiDefs(; kwargs...)
        obj = new(meta(ApiDefs), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ApiDefs
const __meta_ApiDefs = Ref{ProtoMeta}()
function meta(::Type{ApiDefs})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ApiDefs)
            __meta_ApiDefs[] = target = ProtoMeta(ApiDefs)
            allflds = Pair{Symbol,Union{Type,String}}[:op => Base.Vector{ApiDef}]
            meta(target, ApiDefs, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ApiDefs[]
    end
end
function Base.getproperty(obj::ApiDefs, name::Symbol)
    if name === :op
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ApiDef}
    else
        getfield(obj, name)
    end
end

export ApiDef_Visibility, ApiDef_Endpoint, ApiDef_Arg, ApiDef_Attr, ApiDef, ApiDefs


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

module test_withoutMacros

import ModiaMath

# Desired:
#   using Test
#
# In order that Test needs not to be defined in the user environment, it is included via ModiaMath:
using ModiaMath.Test

include(joinpath(ModiaMath.path, "examples", "withoutMacros", "Simulate_Pendulum.jl"))
import .Simulate_Pendulum
pendulumWithoutMacro = ModiaMath.SimulationModel(Simulate_Pendulum.PendulumWithoutMacro(L=0.5, m=1.0, d=0.1, phi0_deg=60.0))

include(joinpath(ModiaMath.path, "examples", "withoutMacros", "Simulate_FreeBodyRotation.jl"))
import .Simulate_FreeBodyRotation
freeBodyRotationWithoutMacro = ModiaMath.SimulationModel(Simulate_FreeBodyRotation.FreeBodyRotationWithoutMacro(), stopTime=5.0, tolerance=1e-8)

@testset "ModiaMath: withoutMacro/Simulate_Pendulum.jl" begin
    result = ModiaMath.simulate!(pendulumWithoutMacro, stopTime=10.0, tolerance=1e-8, interval=0.1, log=true)
    phi = result.series["phi"]
    w   = result.series["w"]

    @test isapprox(phi[end], 0.120289; atol=1e-3 )
    @test isapprox(w[end], 0.273965; atol=1e-2 )
end

@testset "ModiaMath: withoutMacro/Simulate_FreeBodyRotation.jl" begin
    result = ModiaMath.simulate!(freeBodyRotationWithoutMacro, log=true)
    ModiaMath.plot(result, [:q, :w, :derw, :tau])
end

end


text:
# ------------------------------------------------------------------------------------------
# Data
# ------------------------------------------------------------------------------------------
struct Data{N<:Name}
    seed::Int
    shuffle::Bool
    poslabels::Int
    batchsize::Int

    function Data(
        N;
        seed = 1234,
        shuffle = true,
        poslabels = 0,
        batchsize = 0,
    )
        return new{N}(seed, shuffle, poslabels, batchsize)
    end
end

function Base.show(io::IO, d::Data{N}) where {N}
    pars = [:seed, :shuffle,:poslabels]
    d.batchsize != 0 && push!(pars, :batchsize)
    print(io, nameof(N), NamedTuple{(pars...,)}(getfield.(Ref(d), pars)))
    return
end

function todict(d::Data{N}) where {N}
    return Dict(
        :dataset => nameof(N),
        :seed => d.seed,
        :shuffle => d.shuffle,
        :poslabels => d.poslabels,
        :batchsize => d.batchsize,
    )
end

function loaddata(d::Data{N}) where {N}
    data = DatasetProvider.Dataset(
        N;
        seed = d.seed,
        shuffle = d.shuffle,
        asmatrix = true,
        poslabels = d.poslabels,
        binarize = true,
    )

    return  reshape_dataset.(DatasetProvider.load(TrainValidTest(), data))
end

# Reshape to Flux style
reshape_dataset(data::Tuple) = (reshape_samples(data[1]), reshape_labels(data[2]))

reshape_labels(y) = y
reshape_labels(y::AbstractVector) = Array(reshape(y, 1, :))

reshape_samples(x; T = Float32) = T.(x)
reshape_samples(x::AbstractMatrix; T = Float32) = Matrix{T}(x')
reshape_samples(x::AbstractArray; T = Float32) = Array{T}(reshape(x, :, size(x)[end]))

# ------------------------------------------------------------------------------------------
# Train
# ------------------------------------------------------------------------------------------
Base.@kwdef struct Train
    seed::Int = 1
    iters::Int = 1000
    saveat::Int = 200
    optimiser = ADAM
    step::Float64 = 0.01
end

function Base.show(io::IO, m::Train)
    @unpack seed, iters, saveat, optimiser, step = m
    print(io, "Train", (; seed, iters, saveat, optimiser, step))
    return
end

function todict(t::Train)
    return Dict(
        :seed => t.seed,
        :iters => t.iters,
        :saveat => t.saveat,
        :optimiser => nameof(t.optimiser),
        :step => t.step,
    )
end

function build_network(::Data{N}, t::Train; zero::Bool = false) where {N}
    Random.seed!(t.seed)
    model = Dense(prod(DatasetProvider.nattributes(N)), 1)
    pars = params(model)
    delete!(pars, model.b)
    if zero
        for p in pars
            p .= 0
        end
    end
    return model, pars
end

# ------------------------------------------------------------------------------------------
# Model
# ------------------------------------------------------------------------------------------
struct Model{T<:AbstractThreshold}
    pars
    surr
    λ

    Model(T, pars, surr, λ) = new{T}(pars, surr, λ)
end

function Model(T::Type{<:AbstractThreshold}; λ = 0.001, surrogate = hinge, kwargs...)
    d = Dict(:pars => parameters(T; surrogate, kwargs...), :λ => λ, :surr => surrogate)

    return map(dict_list(d)) do dict
        @unpack pars, surr, λ = dict
        return Model(T, pars, surr, λ)
    end
end

function Base.show(io::IO, m::Model{T}) where {T}
    return print(io, nameof(T), (; m.pars..., surrogate = nameof(m.surr), λ = m.λ))
end

function todict(m::Model{T}) where {T}
    return Dict(
        :model => nameof(T),
        :pars => ntuple2dict(m.pars),
        :surr => nameof(m.surr),
        :λ => m.λ,
    )
end

# Thresholds
build(T::Type{<:Union{PatMat, PatMatNP}}, τ, surr, β) = T(τ, x -> surr(x, β))
validkeys(::Type{<:Union{PatMat, PatMatNP}}) = (:τ, :surrogate, :β)
savekeys(::Union{Type{PatMat}, Type{PatMatNP}}) = (:τ, :surrogate, :β, :λ)

build(T::Type{<:Union{Grill, GrillNP, τFPL, TopMean}}, τ) = T(τ)
validkeys(::Type{<:Union{Grill, GrillNP, τFPL, TopMean}}) = (:τ, )
savekeys(::Type{<:Union{Grill, GrillNP, τFPL, TopMean}}) = (:τ, :surrogate, :λ)

build(T::Type{TopPush}) = T()
validkeys(::Type{TopPush}) = tuple()
savekeys(::Type{TopPush}) = (:surrogate, :λ)

build(T::Type{TopPushK}, K) = T(K)
validkeys(::Type{TopPushK}) = (:K, )
savekeys(::Type{TopPushK}) = (:K, :surrogate, :λ)

function parameters(T::Type{<:AbstractThreshold}; kwargs...)
    keys = validkeys(T)
    vals = getindex.(Ref(kwargs), keys)

    return map(dict_list(Dict(tuple.(keys, vals)))) do dict
        return NamedTuple{keys}(getindex.(Ref(dict), keys))
    end
end

# Loss functions
sqsum(x) = sum(abs2, x)

function build_loss(m::Model{T}) where {T}
    @unpack surr, λ = m
    thres = build(T, m.pars...)

    loss(x, y, model, pars) = loss(y, model(x), pars)

    function loss(y, s, pars)
        t = threshold(thres, y, s)
        return fnr(y, s, t, surr) + eltype(s)(λ) * sum(sqsum, pars)
    end
    return loss
end

function build_loss(m::Model{T}) where {T<:Type{<:Union{Grill, GrillNP}}}
    @unpack surr, λ = m
    thres = build(T, m.pars...)

    loss(x, y, model, pars) = loss(y, model(x), pars)

    function loss(y, s, pars)
        t = threshold(thres, y, s)
        return fnr(y, s, t, surr) + fpr(y, s, t, surr) + eltype(s)(λ) * sum(sqsum, pars)
    end
    return loss
end


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

mutable struct IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec <: SwaggerModel
    maxReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: maxReplicas
    metrics::Any # spec type: Union{ Nothing, Vector{IoK8sApiAutoscalingV2beta1MetricSpec} } # spec name: metrics
    minReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: minReplicas
    scaleTargetRef::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1CrossVersionObjectReference } # spec name: scaleTargetRef

    function IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec(;maxReplicas=nothing, metrics=nothing, minReplicas=nothing, scaleTargetRef=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec, Symbol("maxReplicas"), maxReplicas)
        setfield!(o, Symbol("maxReplicas"), maxReplicas)
        validate_property(IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec, Symbol("metrics"), metrics)
        setfield!(o, Symbol("metrics"), metrics)
        validate_property(IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec, Symbol("minReplicas"), minReplicas)
        setfield!(o, Symbol("minReplicas"), minReplicas)
        validate_property(IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec, Symbol("scaleTargetRef"), scaleTargetRef)
        setfield!(o, Symbol("scaleTargetRef"), scaleTargetRef)
        o
    end
end # type IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec

const _property_map_IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec = Dict{Symbol,Symbol}(Symbol("maxReplicas")=>Symbol("maxReplicas"), Symbol("metrics")=>Symbol("metrics"), Symbol("minReplicas")=>Symbol("minReplicas"), Symbol("scaleTargetRef")=>Symbol("scaleTargetRef"))
const _property_types_IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec = Dict{Symbol,String}(Symbol("maxReplicas")=>"Int32", Symbol("metrics")=>"Vector{IoK8sApiAutoscalingV2beta1MetricSpec}", Symbol("minReplicas")=>"Int32", Symbol("scaleTargetRef")=>"IoK8sApiAutoscalingV2beta1CrossVersionObjectReference")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec)
    (getproperty(o, Symbol("maxReplicas")) === nothing) && (return false)
    (getproperty(o, Symbol("scaleTargetRef")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1HorizontalPodAutoscalerSpec }, name::Symbol, val)
end


text:
function reverse_dependency_tree(server_module)
    output_fns = module_functions(server_module)
    # Map from input variables to the outputs they affect
    dependencies = Dict{String, Array{String, 1}}()
    for fn in output_fns
        inputs = argument_names(get_handler(fn, server_module))
        for input_var in string.(inputs)
            if haskey(dependencies, input_var)
                push!(dependencies[input_var], string(fn))
            else
                dependencies[input_var] = [string(fn)]
            end
        end
    end

    filter(k -> first(k) != "session", dependencies)
end

function dependency_tree(server_module)
    output_fns = module_functions(server_module)
    dependencies = Dict{String, Array{String, 1}}()
    for fn in output_fns
        args = argument_names(get_handler(fn, server_module))
        args = filter(x -> x != :session, args)
        dependencies[string(fn)] = string.(args)
    end
    dependencies
end

function is_user_defined(symbol)
    string(symbol)[1] != '#' && symbol != :eval && symbol != :include
end

function module_functions(server_module)
    output_fns = names(server_module, all = true)
    output_fns = filter(is_user_defined, output_fns)[2:end]
    filter(f -> length(methods(getfield(server_module, f))) > 0, output_fns) # Drop functions with no methods. This happens when using Revise if functions are deleted
end

# this works for closures
function argument_names(fn::Function)
    m = methods(fn)
    length(m) > 1 && error("You cannot define more than one method for each output (for output `$id`)")
    m, _ = iterate(m)
    Base.method_argnames(m)[2:end]
end

function get_handler(id, server_module)
    try
        getfield(server_module, Symbol(id))
    catch e
        if typeof(e) <: UndefVarError
            rethrow(ErrorException("You have not defined output for `$id`"))
        else
            rethrow(e)
        end
    end
end


text:
using Grist, MarketData, FactCheck

include("transaction.jl")
include("orderbook.jl")
include("blotter.jl")
include("portfolio.jl")
include("account.jl")
include("operators.jl")

exitstatus()


text:
module TIMIT

import ..LABELS
import ..LOCATIONS

using Glob

const URL_PHONE_MAP = "https://raw.githubusercontent.com/kaldi-asr/kaldi/master/egs/timit/s5/conf/phones.60-48-39.map"
const URL_DEV_SPKS = "https://raw.githubusercontent.com/kaldi-asr/kaldi/master/egs/timit/s5/conf/dev_spk.list"
const URL_TEST_SPKS = "https://raw.githubusercontent.com/kaldi-asr/kaldi/master/egs/timit/s5/conf/test_spk.list"

isspeechunit(u) = u ≠ LABELS[:sil] && u ≠ "h#" && u ≠ "epi" && u ≠ "pau"

function filter_phone(p)
    if p == "sil" return LABELS[:sil] end
    p
end

function prepare(datadir, rootdir)
    localdir = mkpath(joinpath(datadir, LOCATIONS[:local]))

    # phone mapping
    if ! ispath(joinpath(localdir, "phones.60-48-39.map"))
        @info "downloading phone mapping..."
        run(`wget -P $localdir $URL_PHONE_MAP`)
    end

    # Load the the phone mapping
    map_60_48 = Dict()
    map_60_39 = Dict()
    open(joinpath(localdir, "phones.60-48-39.map"), "r") do f
        for line in eachline(f)
            tokens = split(line)
            if length(tokens) == 3
                map_60_48[tokens[1]] = filter_phone(tokens[2])
                map_60_39[tokens[1]] = filter_phone(tokens[3])
            else
                map_60_48[tokens[1]] = nothing
                map_60_39[tokens[1]] = nothing
            end

        end
    end

    # speakers (dev)
    if ! ispath(joinpath(localdir, "dev_spk.list"))
        @info "downloading dev speaker list..."
        run(`wget -P $localdir $URL_DEV_SPKS`)
    end

    # speakers (test)
    if ! ispath(joinpath(localdir, "test_spk.list"))
        @info "downloading test speaker list..."
        run(`wget -P $localdir $URL_TEST_SPKS`)
    end

    # speakers (train)
    run(pipeline(`ls -d $(glob("dr*/*", "$rootdir/train"))`, `sed -e "s:^.*/::"`,
                 joinpath(localdir, "train_spk.list")))

    for dataset in ["train", "test", "dev"]
        datasetdir = mkpath(joinpath(datadir, dataset))

        @info "preparing $dataset"

        # Load the speakers
        spks = Set()
        open(joinpath(localdir, "$(dataset)_spk.list"), "r") do f
            for line in eachline(f)
                push!(spks, line)
            end
        end

        # Extract the utterance ids
        uttids = Dict()
        for path in glob("*/dr*/*", rootdir)
            spk = basename(path)
            spk in spks || continue

            for wavpath in glob("*wav", path)
                utt = splitext(basename(wavpath))[1]
                ! startswith(utt, "sa") || continue
                uttids["$(spk)_$(utt)"] = path
            end
        end

        # wav.scp
        scp = joinpath(datasetdir, LOCATIONS[:wavs])
        open(scp, "w") do f
            for (uttid, path) in uttids
                uttname = split(uttid, "_")[2]
                println(f, "$uttid\tsph2pipe -f wav $(joinpath(path, "$(uttname).wav")) |")
            end
        end

        # uttids
        uttidsfile = joinpath(datasetdir, LOCATIONS[:uttids])
        open(uttidsfile, "w") do f
            for uttid in keys(uttids) println(f, uttid) end
        end

        # uttids_speakers
        uttids_spk = joinpath(datasetdir, LOCATIONS[:utt2spk])
        open(uttids_spk, "w") do f
            for uttid in keys(uttids)
                spk = split(uttid, "_")[1]
                println(f, "$uttid\t$spk")
            end
        end

        # speakers
        spkfile = joinpath(datasetdir, LOCATIONS[:speakers])
        open(spkfile, "w") do f
            for spk in sort(collect(spks)) println(f, spk) end
        end

        # trans.phn.60
        open(joinpath(datasetdir, LOCATIONS[:trans]*".phn.60"), "w") do f
            for (uttid, path) in uttids
                print(f, uttid, "\t")
                uttname = split(uttid, "_")[2]
                open(joinpath(path, "$(uttname).phn"), "r") do f2
                    for line in eachline(f2)
                        print(f, split(line)[3], " ")
                    end
                end
                println(f)
            end
        end

        # trans.phn.48
        open(joinpath(datasetdir, LOCATIONS[:trans]*".phn.48"), "w") do f
            for (uttid, path) in uttids
                print(f, uttid, "\t")
                uttname = split(uttid, "_")[2]
                open(joinpath(path, "$(uttname).phn"), "r") do f2
                    for line in eachline(f2)
                        token = split(line)[3]
                        if token != "q"
                            print(f, map_60_48[token], " ")
                        end
                    end
                end
                println(f)
            end
        end

        # trans.phn.39
        open(joinpath(datasetdir, LOCATIONS[:trans]*".phn.39"), "w") do f
            for (uttid, path) in uttids
                print(f, uttid, " \t")
                uttname = split(uttid, "_")[2]
                open(joinpath(path, "$(uttname).phn"), "r") do f2
                    for line in eachline(f2)
                        token = split(line)[3]
                        if token != "q"
                            print(f, map_60_39[token], " ")
                        end
                    end
                end
                println(f)
            end
        end

        # ali.phn
        open(joinpath(datasetdir, LOCATIONS[:ali]*".phn.60"), "w") do f
            for (uttid, path) in uttids
                print(f, uttid, "\t")
                uttname = split(uttid, "_")[2]
                open(joinpath(path, "$(uttname).phn"), "r") do f2
                    for line in eachline(f2)
                        s, e, token = split(line)
                        s = Int(floor(parse(Int, s)/160))
                        e = Int(floor(parse(Int, e)/160))
                        print(f, strip("$token "^(e-s)), " ")
                    end
                end
                println(f)
            end
        end
    end

    # lang dir
    fn = sort ∘ collect
    for (ext, phones) in [(".60", fn(keys(map_60_48))),
                          (".48", fn(Set([map_60_48[k] for k in filter(u -> u ≠ "q", keys(map_60_48))]))),
                          (".39", fn(Set([map_60_39[k] for k in filter(u -> u ≠ "q", keys(map_60_39))])))]
        langdir = mkpath(joinpath(datadir, "lang$ext"))
        @info "preparing $langdir"

        open(joinpath(langdir, LOCATIONS[:units]), "w") do f
            for phone in phones
                println(f, phone, "\t", isspeechunit(phone) ? LABELS[:speechunit] : LABELS[:nonspeechunit])
            end
        end

        open(joinpath(langdir, LOCATIONS[:words]), "w") do f
            for phone in phones
                println(f, phone)
            end
        end

        open(joinpath(langdir, LOCATIONS[:lexicon]), "w") do f
            for phone in phones
                println(f, phone, "\t", phone)
            end
        end
    end
end

end

max_issues_repo_name: lucasondel/SpeechDatasets

max_forks_repo_name: lucasondel/SpeechDatasets


text:
struct DBSample
    database::SQLite.DB
    pkey::Int
    parent::Union{DBSample,Missing}
    owner::DBLaboratory
    name::String
    description::String
end

function Base.show(io::IO, samp::DBSample)
    helper(samp) = (ismissing(samp.parent) ? samp.name : helper(samp.parent)*" : "*samp.name)
    print(io, "$(samp.owner)'s $(helper(samp))")
end

function Base.write(db::SQLite.DB, ::Type{DBSample}, parentkey::Int, ownerkey::Int, name::String, desc::String)::Int
    stmt1 = SQLite.Stmt(db, "INSERT INTO SAMPLE ( PARENT, OWNER, NAME, DESCRIPTION ) VALUES (?, ?, ?, ? );")
    results = DBInterface.execute(stmt1, (parentkey, ownerkey, name, desc ))
    return  DBInterface.lastrowid(results)
end

function Base.findall(db::SQLite.DB, ::Type{DBSample}, owner::DBLaboratory; parent=missing, name=missing, desc=missing)
    bs = "SELECT PKEY FROM SAMPLE WHERE OWNER=?"
    args = Any[ owner.pkey ]
    if !ismissing(parent)
        bs *= " AND PARENT=?"
        push!(args, parent)
    end
    if !ismissing(name)
        bs *= " AND NAME=?"
        push!(args, name)
    end
    if !ismissing(desc)
        bs *= " AND DESCRIPTION=?"
        push!(args, desc)
    end
    stmt = SQLite.Stmt(db,bs*";")
    q = DBInterface.execute(stmt, args)
    return [ read(db, DBSample, r[:PKEY]) for r in q ]
end

function find(db::SQLite.DB, ::Type{DBSample}, owner::Int)::Vector{DBSample}
    stmt = SQLite.Stmt(db,"SELECT PKEY FROM SAMPLE WHERE OWNER=?;")
    q = DBInterface.execute(stmt, (owner,))
    return [ read(db, DBSample, r[:PKEY]) for r in q ]
end

function find(db::SQLite.DB, ::Type{DBSample}, owner::Int, parent::Int)::Vector{DBSample}
    stmt = SQLite.Stmt(db,"SELECT PKEY FROM SAMPLE WHERE OWNER=? AND PARENT=?;")
    q = DBInterface.execute(stmt, (owner,parent))
    return [ read(db, DBSample, r[:PKEY]) for r in q ]
end

function find(db::SQLite.DB, ::Type{DBSample}, owner::DBLaboratory, name::String)::Union{DBSample,Nothing}
    stmt = SQLite.Stmt(db,"SELECT PKEY FROM SAMPLE WHERE OWNER=? AND NAME=?;")
    q = DBInterface.execute(stmt, (owner.pkey, name))
    return SQLite.done(q) ? nothing : read(db, DBSample, SQLite.Row(q)[:PKEY])
end

Base.write(db::SQLite.DB, ::Type{DBSample}, ownerkey::Int, name::String, desc::String)::Int =
    write(db, DBSample, -1, ownerkey, name, desc)

Base.write(db::SQLite.DB, ::Type{DBSample}, parent::DBSample, owner::DBLaboratory, name::String, desc::String)::Int =
    write(db, DBSample, parent.pkey, owner.pkey, name, desc)

Base.write(db::SQLite.DB, ::Type{DBSample}, owner::DBLaboratory, name::String, desc::String)::Int =
    write(db, DBSample, owner.pkey, name, desc)

function Base.read(db::SQLite.DB, ::Type{DBSample}, pkey::Int)
    stmt1 = SQLite.Stmt(db, "SELECT * FROM SAMPLE WHERE PKEY=?;")
    q = DBInterface.execute(stmt1, (pkey, ))
    if SQLite.done(q)
        error("No sample found with pkey=$(pkey)")
    end
    r=SQLite.Row(q)
    parent = r[:PARENT] ≠ -1 ? read(db, DBSample, r[:PARENT]) : missing
    owner = read(db, DBLaboratory, r[:OWNER])
    return DBSample(db, r[:PKEY], parent, owner, r[:NAME], r[:DESCRIPTION])
end


text:
using JLD

include("perfplots.jl")

timeshift = 10.0
itershift = 1
nodeshift = 10

# produced by process_csv.jl in perfprofile mode
jldsource = joinpath(dirname(@__FILE__),"..",ARGS[1],"scale_perf.jld")

d = load(jldsource)
solvers = d["solvers"]
@assert solvers == ["PAJ_Gurobi_MOSEK_noscale_subponly_noinit","PAJ_Gurobi_MOSEK_msd_noscale_subponly_noinit","PAJ_Gurobi_MOSEK_scale_subponly_noinit","PAJ_Gurobi_MOSEK_msd_scale_subponly_noinit"]
instance_names = d["instance_names"]

iter_solvers = solvers[[1,3]]
iter_time_table = d["time_table"][:,[1,3]] + timeshift
iter_iter_table = d["itndcount_table"][:,[1,3]] + itershift

msd_solvers = solvers[[2,4]]
msd_time_table = d["time_table"][:,[2,4]] + timeshift
msd_node_table = d["itndcount_table"][:,[2,4]] + nodeshift

performance_profile(iter_time_table, iter_solvers, logscale=false, ymax=1.0, xmax=2, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"..",ARGS[1],"scale_iter_time.tex"))

performance_profile(iter_iter_table, iter_solvers, logscale=false, ymax=1.0, xmax=2, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"..",ARGS[1],"scale_iter_iters.tex"))

performance_profile(msd_time_table, msd_solvers, logscale=false, ymax=1.0, xmax=2, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"..",ARGS[1],"scale_msd_time.tex"))

performance_profile(msd_node_table, msd_solvers, logscale=false, ymax=1.0, xmax=6, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"..",ARGS[1],"scale_msd_nodes.tex"))

max_stars_repo_name: chriscoey/PajaritoSupplement

max_issues_repo_name: chriscoey/PajaritoSupplement

max_forks_repo_name: chriscoey/PajaritoSupplement


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("World of machine learning")
cd(datapath("mnist"))

# # Load and minibatch MNIST data
# (c) Deniz Yuret, 2018

#-

# * Objective: Learning the structure of the MNIST dataset and usage of the Knet.Data struct.
# * Prerequisites: [The iteration interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1)
# * Knet: mnistdata, mnistview (used and explained)
# * Knet: dir, minibatch (used by mnist.jl)

#-

## This loads MNIST handwritten digit recognition dataset.
## Traininig and test data go to variables dtrn and dtst
using Knet: Knet, minibatch
include(Knet.dir("data","mnist.jl"))  # defines mnistdata and mnistview
dtrn,dtst = mnistdata(xtype=Array{Float32});

#-

## dtrn contains 600 minibatches of 100 images (total 60000)
## dtst contains 100 minibatches of 100 images (total 10000)
length.((dtrn,dtst))

#-

## Each minibatch is an (x,y) pair where x is 100 28x28x1 images and y contains their labels.
## Here is the first minibatch in the test set:
(x,y) = first(dtst)
summary.((x,y))

#-

## Here is the first five images from the test set:
using Images, ImageMagick
hcat([mnistview(x,i) for i=1:5]...)

#-

## Here are their labels (0x0a=10 is used to represent 0)
@show y[1:5];

#-

## dtrn and dtst are implemented as Julia iterators (see https://docs.julialang.org/en/v1/manual/interfaces)
## This means they can be used in for loops, i.e. `for (x,y) in dtst`
cnt = zeros(Int,10)
for (x,y) in dtst
    for label in y
        cnt[label] += 1
    end
end
@show cnt;

# # Convolutional Neural Networks
# (c) Deniz Yuret, 2018

#-

# * Objectives: See the effect of sparse and shared weights implemented by convolutional networks.
# * Prerequisites: MLP models (04.mlp.ipynb), KnetArray, param, param0, dropout, relu, nll
# * Knet: conv4, pool, mat (explained)
# * Knet: dir, gpu, minibatch, KnetArray (used by mnist.jl)
# * Knet: SGD, train!, Train, load, save (used by trainresults)

# ## Introduction to convolution

## Convolution operator in Knet
using Knet: conv4
@doc conv4

#-

## Convolution in 1-D
@show w = reshape([1.0,2.0,3.0], (3,1,1,1))
@show x = reshape([1.0:7.0...], (7,1,1,1))
@show y = conv4(w, x);  # size Y = X - W + 1 = 5 by default

#-

## Padding
@show y2 = conv4(w, x, padding=(1,0));  # size Y = X + 2P - W + 1 = 7 with padding=1
## To preserve input size (Y=X) for a given W, what padding P should we use?

#-

## Stride
@show y3 = conv4(w, x; padding=(1,0), stride=3);  # size Y = 1 + floor((X+2P-W)/S)

#-

## Mode
@show y4 = conv4(w, x, mode=0);  # Default mode (convolution) inverts w
@show y5 = conv4(w, x, mode=1);  # mode=1 (cross-correlation) does not invert w

#-

## Convolution in more dimensions
x = reshape([1.0:9.0...], (3,3,1,1))

#-

w = reshape([1.0:4.0...], (2,2,1,1))

#-

y = conv4(w, x)

#-

## Convolution with multiple channels, filters, and instances
## size X = [X1,X2,...,Xd,Cx,N] where d is the number of dimensions, Cx is channels, N is instances
x = reshape([1.0:18.0...], (3,3,2,1))

#-

## size W = [W1,W2,...,Wd,Cx,Cy] where d is the number of dimensions, Cx is input channels, Cy is output channels
w = reshape([1.0:24.0...], (2,2,2,3));

#-

## size Y = [Y1,Y2,...,Yd,Cy,N]  where Yi = 1 + floor((Xi+2Pi-Wi)/Si), Cy is channels, N is instances
y = conv4(w,x)

# See http://cs231n.github.io/assets/conv-demo/index.html for an animated example.

#-

# ## Introduction to Pooling

## Pooling operator in Knet
using Knet: pool
@doc pool

#-

## 1-D pooling example
@show x = reshape([1.0:6.0...], (6,1,1,1))
@show pool(x);

#-

## Window size
@show pool(x; window=3);  # size Y = floor(X/W)

#-

## Padding
@show pool(x; padding=(1,0));  # size Y = floor((X+2P)/W)

#-

## Stride
@show x = reshape([1.0:10.0...], (10,1,1,1));
@show pool(x; stride=4);  # size Y = 1 + floor((X+2P-W)/S)

#-

## Mode (using Array here; not all modes are implemented on the CPU)
using Knet: KnetArray
x = Array(reshape([1.0:6.0...], (6,1,1,1)))
@show x
@show pool(x; padding=(1,0), mode=0)  # max pooling
@show pool(x; padding=(1,0), mode=1)  # avg pooling
## @show pool(x; padding=(1,0), mode=2); # avg pooling excluding padded values (is not implemented on CPU so will error)

#-

## More dimensions
x = reshape([1.0:16.0...], (4,4,1,1))

#-

pool(x)

#-

## Multiple channels and instances
x = reshape([1.0:32.0...], (4,4,2,1))

#-

## each channel and each instance is pooled separately
pool(x)  # size Y = (Y1,...,Yd,Cx,N) where Yi are spatial dims, Cx and N are identical to input X

# ## Experiment setup

## Load data (see 02.mnist.ipynb)
using Knet: Knet, KnetArray, gpu, minibatch
include(Knet.dir("data","mnist.jl"))  # Load data
dtrn,dtst = mnistdata();              # dtrn and dtst = [ (x1,y1), (x2,y2), ... ] where xi,yi are minibatches of 100

#-

(x,y) = first(dtst)
summary.((x,y))

#-

## For running experiments
using Knet: SGD, train!, nll, zeroone
import ProgressMeter

function trainresults(file,model; o...)
    if (print("Train from scratch? ");readline()[1]=='y')
        results = Float64[]; updates = 0; prog = ProgressMeter.Progress(60000)
        function callback(J)
            if updates % 600 == 0
                push!(results, nll(model,dtrn), nll(model,dtst), zeroone(model,dtrn), zeroone(model,dtst))
                ProgressMeter.update!(prog, updates)
            end
            return (updates += 1) <= 60000
        end
        train!(model, dtrn; callback=callback, optimizer=SGD(lr=0.1), o...)
        Knet.save(file,"results",reshape(results, (4,:)))
    end
    isfile(file) || download("http://people.csail.mit.edu/deniz/models/tutorial/$file",file)
    results = Knet.load(file,"results")
    println(minimum(results,dims=2))
    return results
end

# ## A convolutional neural network model for MNIST

## Redefine Linear layer (See 03.lin.ipynb):
using Knet: param, param0
struct Linear; w; b; end
(f::Linear)(x) = (f.w * mat(x) .+ f.b)
mat(x)=reshape(x,:,size(x)[end])  # Reshapes 4-D tensor to 2-D matrix so we can use matmul
Linear(inputsize::Int,outputsize::Int) = Linear(param(outputsize,inputsize),param0(outputsize))

#-

## Define a convolutional layer:
struct Conv; w; b; end
(f::Conv)(x) = pool(conv4(f.w,x) .+ f.b)
Conv(w1,w2,cx,cy) = Conv(param(w1,w2,cx,cy), param0(1,1,cy,1))

#-

## Define a convolutional neural network:
struct CNN; layers; end

#-

## Weight initialization for a multi-layer convolutional neural network
## h[i] is an integer for a fully connected layer, a triple of integers for convolution filters and tensor inputs
## use CNN(x,h1,h2,...,hn,y) for a n hidden layer model
function CNN(h...)
    w = Any[]
    x = h[1]
    for i=2:length(h)
        if isa(h[i],Tuple)
            (x1,x2,cx) = x
            (w1,w2,cy) = h[i]
            push!(w, Conv(w1,w2,cx,cy))
            x = ((x1-w1+1)÷2,(x2-w2+1)÷2,cy) # assuming conv4 with p=0, s=1 and pool with p=0,w=s=2
        elseif isa(h[i],Integer)
            push!(w, Linear(prod(x),h[i]))
            x = h[i]
        else
            error("Unknown layer type: $(h[i])")
        end
    end
    CNN(w)
end;

#-

using Knet: dropout, relu
function (m::CNN)(x; pdrop=0)
    for (i,layer) in enumerate(m.layers)
        p = (i <= length(pdrop) ? pdrop[i] : pdrop[end])
        x = dropout(x, p)
        x = layer(x)
        x = (layer == m.layers[end] ? x : relu.(x))
    end
    return x
end

#-

lenet = CNN((28,28,1), (5,5,20), (5,5,50), 500, 10)
summary.(l.w for l in lenet.layers)

#-

using Knet: nll
(x,y) = first(dtst)
nll(lenet,x,y)

# ## CNN vs MLP

using Plots; default(fmt=:png,ls=:auto)
ENV["COLUMNS"] = 92

#-

@time cnn = trainresults("cnn.jld2", lenet; pdrop=(0,0,.3)); # 406s [8.83583e-5, 0.017289, 0.0, 0.0048]
## Note that training will take a very long time without a GPU

#-

isfile("mlp.jld2") || download("http://people.csail.mit.edu/deniz/models/tutorial/mlp.jld2","mlp.jld2")
mlp = Knet.load("mlp.jld2","results");

#-

## Comparison to MLP shows faster convergence, better generalization
plot([mlp[1,:], mlp[2,:], cnn[1,:], cnn[2,:]],ylim=(0.0,0.1),
     labels=[:trnMLP :tstMLP :trnCNN :tstCNN],xlabel="Epochs",ylabel="Loss")

#-

plot([mlp[3,:], mlp[4,:], cnn[3,:], cnn[4,:]],ylim=(0.0,0.03),
    labels=[:trnMLP :tstMLP :trnCNN :tstCNN],xlabel="Epochs",ylabel="Error")

# ## Convolution vs Matrix Multiplication

## Convolution and matrix multiplication can be implemented in terms of each other.
## Convolutional networks have no additional representational power, only statistical efficiency.
## Our original 1-D example
@show w = reshape([1.0,2.0,3.0], (3,1,1,1))
@show x = reshape([1.0:7.0...], (7,1,1,1))
@show y = conv4(w, x);  # size Y = X - W + 1 = 5 by default

#-

## Convolution as matrix multiplication (1)
## Turn w into a (Y,X) sparse matrix
w2 = Float64[3 2 1 0 0 0 0; 0 3 2 1 0 0 0; 0 0 3 2 1 0 0; 0 0 0 3 2 1 0; 0 0 0 0 3 2 1]

#-

@show y2 = w2 * mat(x);

#-

## Convolution as matrix multiplication (2)
## Turn x into a (W,Y) dense matrix (aka the im2col operation)
## This is used to speed up convolution with known efficient matmul algorithms
x3 = Float64[1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7]

#-

@show w3 = [3.0 2.0 1.0]
@show y3 = w3 * x3;

#-

## Matrix multiplication as convolution
## This could be used to make a fully connected network accept variable sized inputs.
w = reshape([1.0:6.0...], (2,3))

#-

x = reshape([1.0:3.0...], (3,1))

#-

y = w * x

#-

## Consider w with size (Y,X)
## Treat each of the Y rows of w as a convolution filter
w2 = copy(reshape(Array(w)', (3,1,1,2)))

#-

## Reshape x for convolution
x2 = reshape(x, (3,1,1,1))

#-

## Use conv4 for matrix multiplication
y2 = conv4(w2, x2; mode=1)

#-

## So there is no difference between the class of functions representable with an MLP vs CNN.
## Sparse connections and weight sharing give CNNs more generalization power with images.
## Number of parameters in MLP256: (256x784)+256+(10x256)+10 = 203530
## Number of parameters in LeNet: (5*5*1*20)+20+(5*5*20*50)+50+(500*800)+500+(10*500)+10 = 431080

max_stars_repo_path: Courses/World of machine learning/02.mnist.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
export iauP2pv
"""
Extend a p-vector to a pv-vector by appending a zero velocity.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   p        double[3]       p-vector

Returned:
   pv       double[2][3]    pv-vector

Called:
   iauCp        copy p-vector
   iauZp        zero p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauP2pv(double p[3], double pv[2][3])
function iauP2pv(p::AbstractVector{<:Real})
   pv = zeros(Float64, 3, 2)

   ccall((:iauP2pv, libsofa_c), Cvoid,
       (Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 1}, p),
       pv)

   return SMatrix{2,3}(pv')
end


text:
macro qasm_str(s::String)
    s = Meta.parse("\"\"\"\n$s\"\"\"")
    s isa String && return esc(qasm_str_m(__module__, __source__, s))
    args = map(s.args) do x
        x isa String && return x
        Base.eval(__module__, x)
    end
    return esc(qasm_str_m(__module__, __source__, join(args)))
end

function qasm_str_m(m::Module, lino::LineNumberNode, src::String)
    ast = OpenQASM.parse(src)
    return transpile(m, lino, ast)
end

mutable struct VirtualRegister
    type::Symbol
    address::UnitRange{Int}
end

mutable struct RegisterRecord
    map::Dict{String,VirtualRegister}
    nqubits::Int
    ncbits::Int
end

mutable struct GateRegisterRecord
    map::Dict
    total::Int
end

Base.getindex(x::RegisterRecord, key) = x.map[key]
Base.getindex(x::VirtualRegister, xs...) = x.address[xs...]

RegisterRecord() = RegisterRecord(Dict{String,UnitRange{Int}}(), 0, 0)
GateRegisterRecord() = GateRegisterRecord(Dict(), 0)

function scan_registers(ast::MainProgram)
    return scan_registers!(RegisterRecord(), ast)
end

function scan_registers!(record::RegisterRecord, ast::MainProgram)
    for stmt in ast.prog
        scan_registers!(record, stmt)
    end
    return record
end

function scan_registers!(record::RegisterRecord, ast::RegDecl)
    nbits = transpile_token(ast.size)
    nqubits = record.nqubits
    ncbits = record.ncbits

    if ast.type.str == "qreg"
        record.map[ast.name.str] = VirtualRegister(:quantum, (nqubits+1):(nqubits+nbits))
        record.nqubits += nbits
    else # classical
        record.map[ast.name.str] = VirtualRegister(:classical, (ncbits+1):(ncbits+nbits))
        record.ncbits += nbits
    end
    return record
end

scan_registers!(record::RegisterRecord, ast) = record

struct TranspileCtx
    m::Module
    source::LineNumberNode
    record::Any
    returns::Dict{Symbol, Vector{Symbol}}
end

function transpile(m::Module, lino::LineNumberNode, ast::MainProgram)
    @assert v"2.0.0" ≤ ast.version < v"3.0.0"

    code = Expr(:block)
    body = Expr(:block)
    routines = []
    record = scan_registers(ast)
    ctx = TranspileCtx(m, lino, record, Dict{Symbol, Vector{Any}}())

    for stmt in ast.prog
        if stmt isa RegDecl
            continue
        elseif stmt isa Types.Gate
            push!(routines, transpile(ctx, stmt))
        elseif stmt isa Include
            file = stmt.file.str[2:end-1]
            # use relative path to current file if not in REPL
            # isinteractive can be true in IDEs
            if !isnothing(lino.file) && !(isinteractive() && isempty(PROGRAM_FILE))
                file = joinpath(dirname(string(lino.file)), file)
            end
            source = read(file, String)
            ast = OpenQASM.parse(source)
            push!(code.args, transpile(m, lino, ast))
        else
            ex = transpile(ctx, stmt)
            if !isnothing(ex)
                push!(body.args, ex)
            end
        end
    end

    # if there are classical registers
    # return them in a NamedTuple
    ret = Expr(:tuple)
    if !isempty(ctx.returns)
        for (name, list) in ctx.returns
            if length(list) == 1
                push!(ret.args, Expr(:(=), name, list[1]))
            else
                push!(ret.args, Expr(:(=), name, Expr(:tuple, list...)))
            end
        end
    end

    if !isempty(body.args)
        if isempty(ret.args)
            push!(body.args, :(return))
        else
            push!(body.args, :(return $ret))
        end
    end

    # routines
    for each in routines
        push!(code.args, each)
    end

    if !isempty(body.args)
        def = JLFunction(name=gensym(:qasm), body=body)
        push!(code.args, YaoCompiler.codegen_routine(def))
    end
    return code
end

function transpile(ctx::TranspileCtx, stmt::Gate)
    name = transpile(ctx, stmt.decl.name)
    args = transpile_exp(ctx, stmt.decl.cargs)
    record = transpile_gate_registers(stmt.decl.qargs)
    body = Expr(:block)
    new_ctx = TranspileCtx(ctx.m, ctx.source, record, Dict{Symbol, Vector{Any}}())

    for each in stmt.body
        push!(body.args, transpile(new_ctx, each))
    end
    def = JLFunction(;name, args, body)
    return YaoCompiler.codegen_routine(def)
end

function transpile_gate_registers(stmt::Vector)
    record = GateRegisterRecord()
    for each in stmt
        haskey(record.map, each.str) && throw(Meta.ParseError("duplicated register name $(each.str)"))
        record.total += 1
        record.map[each.str] = record.total
    end
    return record
end

function transpile(ctx::TranspileCtx, stmt)
    # NOTE: apply etc. are treated as keyword inside device
    # no need to put Intrinsics infront of them
    @switch stmt begin
        @case UGate(z1, y, z2, qarg)
            locs = transpile(ctx, qarg)
            quote
                apply($Intrinsics.Rz($(transpile_exp(ctx, z1))), $locs)
                apply($Intrinsics.Ry($(transpile_exp(ctx, y))), $locs)
                apply($Intrinsics.Rz($(transpile_exp(ctx, z2))), $locs)
            end

        @case CXGate(ctrl, qarg)
            locs = transpile(ctx, qarg)
            ctrl = transpile(ctx, ctrl)
            :(apply($Intrinsics.X, $locs, $ctrl))

        @case IfStmt(left, right, body)
            quote
                if $(transpile(ctx, left)) == $(transpile(ctx, right))
                    $(transpile(ctx, body))
                end
            end

        # measure qreg -> creg
        @case Measure(qarg, Bit(name, nothing))
            ctx.returns[Symbol(name.str)] = Symbol[Symbol(name.str)]
            :($(Symbol(name.str)) = measure($(transpile(ctx, qarg))))

        @case Measure(qarg, Bit(name, address))
            ret = get!(ctx.returns, Symbol(name.str), Symbol[])
            name = Symbol(name.str, "_", transpile(ctx, address) + 1)
            push!(ret, name)
            :($name = measure($(transpile(ctx, qarg))))

        @case Barrier(qargs)
            :(barrier($(transpile_qargs(ctx, qargs))))

        @case Instruction(op, cargs, qargs)
            args = transpile_exp(ctx, cargs)
            locs = transpile_qargs(ctx, qargs)
            @match op begin
                "x" => :(apply($Intrinsics.X, $locs))
                "y" => :(apply($Intrinsics.Y, $locs))
                "z" => :(apply($Intrinsics.Z, $locs))
                "h" => :(apply($Intrinsics.H, $locs))
                "s" => :(apply($Intrinsics.S, $locs))
                _ => begin
                    gate = Expr(:call, GlobalRef(ctx.m, Symbol(op)), transpile_exp(ctx, cargs)...)
                    :(apply($gate, $locs))
                end
            end

        # tokens
        @case Bit(name, address)
            ctx.record isa RegisterRecord || return Locations(ctx.record.map[name.str])
            r = ctx.record[name.str]
            r.type === :classical && return Symbol(name.str)
            isnothing(address) && return Locations(r[:])
            address = transpile(ctx, address)
            Locations(r[address+1])

        @case ::Token
            transpile_token(stmt)

        @case _
            error("unknown statement: $stmt")
    end
end

function transpile_qargs(ctx::TranspileCtx, stmts::Vector)
    locs = map(stmts) do stmt
        transpile(ctx, stmt)
    end
    return merge_locations(locs...)
end

function transpile_exp(ctx::TranspileCtx, stmt)
    @switch stmt begin
        @case Call(name, args)
            Expr(:call, name, transpile_exp(ctx, args))
        @case Neg(val)
            Expr(:call, :-, transpile_exp(ctx, val))
        @case Token{:reserved}(_, _, _, "pi", _)
            return π
        @case (Token{:reserved}(_, _, _, "pi", _), Token{:reserved}(_, _, _, "/", _), Token{:int}(_, _, _, "2", _))
            Expr(:call, GlobalRef(Base, :(/)), π, 2)
        @case (lhs, op::Token, rhs)
            Expr(:call, Symbol(op.str), transpile_exp(ctx, lhs), transpile_exp(ctx, rhs))
        @case ::Vector
            map(x->transpile_exp(ctx, x), stmt)
        @case ::Token
            transpile_token(stmt)
        @case _
            error("unknown expression: $stmt")
    end
end

function transpile_token(stmt)
    @switch stmt begin
        @case ::Token{:unnamed}
            Symbol(stmt.str)
        @case ::Token{:reserved}
            Symbol(stmt.str)
        @case ::Token{:id}
            stmt.str == "pi" && return Base.pi
            return Symbol(stmt.str)
        @case ::Token{:float64}
            return Base.parse(Float64, stmt.str)
        @case ::Token{:int}
            return Base.parse(Int, stmt.str)
        @case ::Token
            error("unknown token: $stmt")
    end
end


text:
using StaticArrays

x = @SVector rand(3)
typeof(x.data)

struct MySVec{N, T} <: StaticVector{N, T}
   data::NTuple{N, T}
end

@inline Base.getindex(x::MySVec{N, T}, i::Int) where {N, T} = x.data[i]
StaticArrays.similar_type(x::MySVec{N, T}) where {N, T} = zero(MySVec{N, T})
StaticArrays.similar_type(VT::Type{MySVec{N, T}}) where {N, T} = zero(VT)

y = rand(MySVec{3, Float64})
z = rand(MySVec{3, Float64})

y + z

using StaticArrays

ACE.Invariant( (@SVector rand(3)) )

A2B = [ ACE.Invariant( (@SVector rand(3)) )
        for i = 1:10, j = 1:10 ]
AA = rand(10)
A2B * AA

using SparseArrays
A2Bsp = sparse(A2B)
A2Bsp * AA == A2B * AA


text:
"""
    function opf_PTDF(sys::System; kwargs...)

Solves the Optimal Power Flow (OPF) of a power system using the PTDF matrix and returns the
Locational Marginal Price (lmp) and the results of the optimisation (res).

#Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)

#Keywords
- `network::DataType = StandardPTDFModel`:  Network Representation of the Power System
- `solver= optimizer_with_attributes
    (Ipopt.Optimizer)`:                     Solver to be used for OPF
"""
function opf_PTDF(
    sys::System;
    network::DataType = StandardPTDFModel,
    solver = optimizer_with_attributes(Ipopt.Optimizer)
    )
    # Solve the OPF and calculates the lmp

    # Branch representation
    branches = Dict{Symbol,DeviceModel}(
        :L => DeviceModel(Line, StaticLine),
        :ML => DeviceModel(MonitoredLine, StaticLine),
    )
    # Devices representation
    devices = Dict(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
        :HydroROR => DeviceModel(HydroDispatch, FixedOutput),
        :RenFx => DeviceModel(RenewableFix, FixedOutput),
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    )
    # Services representation
    services = Dict(
        :ReserveUp => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :ReserveDown => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
    )

    # Select Template for Optimisation problem
    template_opt = OperationsProblemTemplate(network, devices, branches, services)

    # Calculate PTDF
    PTDF_matrix = PTDF(sys)

    #Save the information to the Operations Problem
    problem = OperationsProblem(
        GenericOpProblem,
        template_opt,
        sys,
        horizon = 1,
        #use_forecast_data = false,
        optimizer = solver,
        balance_slack_variables = true,
        constraint_duals = [:CopperPlateBalance, :network_flow,
        :P_ub__ThermalStandard__RangeConstraint, :P_lb__ThermalStandard__RangeConstraint],
        PTDF = PTDF_matrix,
    )
    # Solve optimisation problem
    res = solve!(problem)

    # Computing final LMP by subtracting the duals (μ) of network_flow constraint ...
    # multiplied by the PTDF matrix from the dual (λ) of CopperPlateBalance constraint.

    # Convert the results from DataFrame to Array for ease of use.
    λ = convert(Array, res.dual_values[:CopperPlateBalance])
    μ = convert(Array, res.dual_values[:network_flow])

    # Calculate congestion component of the LMP as a product of μ and the PTDF matrix.
    buses = get_components(Bus, sys)
    congestion_lmp = Dict()
    for bus in buses
        congestion_lmp[get_name(bus)] = μ * PTDF_matrix[:, get_number(bus)]
    end
    congestion_lmp = DataFrame(congestion_lmp)

    # Get the LMP for each node in a lossless DC-OPF using the PTDF formulation.
    # resize!(lmp,counter);
    lmp = λ .- congestion_lmp

    return lmp, res
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
module LVAChannels

export Stim, run_lva, LVAParam

using Sundials
using SIUnits
using SIUnits.ShortUnits

const mS = Milli*Siemens

immutable LVAParam
    c_m       # membrane spacific capacitance
    gbar      # Ca channel cross-membrane conductivity
    eca       # Ca channel reversal potential
    gl        # leak conductivity
    el        # leak reversal potential
    q10_1     # rate scaling for 'm' activation gate
    q10_2     # rate scaling for 'dsh' inactivation gate
    vrest     # (as hoc) resting potential

    # constructor with default values, with q10 values
    # corresponding to the body temperature current-clamp experiments
    LVAParam(;
        c_m    = 0.01F*m^-2,
        gbar   = 0.2mS*cm^-2,
        eca    = 120mV,
        gl     = .1mS*cm^-2,
        el     = -65mV,
        q10_1  = 5,
        q10_2  = 3,
#        vrest  = -61.47mV
        vrest  = -63mV
    ) = new(c_m, gbar, eca, gl, el, q10_1, q10_2, vrest)

end

immutable Stim
    t0        # start time of stimulus
    t1        # stop time of stimulus
    i_e       # stimulus current density

    Stim() = new(0s, 0s, 0A/m^2)
    Stim(t0, t1, i_e) = new(t0, t1, i_e)
end

# 'm' activation gate
function m_lims(v, q10)
    quotient = 1+exp(-(v+63mV)/7.8mV)
    mtau = 1ms*(1.7+exp(-(v+28.8mV)/13.5mV))/(q10*quotient)
    minf = 1/quotient
    return mtau, minf
end

# 'dsh' inactivation gate:
# d <-> s (alpha2, beta2)
# s <-> h (alpha1, beta1)
# subject to s=1-h-d
function dsh_lims(v, q10)
    k = sqrt(0.25+exp((v+83.5mV)/6.3mV))-0.5
    alpha1 = q10*exp(-(v+160.3mV)/17.8mV)/1ms
    beta1 = alpha1*k

    tau2 = 240ms/(1+exp((v+37.4mV)/30mV))/q10
    alpha2 = 1/tau2/(1+k)
    beta2 = alpha2*k

    hinf = 1/(1+k+k^2)
    dinf = k^2*hinf

    return alpha1, beta1, alpha2, beta2, dinf, hinf
end

# Choose initial conditions for the system such that the gating variables
# are at steady state for the user-specified voltage v.
function initial_conditions(v, q10_1, q10_2)
    mtau, minf = m_lims(v, q10_1)
    alpha1, beta1, alpha2, beta2, dinf, hinf = dsh_lims(v, q10_2)

    return (v, minf, dinf, hinf)
end

# Given time t and state (v, m, d, h),
# return (vdot, mdot, ddot, hdot).
function f(t, state; p=LVAParam(), istim=0A/m^2)
    v, m, d, h = state

    ica = p.gbar*m^3*h*(v - p.eca)
    il = p.gl*(v - p.el)
    itot = ica + il - istim

    # Calculate the voltage dependent rates for the gating variables.
    mtau, minf = m_lims(v, p.q10_1)
    alpha1, beta1, alpha2, beta2, dinf, hinf = dsh_lims(v, p.q10_2)

    mdot = (minf-m)/mtau
    hdot = alpha1*(1-h-d)-beta1*h
    ddot = beta2*(1-h-d)-alpha2*d

    return (-itot/p.c_m, mdot, ddot, hdot)
end

function make_range(t_start, dt, t_end)
    r = collect(t_start: dt: t_end)
    if length(r)>0 && r[length(r)]<t_end
        push!(r, t_end)
    end
    return r
end

function run_lva(t_end; stim=Stim(), param=LVAParam(), sample_dt=0.01ms)
    v_scale = 1V
    t_scale = 1s

    v0, m0, h0, d0 = initial_conditions(param.vrest, param.q10_1, param.q10_2)
    y0 = [ v0/v_scale, m0, h0, d0 ]

    fbis(t, y, ydot, istim) = begin
        vdot, mdot, hdot, ddot =
            f(t*t_scale, (y[1]*v_scale, y[2], y[3], y[4]), istim=istim, p=param)

        ydot[1], ydot[2], ydot[3], ydot[4] =
            vdot*t_scale/v_scale, mdot*t_scale, hdot*t_scale, ddot*t_scale

        return Sundials.CV_SUCCESS
    end

    fbis_nostim(t, y, ydot) = fbis(t, y, ydot, 0A/m^2)
    fbis_stim(t, y, ydot) = fbis(t, y, ydot, stim.i_e)

    # Ideally would run with vector absolute tolerance to account for v_scale,
    # but this would prevent us using the nice cvode wrapper.

    res = []
    samples = []

    t1 = clamp(stim.t0, 0s, t_end)
    if t1>0s
        ts = make_range(0s, sample_dt, t1)
        r = Sundials.cvode(fbis_nostim, y0, map(t->t/t_scale, ts), abstol=1e-6, reltol=5e-10)
        y0 = vec(r[size(r)[1], :])
        push!(res, r)
        push!(samples, ts)
    end
    t2 = clamp(stim.t1, t1, t_end)
    if t2>t1
        ts = make_range(t1, sample_dt, t2)
        r = Sundials.cvode(fbis_stim, y0, map(t->t/t_scale, ts), abstol=1e-6, reltol=5e-10)
        y0 = vec(r[size(r)[1], :])
        push!(res, r)
        push!(samples, ts)
    end
    if t_end>t2
        ts = make_range(t2, sample_dt, t_end)
        r = Sundials.cvode(fbis_nostim, y0, map(t->t/t_scale, ts), abstol=1e-6, reltol=5e-10)
        y0 = vec(r[size(r)[1], :])
        push!(res, r)
        push!(samples, ts)
    end

    res = vcat(res...)
    samples = vcat(samples...)

    # Use map here because of issues with type deduction with arrays and SIUnits.
    return samples, map(v->v*v_scale, res[:, 1]), res[:, 2], res[:, 3], res[:, 4]
end

end # module LVAChannels

max_stars_repo_name: kabicm/arbor

max_issues_repo_name: kabicm/arbor

max_forks_repo_name: kabicm/arbor


text:
@defcomp PermafrostJULES begin
    rt_g = Parameter(index=[time],unit="degreeC")
    y_year_0 = Parameter(unit="year")
    y_year = Parameter(index=[time], unit="year")

    perm_jul_af = Parameter(default=1.9359078717122433)
    perm_jul_limit_t = Parameter(default=18.86448815404626)
    perm_jul_sens_slope_vs_t_co2 = Parameter(default=-0.12186871366342877)
    perm_jul_lag_slope_vs_t_co2 = Parameter(default=-0.6550115433585014)
    perm_jul_pow_slope_vs_t_co2 = Parameter(default=1.618876347278246)
    perm_jul_sens_c_co2 = Parameter(default=61867.779459651385)
    perm_jul_limit_c = Parameter(default=738145.3283931707)
    perm_jul_ce_c_co2_0 = Parameter(default=20040.99686065949)
    perm_jul_pow_c_co2 = Parameter(default=0.457030756829357)
    perm_jul_lag_c_co2 = Parameter(default=543.6163557759747)
    perm_jul_ch4_co2_c_ratio = Parameter(default=6.116162272246744)

    perm_jul_temp = Variable(index=[time], unit="degreeC")
    perm_jul_sens_c_co2_correct = Variable(index=[time])
    perm_jul_lag_c_co2_correct = Variable(index=[time])
    perm_jul_pow_c_co2_correct = Variable(index=[time])
    perm_jul_equilib_c_co2 = Variable(index=[time], unit="Mtonne C")
    perm_jul_ce_c_co2 = Variable(index=[time], unit="Mtonne C")
    perm_jul_ce_co2 = Variable(index=[time], unit="Mtonne CO2")
    perm_jul_e_co2 = Variable(index=[time], unit="MtonCO2/yr")
    perm_jul_ce_ch4 = Variable(index=[time], unit="Mtonne CH4")

    function run_timestep(p, v, d, tt)
        if is_first(tt)
            yp = p.y_year[tt] - p.y_year_0
        else
            yp = p.y_year[tt] - p.y_year[tt-1]
        end

        # Permafrost temperature
        v.perm_jul_temp[tt] = p.rt_g[tt] * p.perm_jul_af
        # Correction for perm sensitivity, CO2
        v.perm_jul_sens_c_co2_correct[tt] = (v.perm_jul_temp[tt]/(0.5*p.perm_jul_limit_t))^p.perm_jul_sens_slope_vs_t_co2
        # Correction for perm lag, CO2
        v.perm_jul_lag_c_co2_correct[tt] = (v.perm_jul_temp[tt]/(0.5*p.perm_jul_limit_t))^p.perm_jul_lag_slope_vs_t_co2
        # Correction for perm power, CO2
        v.perm_jul_pow_c_co2_correct[tt] = ifelse(1 + p.perm_jul_pow_slope_vs_t_co2 * (v.perm_jul_temp[tt] - 0.5*p.perm_jul_limit_t) /p.perm_jul_limit_t <= 0, 0.000001, 1 + p.perm_jul_pow_slope_vs_t_co2 * (v.perm_jul_temp[tt] - 0.5*p.perm_jul_limit_t) /p.perm_jul_limit_t)
        # Equilibrium perm cumul carbon, CO2
        v.perm_jul_equilib_c_co2[tt] = ifelse((p.perm_jul_sens_c_co2 * v.perm_jul_sens_c_co2_correct[tt]) * v.perm_jul_temp[tt] < p.perm_jul_limit_c, (p.perm_jul_sens_c_co2 * v.perm_jul_sens_c_co2_correct[tt]) * v.perm_jul_temp[tt], p.perm_jul_limit_c)
        # Released perm cumul carbon, CO2
        if is_first(tt)
            v.perm_jul_ce_c_co2[tt] = ((v.perm_jul_equilib_c_co2[tt] < p.perm_jul_ce_c_co2_0) ? p.perm_jul_ce_c_co2_0 : ((((v.perm_jul_equilib_c_co2[tt] - p.perm_jul_ce_c_co2_0)/p.perm_jul_limit_c)^(1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) < (1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) * yp / (p.perm_jul_lag_c_co2*v.perm_jul_lag_c_co2_correct[tt])) ? v.perm_jul_equilib_c_co2[tt] : (v.perm_jul_equilib_c_co2[tt] - p.perm_jul_limit_c * ((((v.perm_jul_equilib_c_co2[tt] - p.perm_jul_ce_c_co2_0)/p.perm_jul_limit_c)^(1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) - (1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) * yp / (p.perm_jul_lag_c_co2*v.perm_jul_lag_c_co2_correct[tt]))^(1/(1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]))))))
        else
            v.perm_jul_ce_c_co2[tt] = ((v.perm_jul_equilib_c_co2[tt] < v.perm_jul_ce_c_co2[tt-1]) ? v.perm_jul_ce_c_co2[tt-1] : ((((v.perm_jul_equilib_c_co2[tt] - v.perm_jul_ce_c_co2[tt-1])/p.perm_jul_limit_c)^(1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) < (1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) * yp / (p.perm_jul_lag_c_co2*v.perm_jul_lag_c_co2_correct[tt])) ? v.perm_jul_equilib_c_co2[tt] : (v.perm_jul_equilib_c_co2[tt] - p.perm_jul_limit_c * ((((v.perm_jul_equilib_c_co2[tt] - v.perm_jul_ce_c_co2[tt-1])/p.perm_jul_limit_c)^(1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) - (1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]) * yp / (p.perm_jul_lag_c_co2*v.perm_jul_lag_c_co2_correct[tt]))^(1/(1-(1+p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]))))))
        end
        # Released perm cumulative CO2
        v.perm_jul_ce_co2[tt] = v.perm_jul_ce_c_co2[tt] * (44/12)
        # Annual emissions CO2 perm
        v.perm_jul_e_co2[tt] = (44/12) * (p.perm_jul_limit_c / (p.perm_jul_lag_c_co2 * v.perm_jul_lag_c_co2_correct[tt])) * (ifelse(v.perm_jul_equilib_c_co2[tt] > v.perm_jul_ce_c_co2[tt], ((v.perm_jul_equilib_c_co2[tt] - v.perm_jul_ce_c_co2[tt])/p.perm_jul_limit_c)^((1 + p.perm_jul_pow_c_co2)*v.perm_jul_pow_c_co2_correct[tt]), 0))
        # Released perm cumulative CH4
        v.perm_jul_ce_ch4[tt] = v.perm_jul_ce_c_co2[tt] * (p.perm_jul_ch4_co2_c_ratio/100) * (16/12)
    end
end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
using PyPlot

function plot_stim(settings::RunSettings, layerNames::Vector{ASCIIString})
	dt = settings.dt
	T = settings.T
	time_trace = 0:dt:T
	nsteps = length(time_trace)
	fstim = settings.stim

	figure()
	for (l,lname) in enumerate(layerNames)
	    # show injections
	    stim_trace = [fstim(l,1,i) for i in 1:nsteps]
	    plot(time_trace,stim_trace)
	end
	legend(layerNames)
	ylim([-10, 80])
	show()
end

# raster is a units going down rows, true/false going across columns (time)
function plot_raster(spikes,N,nsteps)
	C = [t for (t,i) in spikes]
	R = [i for (t,i) in spikes]
	vlines(C,R-0.5,R+0.5) # vlines(x,ymin,ymax)
	xlim([0, nsteps])
	ylim([0, N])
end

# generic network plotting, from first layer to last
function plot_all_rasters(results::RunResults)
	T = results.T
	dt = results.dt
	nsteps = length(0:dt:T)
	nL = length(results.names)
	i=0
	for (i,layer_name) in results.names
		subplot(1,nL,i)
		N = results.N[i]
		spikes = results.spikes[i]
		plot_raster(spikes,N,nsteps)
		title(layer_name)
		i+=1
	end
end

# for plotting the vision architecture
function plot_vision_rasters(results::RunResults)
	T = results.T
	dt = results.dt
	nsteps = length(0:dt:T)
	D = (3,13)

	nL = length(results.names)
	for i=1:nL
		layer_name = results.names[i]
		println((i,layer_name))
		N = results.N[i]
		spikes = results.spikes[i]
		if i in 1:3
			subplot2grid(D,(i-1,0),colspan=2)
		elseif i in 4:6
			subplot2grid(D,(i%4,2),colspan=2)
		elseif i in 7:9
			subplot2grid(D,(0,(i-6)*3+1),rowspan=3, colspan=3) # V1_4
		else
			error("plot_vision_rasters only works with the 9-layer network")
		end
		plot_raster(spikes,N,nsteps)
		#xticks(0:dt:T) # in ms
		title(layer_name)
	end
end


text:
function calculateLikelihood(matrix, beta, gammas)
    llh = 0
    D, V = size(matrix) # number of documents

    for d = 1:D
        g = gammas[:, d]  / sum(gammas[:, d]) # normalize
        for j = 1:V
            if matrix[d, j] != 0
                llh += matrix[d, j] * log(dot(g, beta[:, j]))
            end
        end
    end

    return llh
end

function EStep(K, alpha, beta, docVector, maxIter=10)
    # doing E-step
    # K is number of topic
    # alpha, beta are parameters
    # docVector: term vector of doc 1 x N
    # docVector[i] = j means i-th word has j occurrence in corpus

    # init phi and gamma
    wIds = find(docVector) # word id of term in document
    N = length(wIds)
    phi = (ones(N, K) * 1 / K)'
    gammas = copy(alpha) + N / K # size of gamma K x 1

    # loop until convergence
    convergent = false
    threshold = 0.01 # parameter for checking convergence
    counter = 1

    while ~convergent
        oldPhi = copy(phi)
        oldGamma = copy(gammas)

        # update phi
        e_dig = exp(digamma(gammas))
        phi = beta[:, wIds] .* e_dig
        phi = phi ./ sum(phi, 1)

        # update gamma
        gammas = alpha + sum(phi, 2)

        # checking convergence
        m = max(sum(abs(phi - oldPhi)), sum(abs(gammas - oldGamma)))
        if m < threshold || counter == 10
            convergent = true
        end
        counter += 1
    end

    # return result
    return (phi, gammas)
end

function updateBeta(K, phi, matrix)
    D, N = size(matrix)
    beta = zeros(K, N)

    # update beta
    for d=1:D
        wIds = find(matrix[d, :])
        beta[:, wIds] += (diagm(vec(matrix[d, wIds])) * phi[d]')'
    end
    #normalize beta
    beta = broadcast( * , beta, 1./sum(beta, 1))
    return beta
end

function updateAlpha(K, gammas, maxIter=10)
    # K is number of topic
    # D is total number of documents
    # N is total number of words
    # phi is
    # gammas is K x D matrix
    # matrix is (

    #alpha = zeros(K, 1) + 1
    (K, D) = size(gammas)
    alpha = mean(gammas, 2) ./ K
    # update alpha
    convergent = false
    threshold = 0.01
    counter = 1
    pg = sum(digamma(gammas),2) - sum(digamma(sum(gammas,1)));

    while ~convergent
        #L = L_α(alpha, gammas)
        #println("MStep: ", counter, ":L:", L, ":alpha:", alpha)
        oldAlpha = copy(alpha)

        alpha0 = sum(alpha)
        g = D * (digamma(alpha0) - digamma(alpha)) + pg
        h = -1 ./ trigamma(alpha)
        hgz = h' * g / (1 / trigamma(alpha0) + sum(h))

        for i=1:K
            alpha[i] = alpha[i] - h[i] * ( g[i] - hgz[1, 1]) / D
        end

        # check convergence
        if sum(abs(alpha - oldAlpha)) < threshold || counter == maxIter
            convergent = true
        end
        counter +=1
    end

    return alpha
end

function L_α(α, γs)
    (K, D) = size(γs)
    L = D * (  log(gamma(sum(α))) - sum(log(gamma(α))))
    for d=1:D
        for i=1:K
            L += ( α[i] - 1 ) * ( digamma(γs[i, d]) - digamma(sum(γs[:, d])))
        end
    end
    return L
end

function doingEM(K, vocFile="../data/testvocab.txt", dataFile="../data/test.dat")
    # K is number of topic
    # D number of documents
    # N total number of words
    # alpha is parameter of Dirichlet distribution (size K x 1)
    # beta is matrix K x D
    corpus = readFile(dataFile)
    D = length(corpus) # number of doc in dataset
    println(D)
    vocab = readVoc(vocFile)
    N = length(vocab) # number of vocabulary in dataset
    println(N)
    matrix = convertToMatrix(corpus, vocab) # size D x N

    convergent = false
    # create phi and gamma to allocate memory only
    phi = Dict()
    gammas = zeros(K, D)
    # init alpha and beta
    alpha = rand(K, 1)
    alpha = alpha ./ sum(alpha) # normalize alpha
    b = rand(K, N)
    beta = broadcast(*, b, 1 ./ sum(b, 1)) # normalize

    counter = 1

    pllh = 0 # previous likelihood
    threshold = 0.01

    while ~convergent
        # loop until convergence

        # E-step
        #println("E step")
        for d = 1:D
            # doing for each document in corpus
            (subPhi, subGamma) = EStep(K, alpha, beta, matrix[d, :])
            phi[d] = subPhi
            gammas[:, d] = subGamma
        end

        # M-step
        #println("M Step")
        beta = updateBeta(K, phi, matrix)
        alpha = updateAlpha(K, gammas)

        # checking convergence
        llh = calculateLikelihood(matrix, beta, gammas)
        println("iter:", counter,"\tLLH:", llh)
        #if abs(pllh - llh) < threshold
        if counter == 10
            convergent = true
        #else
        #    pllh = llh
        #    println(llh)
        end
        counter += 1
    end
end

max_stars_repo_name: tndoan/SimpleLDA

max_issues_repo_name: tndoan/SimpleLDA

max_forks_repo_name: tndoan/SimpleLDA


text:
struct ConfusionMatrix
  tp::Int
  fp::Int
  fn::Int
  tn::Int
end

function getConfusionMatrix(Δ::AbstractMatrix{T}, hΔ::AbstractMatrix{T}, thr=1e-4) where T
  p = size(Δ, 1)

  tp = fp = fn = tn = 0
  for c=1:p
    for r=c:p
      if Δ[r, c] != zero(T)
        if abs(hΔ[r, c]) > thr
          # not zero estimated as not zero
          tp += 1
        else
          # not zero estimated as zero
          fn += 1
        end
      else
        if abs(hΔ[r, c]) > thr
          # zero estimated as not zero
          fp += 1
        else
          # zero estimated as zero
          tn += 1
        end
      end
    end
  end
  ConfusionMatrix(tp, fp, fn, tn)
end

function getConfusionMatrix(Δ::AbstractMatrix{T}, hΔ::SparseIterate{T}, thr=1e-4) where T
  p = size(Δ, 1)

  tp = fp = fn = tn = 0
  for c=1:p
    for r=c:p
      if Δ[r, c] != zero(T)
        if abs(hΔ[sub2indLowerTriangular(p, r, c)]) > thr
          # not zero estimated as not zero
          tp += 1
        else
          # not zero estimated as zero
          fn += 1
        end
      else
        if abs(hΔ[sub2indLowerTriangular(p, r, c)]) > thr
          # zero estimated as not zero
          fp += 1
        else
          # zero estimated as zero
          tn += 1
        end
      end
    end
  end
  ConfusionMatrix(tp, fp, fn, tn)
end

function tpr(t::ConfusionMatrix)
  dn = t.tp + t.fn
  dn > 0 ? t.tp / dn : 0.
end
function fpr(t::ConfusionMatrix)
  dn = t.fp + t.tn
  dn > 0 ? t.fp / dn : 0.
end
function precision(t::ConfusionMatrix)
  dn = t.tp + t.fp
  dn > 0 ? t.tp / dn : 1.
end
recall(t::ConfusionMatrix) = tpr(t)


text:
# manual fixes

const ptrdiff_t = Cptrdiff_t
const TIFF_STRIPCHOP = 8000
const FILE = 0 # This fixes loading, but breaks TIFFPrintDirectory(arg1, arg2, arg3)


text:
function to_rna(dna::AbstractString)

end


text:
struct Metapopulation{T} <: AbstractEpidemicModel where T
    h::AbstractSimpleGraph # Population network
    D::Vector{Float64}
    dynamics::T
end

# handy little wrapper
Metapopulation(h, D::Real, dynamics) = Metapopulation(
    h,
    fill(float(D),
    num_states(dynamics)),
    dynamics
    )

# Implement the AbstractEpidemicModel interface

function init_rates(mp::Metapopulation, state)
    N = nv(mp.h)
    a =  zeros(N) # Top level reaction rates per Population node
    init_rates!(a, state, mp) # pass it down to the dynamics
    return a
end

output_type(::Metapopulation) = Int

progressbar(::Metapopulation) = true

# Not part of the interface


text:
## Define the model
T = [
    -2 1 1;
    0.5 -1 0.5;
    1 1 -2;
    ]

C = [1; -1; 0]

r₁(x) = ones(size(x))
r₂(x) = -ones(size(x))

R₁(x) = x
R₂(x) = -x

r = (
    r = function (x)
        [r₁(x) r₂(x) r₂(x)]
    end,
    R = function (x)
        [R₁(x) R₂(x) R₂(x)]
    end,
)

## make models
simBounds = [0 Inf; -Inf Inf] # bounds for simulation only
simModel = SFFM.Model(T, C, r, Bounds = simBounds)

approxBounds = [0 48; -Inf Inf] # bounds for approximation only
model = SFFM.Model(T, C, r, Bounds = approxBounds)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


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

mutable struct CreationData <: SwaggerModel
    createOption::Any # spec type: Union{ Nothing, String } # spec name: createOption
    storageAccountId::Any # spec type: Union{ Nothing, String } # spec name: storageAccountId
    imageReference::Any # spec type: Union{ Nothing, ImageDiskReference } # spec name: imageReference
    sourceUri::Any # spec type: Union{ Nothing, String } # spec name: sourceUri
    sourceResourceId::Any # spec type: Union{ Nothing, String } # spec name: sourceResourceId

    function CreationData(;createOption=nothing, storageAccountId=nothing, imageReference=nothing, sourceUri=nothing, sourceResourceId=nothing)
        o = new()
        validate_property(CreationData, Symbol("createOption"), createOption)
        setfield!(o, Symbol("createOption"), createOption)
        validate_property(CreationData, Symbol("storageAccountId"), storageAccountId)
        setfield!(o, Symbol("storageAccountId"), storageAccountId)
        validate_property(CreationData, Symbol("imageReference"), imageReference)
        setfield!(o, Symbol("imageReference"), imageReference)
        validate_property(CreationData, Symbol("sourceUri"), sourceUri)
        setfield!(o, Symbol("sourceUri"), sourceUri)
        validate_property(CreationData, Symbol("sourceResourceId"), sourceResourceId)
        setfield!(o, Symbol("sourceResourceId"), sourceResourceId)
        o
    end
end # type CreationData

const _property_map_CreationData = Dict{Symbol,Symbol}(Symbol("createOption")=>Symbol("createOption"), Symbol("storageAccountId")=>Symbol("storageAccountId"), Symbol("imageReference")=>Symbol("imageReference"), Symbol("sourceUri")=>Symbol("sourceUri"), Symbol("sourceResourceId")=>Symbol("sourceResourceId"))
const _property_types_CreationData = Dict{Symbol,String}(Symbol("createOption")=>"String", Symbol("storageAccountId")=>"String", Symbol("imageReference")=>"ImageDiskReference", Symbol("sourceUri")=>"String", Symbol("sourceResourceId")=>"String")
Base.propertynames(::Type{ CreationData }) = collect(keys(_property_map_CreationData))
Swagger.property_type(::Type{ CreationData }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_CreationData[name]))}
Swagger.field_name(::Type{ CreationData }, property_name::Symbol) =  _property_map_CreationData[property_name]

const _allowed_CreationData_createOption = ["Empty", "Attach", "FromImage", "Import", "Copy"]

function check_required(o::CreationData)
    true
end

function validate_property(::Type{ CreationData }, name::Symbol, val)
    if name === Symbol("createOption")
        Swagger.validate_param(name, "CreationData", :enum, val, _allowed_CreationData_createOption)
    end
end


text:
const act_fn = relu

convx(ni, nf; ks = 3, stride = 1, ndim = 2) = Conv(
    ntuple(_ -> ks, ndim),
    ni => nf,
    stride = stride,
    pad = ks ÷ 2,
    init = Flux.kaiming_normal,
)

function convxlayer(ni, nf; ks = 3, stride = 1, zero_bn = false, act = true, ndim = 2)
    bn = BatchNorm(nf, act ? act_fn : identity)
    fill!(bn.γ, zero_bn ? 0 : 1)
    return Chain(convx(ni, nf; ks = ks, stride = stride, ndim = ndim), bn)
end

struct ResBlock
    convs::Any
    idconv::Any
    pool::Any
end

function ResBlock(expansion::Int, ni::Int, nh::Int; stride::Int = 1, ndim = 2)
    nf, ni = nh * expansion, ni * expansion
    if expansion == 1
        layers = [
            convxlayer(ni, nh; stride = stride, ndim = ndim),
            convxlayer(nh, nf; zero_bn = true, act = false, ndim = ndim),
        ]
    else
        layers = [
            convxlayer(ni, nh; ks = 1, ndim = ndim),
            convxlayer(nh, nh; stride = stride, ndim = ndim),
            convxlayer(nh, nf; ks = 1, zero_bn = true, act = false, ndim = ndim),
        ]
    end

    return ResBlock(
        Chain(layers...),
        ni == nf ? identity : convxlayer(ni, nf; ks = 1, stride = 1, ndim = ndim),
        stride == 1 ? identity : MeanPool(ntuple(_ -> 2, ndim)),
    )
end
Flux.@functor ResBlock

(r::ResBlock)(x) = act_fn.(r.convs(x) .+ r.idconv(r.pool(x)))

function make_layer(expansion, ni, nf, n_blocks, stride; ndim = 2)
    return Chain(
        [
            ResBlock(
                expansion,
                i == 1 ? ni : nf,
                nf;
                stride = (i == 1 ? stride : 1),
                ndim = ndim,
            ) for i = 1:n_blocks
        ]...,
    )
end

"""
    XResNet(expansion, layers; c_in = 3, ndim = 2)

Create an XResNet model backbone following the [implementation in
fastai](https://github.com/fastai/fastai/blob/master/fastai/vision/models/xresnet.py#L22).

- `c_in::Int = 3`: The number of input channels, e.g. `1` for grayscale images and `3` for
    RGB images
- `ndim::Int = 2`: The number of dimensions for the convolutional and pooling layers, e.g.
    `2` for 2D input images and `3` for 3D volumes.
"""
function XResNet(expansion, layers; c_in = 3, ndim = 2)
    nfs = [c_in, (c_in + 1) * 8, 64, 64]
    stem = [convxlayer(nfs[i], nfs[i+1]; stride = i == 1 ? 2 : 1, ndim = ndim) for i = 1:3]

    nfs = [64 ÷ expansion, 64, 128, 256, 512]
    res_layers = [
        make_layer(expansion, nfs[i], nfs[i+1], l, i == 1 ? 1 : 2, ndim = ndim) for
        (i, l) in enumerate(layers)
    ]

    return Chain(stem..., MaxPool(ntuple(_ -> 3, ndim); pad = 1, stride = 2), res_layers...)
end

xresnet18(; kwargs...) = XResNet(1, [2, 2, 2, 2]; kwargs...)
xresnet50(; kwargs...) = XResNet(4, [3, 4, 6, 3]; kwargs...)

@testset "XResNet [model]" begin
    @testset "Basic" begin
        @test_nowarn begin
            model = xresnet18()
            @test Flux.outputsize(model, (128, 128, 3, 1)) == (4, 4, 512, 1)
        end
    end

    @testset "3D" begin
        @test_nowarn begin
            model = xresnet18(ndim=3)
            @test Flux.outputsize(model, (128, 128, 128, 3, 1)) == (4, 4, 4, 512, 1)
        end
    end
end


text:
module ReconBuckets
using FFTW
using Gadgetron
using Gadgetron.MRD
using Gadgetron.Types:AcquisitionBucket,ReconData
using Gadgetron.Stream
using Gadgetron.Default
using PartialFunctions
using LinearAlgebra
using Setfield

fft_scaling(xsize,dims=[1]) = sqrt(reduce((*),map(d->xsize[d],dims)))

ortho_fft(x,dims=[1]) = fft(x,dims) ./ Float32(fft_scaling(size(x),dims))

ortho_ifft(x,dims=[1]) = bfft(x,dims) ./ Float32(fft_scaling(size(x),dims))

cfft(x, dims=[1]) = ifftshift(ortho_fft(fftshift(x,dims),dims),dims)

cifft(x, dims=[1]) = fftshift(ortho_ifft(ifftshift(x,dims),dims),dims)

function make_kspace(header::MRDHeader)

	matrix_size = header.encoding[1].encodedSpace.matrixSize

	function inner(bucket::AcquisitionBucket)
        acqs = bucket.data
		number_of_samples, number_of_channels = size(acqs[1].data)
		buffer = zeros(ComplexF32,number_of_samples,matrix_size.y,matrix_size.z,number_of_channels)
		for acq in acqs
			buffer[:,acq.header.idx.kspace_encode_step_1+1,acq.header.idx.kspace_encode_step_2+1,:] = acq.data
		end
		return [(acqs[1].header,buffer)]

	end

    function inner(buffer::ReconData)
        data= buffer[1].data.data
        dims = size(data )
        data_view = reshape(data,dims[1],dims[2],dims[3],dims[4],:)
        return (( buffer[1].data.header[1], data) for data in eachslice(data_view,dims=5))
    end

	return MapCat(inner)
end

function reconstruct_image(header::MRDHeader)

    field_of_view = header.encoding[1].reconSpace.fieldOfView_mm
    fov = (field_of_view.x,field_of_view.y,field_of_view.z)

	recon(x) = cifft(x,[1,2,3])
    coil_combine(x) = sqrt.(sum(abs2,x,dims=4))

	inner(index, reference::AcquisitionHeader, buffer) =   Image(ImageHeader(reference; image_index=index, image_type = MRD.ImageType.magnitude, field_of_view=fov),  coil_combine(recon(buffer)))
    inner_splat(args ) = inner(args[1],args[2][1],args[2][2])
	return Map(inner_splat)
end

function reconstruct_buckets(connection)
    start = time()
    header = connection.header
    connection |>
    make_kspace(header) |>
    Enumerate() |>
    reconstruct_image(header) .|>
    push! $ connection

    tot_time = time()-start
    @info "Total recon time $tot_time"
end

Base.precompile(Tuple{typeof(reconstruct_buckets),Gadgetron.MRDChannel})   # time: 0.5008968

end


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using PyPlot
using Elefridge
using JLD

path = "/Users/milan/cams/unstructured/"
filelist = filter(x->endswith(x,".grib"),readdir(path))

n = length(filelist)
varnames = fill("",n)

# example histogram for preallocation
Hex = fit(Histogram,rand(Float32,100),[0f0,0.5f0,1f0])
S = fill(Hex,n)

for (i,file) in enumerate(filelist)
    varname = split(split(file,"0_")[end],".")[1]
    varnames[i] = varname
    println("---")
    println("Reading $varname")
    local X = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    if all(X .== zero(eltype(X)))
        println("Only zeros in $varname")
    else
        mi = log(Elefridge.minpos(X))
        ma = log(maximum(X))

        local n = 1000
        Δ = (ma-mi)/n

        bins = vcat([-Inf32,0f0,floatmin(Float32)],exp.(mi:Δ:ma),[prevfloat(Inf32),Inf32])

        S[i] = fit(Histogram,vec(X),bins)
    end
end

@save "/Users/milan/cams/entropy/distribution_all_1k.jld" varnames S