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

const Transaction_Type = (;[
    Symbol("READ") => Int32(0),
    Symbol("WRITE") => Int32(1),
]...)

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

    function Transaction_Open_Req(; kwargs...)
        obj = new(meta(Transaction_Open_Req), 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 Transaction_Open_Req
const __meta_Transaction_Open_Req = Ref{ProtoMeta}()
function meta(::Type{Transaction_Open_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Open_Req)
            __meta_Transaction_Open_Req[] = target = ProtoMeta(Transaction_Open_Req)
            allflds = Pair{Symbol,Union{Type,String}}[:session_id => Vector{UInt8}, :_type => Int32, :options => Options, :network_latency_millis => Int32]
            meta(target, Transaction_Open_Req, 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_Transaction_Open_Req[]
    end
end
function Base.getproperty(obj::Transaction_Open_Req, name::Symbol)
    if name === :session_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Options
    elseif name === :network_latency_millis
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Open_Res(; kwargs...)
        obj = new(meta(Transaction_Open_Res), 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 Transaction_Open_Res
const __meta_Transaction_Open_Res = Ref{ProtoMeta}()
function meta(::Type{Transaction_Open_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Open_Res)
            __meta_Transaction_Open_Res[] = target = ProtoMeta(Transaction_Open_Res)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Open_Res, 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_Transaction_Open_Res[]
    end
end

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

    function Transaction_Open(; kwargs...)
        obj = new(meta(Transaction_Open), 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 Transaction_Open
const __meta_Transaction_Open = Ref{ProtoMeta}()
function meta(::Type{Transaction_Open})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Open)
            __meta_Transaction_Open[] = target = ProtoMeta(Transaction_Open)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Open, 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_Transaction_Open[]
    end
end

const Transaction_Stream_State = (;[
    Symbol("CONTINUE") => Int32(0),
    Symbol("DONE") => Int32(1),
]...)

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

    function Transaction_Stream_Req(; kwargs...)
        obj = new(meta(Transaction_Stream_Req), 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 Transaction_Stream_Req
const __meta_Transaction_Stream_Req = Ref{ProtoMeta}()
function meta(::Type{Transaction_Stream_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Stream_Req)
            __meta_Transaction_Stream_Req[] = target = ProtoMeta(Transaction_Stream_Req)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Stream_Req, 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_Transaction_Stream_Req[]
    end
end

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

    function Transaction_Stream_ResPart(; kwargs...)
        obj = new(meta(Transaction_Stream_ResPart), 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 Transaction_Stream_ResPart
const __meta_Transaction_Stream_ResPart = Ref{ProtoMeta}()
function meta(::Type{Transaction_Stream_ResPart})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Stream_ResPart)
            __meta_Transaction_Stream_ResPart[] = target = ProtoMeta(Transaction_Stream_ResPart)
            allflds = Pair{Symbol,Union{Type,String}}[:state => Int32]
            meta(target, Transaction_Stream_ResPart, 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_Transaction_Stream_ResPart[]
    end
end
function Base.getproperty(obj::Transaction_Stream_ResPart, name::Symbol)
    if name === :state
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

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

    function Transaction_ResPart(; kwargs...)
        obj = new(meta(Transaction_ResPart), 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 Transaction_ResPart
const __meta_Transaction_ResPart = Ref{ProtoMeta}()
function meta(::Type{Transaction_ResPart})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_ResPart)
            __meta_Transaction_ResPart[] = target = ProtoMeta(Transaction_ResPart)
            allflds = Pair{Symbol,Union{Type,String}}[:req_id => Vector{UInt8}, :stream_res_part => Transaction_Stream_ResPart, :query_manager_res_part => QueryManager_ResPart, :logic_manager_res_part => LogicManager_ResPart, :type_res_part => Type_ResPart, :thing_res_part => Thing_ResPart]
            oneofs = Int[0,1,1,1,1,1]
            oneof_names = Symbol[Symbol("res")]
            meta(target, Transaction_ResPart, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Transaction_ResPart[]
    end
end
function Base.getproperty(obj::Transaction_ResPart, name::Symbol)
    if name === :req_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :stream_res_part
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Stream_ResPart
    elseif name === :query_manager_res_part
        return (obj.__protobuf_jl_internal_values[name])::QueryManager_ResPart
    elseif name === :logic_manager_res_part
        return (obj.__protobuf_jl_internal_values[name])::LogicManager_ResPart
    elseif name === :type_res_part
        return (obj.__protobuf_jl_internal_values[name])::Type_ResPart
    elseif name === :thing_res_part
        return (obj.__protobuf_jl_internal_values[name])::Thing_ResPart
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Stream(; kwargs...)
        obj = new(meta(Transaction_Stream), 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 Transaction_Stream
const __meta_Transaction_Stream = Ref{ProtoMeta}()
function meta(::Type{Transaction_Stream})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Stream)
            __meta_Transaction_Stream[] = target = ProtoMeta(Transaction_Stream)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Stream, 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_Transaction_Stream[]
    end
end

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

    function Transaction_Commit_Req(; kwargs...)
        obj = new(meta(Transaction_Commit_Req), 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 Transaction_Commit_Req
const __meta_Transaction_Commit_Req = Ref{ProtoMeta}()
function meta(::Type{Transaction_Commit_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Commit_Req)
            __meta_Transaction_Commit_Req[] = target = ProtoMeta(Transaction_Commit_Req)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Commit_Req, 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_Transaction_Commit_Req[]
    end
end

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

    function Transaction_Commit_Res(; kwargs...)
        obj = new(meta(Transaction_Commit_Res), 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 Transaction_Commit_Res
const __meta_Transaction_Commit_Res = Ref{ProtoMeta}()
function meta(::Type{Transaction_Commit_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Commit_Res)
            __meta_Transaction_Commit_Res[] = target = ProtoMeta(Transaction_Commit_Res)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Commit_Res, 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_Transaction_Commit_Res[]
    end
end

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

    function Transaction_Commit(; kwargs...)
        obj = new(meta(Transaction_Commit), 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 Transaction_Commit
const __meta_Transaction_Commit = Ref{ProtoMeta}()
function meta(::Type{Transaction_Commit})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Commit)
            __meta_Transaction_Commit[] = target = ProtoMeta(Transaction_Commit)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Commit, 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_Transaction_Commit[]
    end
end

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

    function Transaction_Rollback_Req(; kwargs...)
        obj = new(meta(Transaction_Rollback_Req), 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 Transaction_Rollback_Req
const __meta_Transaction_Rollback_Req = Ref{ProtoMeta}()
function meta(::Type{Transaction_Rollback_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Rollback_Req)
            __meta_Transaction_Rollback_Req[] = target = ProtoMeta(Transaction_Rollback_Req)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Rollback_Req, 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_Transaction_Rollback_Req[]
    end
end

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

    function Transaction_Req_MetadataEntry(; kwargs...)
        obj = new(meta(Transaction_Req_MetadataEntry), 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 Transaction_Req_MetadataEntry (mapentry)
const __meta_Transaction_Req_MetadataEntry = Ref{ProtoMeta}()
function meta(::Type{Transaction_Req_MetadataEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Req_MetadataEntry)
            __meta_Transaction_Req_MetadataEntry[] = target = ProtoMeta(Transaction_Req_MetadataEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => AbstractString]
            meta(target, Transaction_Req_MetadataEntry, 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_Transaction_Req_MetadataEntry[]
    end
end
function Base.getproperty(obj::Transaction_Req_MetadataEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Req(; kwargs...)
        obj = new(meta(Transaction_Req), 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 Transaction_Req
const __meta_Transaction_Req = Ref{ProtoMeta}()
function meta(::Type{Transaction_Req})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Req)
            __meta_Transaction_Req[] = target = ProtoMeta(Transaction_Req)
            allflds = Pair{Symbol,Union{Type,String}}[:req_id => Vector{UInt8}, :metadata => Base.Dict{AbstractString,AbstractString}, :open_req => Transaction_Open_Req, :stream_req => Transaction_Stream_Req, :commit_req => Transaction_Commit_Req, :rollback_req => Transaction_Rollback_Req, :query_manager_req => QueryManager_Req, :concept_manager_req => ConceptManager_Req, :logic_manager_req => LogicManager_Req, :rule_req => Rule_Req, :type_req => Type_Req, :thing_req => Thing_Req]
            oneofs = Int[0,0,1,1,1,1,1,1,1,1,1,1]
            oneof_names = Symbol[Symbol("req")]
            meta(target, Transaction_Req, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Transaction_Req[]
    end
end
function Base.getproperty(obj::Transaction_Req, name::Symbol)
    if name === :req_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :metadata
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,AbstractString}
    elseif name === :open_req
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Open_Req
    elseif name === :stream_req
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Stream_Req
    elseif name === :commit_req
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Commit_Req
    elseif name === :rollback_req
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Rollback_Req
    elseif name === :query_manager_req
        return (obj.__protobuf_jl_internal_values[name])::QueryManager_Req
    elseif name === :concept_manager_req
        return (obj.__protobuf_jl_internal_values[name])::ConceptManager_Req
    elseif name === :logic_manager_req
        return (obj.__protobuf_jl_internal_values[name])::LogicManager_Req
    elseif name === :rule_req
        return (obj.__protobuf_jl_internal_values[name])::Rule_Req
    elseif name === :type_req
        return (obj.__protobuf_jl_internal_values[name])::Type_Req
    elseif name === :thing_req
        return (obj.__protobuf_jl_internal_values[name])::Thing_Req
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Client(; kwargs...)
        obj = new(meta(Transaction_Client), 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 Transaction_Client
const __meta_Transaction_Client = Ref{ProtoMeta}()
function meta(::Type{Transaction_Client})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Client)
            __meta_Transaction_Client[] = target = ProtoMeta(Transaction_Client)
            allflds = Pair{Symbol,Union{Type,String}}[:reqs => Base.Vector{Transaction_Req}]
            meta(target, Transaction_Client, 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_Transaction_Client[]
    end
end
function Base.getproperty(obj::Transaction_Client, name::Symbol)
    if name === :reqs
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Transaction_Req}
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Rollback_Res(; kwargs...)
        obj = new(meta(Transaction_Rollback_Res), 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 Transaction_Rollback_Res
const __meta_Transaction_Rollback_Res = Ref{ProtoMeta}()
function meta(::Type{Transaction_Rollback_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Rollback_Res)
            __meta_Transaction_Rollback_Res[] = target = ProtoMeta(Transaction_Rollback_Res)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Rollback_Res, 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_Transaction_Rollback_Res[]
    end
end

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

    function Transaction_Res(; kwargs...)
        obj = new(meta(Transaction_Res), 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 Transaction_Res
const __meta_Transaction_Res = Ref{ProtoMeta}()
function meta(::Type{Transaction_Res})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Res)
            __meta_Transaction_Res[] = target = ProtoMeta(Transaction_Res)
            allflds = Pair{Symbol,Union{Type,String}}[:req_id => Vector{UInt8}, :open_res => Transaction_Open_Res, :commit_res => Transaction_Commit_Res, :rollback_res => Transaction_Rollback_Res, :query_manager_res => QueryManager_Res, :concept_manager_res => ConceptManager_Res, :logic_manager_res => LogicManager_Res, :rule_res => Rule_Res, :type_res => Type_Res, :thing_res => Thing_Res]
            oneofs = Int[0,1,1,1,1,1,1,1,1,1]
            oneof_names = Symbol[Symbol("res")]
            meta(target, Transaction_Res, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Transaction_Res[]
    end
end
function Base.getproperty(obj::Transaction_Res, name::Symbol)
    if name === :req_id
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :open_res
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Open_Res
    elseif name === :commit_res
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Commit_Res
    elseif name === :rollback_res
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Rollback_Res
    elseif name === :query_manager_res
        return (obj.__protobuf_jl_internal_values[name])::QueryManager_Res
    elseif name === :concept_manager_res
        return (obj.__protobuf_jl_internal_values[name])::ConceptManager_Res
    elseif name === :logic_manager_res
        return (obj.__protobuf_jl_internal_values[name])::LogicManager_Res
    elseif name === :rule_res
        return (obj.__protobuf_jl_internal_values[name])::Rule_Res
    elseif name === :type_res
        return (obj.__protobuf_jl_internal_values[name])::Type_Res
    elseif name === :thing_res
        return (obj.__protobuf_jl_internal_values[name])::Thing_Res
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Server(; kwargs...)
        obj = new(meta(Transaction_Server), 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 Transaction_Server
const __meta_Transaction_Server = Ref{ProtoMeta}()
function meta(::Type{Transaction_Server})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Server)
            __meta_Transaction_Server[] = target = ProtoMeta(Transaction_Server)
            fnum = Int[2,3]
            allflds = Pair{Symbol,Union{Type,String}}[:res => Transaction_Res, :res_part => Transaction_ResPart]
            oneofs = Int[1,1]
            oneof_names = Symbol[Symbol("server")]
            meta(target, Transaction_Server, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Transaction_Server[]
    end
end
function Base.getproperty(obj::Transaction_Server, name::Symbol)
    if name === :res
        return (obj.__protobuf_jl_internal_values[name])::Transaction_Res
    elseif name === :res_part
        return (obj.__protobuf_jl_internal_values[name])::Transaction_ResPart
    else
        getfield(obj, name)
    end
end

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

    function Transaction_Rollback(; kwargs...)
        obj = new(meta(Transaction_Rollback), 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 Transaction_Rollback
const __meta_Transaction_Rollback = Ref{ProtoMeta}()
function meta(::Type{Transaction_Rollback})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction_Rollback)
            __meta_Transaction_Rollback[] = target = ProtoMeta(Transaction_Rollback)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction_Rollback, 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_Transaction_Rollback[]
    end
end

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

    function Transaction(; kwargs...)
        obj = new(meta(Transaction), 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 Transaction
const __meta_Transaction = Ref{ProtoMeta}()
function meta(::Type{Transaction})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Transaction)
            __meta_Transaction[] = target = ProtoMeta(Transaction)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(target, Transaction, 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_Transaction[]
    end
end

export Transaction_Type, Transaction_Client, Transaction_Server, Transaction_Req_MetadataEntry, Transaction_Req, Transaction_Res, Transaction_ResPart, Transaction_Open_Req, Transaction_Open_Res, Transaction_Open, Transaction_Stream_State, Transaction_Stream_Req, Transaction_Stream_ResPart, Transaction_Stream, Transaction_Commit_Req, Transaction_Commit_Res, Transaction_Commit, Transaction_Rollback_Req, Transaction_Rollback_Res, Transaction_Rollback, Transaction
# mapentries: "Transaction_Req_MetadataEntry" => ("AbstractString", "AbstractString")


text:
"""

stan_sample()

Draw from a StanJulia SampleModel (<: CmdStanModel.)

## Required argument
```julia
* `m <: CmdStanModels`                 # SampleModel
```

### Most frequently used keyword arguments
```julia
* `data`                               # Observations Dict or NamedTuple.
* `init`                               # Init Dict or NT (default: -2 to +2).
```

See extended help for other keyword arguments ( `??stan_sample` ).

### Returns
```julia
* `rc`                                 # Return code, 0 is success.
```

# Extended help

### Additional configuration keyword arguments
```julia
* `num_threads=4`                      # Update number of c++ threads.
* `check_num_chains=true`              # Check for C++ chains or Julia level chains
* `num_cpp_chains=1`                   # Update number of c++ chains.
* `num_julia_chains=1`                 # Update number of Julia chains.
                                       # Both initialized from num_chains

* `use_cpp_chains=false`               # Run num_chains on c++ level
                                       # Set to false to use Julia processes

* `num_chains=4`                       # Actual number of chains.
* `num_samples=1000`                   # Number of samples.
* `num_warmups=1000`                   # Number of warmup samples.
* `save_warmup=false`                  # Save warmup samples.

* `thin=1`                             # Set thinning value.
* `seed=-1`                            # Set seed value.

* `engaged=true`                       # Adaptation engaged.
* `gamma=0.05`                         # Adaptation regularization scale.
* `delta=0.8`                          # Adaptation target acceptance statistic.
* `kappa=0.75`                         # Adaptation relaxation exponent.
* `t0=10`                              # Adaptation iteration offset.
* `init_buffer=75`                     # Inital adaptation interval.
* `term_buffer=50`                     # Final fast adaptation interval.
* `window=25`                          # Initia; slow adaptation interval.

* `algorithm=:hmc`                     # Sampling algorithm.
* `engine=:nuts`                       # :nuts or :static.
* `max_depth=10`                       # Max tree depth for :nuts engine.
* `int_time=2 * pi`                    # Integration time for :static engine.

* `metric=:diag_e`                     # Geometry of manifold setting:
                                       # :diag_e, :unit_e or :dense_e.
* `metric_file=""`                     # Precompiled Euclidean metric.
* `stepsize=1.0`                       # Step size for discrete evolution
* `stepsize_jitter=0.0`                # Random jitter on step size ( [%] )

* `use_json=true`                      # Set to false for .R data files.

* `summary=true`                       # Create stansummary .csv file
* `print_summary=false`                # Display summary
```

Note: Currently I do not suggest to use both C++ level chains and Julia
level chains. By default, based on  `use_cpp_chains` the `stan_sample()`
method will set either `num_cpp_chains=num_chains; num_julia_chains=1`
(the default) or `num_julia_chains=num_chains;num_cpp_chain=1`. Set the
`check_num_chains` keyword argument in the call to `stan_sample()` to
`false` to prevent this default behavior.

Threads on C++ level can be used in multiple ways, e.g. to run separate
chains and to speed up certain operations. By default StanSample.jl's
SampleModel sets the C++ num_threads to 4. See the `graphs` subdirectory
in the RedCardsStudy in the Examples directory for an example.

"""
function stan_run(m::T; kwargs...) where {T <: CmdStanModels}

    handle_keywords!(m, kwargs)

    # Diagnostics files requested?
    diagnostics = false
    if :diagnostics in keys(kwargs)
        diagnostics = kwargs[:diagnostics]
        setup_diagnostics(m, m.num_julia_chains)
    end

    # Remove existing sample files
    if m.use_cpp_chains
            sfile = m.output_base  * "_chain.csv"
            isfile(sfile) && rm(sfile)
    else
        for id in 1:m.num_julia_chains
            sfile = sample_file_path(m.output_base, id)
            isfile(sfile) && rm(sfile)
        end
    end

    # Create cmdstan data & init input files (.json or .R)
    m.data_file = String[]
    m.init_file = String[]
    if m.use_json
        :init in keys(kwargs) && update_json_files(m, kwargs[:init],
            m.num_julia_chains, "init")
        :data in keys(kwargs) && update_json_files(m, kwargs[:data],
            m.num_julia_chains, "data")
    else
        :init in keys(kwargs) && update_R_files(m, kwargs[:init],
            m.num_julia_chains, "init")
        :data in keys(kwargs) && update_R_files(m, kwargs[:data],
            m.num_julia_chains, "data")
    end

    m.sample_file = String[]
    m.log_file = String[]
    m.diagnostic_file = String[]
    m.cmds = [stan_cmds(m, id; kwargs...) for id in 1:m.num_julia_chains]

    #println(m.cmds)

    run(pipeline(par(m.cmds), stdout=m.log_file[1]))
 end

"""

Generate a cmdstan command line (a run `cmd`).

$(SIGNATURES)

Internal, not exported.
"""
function stan_cmds(m::T, id::Int; kwargs...) where {T <: CmdStanModels}
    if m.use_cpp_chains && m.check_num_chains
        append!(m.sample_file, [m.output_base  * "_chain.csv"])
    else
        for id in 1:m.num_julia_chains
            append!(m.sample_file, [sample_file_path(m.output_base, id)])
        end
    end
    append!(m.log_file, [log_file_path(m.output_base, id)])
    if length(m.diagnostic_file) > 0
      append!(m.diagnostic_file, [diagnostic_file_path(m.output_base, id)])
    end
    cmdline(m, id)
end

ext: jl

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text: using eFEM, Test


text:
using ..Atmos
using ..Atmos: AbstractMoistureModel, PrecipitationModel, recover_thermo_state
using ..Mesh.Topologies
using ..Mesh.Grids
using Thermodynamics
using ..TurbulenceClosures
using LinearAlgebra

"""
    setup_atmos_default_diagnostics(
        ::AtmosLESConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create the "AtmosLESDefault" `DiagnosticsGroup` which contains the following
diagnostic variables, all of which are density-averaged horizontal averages,
variances and co-variances:

- u: x-velocity
- v: y-velocity
- w: z-velocity
- avg_rho: air density (_not_ density-averaged)
- rho: air density
- temp: air temperature
- pres: air pressure
- thd: dry potential temperature
- et: total specific energy
- ei: specific internal energy
- ht: specific enthalpy based on total energy
- hi: specific enthalpy based on internal energy
- w_ht_sgs: vertical sgs flux of total specific enthalpy
- var_u: variance of x-velocity
- var_v: variance of y-velocity
- var_w: variance of z-velocity
- w3: third moment of z-velocity
- tke: turbulent kinetic energy
- var_ei: variance of specific internal energy
- cov_w_u: vertical eddy flux of x-velocity
- cov_w_v: vertical eddy flux of y-velocity
- cov_w_rho: vertical eddy flux of density
- cov_w_thd: vertical eddy flux of dry potential temperature
- cov_w_ei: vertical eddy flux of specific internal energy

When an `EquilMoist` or a `NonEquilMoist` moisture model is used, the following
diagnostic variables are also output, also density-averaged horizontal averages,
variances and co-variances:

- qt: mass fraction of total water in air
- ql: mass fraction of liquid water in air
- qv: mass fraction of water vapor in air
- qi: mass fraction of ice in air
- thv: virtual potential temperature
- thl: liquid-ice potential temperature
- w_qt_sgs: vertical sgs flux of total specific humidity
- var_qt: variance of total specific humidity
- var_thl: variance of liquid-ice potential temperature
- cov_w_qt: vertical eddy flux of total specific humidity
- cov_w_ql: vertical eddy flux of liquid water specific humidity
- cov_w_qi: vertical eddy flux of cloud ice specific humidity
- cov_w_qv: vertical eddy flux of water vapor specific humidity
- cov_w_thv: vertical eddy flux of virtual potential temperature
- cov_w_thl: vertical eddy flux of liquid-ice potential temperature
- cov_qt_thl: covariance of total specific humidity and liquid-ice potential temperature
- cov_qt_ei: covariance of total specific humidity and specific internal energy

All these variables are output with the `z` dimension (`x3id`) on the DG grid
(`interpol` may _not_ be specified) as well as a (unlimited) `time` dimension
at the specified `interval`.
"""
function setup_atmos_default_diagnostics(
    ::AtmosLESConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosLESDefault",
        Diagnostics.atmos_les_default_init,
        Diagnostics.atmos_les_default_fini,
        Diagnostics.atmos_les_default_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

# Simple horizontal averages
function vars_atmos_les_default_simple(m::AtmosModel, FT)
    @vars begin
        u::FT
        v::FT
        w::FT
        avg_rho::FT             # ρ
        rho::FT                 # ρρ
        temp::FT
        pres::FT
        thd::FT                 # θ_dry
        et::FT                  # e_tot
        ei::FT                  # e_int
        ht::FT
        hi::FT
        w_ht_sgs::FT

        moisture::vars_atmos_les_default_simple(moisture_model(m), FT)
        precipitation::vars_atmos_les_default_simple(precipitation_model(m), FT)
    end
end
vars_atmos_les_default_simple(::AbstractMoistureModel, FT) = @vars()
function vars_atmos_les_default_simple(m::Union{EquilMoist, NonEquilMoist}, FT)
    @vars begin
        qt::FT                  # q_tot
        ql::FT                  # q_liq
        qi::FT                  # q_ice
        qv::FT                  # q_vap
        thv::FT                 # θ_vir
        thl::FT                 # θ_liq
        w_qt_sgs::FT
    end
end
vars_atmos_les_default_simple(::PrecipitationModel, FT) = @vars()
function vars_atmos_les_default_simple(::RainModel, FT)
    @vars begin
        qr::FT                  # q_rai
    end
end
function vars_atmos_les_default_simple(::RainSnowModel, FT)
    @vars begin
        qr::FT                  # q_rai
        qs::FT                  # q_sno
    end
end
num_atmos_les_default_simple_vars(m, FT) =
    varsize(vars_atmos_les_default_simple(m, FT))
atmos_les_default_simple_vars(m, array) =
    Vars{vars_atmos_les_default_simple(m, eltype(array))}(array)

function atmos_les_default_simple_sums!(
    atmos::AtmosModel,
    state,
    gradflux,
    aux,
    thermo,
    currtime,
    MH,
    sums,
)
    sums.u += MH * state.ρu[1]
    sums.v += MH * state.ρu[2]
    sums.w += MH * state.ρu[3]
    sums.avg_rho += MH * state.ρ
    sums.rho += MH * state.ρ * state.ρ
    sums.temp += MH * thermo.temp * state.ρ
    sums.pres += MH * thermo.pres * state.ρ
    sums.thd += MH * thermo.θ_dry * state.ρ
    sums.et += MH * state.energy.ρe
    sums.ei += MH * thermo.e_int * state.ρ
    sums.ht += MH * thermo.h_tot * state.ρ
    sums.hi += MH * thermo.h_int * state.ρ

    ν, D_t, _ = turbulence_tensors(atmos, state, gradflux, aux, currtime)
    d_h_tot = -D_t .* gradflux.energy.∇h_tot
    sums.w_ht_sgs += MH * d_h_tot[end] * state.ρ

    atmos_les_default_simple_sums!(
        moisture_model(atmos),
        state,
        gradflux,
        thermo,
        MH,
        D_t,
        sums,
    )
    atmos_les_default_simple_sums!(
        precipitation_model(atmos),
        state,
        gradflux,
        thermo,
        MH,
        D_t,
        sums,
    )
    return nothing
end
function atmos_les_default_simple_sums!(
    ::AbstractMoistureModel,
    state,
    gradflux,
    thermo,
    MH,
    D_t,
    sums,
)
    return nothing
end
function atmos_les_default_simple_sums!(
    moist::Union{EquilMoist, NonEquilMoist},
    state,
    gradflux,
    thermo,
    MH,
    D_t,
    sums,
)
    sums.moisture.qt += MH * state.moisture.ρq_tot
    sums.moisture.ql += MH * thermo.moisture.q_liq * state.ρ
    sums.moisture.qi += MH * thermo.moisture.q_ice * state.ρ
    sums.moisture.qv += MH * thermo.moisture.q_vap * state.ρ
    sums.moisture.thv += MH * thermo.moisture.θ_vir * state.ρ
    sums.moisture.thl += MH * thermo.moisture.θ_liq_ice * state.ρ
    d_q_tot = (-D_t) .* gradflux.moisture.∇q_tot
    sums.moisture.w_qt_sgs += MH * d_q_tot[end] * state.ρ

    return nothing
end
function atmos_les_default_simple_sums!(
    ::PrecipitationModel,
    state,
    gradflux,
    thermo,
    MH,
    D_t,
    sums,
)
    return nothing
end
function atmos_les_default_simple_sums!(
    precipitation::RainModel,
    state,
    gradflux,
    thermo,
    MH,
    D_t,
    sums,
)
    sums.precipitation.qr += MH * state.precipitation.ρq_rai

    return nothing
end
function atmos_les_default_simple_sums!(
    precipitation::RainSnowModel,
    state,
    gradflux,
    thermo,
    MH,
    D_t,
    sums,
)
    sums.precipitation.qr += MH * state.precipitation.ρq_rai
    sums.precipitation.qs += MH * state.precipitation.ρq_sno

    return nothing
end

function atmos_les_default_clouds(
    ::AbstractMoistureModel,
    thermo,
    idx,
    qc_gt_0_z,
    qc_gt_0_full,
    z,
    cld_top,
    cld_base,
)
    return cld_top, cld_base
end
function atmos_les_default_clouds(
    moist::Union{EquilMoist, NonEquilMoist},
    thermo,
    idx,
    qc_gt_0_z,
    qc_gt_0_full,
    z,
    cld_top,
    cld_base,
)
    if thermo.moisture.has_condensate
        FT = eltype(qc_gt_0_z)
        qc_gt_0_z[idx] = one(FT)
        qc_gt_0_full[idx] = one(FT)

        cld_top = max(cld_top, z)
        cld_base = min(cld_base, z)
    end
    return cld_top, cld_base
end

# Variances and covariances
function vars_atmos_les_default_ho(m::AtmosModel, FT)
    @vars begin
        var_u::FT               # u′u′
        var_v::FT               # v′v′
        var_w::FT               # w′w′
        w3::FT                  # w′w′w′
        tke::FT
        var_ei::FT              # e_int′e_int′

        cov_w_u::FT             # w′u′
        cov_w_v::FT             # w′v′
        cov_w_rho::FT           # w′ρ′
        cov_w_thd::FT           # w′θ_dry′
        cov_w_ei::FT            # w′e_int′

        moisture::vars_atmos_les_default_ho(moisture_model(m), FT)
        precipitation::vars_atmos_les_default_ho(precipitation_model(m), FT)
    end
end
vars_atmos_les_default_ho(::AbstractMoistureModel, FT) = @vars()
function vars_atmos_les_default_ho(m::Union{EquilMoist, NonEquilMoist}, FT)
    @vars begin
        var_qt::FT              # q_tot′q_tot′
        var_thl::FT             # θ_liq_ice′θ_liq_ice′

        cov_w_qt::FT            # w′q_tot′
        cov_w_ql::FT            # w′q_liq′
        cov_w_qi::FT            # w′q_ice′
        cov_w_qv::FT            # w′q_vap′
        cov_w_thv::FT           # w′θ_v′
        cov_w_thl::FT           # w′θ_liq_ice′
        cov_qt_thl::FT          # q_tot′θ_liq_ice′
        cov_qt_ei::FT           # q_tot′e_int′
    end
end
vars_atmos_les_default_ho(::PrecipitationModel, FT) = @vars()
function vars_atmos_les_default_ho(m::RainModel, FT)
    @vars begin
        var_qr::FT              # q_rai′q_rai′
        cov_w_qr::FT            # w′q_rai′
    end
end
function vars_atmos_les_default_ho(m::RainSnowModel, FT)
    @vars begin
        var_qr::FT              # q_rai′q_rai′
        var_qs::FT              # q_sno′q_sno′
        cov_w_qr::FT            # w′q_rai′
        cov_w_qs::FT            # w′q_sno′
    end
end
num_atmos_les_default_ho_vars(m, FT) = varsize(vars_atmos_les_default_ho(m, FT))
atmos_les_default_ho_vars(m, array) =
    Vars{vars_atmos_les_default_ho(m, eltype(array))}(array)

function atmos_les_default_ho_sums!(
    atmos::AtmosModel,
    state,
    thermo,
    MH,
    ha,
    sums,
)
    u = state.ρu[1] / state.ρ
    u′ = u - ha.u
    v = state.ρu[2] / state.ρ
    v′ = v - ha.v
    w = state.ρu[3] / state.ρ
    w′ = w - ha.w
    e_int′ = thermo.e_int - ha.ei
    θ_dry′ = thermo.θ_dry - ha.thd

    sums.var_u += MH * u′^2 * state.ρ
    sums.var_v += MH * v′^2 * state.ρ
    sums.var_w += MH * w′^2 * state.ρ
    sums.w3 += MH * w′^3 * state.ρ
    sums.tke +=
        0.5 * (MH * u′^2 * state.ρ + MH * v′^2 * state.ρ + MH * w′^2 * state.ρ)
    sums.var_ei += MH * e_int′^2 * state.ρ

    sums.cov_w_u += MH * w′ * u′ * state.ρ
    sums.cov_w_v += MH * w′ * v′ * state.ρ
    sums.cov_w_rho += MH * w′ * (state.ρ - ha.avg_rho) * state.ρ
    sums.cov_w_thd += MH * w′ * θ_dry′ * state.ρ
    sums.cov_w_ei += MH * w′ * e_int′ * state.ρ

    atmos_les_default_ho_sums!(
        moisture_model(atmos),
        state,
        thermo,
        MH,
        ha,
        w′,
        e_int′,
        sums,
    )
    atmos_les_default_ho_sums!(
        precipitation_model(atmos),
        state,
        thermo,
        MH,
        ha,
        w′,
        e_int′,
        sums,
    )
    return nothing
end
function atmos_les_default_ho_sums!(
    ::AbstractMoistureModel,
    state,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    return nothing
end
function atmos_les_default_ho_sums!(
    moist::Union{EquilMoist, NonEquilMoist},
    state,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    q_tot = state.moisture.ρq_tot / state.ρ
    q_tot′ = q_tot - ha.moisture.qt
    q_liq′ = thermo.moisture.q_liq - ha.moisture.ql
    q_ice′ = thermo.moisture.q_ice - ha.moisture.qi
    q_vap′ = thermo.moisture.q_vap - ha.moisture.qv
    θ_vir′ = thermo.moisture.θ_vir - ha.moisture.thv
    θ_liq_ice′ = thermo.moisture.θ_liq_ice - ha.moisture.thl

    sums.moisture.var_qt += MH * q_tot′^2 * state.ρ
    sums.moisture.var_thl += MH * θ_liq_ice′^2 * state.ρ

    sums.moisture.cov_w_qt += MH * w′ * q_tot′ * state.ρ
    sums.moisture.cov_w_ql += MH * w′ * q_liq′ * state.ρ
    sums.moisture.cov_w_qi += MH * w′ * q_ice′ * state.ρ
    sums.moisture.cov_w_qv += MH * w′ * q_vap′ * state.ρ
    sums.moisture.cov_w_thv += MH * w′ * θ_vir′ * state.ρ
    sums.moisture.cov_w_thl += MH * w′ * θ_liq_ice′ * state.ρ
    sums.moisture.cov_qt_thl += MH * q_tot′ * θ_liq_ice′ * state.ρ
    sums.moisture.cov_qt_ei += MH * q_tot′ * e_int′ * state.ρ

    return nothing
end
function atmos_les_default_ho_sums!(
    ::PrecipitationModel,
    state,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    return nothing
end
function atmos_les_default_ho_sums!(
    moist::RainModel,
    state,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    q_rai = state.precipitation.ρq_rai / state.ρ
    q_rai′ = q_rai - ha.precipitation.qr

    sums.precipitation.var_qr += MH * q_rai′^2 * state.ρ

    sums.precipitation.cov_w_qr += MH * w′ * q_rai′ * state.ρ

    return nothing
end
function atmos_les_default_ho_sums!(
    moist::RainSnowModel,
    state,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    q_rai = state.precipitation.ρq_rai / state.ρ
    q_rai′ = q_rai - ha.precipitation.qr
    q_sno = state.precipitation.ρq_sno / state.ρ
    q_sno′ = q_sno - ha.precipitation.qs

    sums.precipitation.var_qr += MH * q_rai′^2 * state.ρ
    sums.precipitation.var_qs += MH * q_sno′^2 * state.ρ

    sums.precipitation.cov_w_qr += MH * w′ * q_rai′ * state.ρ
    sums.precipitation.cov_w_qs += MH * w′ * q_sno′ * state.ρ

    return nothing
end

function prefix_filter(s)
    if startswith(s, "moisture.")
        s[10:end]
    elseif startswith(s, "precipitation.")
        s[15:end]
    else
        s
    end
end

"""
    atmos_les_default_init(dgngrp, currtime)

Initialize the 'AtmosLESDefault' diagnostics group.
"""
function atmos_les_default_init(dgngrp::DiagnosticsGroup, currtime)
    atmos = Settings.dg.balance_law
    FT = eltype(Settings.Q)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    atmos_collect_onetime(mpicomm, Settings.dg, Settings.Q)

    if mpirank == 0
        dims = OrderedDict("z" => (AtmosCollected.zvals, Dict()))

        # set up the variables we're going to be writing
        vars = OrderedDict()
        varnames = map(
            prefix_filter,
            flattenednames(vars_atmos_les_default_simple(atmos, FT)),
        )
        ho_varnames = map(
            prefix_filter,
            flattenednames(vars_atmos_les_default_ho(atmos, FT)),
        )
        append!(varnames, ho_varnames)
        for varname in varnames
            var = Variables[varname]
            vars[varname] = (("z",), FT, var.attrib)
        end
        vars["cld_frac"] = (("z",), FT, Variables["cld_frac"].attrib)
        vars["cld_top"] = ((), FT, Variables["cld_top"].attrib)
        vars["cld_base"] = ((), FT, Variables["cld_base"].attrib)
        vars["cld_cover"] = ((), FT, Variables["cld_cover"].attrib)
        vars["lwp"] = ((), FT, Variables["lwp"].attrib)
        vars["iwp"] = ((), FT, Variables["iwp"].attrib)
        vars["rwp"] = ((), FT, Variables["rwp"].attrib)
        vars["swp"] = ((), FT, Variables["swp"].attrib)

        # create the output file
        dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
        dfilename = joinpath(Settings.output_dir, dprefix)
        init_data(dgngrp.writer, dfilename, Settings.no_overwrite, dims, vars)
    end

    return nothing
end

"""
    atmos_les_default_collect(dgngrp, currtime)

Collect the various 'AtmosLESDefault' diagnostic variables for the
current timestep and write them into a file.
"""
function atmos_les_default_collect(dgngrp::DiagnosticsGroup, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    atmos = dg.balance_law
    grid = dg.grid
    grid_info = basic_grid_info(dg)
    topl_info = basic_topology_info(grid.topology)
    Nqk = grid_info.Nqk
    Nqh = grid_info.Nqh
    npoints = prod(grid_info.Nq)
    nrealelem = topl_info.nrealelem
    nvertelem = topl_info.nvertelem
    nhorzelem = topl_info.nhorzrealelem

    energy_model(atmos) isa TotalEnergyModel ||
        error("TotalEnergyModel only supported")

    # get needed arrays onto the CPU
    if array_device(Q) isa CPU
        state_data = Q.realdata
        aux_data = dg.state_auxiliary.realdata
        vgeo = grid.vgeo
        # ω is the weight vector for the vertical direction
        ω = grid.ω[end]
        gradflux_data = dg.state_gradient_flux.realdata
    else
        state_data = Array(Q.realdata)
        aux_data = Array(dg.state_auxiliary.realdata)
        vgeo = Array(grid.vgeo)
        # ω is the weight vector for the vertical direction
        ω = Array(grid.ω[end])
        gradflux_data = Array(dg.state_gradient_flux.realdata)
    end
    FT = eltype(state_data)

    zvals = AtmosCollected.zvals
    MH_z = AtmosCollected.MH_z

    # Visit each node of the state variables array and:
    # - generate and store the thermo variables,
    # - accumulate the simple horizontal sums, and
    # - determine the cloud fraction, top and base
    #
    thermo_array =
        [zeros(FT, num_thermo(atmos, FT)) for _ in 1:npoints, _ in 1:nrealelem]
    simple_sums = [
        zeros(FT, num_atmos_les_default_simple_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    # for liquid, ice, rain and snow water paths
    ρq_liq_z = [zero(FT) for _ in 1:(Nqk * nvertelem)]
    ρq_ice_z = [zero(FT) for _ in 1:(Nqk * nvertelem)]
    ρq_rai_z = [zero(FT) for _ in 1:(Nqk * nvertelem)]
    ρq_sno_z = [zero(FT) for _ in 1:(Nqk * nvertelem)]
    # for cld*
    qc_gt_0_z = [zeros(FT, (Nqh * nhorzelem)) for _ in 1:(Nqk * nvertelem)]
    qc_gt_0_full = zeros(FT, (Nqh * nhorzelem))
    # In honor of PyCLES!
    cld_top = FT(-100000)
    cld_base = FT(100000)
    @traverse_dg_grid grid_info topl_info begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        gradflux = extract_state(dg, gradflux_data, ijk, e, GradientFlux())
        aux = extract_state(dg, aux_data, ijk, e, Auxiliary())
        MH = vgeo[ijk, grid.MHid, e]

        thermo = thermo_vars(atmos, thermo_array[ijk, e])
        compute_thermo!(atmos, state, aux, thermo)

        simple = atmos_les_default_simple_vars(atmos, simple_sums[evk])
        atmos_les_default_simple_sums!(
            atmos,
            state,
            gradflux,
            aux,
            thermo,
            currtime,
            MH,
            simple,
        )

        idx = (Nqh * (eh - 1)) + (grid_info.Nq[2] * (j - 1)) + i
        cld_top, cld_base = atmos_les_default_clouds(
            moisture_model(atmos),
            thermo,
            idx,
            qc_gt_0_z[evk],
            qc_gt_0_full,
            zvals[evk],
            cld_top,
            cld_base,
        )

        # FIXME properly
        if isa(moisture_model(atmos), EquilMoist) ||
           isa(moisture_model(atmos), NonEquilMoist)
            # for LWP
            ρq_liq_z[evk] += MH * thermo.moisture.q_liq * state.ρ * state.ρ
            ρq_ice_z[evk] += MH * thermo.moisture.q_ice * state.ρ * state.ρ
        end
        if isa(precipitation_model(atmos), RainModel)
            # for RWP
            ρq_rai_z[evk] += MH * state.precipitation.ρq_rai * state.ρ
        end
        if isa(precipitation_model(atmos), RainSnowModel)
            # for RWP
            ρq_rai_z[evk] += MH * state.precipitation.ρq_rai * state.ρ
            # for SWP
            ρq_sno_z[evk] += MH * state.precipitation.ρq_sno * state.ρ
        end
    end

    # reduce horizontal sums and cloud data across ranks and compute averages
    simple_avgs = [
        zeros(FT, num_atmos_les_default_simple_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    cld_frac = zeros(FT, Nqk * nvertelem)
    for evk in 1:(Nqk * nvertelem)
        MPI.Allreduce!(simple_sums[evk], simple_avgs[evk], +, mpicomm)
        simple_avgs[evk] .= simple_avgs[evk] ./ MH_z[evk]

        # FIXME properly
        if isa(moisture_model(atmos), EquilMoist) ||
           isa(moisture_model(atmos), NonEquilMoist)
            tot_qc_gt_0_z = MPI.Reduce(sum(qc_gt_0_z[evk]), +, 0, mpicomm)
            tot_horz_z = MPI.Reduce(length(qc_gt_0_z[evk]), +, 0, mpicomm)
            if mpirank == 0
                cld_frac[evk] = tot_qc_gt_0_z / tot_horz_z
            end

            # for LWP and IWP
            tot_ρq_liq_z = MPI.Reduce(ρq_liq_z[evk], +, 0, mpicomm)
            tot_ρq_ice_z = MPI.Reduce(ρq_ice_z[evk], +, 0, mpicomm)
            if mpirank == 0
                ρq_liq_z[evk] = tot_ρq_liq_z / MH_z[evk]
                ρq_ice_z[evk] = tot_ρq_ice_z / MH_z[evk]
            end
        end
        if isa(precipitation_model(atmos), RainModel)
            # for RWP
            tot_ρq_rai_z = MPI.Reduce(ρq_rai_z[evk], +, 0, mpicomm)
            if mpirank == 0
                ρq_rai_z[evk] = tot_ρq_rai_z / MH_z[evk]
            end
        end
        if isa(precipitation_model(atmos), RainSnowModel)
            # for RWP and SWP
            tot_ρq_rai_z = MPI.Reduce(ρq_rai_z[evk], +, 0, mpicomm)
            tot_ρq_sno_z = MPI.Reduce(ρq_sno_z[evk], +, 0, mpicomm)
            if mpirank == 0
                ρq_rai_z[evk] = tot_ρq_rai_z / MH_z[evk]
                ρq_sno_z[evk] = tot_ρq_sno_z / MH_z[evk]
            end
        end
    end
    # FIXME properly
    if isa(moisture_model(atmos), EquilMoist) ||
       isa(moisture_model(atmos), NonEquilMoist)
        cld_top = MPI.Reduce(cld_top, max, 0, mpicomm)
        if cld_top == FT(-100000)
            cld_top = NaN
        end
        cld_base = MPI.Reduce(cld_base, min, 0, mpicomm)
        if cld_base == FT(100000)
            cld_base = NaN
        end
        tot_qc_gt_0_full = MPI.Reduce(sum(qc_gt_0_full), +, 0, mpicomm)
        tot_horz_full = MPI.Reduce(length(qc_gt_0_full), +, 0, mpicomm)
        cld_cover = zero(FT)
        if mpirank == 0
            cld_cover = tot_qc_gt_0_full / tot_horz_full
        end
    end

    simple_varnames = map(
        prefix_filter,
        flattenednames(vars_atmos_les_default_simple(atmos, FT)),
    )

    # complete density averaging
    for evk in 1:(Nqk * nvertelem)
        simple_ha = atmos_les_default_simple_vars(atmos, simple_avgs[evk])
        avg_rho = simple_ha.avg_rho
        for vari in 1:length(simple_varnames)
            if simple_varnames[vari] != "avg_rho"
                simple_avgs[evk][vari] /= avg_rho
            end
        end

        # for all the water paths
        # FIXME properly
        if isa(moisture_model(atmos), EquilMoist) ||
           isa(moisture_model(atmos), NonEquilMoist)
            ρq_liq_z[evk] /= avg_rho
            ρq_ice_z[evk] /= avg_rho
        end
        if isa(precipitation_model(atmos), RainModel)
            ρq_rai_z[evk] /= avg_rho
        end
        if isa(precipitation_model(atmos), RainSnowModel)
            ρq_rai_z[evk] /= avg_rho
            ρq_sno_z[evk] /= avg_rho
        end
    end

    # compute all the water paths
    lwp = NaN
    iwp = NaN
    rwp = NaN
    swp = NaN
    if mpirank == 0
        JcV = reshape(
            view(vgeo, :, grid.JcVid, grid.topology.realelems),
            Nqh,
            Nqk,
            nvertelem,
            nhorzelem,
        )
        Mvert = (ω .* JcV[1, :, :, 1])[:]
        lwp = FT(sum(ρq_liq_z .* Mvert))
        iwp = FT(sum(ρq_ice_z .* Mvert))
        rwp = FT(sum(ρq_rai_z .* Mvert))
        swp = FT(sum(ρq_sno_z .* Mvert))
    end

    # compute the variances and covariances
    ho_sums = [
        zeros(FT, num_atmos_les_default_ho_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    @traverse_dg_grid grid_info topl_info begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        thermo = thermo_vars(atmos, thermo_array[ijk, e])
        MH = vgeo[ijk, grid.MHid, e]

        simple_ha = atmos_les_default_simple_vars(atmos, simple_avgs[evk])
        ho = atmos_les_default_ho_vars(atmos, ho_sums[evk])
        atmos_les_default_ho_sums!(atmos, state, thermo, MH, simple_ha, ho)
    end

    # reduce across ranks and compute averages
    ho_avgs = [
        zeros(FT, num_atmos_les_default_ho_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    for evk in 1:(Nqk * nvertelem)
        MPI.Reduce!(ho_sums[evk], ho_avgs[evk], +, 0, mpicomm)
        if mpirank == 0
            ho_avgs[evk] .= ho_avgs[evk] ./ MH_z[evk]
        end
    end

    # complete density averaging and prepare output
    if mpirank == 0
        varvals = OrderedDict()
        for (vari, varname) in enumerate(simple_varnames)
            davg = zeros(FT, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                davg[evk] = simple_avgs[evk][vari]
            end
            varvals[varname] = davg
        end

        ho_varnames = map(
            prefix_filter,
            flattenednames(vars_atmos_les_default_ho(atmos, FT)),
        )
        for (vari, varname) in enumerate(ho_varnames)
            davg = zeros(FT, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                simple_ha =
                    atmos_les_default_simple_vars(atmos, simple_avgs[evk])
                avg_rho = simple_ha.avg_rho
                davg[evk] = ho_avgs[evk][vari] / avg_rho
            end
            varvals[varname] = davg
        end

        if isa(moisture_model(atmos), EquilMoist) ||
           isa(moisture_model(atmos), NonEquilMoist)
            varvals["cld_frac"] = cld_frac
            varvals["cld_top"] = cld_top
            varvals["cld_base"] = cld_base
            varvals["cld_cover"] = cld_cover
            varvals["lwp"] = lwp
            varvals["iwp"] = iwp
        end
        if isa(precipitation_model(atmos), RainModel)
            varvals["rwp"] = rwp
        end
        if isa(precipitation_model(atmos), RainSnowModel)
            varvals["rwp"] = rwp
            varvals["swp"] = swp
        end

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_les_default_fini(dgngrp::DiagnosticsGroup, currtime) end


text:
include("./topoprovision.jl")

ProvisionIPTopology("B4", scaling=5, density=0.98, linknum=52, ilp_planning=false, tofile=true)
ProvisionIPTopology("IBM", scaling=2, density=0.98, linknum=85, ilp_planning=false, tofile=true)

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
#=
Cross validation tuning methods.
Dependencies: sgd.jl, reformulation.jl, spo_util.jl
=#

using Parameters, Gadfly

## Type for storing validation set algorithm parameters
@with_kw struct val_parms
    algorithm_type::Symbol = :spo_plus_sgd
    validation_set_percent::Float64 = 0.2
    validation_loss::Symbol = :spo_loss
    plot_results::Bool = false
    resolve_sgd::Bool = false
    resolve_sgd_accuracy::Float64 = 0.00001
    resolve_iteration_limit::Int = 50000
end

"""
    validation_set_alg(X::Matrix{Float64}, c::Matrix{Float64}, oracle;
        sp_graph = missing, train_ind = missing, validation_ind = missing,
        val_alg_parms::val_parms = val_parms(), path_alg_parms = sgd_path_parms())

Applies the validation set approach in order to tune the regularization parameter
for either least squares or SPO+ approach. Works with either the SGD path algorithm
for SPO+ or LS, or the shoretst path reformulation approach path algorithm. The user is allowed
to prespecify the training and validation sets by setting the optional parameters `train_ind` and `validation_ind`, which
should be mutually exclusvie and exhaustive subsets of `collect(1:n)`. If these are not provided, then
the algorithm will randomly use `validation_set_percent` of the `n` data points as the validation set.

The default values for the algorithm parameters `val_alg_parms` are (explanations are given below):
    algorithm_type::Symbol = :spo_plus_sgd
    validation_set_percent::Float64 = 0.2
    plot_results::Bool = false

Note that the overall default behavior is to run SGD to compute the SPO+ path.

# Arguments
- `X`: p x n training set feature matrix
- `c`: d x n training set matrix of cost vectors
- `oracle` is a function representing the linear optimization oracle for S, which always needs to be provided (even if doing reformulation)
- `sp_graph` is needed in the reformulation approach and is of the type `shortest_path_graph`. It provides information for reformulation approach (sources, destination, start_node, end_node)
- `train_ind` and `validation_ind` allow the user to specify the training and validation sets
- `algorithm_type` specifies which method to use and is one of:
    - `:spo_plus_sgd` = SGD to compute the path for SPO+ loss with ridge penalty
    - `:ls_sgd` = SGD to compute the path for least squares loss with ridge penalty
    - `:sp_spo_plus_reform` = reformulation approach for SPO+ loss on shortest path problem
    - `:ls_jump` = use JuMP and a solver to compute the least squares path
- `validation_set_percent` = percent of data to use for validations set with indicies are not provided
- If `plot_results` are true, then generate plots and save to file
- `path_alg_parms` specify the parameters of the path solver (either SGD or reformulation/JuMP) and should be of the type `sgd_path_parms` or `reformulation_path_parms`
"""
function validation_set_alg(X::Matrix{Float64}, c::Matrix{Float64}, oracle;
    sp_graph = missing, train_ind = missing, validation_ind = missing,
    val_alg_parms::val_parms = val_parms(), path_alg_parms = sgd_path_parms())

    # Get validation set parms
    @unpack algorithm_type,
        validation_set_percent,
        plot_results,
        validation_loss,
        resolve_sgd,
        resolve_sgd_accuracy,
        resolve_iteration_limit = val_alg_parms

    # Check for validitiy of path_alg_parms
    if typeof(path_alg_parms) != sgd_path_parms && typeof(path_alg_parms) != reformulation_path_parms
        error("path_alg_parms is not valid, should be either sgd_path_parms or reformulation_path_parms")
    end

    # Check for sp_graph
    if !ismissing(sp_graph)
        @unpack sources, destinations, start_node, end_node = sp_graph
    end

    # Get dimensions of input
    (p, n) = size(X)
    (d, n2) = size(c)

    # split into train/validation if not given already
    if ismissing(train_ind) && ismissing(validation_ind)
        validation_ind = randsubseq(collect(1:n), validation_set_percent)
        train_ind = setdiff(collect(1:n), validation_ind)
    end

    X_train = X[:, train_ind]
    X_validation = X[:, validation_ind]
    c_train = c[:, train_ind]
    c_validation = c[:, validation_ind]

    # train models
    if algorithm_type == :spo_plus_sgd
        B_soln_list, lambdas = spoPlus_sgd_path(X_train, c_train, oracle; path_alg_parms = path_alg_parms)
    elseif algorithm_type == :ls_sgd
        B_soln_list, lambdas = leastSquares_sgd_path(X_train, c_train, oracle; path_alg_parms = path_alg_parms)
    elseif algorithm_type == :sp_spo_plus_reform
        B_soln_list, lambdas = sp_reformulation_path_jump(X_train, c_train, sp_graph; sp_oracle = oracle, reform_alg_parms = path_alg_parms)
    elseif algorithm_type == :ls_jump
        B_soln_list, lambdas = leastSquares_path_jump(X_train, c_train; reform_alg_parms = path_alg_parms)
    else
        error("Enter a valid algorithm type")
    end

    # get best model
    num_lambda = length(lambdas)
    validation_loss_list = zeros(num_lambda)
    for i = 1:num_lambda
        if validation_loss == :spo_loss
            validation_loss_list[i] = spo_loss(B_soln_list[i], X_validation, c_validation, oracle)
        elseif validation_loss == :spo_plus_loss
            validation_loss_list[i] = spo_plus_loss(B_soln_list[i], X_validation, c_validation, oracle)
        elseif validation_loss == :least_squares_loss
            validation_loss_list[i] = least_squares_loss(B_soln_list[i], X_validation, c_validation)
        elseif validation_loss == :absolute_loss
            validation_loss_list[i] = absolute_loss(B_soln_list[i], X_validation, c_validation)
        elseif validation_loss == :huber_loss
            validation_loss_list[i] = huber_loss(B_soln_list[i], X_validation, c_validation, path_alg_parms.huber_delta)
        elseif validation_loss == :hamming_loss
            (z_validation, w_validation) = oracle_dataset(c_validation, oracle)
            (d_wval, n_wval) = size(w_validation)
            c_ham_validation = ones(d_wval, n_wval) - w_validation
            validation_loss_list[i] = spo_loss(B_soln_list[i], X_validation, c_ham_validation, oracle)
        else
            error("Enter a valid validation set loss function.")
        end
    end

    if plot_results
        val_plot = plot(x = log.(lambdas), y = validation_loss_list, Geom.line,
        Guide.xlabel("Log(Lambda)"), Guide.ylabel("Validation Set SPO Loss"))
        img = SVG("plots/validation_set_plot_spo_loss.svg", 6inch, 4inch)
        draw(img, val_plot)

        # Plot SPO+ or least squares loss as well
        surrogate_loss = zeros(num_lambda)
        if algorithm_type == :spo_plus_sgd || algorithm_type == :sp_spo_plus_reform
            for i = 1:num_lambda
                surrogate_loss[i] = spo_plus_loss(B_soln_list[i], X_validation, c_validation, oracle)
            end
            val_plot2 = plot(x = log.(lambdas), y = surrogate_loss, Geom.line,
            Guide.xlabel("Log(Lambda)"), Guide.ylabel("Validation Set SPO+ Loss"))
            img = SVG("plots/validation_set_plot_spoplus_loss.svg", 6inch, 4inch)
            draw(img, val_plot2)
        elseif algorithm_type == :ls_sgd || algorithm_type == :ls_jump
            for i = 1:num_lambda
                surrogate_loss[i] = least_squares_loss(B_soln_list[i], X_validation, c_validation)
            end
            val_plot2 = plot(x = log.(lambdas), y = surrogate_loss, Geom.line,
            Guide.xlabel("Log(Lambda)"), Guide.ylabel("Validation Set Least Squares Loss"))
            img = SVG("plots/validation_set_plot_leastsquares_loss.svg", 6inch, 4inch)
            draw(img, val_plot2)
        end
    end

    best_ind = argmin(validation_loss_list)
    best_lambda = lambdas[best_ind]
    best_B_matrix = B_soln_list[best_ind]

    if resolve_sgd && algorithm_type == :spo_plus_sgd
        # Make SGD Alg parms
        n_train = length(train_ind)
        second_moment_bound = (1/n_train)*(norm(X_train)^2)
        num_iters = ceil(Int, (2*second_moment_bound)/(best_lambda*resolve_sgd_accuracy)) + 1
        num_iters = min(num_iters, resolve_iteration_limit)
        sgd_alg_parms = sgd_parms(grad_type = path_alg_parms.grad_type, lambda = best_lambda,
            numiter = num_iters, batchsize = path_alg_parms.batchsize, step_type = :short)

        B_spo_plus_final, spo_sgd_history_trace_final = spoPlus_sgd(X_train, c_train, oracle;
            alg_parms = sgd_alg_parms, B_init = best_B_matrix)
    else
        B_spo_plus_final = best_B_matrix
    end

    return B_spo_plus_final, best_lambda
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
using ESDLPlots
using Test

# write your own tests here
using ESDL
using Widgets
using Compose
import Plots

c = Cube()
d = subsetcube(c,variable = ["gross_primary_productivity","net_ecosystem_exchange","terrestrial_ecosystem_respiration]"], time=Date(2001)..Date(2001,1,15))
pm = plotMAP(d)
pm1 = plotMAP(d,variable="gross",time=Date(2001))
@test isa(pm,Widget)
@test isa(pm1,Context)

pts = plotTS(d);

pts1 = plotTS(d,lon=35,lat=55,var="net");

@test isa(pts,Widget)
@test isa(pts1,Plots.Plot)

pxy = plotXY(d,xaxis="var")

pscatter = plotScatter(d)

pimage = plotMAPRGB(d)

@test isa(pimage,Widget)
@test isa(pxy, Widget)
@test isa(pscatter, Widget)


text:
__precompile__()
module BinaryProvider

using Compat
using Compat.Libdl

# Utilities for controlling verbosity
include("LoggingUtils.jl")
# Include our subprocess running functionality
include("OutputCollector.jl")
# External utilities such as downloading/decompressing tarballs
include("PlatformEngines.jl")
# Platform naming
include("PlatformNames.jl")
# Everything related to file/path management
include("Prefix.jl")
# Abstraction of "needing" a file, that would trigger an install
include("Products.jl")

function __init__()
    global global_prefix

    # Initialize our global_prefix
    global_prefix = Prefix(joinpath(dirname(@__FILE__), "../", "global_prefix"))
    activate(global_prefix)

    # Find the right download/compression engines for this platform
    probe_platform_engines!()

    # If we're on a julia that's too old, then fixup the color mappings
    if !haskey(Base.text_colors, :default)
        Base.text_colors[:default] = Base.color_normal
    end
end

end # module


text:
using JetPackTransforms, Jets, Test

@testset "JetSlantStack, dot product test" for T in (Float32, Float64)
    A = JopSlantStack(JetSpace(T, 64, 128); dz=10.0, dh=10.0, h0=-1000.0)

    m = rand(domain(A))
    d = A*m

    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
    @test isapprox(lhs,rhs,rtol=1e-4)

    A = JopSlantStack(JetSpace(T, 64, 128); dz=10.0, dh=10.0, h0=-1000.0, taperz=(0.3,0.3), taperh=(0.3,0.3), taperkz=(0.3,0.3), taperkh=(0.3,0.3))
    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
    @test isapprox(lhs,rhs,rtol=1e-4)
end

@testset "JetSlantStack, correctness" begin
    A = JopSlantStack(JetSpace(Float64, 64, 128); dz=10.0, dh=10.0, h0=-1000.0)
    m = zeros(domain(A))
    m[32,:] .= 1
    d = A*m
    v,i = findmax(d)
    @test i[1] == 32
    @test i[2] == findfirst(x->x≈0, state(A).cp)
end


text:
struct Staged{T}
    x::T
    k::Int
end

Base.:+(a::Staged, b::Staged) = Staged(a.x+b.x, a.k)
Base.:*(a::Staged, b::Number) = Staged(a.x*b, a.k)
Base.:/(a::Staged, b::Number) = Staged(a.x/b, a.k)

struct LQG1D <: POMDP{Staged{Float64}, Float64, Float64}
    a::Float64 # state transition
    b::Float64 # control
    c::Float64 # observation
    q::Float64 # state cost
    r::Float64 # control cost
    m::Float64 # process noise standard deviation
    n::Float64 # process noise standard deviation
    p::Float64 # initial distribution standard deviation
    k::Int     # number of stages
end

discount(::LQG1D) = 1.0
isterminal(m::LQG1D, s::Staged{Float64}) = s.k == m.k

rand(rng::AbstractRNG, d::Staged{D}) where D = Staged(rand(rng, d.x), d.k)
mean(d::Staged) = Staged(mean(d.x), d.k)

transition(m::LQG1D, s::Staged{Float64}, u::Float64) = Staged(Normal(m.a*s.x + m.b*u, m.m), s.k+1)
observation(m::LQG1D, s::Staged{Float64}) = Normal(m.c*s.x, m.n)
function reward(m::LQG1D, s::Staged{Float64}, u::Float64, sp::Staged{Float64})
    c = m.q*s.x^2 + m.r*u^2
    if isterminal(m, sp)
        c += m.q*sp.x^2
    end
    return -c
end
initial_state_distribution(m::LQG1D) = Staged(Normal(0.0, m.p), 1)

# function MCTS.next_action(gen::RandomActionGenerator, m::LQG1D, b, snode::MCTS.AbstractStateNode)
#     rs = rand(gen.rng, b)
#     x = rs.x
#     return rand(gen.rng, Normal(-1/2*x, 1.0))
# end

function MCTS.next_action(gen::RandomActionGenerator, m::LQG1D, b, snode::MCTS.AbstractStateNode)
    ms = mean(b)
    if ms.k == 1
        return 0.0
    else
        x = ms.x
        return rand(gen.rng, Normal(-1/2*x, 0.01))
    end
end

struct LinearFeedback <: Policy
    k::Float64
end

function action(p::LinearFeedback, b::Staged)
    return -p.k*mean(b.x)
end

struct Kalman1D <: Updater
    m::LQG1D
end

function update(kf::Kalman1D, b::Staged{D}, a::Float64, o::Float64) where D <: Normal
    m = kf.m
    propped = m.a*mean(b.x) + m.b*a
    propped_var = m.a^2*var(b.x) + m.m^2
    variance = propped_var - propped_var^2*m.c^2/(m.c^2*propped_var + m.n)
    xhat = propped + propped + variance*m.c/m.n*(o - m.c*propped)
    return Staged(Normal(xhat, sqrt(variance)), b.k+1)
end


text:
#
# Provides common DocStringExtensions abbreviations
# for UnitfulAstrodynamics.jl sub-modules.
#

@template (FUNCTIONS, METHODS, MACROS) =
"""
$(SIGNATURES)
$(DOCSTRING)
"""

@template (DEFAULT) =
"""
$(DOCSTRING)
"""


text:
export PolyForm, simplify_fractions, quick_cancel, flatten_fractions
using Bijections
using DynamicPolynomials: PolyVar

"""
    PolyForm{T} <: Symbolic{T}

Abstracts a [MultivariatePolynomials.jl](https://juliaalgebra.github.io/MultivariatePolynomials.jl/stable/) as a SymbolicUtils expression and vice-versa.

The SymbolicUtils term interface (`istree`, `operation, and `arguments`) works on PolyForm lazily:
the `operation` and `arguments` are created by converting one level of arguments into SymbolicUtils expressions. They may further contain PolyForm within them.
We use this to hold polynomials in memory while doing `simplify_fractions`.

    PolyForm{T}(x; Fs=Union{typeof(*),typeof(+),typeof(^)}, recurse=false)

Turn a Symbolic expression `x` into a polynomial and return a PolyForm that abstracts it.

`Fs` are the types of functions which should be applied if arguments are themselves
polynomialized. For example, if you only want to polynomialize the base of power
expressions, you would  leave out `typeof(^)` from the union. In this case `^`
is not called, but maintained as a `Pow` term.

`recurse` is a flag which calls `PolyForm` recursively on subexpressions. For example:

```julia
PolyForm(sin((x+y)^2))               #=> sin((x+y)^2)
PolyForm(sin((x+y)^2), recurse=true) #=> sin((x^2 + (2x)y + y^2))
```
"""
struct PolyForm{T, M} <: Symbolic{T}
    p::MP.AbstractPolynomialLike
    pvar2sym::Bijection{Any,Any}   # @polyvar x --> @sym x  etc.
    sym2term::Dict{Sym,Any}        # Symbol("sin-$hash(sin(x+y))") --> sin(x+y) => sin(PolyForm(...))
    metadata::M
end

function (::Type{PolyForm{T}})(p, d1, d2, m=nothing) where {T}
    p isa Number && return p
    p isa MP.AbstractPolynomialLike && MP.isconstant(p) && return convert(Number, p)
    PolyForm{T, typeof(m)}(p, d1, d2, m)
end

Base.hash(p::PolyForm, u::UInt64) = xor(hash(p.p, u),  trunc(UInt, 0xbabacacababacaca))
Base.isequal(x::PolyForm, y::PolyForm) = isequal(x.p, y.p)

# We use the same PVAR2SYM bijection to maintain the PolyVar <-> Sym mapping,
# When all PolyForms go out of scope in a session, we allow it to free up memory and
# start over if necessary
const PVAR2SYM = Ref(WeakRef())
const SYM2TERM = Ref(WeakRef())
clear_dicts() = (PVAR2SYM[] = WeakRef(nothing); SYM2TERM[] = WeakRef(nothing); nothing)
function get_pvar2sym()
    v = PVAR2SYM[].value
    if v === nothing
        d = Bijection{Any,Any}()
        PVAR2SYM[] = WeakRef(d)
        return d
    else
        return v
    end
end

function get_sym2term()
    v = SYM2TERM[].value
    if v === nothing
        d = Dict{Sym,Any}()
        SYM2TERM[] = WeakRef(d)
        return d
    else
        return v
    end
end

function mix_dicts(p, q)
    p.pvar2sym !== q.pvar2sym && error("pvar2sym mappings don't match for $p and $q")
    p.sym2term !== q.sym2term && error("sym2term mappings don't match for $p and $q")

    p.pvar2sym, p.sym2term
end

# forward gcd

PF = :(PolyForm{promote_symtype(/, symtype(x), symtype(y))})
@eval begin
    Base.div(x::PolyForm, y::PolyForm) = $PF(div(x.p, y.p), mix_dicts(x, y)...)
    Base.div(x::Integer, y::PolyForm)  = $PF(div(x, y.p), y.pvar2sym, y.sym2term)
    Base.div(x::PolyForm, y::Integer)  = $PF(div(x.p, y), x.pvar2sym, x.sym2term)

    Base.gcd(x::PolyForm, y::PolyForm) = $PF(_gcd(x.p, y.p), mix_dicts(x, y)...)
    Base.gcd(x::Integer, y::PolyForm)  = $PF(_gcd(x, y.p), y.pvar2sym, y.sym2term)
    Base.gcd(x::PolyForm, y::Integer)  = $PF(_gcd(x.p, y), x.pvar2sym, x.sym2term)
end

_isone(p::PolyForm) = isone(p.p)

function polyize(x, pvar2sym, sym2term, vtype, pow, Fs, recurse)
    if x isa Number
        return x
    elseif istree(x)
        if !(symtype(x) <: Number)
            error("Cannot convert $x of symtype $(symtype(x)) into a PolyForm")
        end

        op = operation(x)
        args = arguments(x)

        local_polyize(y) = polyize(y, pvar2sym, sym2term, vtype, pow, Fs, recurse)

        if typeof(+) <: Fs && op == (+)
            return sum(local_polyize, args)
        elseif typeof(*) <: Fs && op == (*)
            return prod(local_polyize, args)
        elseif typeof(^) <: Fs && op == (^) && args[2] isa Integer && args[2] > 0
            @assert length(args) == 2
            return local_polyize(args[1])^(args[2])
        else
            # create a new symbol to store this

            y = if recurse
                similarterm(x,
                            op,
                            map(a->PolyForm(a, pvar2sym, sym2term, vtype; Fs, recurse),
                                args), symtype(x))
            else
                x
            end

            name = Symbol(string(op), "_", hash(y))

            @label lookup
            sym = Sym{symtype(x)}(name)
            if haskey(sym2term, sym)
                if isequal(sym2term[sym][1], x)
                    return local_polyize(sym)
                else # hash collision
                    name = Symbol(name, "_")
                    @goto lookup
                end
            end

            sym2term[sym] = (x => y)

            return local_polyize(sym)
        end
    elseif issym(x)
        if haskey(active_inv(pvar2sym), x)
            return pvar2sym(x)
        end
        pvar = MP.similarvariable(vtype, nameof(x))
        pvar2sym[pvar] = x
        return pvar
    end
end

function PolyForm(x,
        pvar2sym=get_pvar2sym(),
        sym2term=get_sym2term(),
        vtype=DynamicPolynomials.PolyVar{true};
        Fs = Union{typeof(+), typeof(*), typeof(^)},
        recurse=false,
        metadata=metadata(x))

    if !(symtype(x) <: Number)
        return x
    end

    # Polyize and return a PolyForm
    p = polyize(x, pvar2sym, sym2term, vtype, pow, Fs, recurse)
    PolyForm{symtype(x)}(p, pvar2sym, sym2term, metadata)
end

TermInterface.istree(x::Type{PolyForm}) = true

TermInterface.operation(x::PolyForm) = MP.nterms(x.p) == 1 ? (*) : (+)

function TermInterface.arguments(x::PolyForm{T}) where {T}

    function is_var(v)
        MP.nterms(v) == 1 &&
        isone(MP.coefficient(MP.terms(v)[1])) &&
        MP.degree(MP.monomial(v)) == 1
    end

    function get_var(v)
        # must be called only after a is_var check
        MP.variable(MP.monomial(v))
    end

    function resolve(p)
        !is_var(p) && return p
        pvar = get_var(p)
        s = x.pvar2sym[pvar]
        haskey(x.sym2term, s) ? x.sym2term[s][2] : s
    end

    if MP.nterms(x.p) == 1
        MP.isconstant(x.p) && return [convert(Number, x.p)]
        t = MP.term(x.p)
        c = MP.coefficient(t)
        m = MP.monomial(t)

        if !isone(c)
            [c, (unstable_pow(resolve(v), pow)
                        for (v, pow) in MP.powers(m) if !iszero(pow))...]
        else
            [unstable_pow(resolve(v), pow)
                    for (v, pow) in MP.powers(m) if !iszero(pow)]
        end
    else
        ts = MP.terms(x.p)
        return [MP.isconstant(t) ?
                convert(Number, t) :
                (is_var(t) ?
                 resolve(t) :
                 PolyForm{T, Nothing}(t, x.pvar2sym, x.sym2term, nothing)) for t in ts]
    end
end

Base.show(io::IO, x::PolyForm) = show_term(io, x)

"""
    expand(expr)

Expand expressions by distributing multiplication over addition, e.g.,
`a*(b+c)` becomes `ab+ac`.

`expand` uses replace symbols and non-algebraic expressions by variables of type
`variable_type` to compute the distribution using a specialized sparse
multivariate polynomials implementation.
`variable_type` can be any subtype of `MultivariatePolynomials.AbstractVariable`.
"""
expand(expr) = unpolyize(PolyForm(expr, Fs=Union{typeof(+), typeof(*), typeof(^)}, recurse=true))

function unpolyize(x)
    simterm(x, f, args; kw...) = similarterm(x, f, args, symtype(x); kw...)
    Postwalk(identity, similarterm=simterm)(x)
end

## Rational Polynomial form with Div

function polyform_factors(d, pvar2sym, sym2term)
    make(xs) = map(xs) do x
        if x isa Pow && x.base isa Integer && x.exp > 0
            # here we do want to recurse one level, that's why it's wrong to just
            # use Fs = Union{typeof(+), typeof(*)} here.
            Pow(PolyForm(x.base, pvar2sym, sym2term), x.exp)
        else
            PolyForm(x, pvar2sym, sym2term)
        end
    end

    return make(numerators(d)), make(denominators(d))
end

_mul(xs...) = all(isempty, xs) ? 1 : *(Iterators.flatten(xs)...)

function simplify_div(d::Div)
    d.simplified && return d
    ns, ds = polyform_factors(d, get_pvar2sym(), get_sym2term())
    ns, ds = rm_gcds(ns, ds)
    if all(_isone, ds)
        return isempty(ns) ? 1 : simplify_fractions(_mul(ns))
    else
        Div(simplify_fractions(_mul(ns)), simplify_fractions(_mul(ds)), true)
    end
end

add_divs(x::Div, y::Div) = (x.num * y.den + y.num * x.den) / (x.den * y.den)
add_divs(x::Div, y) = (x.num + y * x.den) / x.den
add_divs(x, y::Div) = (x * y.den + y.num) / y.den
add_divs(x, y) = x + y

"""
    simplify_fractions(x; polyform=false)

Find `Div` nodes and simplify them by cancelling a set of factors of numerators
and denominators. If `polyform=true` the factors which were converted into PolyForm
for the purpose of finding polynomial GCDs will be left as they are.
Note that since PolyForms have different `hash`es than SymbolicUtils expressions,
`substitute` may not work if `polyform=true`
"""
function simplify_fractions(x; polyform=false)
    x = Postwalk(quick_cancel)(x)

    !needs_div_rules(x) && return x

    sdiv(a) = a isa Div ? simplify_div(a) : a

    expr = Postwalk(sdiv ∘ quick_cancel)(Postwalk(add_with_div)(x))

    polyform ? expr : unpolyize(expr)
end

function add_with_div(x, flatten=true)
    (!istree(x) || operation(x) != (+)) && return x
    aa = unsorted_arguments(x)
    !any(a->a isa Div, aa) && return x # no rewrite necessary

    divs = filter(a->a isa Div, aa)
    nondivs = filter(a->!(a isa Div), aa)
    nds = isempty(nondivs) ? 0 : +(nondivs...)
    d = reduce(quick_cancel∘add_divs, divs)
    flatten ? quick_cancel(add_divs(d, nds)) : d + nds
end
"""
    flatten_fractions(x)

Flatten nested fractions that are added together.

```julia
julia> flatten_fractions((1+(1+1/a)/a)/a)
(1 + a + a^2) / (a^3)
```
"""
function flatten_fractions(x)
    Fixpoint(Postwalk(add_with_div))(x)
end

function fraction_iszero(x)
    !istree(x) && return _iszero(x)
    ff = flatten_fractions(x)
    # fast path and then slow path
    any(_iszero, numerators(ff)) ||
    any(_iszero∘expand, numerators(ff))
end

function fraction_isone(x)
    !istree(x) && return _isone(x)
    _isone(simplify_fractions(flatten_fractions(x)))
end

function needs_div_rules(x)
    (x isa Div && !(x.num isa Number) && !(x.den isa Number)) ||
    (istree(x) && operation(x) === (+) && count(has_div, unsorted_arguments(x)) > 1) ||
    (istree(x) && any(needs_div_rules, unsorted_arguments(x)))
end

function has_div(x)
    return x isa Div || (istree(x) && any(has_div, unsorted_arguments(x)))
end

flatten_pows(xs) = map(xs) do x
    x isa Pow ? Iterators.repeated(arguments(x)...) : (x,)
end |> Iterators.flatten |> a->collect(Any,a)

coefftype(x::PolyForm) = coefftype(x.p)
coefftype(x::MP.AbstractPolynomialLike{T}) where {T} = T
coefftype(x) = typeof(x)

const MaybeGcd = Union{PolyForm, MP.AbstractPolynomialLike, Integer}
_gcd(x::MaybeGcd, y::MaybeGcd) = (coefftype(x) <: Complex || coefftype(y) <: Complex) ? 1 : gcd(x, y)
_gcd(x, y) = 1

"""
    quick_cancel(d::Div)

Cancel out matching factors from numerator and denominator.
This is not as effective as `simplify_fractions`, for example,
it wouldn't simplify `(x^2 + 15 -  8x)  / (x - 5)` to `(x - 3)`.
But it will simplify `(x - 5)^2*(x - 3) / (x - 5)` to `(x - 5)*(x - 3)`.
Has optimized processes for `Mul` and `Pow` terms.
"""
quick_cancel(d::Div) = Div{symtype(d)}(quick_cancel(d.num, d.den)...)

quick_cancel(x) = x

quick_cancel(x, y) = isequal(x, y) ? (1,1) : (x, y)

function quick_cancel(x::Pow, y)
    x.exp isa Number || return (x, y)
    isequal(x.base, y) && x.exp >= 1 ? (Pow{symtype(x)}(x.base, x.exp - 1),1) : (x, y)
end

quick_cancel(y, x::Pow) = reverse(quick_cancel(x,y))

function quick_cancel(x::Pow, y::Pow)
    if isequal(x.base, y.base)
        !(x.exp isa Number && y.exp isa Number) && return (x, y)
        if x.exp > y.exp
            return Pow{symtype(x)}(x.base, x.exp-y.exp), 1
        elseif x.exp == y.exp
            return 1, 1
        else # x.exp < y.exp
            return 1, Pow{symtype(y)}(y.base, y.exp-x.exp)
        end
    end
    return x, y
end

function quick_cancel(x::Mul, y)
    if haskey(x.dict, y) && x.dict[y] >= 1
        d = copy(x.dict)
        if d[y] > 1
            d[y] -= 1
        elseif d[y] == 1
            delete!(d, y)
        else
            error("Can't reach")
        end

        return Mul(symtype(x), x.coeff, d), 1
    else
        return x, y
    end
end

function quick_cancel(x::Mul, y::Pow)
    y.exp isa Number || return (x, y)
    if haskey(x.dict, y.base)
        d = copy(x.dict)
        if x.dict[y.base] > y.exp
            d[y.base] -= y.exp
            den = 1
        elseif x.dict[y.base] == y.exp
            delete!(d, y.base)
            den = 1
        else
            den = Pow{symtype(y)}(y.base, y.exp-d[y.base])
            delete!(d, y.base)
        end
        return Mul(symtype(x), x.coeff, d), den
    else
        return x, y
    end
end

quick_cancel(x::Pow, y::Mul) = reverse(quick_cancel(y,x))

quick_cancel(y, x::Mul) = reverse(quick_cancel(x,y))

function quick_cancel(x::Mul, y::Mul)
    num_dict, den_dict = _merge_div(x.dict, y.dict)
    Mul(symtype(x), x.coeff, num_dict), Mul(symtype(y), y.coeff, den_dict)
end

function _merge_div(ndict, ddict)
    num = copy(ndict)
    den = copy(ddict)
    for (k, v) in den
        if haskey(num, k)
            nk = num[k]
            if nk > v
                num[k] -= v
                delete!(den, k)
            elseif nk == v
                delete!(num, k)
                delete!(den, k)
            else
                den[k] -= nk
                delete!(num, k)
            end
        end
    end
    num, den
end

function rm_gcds(ns, ds)
    ns = flatten_pows(ns)
    ds = flatten_pows(ds)

    for i = 1:length(ns)
        for j = 1:length(ds)
            g = _gcd(ns[i], ds[j])
            if !_isone(g)
                ns[i] = div(ns[i], g)
                ds[j] = div(ds[j], g)
            end
        end
    end

    filter!(!_isone, ns)
    filter!(!_isone, ds)

    ns,ds
end


text:
#=
   Lotto Problem

   From some unknown source (but I quoted it here: http://www.hakank.org/sims/simulering.html)
   """
   One day I was asked

   What is the probability that no two numbers are consecutive in
   a lottery draw?

   Information:

   National Lottery in UK is to draw 6 different numbers from 1 to 49.
   so, for example,

    (a)   1  4  7   12  19  44    - no consecutive numbers
    (b)   3  6  17  18  44  46    - 17 and 18 are consecutive
    (c)   1  2  3   17  29  49    - 1, 2 and 3 are consecutive

   We are asking the probability that class (a) occurs. Hope this is clear.

   Observation shows that it is NOT a small number (actually near half-and-half).
   ......

   A Monte Carlo simulation experiment produced 50,558 sets of six numbers
   between 1 and 49 with none consecutive out of 100,000 trials giving an
   estimated probability of .50558.
   """

   Summary Statistics
   parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

        prob    0.5054    0.5000     0.0025    0.0034   22057.4684    1.0000     1884.6094

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function lotto_not_two_consecutive_numbers(n=49,m=6)

    lotto ~ filldist(DiscreteUniform(1,n),m)
    # Ensure unicity
    true ~ Dirac(sum([lotto[i] == lotto[j] for i in 1:m for j in i+1:m]) == 0)
    # Are there any consecutive numbers?
    prob ~ Dirac(sum([ abs(lotto[i]-lotto[j]) == 1 for i in 1:m for j in i+1:m]) == 0)
end

n = 49
m = 6
model = lotto_not_two_consecutive_numbers(n,m)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 40_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns[[:prob]])
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module MatrixMarket

export mmread

"""
### mmread(filename, infoonly::Bool=false)

Read the contents of the Matrix Market file 'filename' into a matrix,
which will be either sparse or dense, depending on the Matrix Market format
indicated by 'coordinate' (coordinate sparse storage), or 'array' (dense
array storage).

If infoonly is true (default: false), only information on the size and
structure is returned from reading the header. The actual data for the
matrix elements are not parsed.
"""
function mmread(filename, infoonly::Bool=false)
    mmfile = open(filename,"r")
    # Read first line
    firstline = chomp(readline(mmfile))
    tokens = split(firstline)
    if length(tokens) != 5
        throw(ParseError(string("Not enough words on first line: ", ll)))
    end
    if tokens[1] != "%%MatrixMarket"
        throw(ParseError(string("Not a valid MatrixMarket header:", ll)))
    end
    (head1, rep, field, symm) = map(lowercase, tokens[2:5])
    if head1 != "matrix"
        throw(ParseError("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)"))
    end

    eltype = field == "real" ? Float64 :
             field == "complex" ? Complex128 :
             field == "pattern" ? Bool :
             throw(ParseError("Unsupported field $field (only real and complex are supported)"))

    symlabel = symm == "general" ? identity :
               symm == "symmetric" ? symmetric! :
               symm == "hermitian" ? hermitian! :
               symm == "skew-symmetric" ? skewsymmetric! :
               throw(ParseError("Unknown matrix symmetry: $symm (only general, symmetric, skew-symmetric and hermitian are supported)"))

    # Skip all comments and empty lines
    ll   = readline(mmfile)
    while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
        ll = readline(mmfile)
    end
    # Read matrix dimensions (and number of entries) from first non-comment line
    dd = map(parseint, split(ll))
    if length(dd) < (rep == "coordinate" ? 3 : 2)
        throw(ParseError(string("Could not read in matrix dimensions from line: ", ll)))
    end
    rows = dd[1]
    cols = dd[2]
    entries = (rep == "coordinate") ? dd[3] : (rows * cols)
    if infoonly
        return (rows, cols, entries, rep, field, symm)
    end
    if rep == "coordinate"
        rr = Array(Int, entries)
        cc = Array(Int, entries)
        xx = Array(eltype, entries)
        for i in 1:entries
            flds = split(readline(mmfile))
            rr[i] = parseint(flds[1])
            cc[i] = parsefloat(flds[2])
            if eltype == Complex128
                xx[i] = Complex128(parsefloat(flds[3]), parsefloat(flds[4]))
            elseif eltype == Float64
                xx[i] = parsefloat(flds[3])
            else
                xx[i] = true
            end
        end
        return symlabel(sparse(rr, cc, xx, rows, cols))
    end
    return symlabel(reshape([parsefloat(readline(mmfile)) for i in 1:entries], (rows,cols)))
end

# Hack to represent skew-symmetric matrix as an ordinary matrix with duplicated elements
function skewsymmetric!(M::AbstractMatrix)
    m,n = size(M)
    m == n || throw(DimensionMismatch())
    for i=1:n, j=1:n
        if M[i,j] != 0
            M[j,i] = -M[i,j]
        end
    end
    return M
end

function symmetric!(M::AbstractMatrix)
    m,n = size(M)
    m == n || throw(DimensionMismatch())
    for i=1:n, j=1:n
        if M[i,j] != 0
            M[j,i] = M[i,j]
        end
    end
    return M
end

function hermitian!(M::AbstractMatrix)
    m,n = size(M)
    m == n || throw(DimensionMismatch())
    for i=1:n, j=1:n
        if M[i,j] != 0
            M[j,i] = conj(M[i,j])
        end
    end
    return M
end

end # module