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

# highlighting settings
highlighting = Dict{Symbol, Bool}(
    :warntype => true,
    :llvm => true,
    :native => true,
)

llstyle = Dict{Symbol, Tuple{Bool, Union{Symbol, Int}}}(
    :default     => (false, :normal), # e.g. comma, equal sign, unknown token
    :comment     => (false, :light_black),
    :label       => (false, :light_red),
    :instruction => ( true, :light_cyan),
    :type        => (false, :cyan),
    :number      => (false, :yellow),
    :bracket     => (false, :yellow),
    :variable    => (false, :normal), # e.g. variable, register
    :keyword     => (false, :light_magenta),
    :funcname    => (false, :light_yellow),
)

function printstyled_ll(io::IO, x, s::Symbol, trailing_spaces="")
    printstyled(io, x, bold=llstyle[s][1], color=llstyle[s][2])
    print(io, trailing_spaces)
end

# displaying type warnings

function warntype_type_printer(io::IO, @nospecialize(ty), used::Bool)
    used || return
    str = "::$ty"
    if !highlighting[:warntype]
        print(io, str)
    elseif ty isa Union && Base.is_expected_union(ty)
        Base.emphasize(io, str, Base.warn_color()) # more mild user notification
    elseif ty isa Type && (!Base.isdispatchelem(ty) || ty == Core.Box)
        Base.emphasize(io, str)
    else
        Base.printstyled(io, str, color=:cyan) # show the "good" type
    end
    nothing
end

"""
    code_warntype([io::IO], f, types; debuginfo=:default)

Prints lowered and type-inferred ASTs for the methods matching the given generic function
and type signature to `io` which defaults to `stdout`. The ASTs are annotated in such a way
as to cause "non-leaf" types to be emphasized (if color is available, displayed in red).
This serves as a warning of potential type instability. Not all non-leaf types are particularly
problematic for performance, so the results need to be used judiciously.
In particular, unions containing either [`missing`](@ref) or [`nothing`](@ref) are displayed in yellow, since
these are often intentional.

Keyword argument `debuginfo` may be one of `:source` or `:none` (default), to specify the verbosity of code comments.

See [`@code_warntype`](@ref man-code-warntype) for more information.
"""
function code_warntype(io::IO, @nospecialize(f), @nospecialize(t=Base.default_tt(f));
                       debuginfo::Symbol=:default, optimize::Bool=false, kwargs...)
    debuginfo = Base.IRShow.debuginfo(debuginfo)
    lineprinter = Base.IRShow.__debuginfo[debuginfo]
    for (src, rettype) in code_typed(f, t; optimize, kwargs...)
        lambda_io::IOContext = io
        p = src.parent
        nargs::Int = 0
        if p isa Core.MethodInstance
            println(io, p)
            print(io, "  from ")
            println(io, p.def)
            p.def isa Method && (nargs = p.def.nargs)
            if !isempty(p.sparam_vals)
                println(io, "Static Parameters")
                sig = p.def.sig
                warn_color = Base.warn_color() # more mild user notification
                for i = 1:length(p.sparam_vals)
                    sig = sig::UnionAll
                    name = sig.var.name
                    val = p.sparam_vals[i]
                    print_highlighted(io::IO, v::String, color::Symbol) =
                        if highlighting[:warntype]
                            Base.printstyled(io, v; color)
                        else
                            Base.print(io, v)
                        end
                    if val isa TypeVar
                        if val.lb === Union{}
                            print(io, "  ", name, " <: ")
                            print_highlighted(io, "$(val.ub)", warn_color)
                        elseif val.ub === Any
                            print(io, "  ", sig.var.name, " >: ")
                            print_highlighted(io, "$(val.lb)", warn_color)
                        else
                            print(io, "  ")
                            print_highlighted(io, "$(val.lb)", warn_color)
                            print(io, " <: ", sig.var.name, " <: ")
                            print_highlighted(io, "$(val.ub)", warn_color)
                        end
                    elseif val isa typeof(Vararg)
                        print(io, "  ", name, "::")
                        print_highlighted(io, "Int", warn_color)
                    else
                        print(io, "  ", sig.var.name, " = ")
                        print_highlighted(io, "$(val)", :cyan) # show the "good" type
                    end
                    println(io)
                    sig = sig.body
                end
            end
        end
        if src.slotnames !== nothing
            slotnames = Base.sourceinfo_slotnames(src)
            lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES => slotnames)
            slottypes = src.slottypes
            nargs > 0 && println(io, "Arguments")
            for i = 1:length(slotnames)
                if i == nargs + 1
                    println(io, "Locals")
                end
                print(io, "  ", slotnames[i])
                if isa(slottypes, Vector{Any})
                    warntype_type_printer(io, slottypes[i], true)
                end
                println(io)
            end
        end
        print(io, "Body")
        warntype_type_printer(io, rettype, true)
        println(io)
        irshow_config = Base.IRShow.IRShowConfig(lineprinter(src), warntype_type_printer)
        Base.IRShow.show_ir(lambda_io, src, irshow_config)
        println(io)
    end
    nothing
end
code_warntype(@nospecialize(f), @nospecialize(t=Base.default_tt(f)); kwargs...) =
    code_warntype(stdout, f, t; kwargs...)

import Base.CodegenParams

# Printing code representations in IR and assembly
function _dump_function(@nospecialize(f), @nospecialize(t), native::Bool, wrapper::Bool,
                        strip_ir_metadata::Bool, dump_module::Bool, syntax::Symbol,
                        optimize::Bool, debuginfo::Symbol, binary::Bool,
                        params::CodegenParams=CodegenParams(debug_info_kind=Cint(0)))
    ccall(:jl_is_in_pure_context, Bool, ()) && error("code reflection cannot be used from generated functions")
    if isa(f, Core.Builtin)
        throw(ArgumentError("argument is not a generic function"))
    end
    # get the MethodInstance for the method match
    world = typemax(UInt)
    match = Base._which(signature_type(f, t), world)
    linfo = Core.Compiler.specialize_method(match)
    # get the code for it
    if debuginfo === :default
        debuginfo = :source
    elseif debuginfo !== :source && debuginfo !== :none
        throw(ArgumentError("'debuginfo' must be either :source or :none"))
    end
    if native
        if syntax !== :att && syntax !== :intel
            throw(ArgumentError("'syntax' must be either :intel or :att"))
        end
        if dump_module
            str = _dump_function_linfo_native(linfo, world, wrapper, syntax, debuginfo, binary, params)
        else
            str = _dump_function_linfo_native(linfo, world, wrapper, syntax, debuginfo, binary)
        end
    else
        str = _dump_function_linfo_llvm(linfo, world, wrapper, strip_ir_metadata, dump_module, optimize, debuginfo, params)
    end
    # TODO: use jl_is_cacheable_sig instead of isdispatchtuple
    isdispatchtuple(linfo.specTypes) || (str = "; WARNING: This code may not match what actually runs.\n" * str)
    return str
end

function _dump_function_linfo_native(linfo::Core.MethodInstance, world::UInt, wrapper::Bool, syntax::Symbol, debuginfo::Symbol, binary::Bool)
    str = ccall(:jl_dump_method_asm, Ref{String},
                (Any, UInt, Bool, Bool, Ptr{UInt8}, Ptr{UInt8}, Bool),
                linfo, world, false, wrapper, syntax, debuginfo, binary)
    return str
end

function _dump_function_linfo_native(linfo::Core.MethodInstance, world::UInt, wrapper::Bool, syntax::Symbol, debuginfo::Symbol, binary::Bool, params::CodegenParams)
    llvmf = ccall(:jl_get_llvmf_defn, Ptr{Cvoid}, (Any, UInt, Bool, Bool, CodegenParams), linfo, world, wrapper, true, params)
    llvmf == C_NULL && error("could not compile the specified method")
    str = ccall(:jl_dump_function_asm, Ref{String},
                (Ptr{Cvoid}, Bool, Ptr{UInt8}, Ptr{UInt8}, Bool),
                llvmf, false, syntax, debuginfo, binary)
    return str
end

function _dump_function_linfo_llvm(
        linfo::Core.MethodInstance, world::UInt, wrapper::Bool,
        strip_ir_metadata::Bool, dump_module::Bool,
        optimize::Bool, debuginfo::Symbol,
        params::CodegenParams)
    llvmf = ccall(:jl_get_llvmf_defn, Ptr{Cvoid}, (Any, UInt, Bool, Bool, CodegenParams), linfo, world, wrapper, optimize, params)
    llvmf == C_NULL && error("could not compile the specified method")
    str = ccall(:jl_dump_function_ir, Ref{String},
                (Ptr{Cvoid}, Bool, Bool, Ptr{UInt8}),
                llvmf, strip_ir_metadata, dump_module, debuginfo)
    return str
end

"""
    code_llvm([io=stdout,], f, types; raw=false, dump_module=false, optimize=true, debuginfo=:default)

Prints the LLVM bitcodes generated for running the method matching the given generic
function and type signature to `io`.

If the `optimize` keyword is unset, the code will be shown before LLVM optimizations.
All metadata and dbg.* calls are removed from the printed bitcode. For the full IR, set the `raw` keyword to true.
To dump the entire module that encapsulates the function (with declarations), set the `dump_module` keyword to true.
Keyword argument `debuginfo` may be one of source (default) or none, to specify the verbosity of code comments.
"""
function code_llvm(io::IO, @nospecialize(f), @nospecialize(types), raw::Bool,
                   dump_module::Bool=false, optimize::Bool=true, debuginfo::Symbol=:default)
    d = _dump_function(f, types, false, false, !raw, dump_module, :att, optimize, debuginfo, false)
    if highlighting[:llvm] && get(io, :color, false)
        print_llvm(io, d)
    else
        print(io, d)
    end
end
code_llvm(io::IO, @nospecialize(f), @nospecialize(types=Base.default_tt(f)); raw::Bool=false, dump_module::Bool=false, optimize::Bool=true, debuginfo::Symbol=:default) =
    code_llvm(io, f, types, raw, dump_module, optimize, debuginfo)
code_llvm(@nospecialize(f), @nospecialize(types=Base.default_tt(f)); raw=false, dump_module=false, optimize=true, debuginfo::Symbol=:default) =
    code_llvm(stdout, f, types; raw, dump_module, optimize, debuginfo)

"""
    code_native([io=stdout,], f, types; syntax=:att, debuginfo=:default, binary=false, dump_module=true)

Prints the native assembly instructions generated for running the method matching the given
generic function and type signature to `io`.
Switch assembly syntax using `syntax` symbol parameter set to `:att` for AT&T syntax or `:intel` for Intel syntax.
Keyword argument `debuginfo` may be one of source (default) or none, to specify the verbosity of code comments.
If `binary` is `true`, it also prints the binary machine code for each instruction precedented by an abbreviated address.
"""
function code_native(io::IO, @nospecialize(f), @nospecialize(types=Base.default_tt(f));
                     dump_module::Bool=true, syntax::Symbol=:att, debuginfo::Symbol=:default, binary::Bool=false)
    d = _dump_function(f, types, true, false, false, dump_module, syntax, true, debuginfo, binary)
    if highlighting[:native] && get(io, :color, false)
        print_native(io, d)
    else
        print(io, d)
    end
end
code_native(@nospecialize(f), @nospecialize(types=Base.default_tt(f)); dump_module::Bool=true, syntax::Symbol=:att, debuginfo::Symbol=:default, binary::Bool=false) =
    code_native(stdout, f, types; dump_module, syntax, debuginfo, binary)
code_native(::IO, ::Any, ::Symbol) = error("invalid code_native call") # resolve ambiguous call

## colorized IR and assembly printing

const num_regex = r"^(?:\$?-?\d+|0x[0-9A-Fa-f]+|-?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?)$"

function print_llvm(io::IO, code::String)
    buf = IOBuffer(code)
    for line in eachline(buf)
        m = match(r"^(\s*)((?:[^;]|;\")*)(.*)$", line)
        m === nothing && continue
        indent, tokens, comment = m.captures
        print(io, indent)
        print_llvm_tokens(io, tokens)
        printstyled_ll(io, comment, :comment)
        println(io)
    end
end

const llvm_types =
    r"^(?:void|half|float|double|x86_\w+|ppc_\w+|label|metadata|type|opaque|token|i\d+)$"
const llvm_cond = r"^(?:[ou]?eq|[ou]?ne|[uso][gl][te]|ord|uno)$" # true|false

function print_llvm_tokens(io, tokens)
    m = match(r"^((?:[^\s:]+:)?)(\s*)(.*)", tokens)
    if m !== nothing
        label, spaces, tokens = m.captures
        printstyled_ll(io, label, :label, spaces)
    end
    m = match(r"^(%[^\s=]+)(\s*)=(\s*)(.*)", tokens)
    if m !== nothing
        result, spaces, spaces2, tokens = m.captures
        printstyled_ll(io, result, :variable, spaces)
        printstyled_ll(io, '=', :default, spaces2)
    end
    m = match(r"^([a-z]\w*)(\s*)(.*)", tokens)
    if m !== nothing
        inst, spaces, tokens = m.captures
        iskeyword = occursin(r"^(?:define|declare|type)$", inst) || occursin("=", tokens)
        printstyled_ll(io, inst, iskeyword ? :keyword : :instruction, spaces)
    end

    print_llvm_operands(io, tokens)
end

function print_llvm_operands(io, tokens)
    while !isempty(tokens)
        tokens = print_llvm_operand(io, tokens)
    end
    return tokens
end

function print_llvm_operand(io, tokens)
    islabel = false
    while !isempty(tokens)
        m = match(r"^,(\s*)(.*)", tokens)
        if m !== nothing
            spaces, tokens = m.captures
            printstyled_ll(io, ',', :default, spaces)
            break
        end
        m = match(r"^(\*+|=)(\s*)(.*)", tokens)
        if m !== nothing
            sym, spaces, tokens = m.captures
            printstyled_ll(io, sym, :default, spaces)
            continue
        end
        m = match(r"^(\"[^\"]*\")(\s*)(.*)", tokens)
        if m !== nothing
            str, spaces, tokens = m.captures
            printstyled_ll(io, str, :variable, spaces)
            continue
        end
        m = match(r"^([({\[<])(\s*)(.*)", tokens)
        if m !== nothing
            bracket, spaces, tokens = m.captures
            printstyled_ll(io, bracket, :bracket, spaces)
            tokens = print_llvm_operands(io, tokens) # enter
            continue
        end
        m = match(r"^([)}\]>])(\s*)(.*)", tokens)
        if m !== nothing
            bracket, spaces, tokens = m.captures
            printstyled_ll(io, bracket, :bracket, spaces)
            break # leave
        end

        m = match(r"^([^\s,*=(){}\[\]<>]+)(\s*)(.*)", tokens)
        m === nothing && break
        token, spaces, tokens = m.captures
        if occursin(llvm_types, token)
            printstyled_ll(io, token, :type)
            islabel = token == "label"
        elseif occursin(llvm_cond, token) # condition code is instruction-level
            printstyled_ll(io, token, :instruction)
        elseif occursin(num_regex, token)
            printstyled_ll(io, token, :number)
        elseif occursin(r"^@.+$", token)
            printstyled_ll(io, token, :funcname)
        elseif occursin(r"^%.+$", token)
            islabel |= occursin(r"^%[^\d].*$", token) & occursin(r"^\]", tokens)
            printstyled_ll(io, token, islabel ? :label : :variable)
            islabel = false
        elseif occursin(r"^[a-z]\w+$", token)
            printstyled_ll(io, token, :keyword)
        else
            printstyled_ll(io, token, :default)
        end
        print(io, spaces)
    end
    return tokens
end

function print_native(io::IO, code::String, arch::Symbol=sys_arch_category())
    archv = Val(arch)
    buf = IOBuffer(code)
    for line in eachline(buf)
        m = match(r"^(\s*)((?:[^;#/]|#\S|;\"|/[^/])*)(.*)$", line)
        m === nothing && continue
        indent, tokens, comment = m.captures
        print(io, indent)
        print_native_tokens(io, tokens, archv)
        printstyled_ll(io, comment, :comment)
        println(io)
    end
end

function sys_arch_category()
    if Sys.ARCH === :x86_64 || Sys.ARCH === :i686
        :x86
    elseif Sys.ARCH === :aarch64 || startswith(string(Sys.ARCH), "arm")
        :arm
    else
        :unsupported
    end
end

print_native_tokens(io, line, ::Val) = print(io, line)

const x86_ptr = r"^(?:(?:[xyz]mm|[dq])?word|byte|ptr|offset)$"
const avx512flags = r"^(?:z|r[nduz]-sae|sae|1to1?\d)$"
const arm_cond = r"^(?:eq|ne|cs|ho|cc|lo|mi|pl|vs|vc|hi|ls|[lg][te]|al|nv)$"
const arm_keywords = r"^(?:lsl|lsr|asr|ror|rrx|!|/[zm])$"

function print_native_tokens(io, tokens, arch::Union{Val{:x86}, Val{:arm}})
    x86 = arch isa Val{:x86}
    m = match(r"^((?:[^\s:]+:|\"[^\"]+\":)?)(\s*)(.*)", tokens)
    if m !== nothing
        label, spaces, tokens = m.captures
        printstyled_ll(io, label, :label, spaces)
    end
    haslabel = false
    m = match(r"^([a-z][\w.]*)(\s*)(.*)", tokens)
    if m !== nothing
        instruction, spaces, tokens = m.captures
        printstyled_ll(io, instruction, :instruction, spaces)
        haslabel = occursin(r"^(?:bl?|bl?\.\w{2,5}|[ct]bn?z)?$", instruction)
    end

    isfuncname = false
    while !isempty(tokens)
        m = match(r"^([,:*])(\s*)(.*)", tokens)
        if m !== nothing
            sym, spaces, tokens = m.captures
            printstyled_ll(io, sym, :default, spaces)
            isfuncname = false
            continue
        end
        m = match(r"^([(){}\[\]])(\s*)(.*)", tokens)
        if m !== nothing
            bracket, spaces, tokens = m.captures
            printstyled_ll(io, bracket, :bracket, spaces)
            continue
        end
        m = match(r"^#([0-9a-fx.-]+)(\s*)(.*)", tokens)
        if !x86 && m !== nothing && occursin(num_regex, m.captures[1])
            num, spaces, tokens = m.captures
            printstyled_ll(io, "#" * num, :number, spaces)
            continue
        end

        m = match(r"^([^\s,:*(){}\[\]][^\s,:*/(){}\[\]]*)(\s*)(.*)", tokens)
        m === nothing && break
        token, spaces, tokens = m.captures
        if occursin(num_regex, token)
            printstyled_ll(io, token, :number)
        elseif x86 && occursin(x86_ptr, token) || occursin(avx512flags, token)
            printstyled_ll(io, token, :keyword)
            isfuncname = token == "offset"
        elseif !x86 && (occursin(arm_keywords, token) || occursin(arm_cond, token))
            printstyled_ll(io, token, :keyword)
        elseif occursin(r"^L.+$", token)
            printstyled_ll(io, token, :label)
        elseif occursin(r"^\$.+$", token)
            printstyled_ll(io, token, :funcname)
        elseif occursin(r"^%?(?:[a-z][\w.]+|\"[^\"]+\")$", token)
            islabel = haslabel & !occursin(',', tokens)
            printstyled_ll(io, token, islabel ? :label : isfuncname ? :funcname : :variable)
            isfuncname = false
        else
            printstyled_ll(io, token, :default)
        end
        print(io, spaces)
    end
end

max_stars_repo_name: wentasah/julia

max_issues_repo_name: wentasah/julia

max_forks_repo_name: wentasah/julia


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

struct DiskEncryptionSetsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a disk encryption set
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskEncryptionSetName::String (required)
Param: api_version::String (required)
Param: diskEncryptionSet::DiskEncryptionSet (required)
Return: DiskEncryptionSet
"""
function _swaggerinternal_diskEncryptionSetsCreateOrUpdate(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String, diskEncryptionSet; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", DiskEncryptionSet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/diskEncryptionSets/{diskEncryptionSetName}", ["azure_auth"], diskEncryptionSet)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskEncryptionSetName", diskEncryptionSetName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function diskEncryptionSetsCreateOrUpdate(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String, diskEncryptionSet; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsCreateOrUpdate(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version, diskEncryptionSet; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function diskEncryptionSetsCreateOrUpdate(_api::DiskEncryptionSetsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String, diskEncryptionSet; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsCreateOrUpdate(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version, diskEncryptionSet; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes a disk encryption set.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskEncryptionSetName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_diskEncryptionSetsDelete(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/diskEncryptionSets/{diskEncryptionSetName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskEncryptionSetName", diskEncryptionSetName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function diskEncryptionSetsDelete(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsDelete(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function diskEncryptionSetsDelete(_api::DiskEncryptionSetsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsDelete(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets information about a disk encryption set.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskEncryptionSetName::String (required)
Param: api_version::String (required)
Return: DiskEncryptionSet
"""
function _swaggerinternal_diskEncryptionSetsGet(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DiskEncryptionSet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/diskEncryptionSets/{diskEncryptionSetName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskEncryptionSetName", diskEncryptionSetName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function diskEncryptionSetsGet(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsGet(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function diskEncryptionSetsGet(_api::DiskEncryptionSetsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsGet(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all the disk encryption sets under a subscription.
Param: subscriptionId::String (required)
Param: api_version::String (required)
Return: DiskEncryptionSetList
"""
function _swaggerinternal_diskEncryptionSetsList(_api::DiskEncryptionSetsApi, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DiskEncryptionSetList, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/diskEncryptionSets", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function diskEncryptionSetsList(_api::DiskEncryptionSetsApi, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsList(_api, subscriptionId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function diskEncryptionSetsList(_api::DiskEncryptionSetsApi, response_stream::Channel, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsList(_api, subscriptionId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all the disk encryption sets under a resource group.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Return: DiskEncryptionSetList
"""
function _swaggerinternal_diskEncryptionSetsListByResourceGroup(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DiskEncryptionSetList, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/diskEncryptionSets", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function diskEncryptionSetsListByResourceGroup(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function diskEncryptionSetsListByResourceGroup(_api::DiskEncryptionSetsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates (patches) a disk encryption set.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskEncryptionSetName::String (required)
Param: api_version::String (required)
Param: diskEncryptionSet::DiskEncryptionSetUpdate (required)
Return: DiskEncryptionSet
"""
function _swaggerinternal_diskEncryptionSetsUpdate(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String, diskEncryptionSet; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", DiskEncryptionSet, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/diskEncryptionSets/{diskEncryptionSetName}", ["azure_auth"], diskEncryptionSet)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskEncryptionSetName", diskEncryptionSetName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function diskEncryptionSetsUpdate(_api::DiskEncryptionSetsApi, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String, diskEncryptionSet; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsUpdate(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version, diskEncryptionSet; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function diskEncryptionSetsUpdate(_api::DiskEncryptionSetsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskEncryptionSetName::String, api_version::String, diskEncryptionSet; _mediaType=nothing)
    _ctx = _swaggerinternal_diskEncryptionSetsUpdate(_api, subscriptionId, resourceGroupName, diskEncryptionSetName, api_version, diskEncryptionSet; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export diskEncryptionSetsCreateOrUpdate, diskEncryptionSetsDelete, diskEncryptionSetsGet, diskEncryptionSetsList, diskEncryptionSetsListByResourceGroup, diskEncryptionSetsUpdate

ext: jl


text:
tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_sbm_gwf_model(config)
@unpack network = model

model = Wflow.update_sbm_gwf(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    csv = CSV.File(model.writer.csv_path)
    row = csv[1]

    @test row.time == DateTime("2000-06-01T00:00:00")
    @test row.Q_av ≈ 0.016260519762890745
    @test row.head ≈ 1.8416896245327632
end

@testset "first timestep" begin
    sbm = model.vertical

    @test model.clock.iteration == 2
    @test sbm.altitude[1] == 1.937000036239624
    @test sbm.θₛ[1] == 0.44999998807907104
    @test sbm.runoff[1] == 0.0
    @test sbm.soilevap[1] == 0.0
    @test sbm.transpiration[1] == 0.6117526566330049
end

# run the second timestep
model = Wflow.update_sbm_gwf(model)

@testset "second timestep" begin
    sbm = model.vertical
    @test sbm.altitude[1] == 1.937000036239624
    @test sbm.θₛ[1] == 0.44999998807907104
    @test sbm.runoff[1] == 0.0
    @test sbm.soilevap[1] == 0.0
    @test sbm.transpiration[1] == 1.0122634204681036
end

@testset "overland flow" begin
    q = model.lateral.land.q_av
    @test sum(q) ≈ 0.0
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 0.03515256649675796
    @test q[6] ≈ 0.00795266793733648
    @test q[13] ≈  0.0005980549864213564
    @test q[network.river.order[end]] ≈  0.008482646477950898
end

@testset "groundwater" begin
    gw = model.lateral.subsurface
    @test gw.river.stage[1] ≈ 1.21057153442702
    @test gw.flow.aquifer.head[19] ≈ 1.7999999523162842
    @test gw.river.flux[1] ≈ -50.751898489778426
    @test gw.drain.flux[1] ≈  0.0
    @test gw.recharge.rate[19] ≈  -0.0014241196552847502
end

Wflow.close_files(model, delete_output = false)


text:
import LPDM: bounds, best_lb_action, best_ub_action

mutable struct LDBounds1d{S,A,O}
    lb_    ::Float64
    ub_    ::Float64
    best_lb_action_::Float64
    best_ub_action_::Float64
    rng::RNGVector # dummy RNG, not used

    function LDBounds1d{S,A,O}(::POMDP{S,A,O}) where {S,A,O}
        this = new()
        this.lb_ = +Inf
        this.ub_ = -Inf
        this.best_lb_action_ = NaN
        this.best_ub_action_ = NaN
        this.rng = RNGVector(1, UInt32(42))
        return this
    end
end

function LPDM.bounds(b::LDBounds1d{S,A,O},
                     pomdp::AbstractLD1,
                     particles::Vector{LPDMParticle{S}},
                     config::LPDMConfig) where {S,A,O}

    # reset every time bounds are about to be recomputed
    b.lb_ = +Inf
    b.ub_ = -Inf
    b.best_lb_action_ = NaN
    b.best_ub_action_ = NaN

    tmp_lb = NaN
    tmp_ub = NaN
    tmp_lb_action = NaN
    tmp_ub_action = NaN

    for p in particles
        tmp_ub, tmp_ub_action = upper_bound(pomdp,p)
        if p.state >= pomdp.min_noise_loc # assume min_noise_loc > 0
            # both will have to do roughly the same thing (+/- discretization differences),
            # so make them the same
            #TODO: review for lb > ub
            tmp_lb, tmp_lb_action = tmp_ub, tmp_ub_action
        else
            tmp_lb, tmp_lb_action = lower_bound(pomdp,p)
        end

        if tmp_lb < b.lb_
            b.lb_             = tmp_lb
            b.best_lb_action_ = tmp_lb_action
            # if abs(tmp_lb + 30.0) < 0.1 # DEBUG: remove
            #     println("particles: $particles")
            #     println("s: $(p.state), b.lb_=$(b.lb_), tmp_lb_action = $(tmp_lb_action), b.lba_ = $(b.best_lb_action_)")
            # end
        end
        if tmp_ub > b.ub_
            b.ub_             = tmp_ub
            b.best_ub_action_ = tmp_ub_action
        end
    end

    # Sanity check
    if b.ub_ < b.lb_
        # show(particles); println("")
        # error("BOUNDS: ub=$(b.ub_) < lb=$(b.lb_), lba = $(b.best_lb_action_), uba = $(b.best_ub_action_), tmp_lb = $(tmp_lb), tmp_ub = $(tmp_ub), tmp_lb_action = $(tmp_lb_action), tmp_ub_action = $(tmp_ub_action)")
        error("BOUNDS: ub=$(b.ub_) < lb=$(b.lb_)")
    end
    return b.lb_, b.ub_
end

LPDM.best_lb_action(b::LDBounds1d) = isnan(b.best_lb_action_) ? error("best_lb_action undefined. Call bounds() first") : b.best_lb_action_
LPDM.best_ub_action(b::LDBounds1d) = isnan(b.best_ub_action_) ? error("best_ub_action undefined. Call bounds() first") : b.best_ub_action_

function move(p::AbstractLD1, x1::Float64, x2::Float64)#::(Float64,Float64)
    direction = x2 > x1 ? 1 : -1
    all_actions = POMDPs.actions(p)
    # I = findall(!iszero, all_actions)
    pos_actions = all_actions[all_actions .> 0]
    # nz_actions = all_actions[I]
    min_a = minimum(pos_actions)

    r = 0.0
    x = x1
    first_a = NaN
    a_dir = NaN

    if abs(x1) <= p.term_radius
        return reward(p,x1,0.0), 0.0
    end

    if abs(x2-x1) < min_a #too close to take any action
        return reward(p,x1,0.0), 0.0
    end

    while (abs(x2-x) > min_a) && (abs(x) >= p.term_radius)
        a = maximum(pos_actions[pos_actions .<= abs(x2-x)]) # maximum action not exceeding Δx
        a_dir = direction * a
        if isnan(first_a)
            first_a = a_dir # assign first action (directional)
        end
        r += reward(p, x, a_dir) # use current state for computing the reward
        # println("BOUNDS: x=$x, x1=$x1, x2=$x2, min_a = $min_a, all_actions=$all_actions, pos_actions=$pos_actions, av_actions=$(all_actions[all_actions .< abs(x2-x)]), a_dir=$a_dir,  r=$r ")
        x += a_dir # take the step
    end
    # error("review below, might be a bug with assigning first_a = 0.0")
    if abs(x2) <= p.term_radius
        # x1 < 0.7 && println("TERMINATION REWARD #2 FOR x2=$x")
        r += reward(p,x2,0.0) #termination reward
        first_a = 0.0
    end
    # x1 < 0.7 && println("exiting move $x1 -> $x2")
    # if isnan(first_a) #DEBUG: remove
    #     error("move: first_a = $first_a, x1=$x1, x2=$x2")
    # end
    return r, first_a #action sign depends on the direction
end

# NOTE: not for direct calling
# computes the cost of traveling to the low noise region and only then towards target. i.e. a slow approach
function lower_bound(p::AbstractLD1, particle::LPDMParticle{Float64})

    if abs(particle.state) < p.term_radius
        return reward(p, particle.state, 0.0), 0.0
    else
        r1,a1 = move(p, particle.state, p.min_noise_loc)         # estimate cost of moving x to low noise x
        # estimate cost of moving x to the target from the low noise region (will terminate earlier)
        r2,a2 = move(p, p.min_noise_loc, 0.0)
    end
    return r1+r2, a1
end
# lower_bound(p::LightDarkLPDM.LightDark1DDespot, particle::LPDM.LPDMParticle{Float64}) = lower_bound(p, POMDPToolbox.Particle{Float64}(particle.state, particle.weight))

# NOTE: not for direct calling
# computes the reward for the straight-line path to target
function upper_bound(p::AbstractLD1, particle::LPDMParticle{Float64})
    if abs(particle.state) < p.term_radius
        return reward(p, particle.state, 0.0), 0.0
    else
        return move(p, particle.state, 0.0)
    end
end

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
function _readgml(io::IO, line, ::Type{G}) where G
    mapping = Dict{Int,Int}()
    i = 0
    while startswith(line, "node")
        i += 1
        line = readline(io) |> strip
        line == "[" && (line = readline(io) |> strip)
        while !startswith(line, "]")
            name, valstr = splitgml(line)
            if name == "id"
                mapping[parse(Int, valstr)] = i
            end
            line = readline(io) |> strip
        end
        line = readline(io) |> strip #skip ]
    end
    g = G(i)
    while startswith(line, "edge")
        line = readline(io) |> strip
        line == "[" && (line = readline(io) |> strip)
        u = -1
        v = -1
        while !startswith(line, "]")
            name, valstr = splitgml(line)
            if name == "source"
                id = parse(Int, valstr)
                u = mapping[id]
            elseif name == "target"
                id = parse(Int, valstr)
                v = mapping[id]
            end
            line = readline(io) |> strip
        end
        @assert u > 0 && v > 0
        add_edge!(g, u, v)
        line = readline(io) |> strip #skip ]
    end
    return g
end

function readgml(io::IO, ::Type{G}) where G
    H = G
    line = readline(io) |> strip
    @assert startswith(line, "graph")
    line = readline(io) |> strip
    line == "[" && (line = readline(io) |> strip)
    while !startswith(line, "node") && !isempty(line)
        if startswith(line, "directed")
            H = parse(Int, line[10:end]) == 1 ? digraphtype(G) : graphtype(G)
        end
        line = readline(io) |> strip
    end
    return _readgml(io, line, H)
end

function gmltypeof(valstr)
    str = strip(valstr, '\"')
    if length(str) < length(valstr)
        return String
    else
        return Float64
    end
end

gmlval(::Type{T}, x) where {T} = parse(T, x)
gmlval(::Type{String}, x) = strip(x, '\"')
function splitgml(s::AbstractString)
    i = findfirst(isequal(' '), s)
    return SubString(s, 1, i-1), SubString(s, i+1, length(s))
end

gmlprintval(x) = x
gmlprintval(x::String) = "\"" * x * "\""

function _readnetgml(io::IO, line, g)
    mapping = Dict{Int,Int}()
    i = 0
    while startswith(line, "node")
        i += 1
        line = readline(io) |> strip
        line == "[" && (line = readline(io) |> strip)
        @assert startswith(line, "id")
        name, valstr = splitgml(line)
        mapping[parse(Int, valstr)] = i
        add_vertex!(g)
        line = readline(io) |> strip
        while !startswith(line, "]")
            name, valstr = splitgml(line)
            !has_vprop(g, name) && vprop!(g, name, String)
            T = valtype(vprop(g, name))
            vprop(g, name)[i] = gmlval(T, valstr)
            line = readline(io) |> strip
        end

        line = readline(io) |> strip #skip ]
    end

    while startswith(line, "edge")
        line = readline(io) |> strip
        line == "[" && (line = readline(io) |> strip)
        startswith(line, "id") && (line = readline(io) |> strip)

        @assert startswith(line, "source")
        name, valstr = splitgml(line)
        u = mapping[parse(Int, valstr)]
        line = readline(io) |> strip
        @assert startswith(line, "target")
        name, valstr = splitgml(line)
        v = mapping[parse(Int, valstr)]
        line = readline(io) |> strip
        ok, e = add_edge!(g, u, v)

        while !startswith(line, "]")
            name, valstr = splitgml(line)
            !has_eprop(g, name) && eprop!(g, name, gmltypeof(valstr))
            T = valtype(eprop(g, name))
            eprop(g, name)[e] = gmlval(T, valstr)
            line = readline(io) |> strip
        end

        line = readline(io) |> strip #skip ]
    end
    return g
end

function readnetgml(io::IO, ::Type{G}) where G<:ANetOrDiNet
    H = G
    line = readline(io) |> strip
    @assert startswith(line, "graph")
    line = readline(io) |> strip
    line == "[" && (line = readline(io) |> strip)
    gdict = Dict{String, String}()
    while !startswith(line, "node") && !isempty(line)
        if startswith(line, "directed")
            H = parse(Int, line[10:end]) == 1 ? digraphtype(G) : graphtype(G)
        else
            name, valstr = splitgml(line)
            gdict[name] = gmlval(String, valstr)
        end
        line = readline(io) |> strip
    end
    g = H()
    for (name, val) in gdict
        gprop!(g, name, val)
    end
    return _readnetgml(io, line, g)
end

"""
    writegml(f, g)

Writes a graph `g` to a file `f` in the
[GML](https://en.wikipedia.org/wiki/Graph_Modelling_Language) format.
"""
function writegml(io::IO, g::AGraphOrDiGraph)
    println(io, "graph [")
    # length(gname) > 0 && println(io, "label \"$gname\"")
    is_directed(g) && println(io, "\tdirected 1")
    for i=1:nv(g)
        println(io,"\tnode [")
        println(io,"\t\tid $i")
        println(io,"\t]")
    end
    for (s,t) in edges(g)
        println(io,"\tedge [")
        println(io,"\t\tsource $s")
        println(io,"\t\ttarget $t")
        println(io,"\t]")
    end
    println(io, "]")
    return 1
end

function writenetgml(io::IO, g::ANetOrDiNet)
    println(io, "graph [")
    is_directed(g) && println(io, "\tdirected 1")
    for (pname, p) in gprop(g)
        println(io,"\t$pname $(gmlprintval(p))")
    end
    for i=1:nv(g)
        println(io,"\tnode [")
        println(io,"\t\tid $i")
        for (name, val) in vprop(g, i)
            println(io,"\t\t$name $(gmlprintval(val))")
        end
        println(io,"\t]")
    end
    for e in edges(g)
        println(io,"\tedge [")
        println(io,"\t\tsource $(src(e))")
        println(io,"\t\ttarget $(dst(e))")
        for (name, val) in eprop(g, e)
            println(io,"\t\t$name $(gmlprintval(val))")
        end
        println(io,"\t]")
    end
    println(io, "]")
    return 1
end

filemap[:gml] = (readgml, writegml, readnetgml, writenetgml)


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ 9ec6149e-6acd-414c-8902-798771573672
function is_notebook(p)
	startswith(read(p, String), "### A Pluto")
end;

# ╔═╡ bdf53020-635d-4cc9-8a0e-a612ca470e85
hook_link(p) = replace(p, ".jl" => ".html");

# ╔═╡ 7c852ed0-36b8-43ea-8ee8-c4bd45e9c8a9
readme = read(joinpath(@__DIR__, "../README.md"), String) |> Markdown.parse;

# ╔═╡ d974bcaa-935e-4a8c-bfc8-827034e79efe
hook_name(file) = file == "notebook.jl" ? "Other hooks" : "@" * replace(basename(file), ".jl" => "");

# ╔═╡ 11c144eb-60b4-4a64-8077-29e97abedce8
docs = let
	files = map(
		f -> " - [`$(hook_name(f))`](./$(hook_link(f)))",
		filter(
			name -> basename(name) != "index.jl" && endswith(name, ".jl") && is_notebook(name),
			readdir("./")
		)
	) |> s -> join(s, "\n") |> Markdown.parse

	md"""
	#### Documentation

	 $(files)
	"""
end;

# ╔═╡ 3e715447-312c-4c5f-85a3-96abbb91fc0c
sources = let
	files = map(
		f -> " - [`$(hook_name(f))`](../src/$(hook_link(f)))",
		filter(
			path -> basename(path) != "PlutoLinks.jl" && endswith(path, ".jl"),
			readdir("../src")
		)
	) |> s -> join(s, "\n") |> Markdown.parse

	md"""
	#### Source code

	$(files)
	"""
end;

# ╔═╡ 2175e675-df15-4934-b410-d0def1decf04
md"""
$(readme)
$(docs)
$(sources)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─2175e675-df15-4934-b410-d0def1decf04
# ╟─9ec6149e-6acd-414c-8902-798771573672
# ╟─11c144eb-60b4-4a64-8077-29e97abedce8
# ╟─3e715447-312c-4c5f-85a3-96abbb91fc0c
# ╟─bdf53020-635d-4cc9-8a0e-a612ca470e85
# ╟─7c852ed0-36b8-43ea-8ee8-c4bd45e9c8a9
# ╟─d974bcaa-935e-4a8c-bfc8-827034e79efe
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
@doc "G_M_1" ->
G_M_1 = 1

@doc """G_M_2""" ->
G_M_2 = 1

# Constant

@doc "G_C_1" ->
const G_C_1 = 1

@doc """G_C_2""" ->
const G_C_2 = 1


text:
# Note all the graphical stuff and IBM specific stuff has been moved into a different file
# marginalDrawing.jl which should be included if you want to use it.

# Copyright Robin Harper 2016-2020

using LsqFit, Hadamard, DelimitedFiles, LinearAlgebra
using Statistics
import Base.show

export readInCSVFile,transformToFidelity,fitTheFidelities,convertAndProject,gibbsRandomField
export getGrainedP,mutualInformation,relativeEntropy,conditionalMutualInfo,covarianceMatrix,JSD
export correlationMatrix

⊗ = kron

"""
readInCSVFile(filename)

## Arguments
    `Filename`: full name of a csv delimited Filename

## File Format
Assumes that each sequence appears on a new line and each line has the same number of entries.
Assumes that it is a raw count for each of the possible measurement outcomes.
Note it is important to know the order of the outcomes (this is significant).
Although this can be overriden, it is a assumed that they are stored as:

```
    00000000
    00000001
    00000010
    00000011

    ...
    11111110
    11111111

```
i.e. lsb to the right.

Returns a full rows x col matrix depending on the number of entries of the file.
for a 16 qubit machine, where the data was taken over 6 sequences we will return a 6x655536 array of Ints.

"""
function readInCSVFile(filename)
    return readdlm(filename,',',Int64)
end

"""
transformToFidelity(data)

## Arguments:
    `Data`, as returned by readInCSVFile - see help for that function.

## Returns:
The Hadamard transformed fidelities for each sequence as a list of lists.

"""
function transformToFidelity(data)
    splitMatrix = [data[i,:]/sum(data[i,:]) for i in 1:size(data,1)]
    return [ifwht_natural(x) for x in splitMatrix]
end

modelF(x, p) = p[1]*(p[2].^x)

"""
fitTheFidelties(lengths, data; no =0)

## Arguments

   `lengths: Array{Int64,1}`, e.g. collect(1:2:24) - the length at which each of the above sets of data was gathered.

   `data: Array{Array{Float64,1},1}`, this is a list of observed probabilities at different sequences.

   `no:` printed with the warning message if we fail to fit. Allowing identification of problem in "batch" fits.

Take the data, transform it  and see if we can fit it to an exponential decay.
This proves problematic where the data is a bit sparse, as the decays can be all over the place.
In general it is *extremely* helpful to have good 'low' gate sequence as that helps to anchor the fits. If we have an extremely high decay rate, the intercept parameter is difficult to find.

Here we check for convergence, and warn if it doesn't converge.

This enforces a cut-off to the data used for the fit. Specifically it gets rid of the 'tail' of data.

## Returns

three things:
-   the fitting parameters, in the form [A,p], where A is the y-axis intercept and p is the decay probability.
-   the length of the sequence that was used for each parameter.
-   the indices of those that still failed to converge.

## Typical usage
```
    params,_ = fitTheFidelities(collect(1:2:22),transformToFidelity(data))
```
"""
function fitTheFidelities(lengths,matrix;no=0)
    revh = [ifwht_natural(x) for x in matrix]
    params = Array{Float64,1}[]
    dataCounted =[]
    failedToC = []
    _negCount =0
    for idx = 2:length(revh[1])
        extract = map(x->x[idx],revh)
        # The idea here is as per https://arxiv.org/abs/1901.00535 i.e. don't let the tail wag.
        # It doesn't make a large difference, but it does make a difference.
        p1 = (extract[1]*(1+1/16))/4
        lastData = findfirst(x->p1>x,extract)
        if lastData == nothing
            lastData = length(extract)
        end
        # two is probably sufficient, but at 3 for now.
        if lastData < 3
            lastData = 3
        end
        fit =[]
        try
            fit = curve_fit(modelF,lengths[1:lastData],extract[1:lastData],[0.8,0.8],upper=[1.0,1.0],lower=[0.01,0.01])
        catch
            print("Fall back on $no, $idx\n")
            fit = curve_fit(modelF,lengths[1:lastData],extract[1:lastData],[0.8,0.4],upper=[1.0,1.0],lower=[0.01,0.01])

        end
        if !fit.converged
            push!(failedToC,idx)
        end
        push!(params,fit.param)
        push!(dataCounted,lastData)

    end
    if (_negCount > 0) print("$no: There were $(_negCount) with a negative first number.\n") end
    if (length(failedToC) > 0) print("$no: $(length(failedToC)) failed to converge.\n") end
    return params,dataCounted,failedToC
end

"""
projectSimplex(probs)

## Arguments

   probs: a probabiliity distribution, that isn't quite

Takes the proffered probability distribution and projects it onto the closest probability distribution (where they are all >=0). Uses an algorithm similar to
[Projection onto the probability simplex:
An efficient algorithm with a simple proof, and an application)[https://arxiv.org/abs/1309.1541]

## Returns:
   the projected probability distrubtion.
"""
function projectSimplex(ps)
    sortedSeq = reverse(sort(ps))
    # This is a bit crytic.
    #   a) Sort (greatest first)
    #   b) create a culmlative sum (cumsum)
    #   c) we want the last one where (Cum[k]-1)/k < sorted[k]
    #   d) we then reduce all probs by the last value we had in c, setting them to zero if they are negative.
    # (The reason why we divide by the length in c is because we subtract it more times the further into the index we are)
    τ=filter(x->x[1]<x[2],collect(zip((cumsum(sortedSeq).-1)./(1:length(ps)),sortedSeq)))[end][1]
    projected = map(x->max(x-τ,0),ps)
    return projected
end

"""
ConvertAndProject(params)

## Arguments
    params: parameters that come from the fittingTheFidelities method (list of A and p)

    Convert and project the fitted fidelities, to rerieve the 'p's representing the joint probabilities.
## Returns:
   the projected probability distribuion.
"""
function convertAndProject(params)
    ps = fwht_natural(map(x->x[2],params))
    pps = projectSimplex(ps)
    if !(isapprox(sum(pps),1))
    @warn "There may be a problem, the probabilities do not add to 1."
    end
    return pps
end

"""
    Returns the tensor indices.
"""
function getIndices(qs;dimension=16)
    sqs = sort(qs)
    indices = [2^(sqs[1]-1),2]
    for i in 2:length(sqs)
            push!(indices,2^(sqs[i]-sqs[i-1]-1))
            push!(indices,2)
    end
    push!(indices,2^(dimension-sqs[end]))
    return indices
end

"""
marginalise(q,pps)

## Arguments
   q the qubits you want to marginalise to, expressed as a list e.g. [1,2] marginalise away everything else
   pps the probability distribution you want.

   Example
   ```
      marginalise([1,3,5],pps)
   ```
   returns the joint probability distribution of 1,3 and 5 with all other qubits marginalised out.

## Returns
    the marginalised joint probability distribution. You may need to vec this.
"""
function marginalise(q,pps)
    # Get indices sorts the entries
    dimension = 0
    try
        dimension = Integer(log(2,length(pps)))
    catch
        @warn "Error, the size of pps needs to be an integer power of 2"
        return
    end
    if 0 in q
        print("Please index off 1, a 0 in the bits to marginalise over will cause grief")
        return nothing
    end
    x = reshape(pps,tuple((getIndices(q,dimension=dimension))...))
    # Below the permutation that will get us back.
    sorted = sort(q)
    # This is the reverse of sortperm is from sorted -> original
    permute = [findfirst(isequal(x),sorted) for x in q]

    # print(permute)
    # Sum over all the odd ones.
    marginalised = sum(x,dims=1:2:length(size(x)))
    # Work out how many 2 variables we have
    indices = length(q)
    # Shove them in a tuple (note splat so it works well with reshape)
    fullIndices = tuple([2 for i in 1:indices]...)
    # Might as well get it in a 2 dim array to return
    # Rows has to be even of course.
    rows = 2^(floor(Int,indices/2))
    if rows == 0
        rows = 1
    end
    cols = round.(Int,2^indices/rows)
    return reshape(permutedims(reshape(marginalised,fullIndices),permute),rows,cols)
end

"""
gibbsRandomField(pps,constraints)

## Arguments
-   `pps: Array{Float64,1}` For example the output of convertAndProject
-   `constraints: Array{Int64,1}`: The division of the gibbs field. See below for example.

    Takes a joint probability and the gibbs variable constraints and returns a vector of reduced probability distributions.
    For example the constraints might be [[1,2,3,4],[3,4,5,6],[5,6,7,8]] over a field of 8 qubits
    This would return an array of 3 ϕ's each 16 long (2^4), from which the joint probability can be extracted on the assumption
    That, say, qubits 1,2 are independent of qubits 5-8.
    Note: the assumption inbuilt here is that the constraints end and start with qubits in common e.g. [1,2,->3,4],[3,4<-,5,6]

Note the other version of is one gibbsRandomField, where the constraints are constraints::Array{Array{Array{Any,1},1},1}
is probably (now) preferred - there you can specify the constraints more generally.
In particular as tuples of qubits and the qubits they are conditioned on.

## Returns
   ϕ as a series of factors, depending on the constraints.
"""
function gibbsRandomField(pps,constraints)
    overlaps = []
    for (idx,i) in enumerate(constraints[1:end-1])
        x = findfirst(x->x==i[end],constraints[idx+1])
        if x == nothing
            print("Error in finding overlaps the constraints need to overlap at least slightly\n")
            return nothing
        end
        push!(overlaps,x)
    end
    pxx = [marginalise(x,pps) for x in constraints]
    for i in 1:length(overlaps)
        pxx[i] = reshape(pxx[i],:,2^overlaps[i])
    end # reshape to correspand to the overlap.
    px =  [vec(marginalise(x,pps))' for x in [i[end-(overlaps[idx]-1):end] for (idx,i) in enumerate(constraints[1:end-1])]]
    # Map any zero entries in px to something non-zero or division will fail.
    # Not an issue as if px is zero, the corresponding pxx has to be zero. (and in this case by definition the probability is zero)
    pxNz = [map(x-> x == 0 ? 1e-8 : x,pxidx) for pxidx in px]
    ϕ = [vec(pxx[i]./pxNz[i]) for i in 1:(length(constraints)-1)]
    push!(ϕ, vec(marginalise(constraints[end],pps)))
    return ϕ
end

"""
    Helper function for the second form of gibbsRandomField. Takes a list of probability tuples and flattens it out.
    So [(1,2),(3,4)],[(5,),(6,7)]] becomes [[1,2,3,4],[5,6,7]] (read the first as p(1,2|3,4)p(5|6,7))
"""
function jointIfy(x)
    return [collect(Iterators.Flatten(y)) for y in x]
end

"""
gibbsRandomField(pps,gen_constraints)

## Arguments
-   `pps: Array{Float64,1}` For example the output of convertAndProject
-   `constraints: Array{Array{Array{Any,1},1},1}`: The division of the gibbs field. See below for example.

    Takes a joint probability and the variable constraints that are conditioned and returns a vector of reduced probability distributions.
    It is up to the user to make sure that the factorization is complete and makes sense.
    For example the constraints might be [[(1,2),(3,4)],[(3,),(4,5,6)],[(4,5,6,7,8),()]] over a field of 8 qubits
    This would return an array of 3 ϕ's the first 16 long (2^4), the second 2^4 long and the third 2^5 long from which the joint probability can be extracted on the assumption
    That, say, qubits 1,2 are independent of qubits 5-8.

## Returns
   ϕ as a series of factors, depending on the constraints.

"""
function gibbsRandomField(pps,constraints::Array{Array{Array{Any,1},1},1})
 # This is more flexible than the previous one, which took overlapping contraints, to try and work out
 # what the qubits were conditioned on, e.g. if give [[1,2,3,4],[3,4,5,6]], it would assume that you wanted:
 # p(1,2|3,4)p(3,4,5,6)

# Here we can expressly state the dependancies.
# It is expecting contraints as follows a list of [ [[bits],[conditioned on]],[[bits],[conditioned on]]]
# eg [[[1,2],[3,5]],[[3,4],[5]],[[5],[]]]
# means that we want p(1,2|3,5)p(3,4|5)p(5)
# Which will be returned as a 2^4, 2^2, 2^1 set of arrays.

    overlaps = [length(i[2]) for i in constraints]
    jointMarginals = jointIfy(constraints)
    pxx = [marginalise(x,pps) for x in jointMarginals]
    for i in 1:length(overlaps)
        pxx[i] = reshape(pxx[i],:,2^overlaps[i])
    end # reshape to correspand to the overlap.

    px =  [x == [] ? 1 : vec(marginalise([x...],pps))' for x in [i[2] for i in constraints]]
    # Map any zero entries in px to something non-zero or division will fail.
    # Not an issue as if px is zero, the corresponding pxx has to be zero. (and in this case by definition the probability is zero)
    pxNz = [map(x-> x == 0 ? 1e-8 : x,pxidx) for pxidx in px]
    ϕ = [vec(pxx[i]./pxNz[i]) for i in 1:(length(constraints))]

    # If we did just allow them to be zero, then it causes problems reconstructing the probability distribution.
    # Here we find where the reshaped vector e.g. a conditioned on b,c,d
    # is equal to zero (it should be 1, i.e.the culmlative percentage chance of a being 0 or 1 must be 1.)
    # We set the chance of it being 0 to be the average of the non 'all zero' chances. Same with 1  etc.
    # Generalised where we have a joint distribution on the left hand size (reshapeFactor)
    for (cidx,toC) in enumerate(ϕ)
        reshapeFactor = 2^length(constraints[cidx][1])
        for (idx,p) in enumerate(sum(reshape(toC,reshapeFactor,:),dims=1))
            if p==0
                newVals = zeros(reshapeFactor)
                #display(reshape(toC,reshapeFactor,:))
                for i in 1:reshapeFactor
                    newVals[i] =  (mean(reshape(toC,reshapeFactor,:)[i,1:idx-1]))
                end
                for i in 1:reshapeFactor
                    reshape(ϕ[cidx],reshapeFactor,:)[i,idx]=newVals[i]
                end
            end
        end
    end
    return ϕ
end

# The function below is probably worth a bit of comment.
# uses the information in 'graining' to tells us which index the pattern we want resides in
# For instance say we want all bits zero apart from qubit 3 (which = 1)
# Then for ϕ[1], which represents bits [2,1,3,6] - we need the entry corresponding to decimal equiv of reverse(0,0,1,0)
# plus one, because 0 0 0 0 = first entry i.e. index 1.
# We will have length(ϕ) of these entries, this gives us the perVecIndex.
# Then we just read it out for ϕ, with a map, using foldl to multiply them for us.

"""
getGrainedP(ϕ,tomatch,graining)
## Arguments
    ϕ: the factor graph (possibly obtained through gibbsRandomField)
    tomatch: the bit pattern we want
    graining: the split for the factor graph.

## Description
   For instance say we want all bits zero apart from qubit 3 (which = 1)
   Then for ϕ[1], which represents bits [2,1,3,6] - we need the entry corresponding to decimal equiv of reverse(0,0,1,0)
   plus one, because 0 0 0 0 = first entry i.e. index 1.
   We will have length(ϕ) of these entries, this gives us the perVecIndex.
   Then we just read it out for ϕ, with a map, using foldl to multiply them for us.

## Example use
```
   generalisedConstraints =[
                        [[1,],[2,14]],
                        [[2,14],[13,3]],
                        [[13,3],[4,12]],
                        [[4,12],[5,11]],
                        [[5,11],[6,10]],
                        [[6,10],[7,9]],
                        [[7,9,8],[]]]
    ϕ=Marginal.gibbsRandomField(singlePps,generalisedConstraints)
    reconstructed = [Marginal.getGrainedP(ϕ,tomatch,[vcat(x[1],x[2]) for x in generalisedConstraints]) for tomatch =0:(2^14-1)]

```
    Takes a global probability distribution, uses it to fill in the factor graph variables spedified in generalisedConstraints
    and the reconstructs the probability distribution as given by the factor graph.

## Returns:
     for the given tomatch bitpattern the extracted probability.
"""
function getGrainedP(ϕ,tomatch,graining)
    vectorIndex = [[x&tomatch>0 ? 1 : 0 for x in (map(x->2^(x-1),i))] for i in graining]
    perVecIndex = map(x->parse(Int64,join(reverse(x)),base=2)+1,vectorIndex)
    return foldl(*,map(x->ϕ[x][perVecIndex[x]],1:length(ϕ)))
end

"""
mutualInformation(p1,p2,p)

## Arguments
   p1 a qubit position
   p2 a qubit position
   p a global probability distribution.

## Reutrns
   the mutualInformation between probabilities at p1, p2 given probability distribution p.

## Comment
    Makes use of the marginalise function to marginalise over probabilities

    \$MI(p1,p2) =\\sum\\limit_{x\\in p1,y\\in p2}\\left[p(x,y)\\log{\\frac{p(x,y)}{p(x),p(y)})\\right]\$

    Some gotchas:
        I return -1 if p1==p2 , this allows an identification of the qubit where I loop.
        if p(x,y) is 0, p(x) or p(y) is 0 then this is defined as a 'zero' part of the sum - even though the log is undefined.

"""
function mutualInformation(p1,p2,p)
    # Allow this for easy looping and identificaton of controlling qubit.
    #print("p1 in $p1\n")
    #print("p2 in $p2\n")
    for x in [p1...]
        if x in [p2...]
            return -1
        end
    end
    p12 = marginalise([p1...,p2...],p)
    #println(p12)
    p1 = sum(p12,dims=2)
    #print("p1 $p1\n")
    p2 = sum(p12,dims=1)
    #print("p2 $p2\n")
    logPart = (log2.(vec(p12)) - log2.(vec(p1⊗p2)))
    # Note that the log of 0 will give -Inf or NAN (if p1 or p2 is zero as well).
    # The log can only be zero (or NAN) if the relevant
    # p12 entry is 0, which by definition means the MI contribution is 0, but of course -Inf*0 = Nan
    # So we zap the -Inf's (and Nan's) to zero which gives us the result we want.
    toReturn = vec(p12)'*map(x->x==-Inf||isnan(x) ? 0 : x,log2.(vec(p12)) - log2.(vec(p1⊗p2)))
    return toReturn
end

"""
    relativeEntropy(P, Q)

    Calculates the relative entropy between two joint probability distributions.
    D(P||Q) = Sum p_j*log(p_j/q_j)

    Undefined if any of the q_js are zero, unless the corresponding p_j is zero.
"""
function relativeEntropy(P,Q)
    # First of all assert that where p̃ is zero pps is also zero
    if !(length(filter(x->x!=0,P.*Q)) == length(filter(x->x!=0,P)))
            @warn "D(P||Q) - the distribution (Q) has zero values not matched by by zero values in (P) - this is undefined."
            return NaN
    end

    # We only need the values of P and Q where P is non zero
    nonZero = vec(map(x->x!=0,P))
    xp = P[nonZero]
    xq = Q[nonZero]
    pDq = xp./xq
    return xp'*log2.(pDq)
end

"""
JSD(dist1,dist2)

## Arguments
    dist1, probability distribuion
    dist2, probability distribution

## Discussion
   calcluates the Jensen-Shannon Divergence between the two distributions. This is symmetric and well defined, even if the probabilities are not in each other's support. The square root is a metric.

   \$JSD(P||Q) = \\frac{1}{2}D(P||M)+\\frac{1}{2}D(Q||M)\$, where

   \$M =\\frac{1}{2}(P+Q)\$

   Returns \$JSD(P||Q)\$
"""
function JSD(dist1,dist2)
    m = 0.5.*(dist1.+dist2)
    return 0.5*(relativeEntropy(dist1,m)+relativeEntropy(dist2,m))
end

"""
reconstructedJS
## Arguments
    distribution: a probability distribution
    constraints: a set of constraints (see first form of gibbsRandomField)
## Discusion
    Given a distribution and a set of constraints.
    Reconstruct the distribution using the constraints.
## Return
    The JSD of the reconstructed distribution and the original distribution.
"""
function reconstructedJS(distribution,constraints)
    return JSD(reconstruct(distribution,constraints),distribution)
end

"""
    Given a distribution and a set of constraints.
    Return the reconstruction of a distribution parameterized by the constraints.
"""
function reconstruct(dist,constraints)
    gibbs2ϕ = gibbsRandomField(dist,constraints);
    reconstructedPps2 = [getGrainedP(gibbs2ϕ,tomatch,constraints) for tomatch=0:(2^14-1)];
    @assert(isapprox(sum(reconstructedPps2),1))
    return reconstructedPps2
end

"""
    Note to self this will often be undefined for the types of distributions I am currently looking at
"""
function relEntropy(dist,constraints)
    return relativeEntropy(reconstruct(dist,constraints),dist)
end

"""
    Helper function that given an xyz correctly indexes and pulls out the
    xyz from the joint probability distribution, jpXYZ
    and the xz and yz from jpXZ and jpYZ
    then does the calculation, inside the sum for that x,y,z
"""
function getSummand(x,y,z,jpXYZ,jpXZ,jpYZ,jpZ)
    multipland = size(jpXZ)[1]
    pxyz = jpXYZ[multipland*(y-1)+x,z]
    if pxyz == 0
        return 0
    end
    pz = jpZ[z]
    pxz = jpXZ[x,z]
    pyz = jpYZ[y,z]
    return pxyz*log2((pz*pxyz)/(pxz*pyz))
end

"""
conditionalMutualInfo
    Pass in vectors of qubits for X Y and Z and the distribution to analyse.

## Discussion
    Gives the conditional mutual information for the qubits
    supplied in X,Y and Z where we want I(XY|Z)

    \$I(XY|Z) = \\Sum( P(x,y,z)\\log(2, P(Z)P(X,Y,Z)/P(X,Z)P(Y,Z) )\$

## Returns
    \$I(XY|Z)\$

"""
function conditionalMutualInfo(X,Y,Z,p)
    XYZ = vcat(X,Y,Z)
    XZ = vcat(X,Z)
    YZ = vcat(Y,Z)
    xSize = 2^length(X)
    ySize = 2^length(Y)
    zSize = 2^length(Z)
    jpXYZ = reshape(marginalise(XYZ,p),:,zSize)
    jpXZ = reshape(marginalise(XZ,p),:,zSize)
    jpYZ = reshape(marginalise(YZ,p),:,zSize)
    jpZ = reshape(marginalise(Z,p),1,zSize)
    return sum([getSummand(x,y,z,jpXYZ,jpXZ,jpYZ,jpZ) for x=1:xSize,y=1:ySize,z=1:zSize])
end

"""
    covarianceMatrix(p)

## Arguments
    p a probability vector

## Discussion
    Computes the covariance matrix between the 0/1 random variables representing no error / error.
    If x is a column vector of bits representing an error pattern, then we can compute the matrix
    Expect_p[(x-μ) (x-μ)^T]
    where μ = Expect_p[x].

    If we set reverse to true, then we reverse the digits.
    The way IBM stores its digits we wouldn't normally want to reverse them, but in general
    this might not be true
## Returns
   the matrix Expect_p[(x-μ) (x-μ)^T]
"""
function covarianceMatrix(p;reverseDigits = false)
# Take a probability distribution on n qubits and compute the covariance matrix as above.
    n = convert(Int,log2(length(p)));
    v = digits.(0:(2^n-1),base=2, pad=n) # the 0/1 random variable of errors
    if reverseDigits
        v = reverse.(v)
    end
    v = convert(Array{Array{Int,1},1}, v);
    mu= sum(p.*v); # compute the mean
    v = hcat(v...).-mu; # center the random variable v
    covmat = p[1]*v[:,1]*v[:,1]';
    for i in 2:2^n
        covmat+=p[i]*v[:,i]*v[:,i]';
    end
    return covmat
end

"""
    marginaliseFromRawData

    rawData - all of the observations ie. an array of 1..|lengths] each being a \$2^n\$ vector.
    constraints - currently a list of constraints in the same from as [`gibbsRandomField(pps,constraints)`][@ref]
    lengths - the lengths the data sets were taken at

    Takes in the raw observations, for each of the lengths specified in length.
    Extracts the joint probabilities needed from constraints.
    Marginalises the observations, transforms, fits, transforms and fills in the joint probabilities.
    Returns the completed gibbs factors.

    Note if your \$2^n\$  vectors are too big, simply re-write to pass in the actual observations and work out the marginalised observations as you scan through them,
    potentially sparse arrays will help (although not sure if marginalise will work)

"""
function marginaliseFromRawData(rawData,constraints,lengths)
    overlaps = []
    for (idx,i) in enumerate(constraints[1:end-1])
        x = findfirst(x->x==i[end],constraints[idx+1])
        if x == nothing
            print("Error in finding overlaps the constraints need to overlap at least slightly\n")
            return nothing
        end
        push!(overlaps,x)
    end
    pxx=[]
    # Marginalise over measurements to get probabilities
    for x in constraints
        marginalMatrix = [marginalise(x,ms) for ms in rawData]
        paramsM,_ = fitTheFidelities(lengths,marginalMatrix)
        pm =  fwht_natural(vcat([1],map(x->x[2],paramsM)))
        push!(pxx,vec(projectSimplex(pm))');
    end
    for i in 1:length(overlaps)
        pxx[i] = reshape(pxx[i],:,2^overlaps[i])
    end # reshape to correspand to the overlap.
    px=[]
    # Here we can use our previous marginalisation to consistently marginalise the constraints
    for (idx,x) in enumerate(constraints[1:end-1])
        lx=convert(Int,log2(length(pxx[idx])))
        toM = [y for y in lx-overlaps[idx]+1:lx]
        push!(px,vec(marginalise(toM,pxx[idx]))');
    end
    # Map any zero entries in px to something non-zero or division will fail.
    # Not an issue as if px is zero, the corresponding pxx has to be zero. (and in this case by definition the probability is zero)
    pxNz = [map(x-> x == 0 ? 1e-8 : x,pxidx) for pxidx in px]
    ϕ = [vec(pxx[i]./pxNz[i]) for i in 1:(length(constraints)-1)]
    push!(ϕ, vec(pxx[end]));
    return ϕ
end

"""
    correlationMatrix(p)

## Correlation

    Computes the covariance matrix \$\\Sigma\$ and then let:
    \$D = \\sqrt{diag{\\Sigma}\$

    Return \$D^{-1}\\Sigma D^{-1}\$
    See also See also: [`covarianceMatrix`](@ref)
"""
function correlationMatrix(p,reverseDigits=false)
    M = covarianceMatrix(p,reverseDigits=reverseDigits)
    d = sqrt(inv(LinearAlgebra.Diagonal(M)));
    return d*M*d
end


text:
# DAY 3 PROBLEM 2 ADVENT OF CODE
# Get the life support rating from oxygen generator rating and CO2 scrubber rating.

diags = """011110011110
101101001111
000000010101
100111001010
110000011010
011101010101
011110001100
010111111001
111011100101
011110000110
010100001111
000001001100
110001001000
001000100101
010100011110
010111011110
110101001111
011010111011
000011010101
100111011100
010010100011
010010011001
000111010011
111001000111
100101111111
011001101110
110010101100
001000010001
110101100011
101110001011
010100000110
001110011110
011101100110
001111101101
000101101000
100110000111
000001010100
111100010110
000100100100
001000110001
101000001011
101000100110
110000101010
010111110000
111110010110
011011111111
111010100000
100001000111
011000001010
011000010001
101101100001
011001101011
111001001010
000010110101
101110100010
001111111101
011110001001
000111100011
010101001001
111110111111
111000111010
001001111100
000011000001
011100010110
110010000111
000010010101
110000001110
010011111011
001010100100
111001010101
001101000011
110100110011
110010100101
001000001011
111110011000
110111011011
110110101001
010010100100
111100001100
100100001100
110001111010
111011001011
010000010001
110010100000
010100011100
000110001001
010011000001
111011010001
010011111101
111101110100
110110101100
010100111110
100011001101
000010000010
011111100001
100000101100
100001110011
100101100100
000100001001
101000000011
110111101100
100100101001
001001101001
000001010000
001001111110
111011111010
001000011110
110000100111
010010010101
011010001011
100010010101
010110101010
011111001011
110010011100
111110101010
000001111010
111101110101
101010110111
100010101111
000110101110
000000010000
000010100110
111110010001
000100110000
110010011011
101001000010
010110000101
100011100100
110000110100
110011000001
101110001000
010100011111
111000100111
110001100010
100001101000
101001000111
111110100000
000011000010
110000000001
111010000111
100100111011
101101110001
111100110100
100110010000
001001001110
100010001001
101000110100
011100001111
001110100110
110100001100
101101001011
000011111110
101101010100
000110110001
011100011001
111000111111
011010101100
111000100101
100111011110
011010100110
110000101111
110010000001
100001110101
110101011111
111111000011
101010010010
000111111111
000001001111
101101010010
111111001000
001001110000
110010011110
110011001111
100111101011
011111000111
110010111100
010110000000
110001100011
111000100001
011000101000
100110001100
001101011011
011001111110
011101000110
110110111001
110101001010
111000111011
100111010001
100110000011
010001000100
001001001000
010010001010
010101101101
001110001001
101110011111
010111001111
110111001110
100000101011
010000111111
101000101011
101011010111
000101110001
100100010010
100111011101
111100000100
100110000001
000001000101
101100001111
100000101001
010101010110
000110010111
011101001100
000111111011
000001100101
100001011011
111000101110
111101100111
011010010100
101001001010
110100110001
111101011101
000011101110
011000100100
101110000111
011101110110
111001010001
011110000100
011110001101
100111100001
010011100001
001001101101
001111111010
010101000101
011100000110
100010001101
011011101011
110111001010
101000010010
111111100011
000101111011
001101010011
100101000001
101101100011
101110010011
011001100000
110001111111
011100101100
010010001110
111110101101
010110001000
000011000110
101010011011
000100111011
001000111101
010011111111
011101100111
111111011111
001100000000
101100001110
101100001000
100001011111
100001000101
000101011001
010100111100
110000111011
010000111100
011100110000
101111110010
100011110101
100100110101
100001100000
101000100101
111000001100
010000001011
110101010111
011011001000
100111000101
011000000100
110110011000
101101101101
100101110010
001100001010
010111000011
111111100001
100010111011
110011001001
000011011101
001100101101
010110011001
111100110010
011111001110
100110001111
000111101111
011000110001
001010001010
000111010100
010011110101
011111000010
000101000010
111101011010
110011111100
100000001000
100001000110
101010100001
000101000101
100100011111
110000101110
000110000101
001001011111
100010011100
000001110001
110000010100
011111010100
001000001111
000111000001
110010001000
001001111101
111000001001
001011000001
111110010100
110010001100
010110110000
011100100010
011000101011
001001000001
011001001110
111000010100
000110000110
011011000100
011000010111
110001011000
010111011000
000101110100
000000011111
001100111100
000010100101
000111110000
000110000001
000100101001
101111110100
010010000001
001100010001
101111001101
111011111101
000011110011
111001011100
010110110110
000000001101
010010000000
000001010011
101000010000
101000001110
010001000111
111010011011
101001111110
000110101100
100110101001
100111111111
101001101001
110010010011
011001011111
111111110110
101110010110
100100001101
010001111110
111111101000
100010100001
011110010100
101000110011
010110110011
011111001000
100010100101
110010010110
111011001000
101110001101
000010110011
001110100100
100000011101
011001100001
011001101010
001100011101
000111101100
001010001001
111001101110
110110001111
001111001010
001001000011
001000010000
000001111110
000101101111
001001000110
100001101110
101111011100
111110010010
010001100011
101000100111
001101000111
111011110011
000011000101
111101001110
110111111111
011100010100
110011110001
001011011001
000011010111
011000001110
000101010111
110101110001
001100001111
110111010111
110010010101
100111100111
101000000111
101111000100
100011100010
000100100101
111011110110
001111100110
011011000010
100010010100
010111110010
100010110011
111111111000
011010101111
101110111011
101110101010
001011001000
110110001011
110100101111
001011001011
110011010000
000000100000
100000110000
101110011011
100101101111
010110100010
000010100011
000110110000
111101000001
011001010011
000000001110
000100110101
111111111110
001010101001
001101100111
010100000101
000111001001
101100101110
011011110111
011010000011
000100101011
100001100110
000111011101
000100011100
000101100101
011111111110
111010111011
100000100101
011001011110
010101100010
010101111101
100001000001
001000100001
010010111110
100011001011
011011000111
110011010101
111101100100
001110001010
110011111101
111011011101
001101000001
101011011011
100000111010
010101110101
000101011100
011000011001
111001010011
100100010110
111101010001
110110011011
011110110100
101000110101
101110111101
011011011110
000100110110
100011011111
111010101001
011111011000
101000101110
110010111001
111010000010
110011011001
110110010100
111010011010
000010110001
011101100010
110010110111
010101001111
011011001101
101111000110
010111010010
101111011011
100111101100
001000001101
010110100000
001010100110
111111000111
100101000100
100101001100
110111010110
101000101000
100001001111
110011110010
111111101111
110001111000
010001110110
101100000010
000110000111
001011100101
101010100101
110110001001
000011011111
110110101110
101100110010
101110101110
101010110101
000111000011
101000011111
001100000111
101001111111
100100011101
111100011011
000010111010
111010101010
101001010010
101100010011
010000110110
011110110011
001101110111
110111110110
001011000100
101111100010
011010011100
110100001010
100100110111
000011010011
000110111010
001100100000
100111011011
001101101011
101011101111
001101100100
111010010111
100111111100
101101000100
000011111100
011111001101
110101001101
011100100110
101011000000
100001010011
011110111110
011011100110
110001011010
111111110101
001111111000
100010011011
100001101100
110011010100
100011000001
011110111001
010111110110
100110001010
100000000011
001001101111
101111011010
010000110100
101001011110
100100011000
101101000111
011101001001
010101110010
000110000000
000100111000
010110111010
110000011111
011001011100
010100110110
011111110011
111011011010
001000101011
000010110100
011100101000
111111100110
011001110100
000100011000
010101010111
001011101000
001000011100
110101101100
010001000010
001010010000
111000110001
101111111111
011101011110
001001111010
001000100110
000011011001
100101100001
110111000000
100111011111
111000110000
110010111111
000010010011
100001100111
010110100100
001101001010
101001100110
010001110101
001011000010
010111110101
000001001110
010010100110
111011001110
101011111001
000100000000
111101001100
100111011010
000010010001
100000010000
101101000010
110100010100
111011111000
101001100100
001110101011
000111001010
111110000011
001110100101
101010000101
010101010011
110000111100
101011001010
010001011011
010000011010
000010101011
100110110101
100111011001
100111010111
000101100000
101110000001
100100001000
010111000010
000101001010
001110010010
111011000111
001110110010
111010000001
011000100101
000111100101
100000100111
100110011011
001011100001
101010100111
011001111000
000011011010
001011101100
010011011001
000111111100
010001011110
111111011010
111010010101
110000010011
001001101010
000101011011
100111111110
110010101001
111000011000
011100101110
011100010001
001101110010
101000111111
101011011001
110111000110
001000010100
000111000100
010000101110
001100100001
111100000011
111011101010
000001111000
100100010000
110110000110
000010111100
110110110111
010000011011
101110000100
011100100100
010111001011
001011010000
000111010001
110010111000
001011011110
000001100111
010000101111
110001100100
111000001101
000110011001
010000010110
101010001110
111000111000
011111011100
111001110110
111100010001
110101111001
001001010011
010001111001
001101101110
010011000010
000100011011
000001101000
001110100001
100011010011
100000011100
011001010111
110011011010
100110001011
011001110111
001010001100
000111010111
110111111101
111001101010
001110001000
011100001010
010011101010
111011010011
000011001010
100101000110
100101010111
110010100010
010100101111
111011011100
100011100101
101111111101
001001110011
110111101001
100100010101
010101100001
100011000011
001010010111
011001001001
110001101001
010011100010
111000001111
010101111111
001001001011
000111110111
100110111101
011011110000
001001011110
011010101101
111011101011
010001100111
110100001111
100010110111
000000010100
100000010001
010011010101
010011100110
100011111011
110110100001
110000100110
110101100000
110001101011
110110100110
101111101110
011010101011
101001010101
110110111101
100111100010
101100100010
101011101101
010011010011
001111011100
100101010100
101010101010
010000011101
011100000010
100111101000
111111010101
011010101000
001010110010
111110010011
110111011010
010010111111
111001100111
010101111000
000111110011
111101000110
001011000000
110100100110
101001110010
101111011001
111110000001
011101110010
011000010110
101011010010
100110011101
000010011101
001111110100
110011001100
111110111100
011001011000
001110000011
110001001101
111110111011
110000000000
010000110001
000100110001
001000101101
001100110111
101011001111
110011110011
010011100101
001001000000
011110100000
001010110111
110101111101
001111000000
100011011001
100000001001
000010101010
110000100000
111111101001
111010101100
101100110111
110100011000
010011000111
010010001000
001011001111
001001001100
100010110110
111100111011
111100101100
000111111101
101110000011
100001011101
111101110011
100100000011
111110011001
100010111101
110001001010
111111100100
010000110101
001100001101
001001110111
110110101010
011110010010
101011010000
101101110010
010010110000
011110100111
010111000100
010100000010
000011110100
111001011101
010111100101
110010101101
011100010010
001110111100
011111110110
001110100010
000001100001
011011110010
000010111101
001011111110
111110111010
001010101101
100011110010
001111110010
011001110010
100110101100
010111010011
001000010011
111110110010
011000000110
011111111100
000111110110
110001110100
101110010101
010000111101
011011101110
011010100111
001010000010
100001100011
110001000000
101111100111
101111110011
001101110011
110001110101
100101010010
010101011110
000011111000
000100001000
101001110011
011100110001
000000011100
111110010111
101011101100
110011101100
110101011000
010111001000
111100001010
001000011111
010100111011
101010101011
001110001110
100011000110
000101101011
110011010111
111001011110
111010100110
010001110001
011010101110
101011110000
001101010010
001010000011
011001000111
010001000110
010110011111
001111011000
000010010111
100011110111
111100111100
010101010001
100011101111
011101011001
001011001100
001111000101
110111011000
011001001010
010000000100
100110111010
000000111110
011100100011
100100111110
000010011010
111000000100
001110001111
011000101010
011100110010
110101000111
110110111111
101110110010
000110100000
110100001101
000111101110
111000111001
000111111010
010110111011
000101110111
111000011011
010110001010
001011111011
100100100011
011110110010
000111110101
110110011101
000011000011
110011001010
000100101100
101100100011
010111000001
110101100001
011110111010
111100101000
010000100010
110011000111
010011011100
100110101011
000100001011
101111000010
010000010000
111110000100
111100011110
010100111001
100000011111
010111110100
100100110011
111101111001
100010011000
111110101111
010111001010
110000101000
001111111111
000111100110
000011101100
011010001000
001011010011
110011111011"""

"""
    parse_diagnostic_codes(codes::AbstractString) ::Vector{AbstractString}

Get the individual diagnostic codes.
"""
function parse_diagnostic_codes(codes::AbstractString) ::Vector{AbstractString}
    return split(codes, "\n")
end

"""
flip_bit(char::AbstractChar) ::AbstractChar

Flip the string bit.
"""
function flip_bit(char::AbstractChar) ::AbstractChar
if char == '1'
    return '0'
else
    return '1'
end
end

"""
    most_common_bit(codes::Vector{AbstractString}, position::Integer, default::AbstractChar='0') ::AbstractChar

Get the most common bit in position i.
"""
function most_common_bit(codes::Vector{<:AbstractString}, position::Integer; default::AbstractChar='0') ::AbstractChar
    ones::Integer = 0
    zeroes::Integer = 0
    # calculate number of ones and zeroes in i'th spot
    for code in codes
        if code[position] == '1'
            ones += 1
        else
            zeroes += 1
        end
    end
    # return the most common.
    if ones > zeroes
        return '1'
    elseif zeroes > ones
        return '0'
    else
        return default
    end
end

"""
    get_oxygen_rate(codes::Vector{AbstractString}) ::AbstractString

Get the oxygen rate from the gamma. This involves getting the entry that matches the bit in
each position of gamma as close as possible.
"""
function get_oxygen_rate(codes::Vector{AbstractString}; flip::Bool=false) ::AbstractString
    # get the length of the binary number.
    len::Int8 = length(codes[1])
    # for each index, calculate the count of each binary digit and use the most common.
    for i=1:len
        # get most common bit in position i.
        most_common = most_common_bit(codes, i, default='1')
        if flip
            most_common = flip_bit(most_common)
        end
        # get only codes that match gamma in this bit.
        codes = [code for code in codes if code[i] == most_common]
        # if we only have one code left, break.
        if length(codes) == 1
            return only(codes)
        end
    end
    return returnval
end

"""
    get_scrubber_rate(codes::Vector{AbstractString}) ::AbstractString

Get the scrubber rate from the gamma. This involves getting the entry that does not match the bit in
each position of gamma as close as possible.
"""
function get_scrubber_rate(codes::Vector{AbstractString}) ::AbstractString
    return get_oxygen_rate(codes, flip=true)
end

# Run the calculation!
diagnostics = parse_diagnostic_codes(diags)
# Our gamma is just the most common digit at each position, with one if they are both equal
# Thus we can calculate from this the oxygen rate by filtering on matches with these
#   digits.
# we can also calculate scrubber by checking it doesn't match these digits, since we
#   want least common (inverse of gamma) and 0 instead of 1 when equal.
gamma = get_gamma_rate(diagnostics)
oxygen = get_oxygen_rate(diagnostics)
scrubber = get_scrubber_rate(diagnostics)
parse(Int64, oxygen, base=2) * parse(Int64, scrubber, base=2)


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.6.0
#   kernelspec:
#     display_name: Julia 1.5.2
#     language: julia
#     name: julia-1.5
# ---

# # WebIO.jl Example

using WebIO
using JSExpr

# +
scope = Scope()
obs = Observable(scope, "rand-value", 0.0)

on(obs) do x
    println("JavaScript sent $(x)!")
end

scope(
    dom"button"(
        "Generate Random Number",
        events=Dict(
            "click" => @js () -> $obs[] = Math.random() + $obs[]
        ),
    ),
)

# -

# # Reference
#
# https://juliagizmos.github.io/WebIO.jl/latest/gettingstarted/#Sending-values-from-JavaScript-to-Julia

max_stars_repo_name: terasakisatoshi/jsexample

max_issues_repo_name: terasakisatoshi/jsexample

max_forks_repo_name: terasakisatoshi/jsexample