text:
const rebug_prompt_string = "rebug> "
const rebug_prompt_ref = Ref{Union{LineEdit.Prompt,Nothing}}(nothing)   # set by __init__

dummy() = nothing
const dummymethod = first(methods(dummy))
const dummyuuid   = UUID(UInt128(0))

uuidextractor(str) = match(r"getstored\(\"([a-z0-9\-]+)\"\)", str)

mutable struct RebugHeader <: AbstractHeader
    warnmsg::String
    errmsg::String
    uuid::UUID
    current_method::Method
    nlines::Int   # size of the printed header
end
RebugHeader() = RebugHeader("", "", dummyuuid, dummymethod, 0)

function header(s::LineEdit.MIState)
    rebug_prompt = rebug_prompt_ref[]
    rebug_prompt.repl.header
end

# Custom methods
set_method!(header::RebugHeader, method::Method) = header.current_method = method

function set_uuid!(header::RebugHeader, uuid::UUID)
    if haskey(stored, uuid)
        header.uuid = uuid
        header.current_method = stored[uuid].method
    else
        header.uuid = dummyuuid
        header.current_method = dummymethod
    end
    uuid
end

function set_uuid!(header::RebugHeader, str::AbstractString)
    m = uuidextractor(str)
    uuid = if m isa RegexMatch && length(m.captures) == 1
        UUID(m.captures[1])
    else
        dummyuuid
    end
    set_uuid!(header, uuid)
end

struct FakePrompt{Buf<:IO}
    input_buffer::Buf
end
LineEdit.mode(p::FakePrompt) = rebug_prompt_ref[]

"""
    stepin(s)

Given a buffer `s` representing a string and "point" (the seek position) set at a call expression,
replace the contents of the buffer with a `let` expression that wraps the *body* of the callee.

For example, if `s` has contents

    <some code>
    if x > 0.5
        ^fcomplex(x)
        <more code>

where in the above `^` indicates `position(s)` ("point"), and if the definition of `fcomplex` is

    function fcomplex(x::A, y=1, z=""; kw1=3.2) where A<:AbstractArray{T} where T
        <body>
    end

rewrite `s` so that its contents are

    @eval ModuleOf_fcomplex let (x, y, z, kw1, A, T) = Main.Rebugger.getstored(id)
        <body>
    end

where `Rebugger.getstored` returns has been pre-loaded with the values that would have been
set when you called `fcomplex(x)` in `s` above.
This line can be edited and `eval`ed at the REPL to analyze or improve `fcomplex`,
or can be used for further `stepin` calls.
"""
function stepin(s::LineEdit.MIState)
    # Add the command we're tracing to the history. That way we can go "up the call stack".
    pos = position(s)
    cmd = String(take!(copy(LineEdit.buffer(s))))
    add_history(s, cmd)
    # Analyze the command string and step in
    local uuid, letcmd
    try
        uuid, letcmd = stepin(LineEdit.buffer(s))
    catch err
        repl = rebug_prompt_ref[].repl
        handled = false
        if err isa StashingFailed
            repl.header.warnmsg = "Execution did not reach point"
            handled = true
        elseif err isa Meta.ParseError || err isa StepException
            repl.header.warnmsg = "Expression at point is not a call expression"
            handled = true
        elseif err isa EvalException
            repl.header.errmsg = "$(typeof(err.exception)) exception while evaluating $(err.exprstring)"
            handled = true
        elseif err isa DefMissing
            repl.header.errmsg = "The expression for method $(err.method) was unavailable. Perhaps it was untracked or generated by code."
            handled = true
        end
        if handled
            buf = LineEdit.buffer(s)
            LineEdit.edit_clear(buf)
            write(buf, cmd)
            seek(buf, pos)
            return nothing
        end
        rethrow(err)
    end
    set_uuid!(header(s), uuid)
    LineEdit.edit_clear(s)
    LineEdit.edit_insert(s, letcmd)
    return nothing
end

function capture_stacktrace(s)
    if mode(s) isa LineEdit.PrefixHistoryPrompt
        # history search, re-enter with the corresponding mode
        LineEdit.accept_result(s, mode(s))
        return capture_stacktrace(s)
    end
    cmdstring = LineEdit.content(s)
    add_history(s, cmdstring)
    print(REPL.terminal(s), '\n')
    expr = Meta.parse(cmdstring)
    uuids = capture_stacktrace(expr)
    io = IOBuffer()
    buf = FakePrompt(io)
    hp = mode(s).hist
    for uuid in uuids
        println(io, generate_let_command(uuid))
        REPL.add_history(hp, buf)
        take!(io)
    end
    hp.cur_idx = length(hp.history) + 1
    if !isempty(uuids)
        set_uuid!(header(s), uuids[end])
        print(REPL.terminal(s), '\n')
        LineEdit.edit_clear(s)
        LineEdit.enter_prefix_search(s, find_prompt(s, LineEdit.PrefixHistoryPrompt), true)
    end
    return nothing
end

function add_history(s, str::AbstractString)
    io = IOBuffer()
    buf = FakePrompt(io)
    hp = mode(s).hist
    println(io, str)
    REPL.add_history(hp, buf)
end

### REPL mode

function HeaderREPLs.print_header(io::IO, header::RebugHeader)
    if header.nlines != 0
        HeaderREPLs.clear_header_area(io, header)
    end
    iocount = IOBuffer()  # for counting lines
    for s in (io, iocount)
        if !isempty(header.warnmsg)
            printstyled(s, header.warnmsg, '\n'; color=Base.warn_color())
        end
        if !isempty(header.errmsg)
            printstyled(s, header.errmsg, '\n'; color=Base.error_color())
        end
        if header.current_method != dummymethod
            printstyled(s, header.current_method, '\n'; color=:light_magenta)
        end
        if header.uuid != dummyuuid
            data = stored[header.uuid]
            ds = displaysize(io)
            printer(args...) = printstyled(args..., '\n'; color=:light_blue)
            for (name, val) in zip(data.varnames, data.varvals)
                # Make sure each only spans one line
                try
                    Revise.printf_maxsize(printer, s, "  ", name, " = ", val; maxlines=1, maxchars=ds[2]-1)
                catch # don't error just because a print method is borked
                    printstyled(s, "  ", name, " errors in its show method"; color=:red)
                end
            end
        end
    end
    header.nlines = count_display_lines(iocount, displaysize(io))
    header.warnmsg = ""
    header.errmsg = ""
    return nothing
end

function HeaderREPLs.setup_prompt(repl::HeaderREPL{RebugHeader}, hascolor::Bool)
    julia_prompt = find_prompt(repl.interface, "julia")

    prompt = REPL.LineEdit.Prompt(
        rebug_prompt_string;
        prompt_prefix = hascolor ? repl.prompt_color : "",
        prompt_suffix = hascolor ?
            (repl.envcolors ? Base.input_color : repl.input_color) : "",
        complete = julia_prompt.complete,
        on_enter = REPL.return_callback)

    prompt.on_done = HeaderREPLs.respond(repl, julia_prompt) do str
        Base.parse_input_line(str; filename="REBUG")
    end
    # hist will be handled automatically if repl.history_file is true
    # keymap_dict is separate
    return prompt, :rebug
end

function HeaderREPLs.append_keymaps!(keymaps, repl::HeaderREPL{RebugHeader})
    julia_prompt = find_prompt(repl.interface, "julia")
    kms = [
        trigger_search_keymap(repl),
        mode_termination_keymap(repl, julia_prompt),
        trigger_prefix_keymap(repl),
        REPL.LineEdit.history_keymap,
        REPL.LineEdit.default_keymap,
        REPL.LineEdit.escape_defaults,
    ]
    append!(keymaps, kms)
end

# To get it to parse the UUID whenever we move through the history, we have to specialize
# this method
function HeaderREPLs.activate_header(header::RebugHeader, p, s, termbuf, term)
    str = String(take!(copy(LineEdit.buffer(s))))
    set_uuid!(header, str)
end

const keybindings = Dict{Symbol,String}(
    :stacktrace => "\es",                        # Alt-s ("[s]tacktrace")
    :stepin => "\ee",                            # Alt-e ("[e]nter")
)
const deprecated_keybindings = Dict{Symbol,Vector{String}}(
    :deprecated_stacktrace => ["\e[15~"],        # F5
    :deprecated_stepin => ["\e\eOM", "\e[23~"],  # Alt-Shift-Enter, F11
)

const modeswitches = Dict{Any,Any}(
    :stacktrace => (s, o...) -> capture_stacktrace(s),
    :stepin => (s, o...) -> (stepin(s); enter_rebug(s)),
    :deprecated_stacktrace => (s, o...) -> (rebug_prompt_ref[].repl.header.warnmsg="stacktrace is now Alt-s (see docs for customization options)"; capture_stacktrace(s)),
    :deprecated_stepin => (s, o...) -> (stepin(s); rebug_prompt_ref[].repl.header.warnmsg="stepin is now Alt-e (see docs for customization options)"; enter_rebug(s)),
)

function get_rebugger_modeswitch_dict()
    rebugger_modeswitch = Dict()
    # Process the deprecated bindings first; that way, if the user customizes it back
    # to one of the deprecated choices, it will override without issuing the warning.
    for (action, list) in deprecated_keybindings
        for keybinding in list
            rebugger_modeswitch[keybinding] = modeswitches[action]
        end
    end
    for (action, keybinding) in keybindings
        rebugger_modeswitch[keybinding] = modeswitches[action]
    end
    rebugger_modeswitch
end

function add_keybindings(; override::Bool=false, kwargs...)
    main_repl = Base.active_repl
    history_prompt = find_prompt(main_repl.interface, LineEdit.PrefixHistoryPrompt)
    julia_prompt = find_prompt(main_repl.interface, "julia")
    rebug_prompt = find_prompt(main_repl.interface, "rebug")
    for (action, keybinding) in kwargs
        if !(action in keys(keybindings)) && !(action in keys(deprecated_keybindings))
            error("$action is not a supported action.")
        end
        if !(keybinding isa Union{String,Vector{String}})
            error("Expected the value for $action to be a String or Vector{String}, got $keybinding instead")
        end
        if haskey(keybindings, action)
            keybindings[action] = keybinding
        end
        # We need Any here because "cannot convert REPL.LineEdit.PrefixHistoryPrompt to an object of type REPL.LineEdit.Prompt"
        prompts = Any[julia_prompt, rebug_prompt]
        if action ∈ (:stacktrace, :deprecated_stacktrace) push!(prompts, history_prompt) end
        for prompt in prompts
            if keybinding isa Vector
                for kb in keybinding
                    LineEdit.add_nested_key!(prompt.keymap_dict, kb, modeswitches[action], override=override)
                end
            else
                LineEdit.add_nested_key!(prompt.keymap_dict, keybinding, modeswitches[action], override=override)
            end
        end
    end
end

# These work only at the `rebug>` prompt
const rebugger_keys = Dict{Any,Any}(
)

## REPL commands TODO?:
## "\em" (meta-m): create REPL line that populates Main with arguments to current method
## "\eS" (meta-S): save version at REPL to file? (a little dangerous, perhaps make it configurable as to whether this is on)
## F1 is "^[OP" (showvalues?), F4 is "^[OS" (showinputs?)

enter_rebug(s) = mode_switch(s, rebug_prompt_ref[])

function mode_switch(s, other_prompt)
    buf = copy(LineEdit.buffer(s))
    LineEdit.edit_clear(s)
    LineEdit.transition(s, other_prompt) do
        LineEdit.state(s, other_prompt).input_buffer = buf
    end
end

# julia_prompt = find_prompt(main_repl.interface, "julia")
# julia_prompt.keymap_dict['|'] = (s, o...) -> enter_count(s)


text:
# This is a collection of toys under the Yggdrasil tree for the good <s>kids</s>
# developers.  These utilities can be employed in builder files.

using BinaryBuilder, Pkg, LibGit2

"""
    should_build_platform(platform) -> Bool

Return whether the tarballs for the given `platform` should be built.

This is useful when the builder has different platform-dependent elements
(sources, script, products, etc...) that make it hard to have a single
`build_tarballs` call.
"""
function should_build_platform(platform)
    # If you need inspiration for how to use this function, look at the builder
    # for Git:
    # https://github.com/JuliaPackaging/Yggdrasil/blob/c3e3c4a96c723306b4da23fc6d05f12995b21ed8/G/Git/build_tarballs.jl#L76-L93

    # Get the list of platforms requested from the command line.  This should be
    # the only argument not prefixed with "--".
    requested_platforms = filter(arg -> !occursin(r"^--.*", arg), ARGS)

    if isone(length(requested_platforms))
        # `requested_platforms` has only one element: the comma-separated list
        # of platform.  We'll run the platform only if it's in the list
        return any(platforms_match.(Ref(platform), split(requested_platforms[1], ",")))
    else
        # `requested_platforms` doesn't have only one element: if its length is
        # zero, no platform has been explicitely passed from the command line
        # and we we'll run all platforms, otherwise we don't know what to do, so
        # let's return false to be safe.
        return iszero(length(requested_platforms))
    end
end

# compatibility for Julia 1.3-
if VERSION < v"1.4"
    Pkg.Types.registry_resolve!(ctx::Pkg.Types.Context, deps) = Pkg.Types.registry_resolve!(ctx.env, deps)
end

"""
    get_tree_hash(tree::LibGit2.GitTree)

Given a `GitTree`, get the `GitHash` that identifies it.
"""
function get_tree_hash(tree::LibGit2.GitTree)
    oid_ptr = Ref(LibGit2.GitHash())
    oid_ptr = ccall((:git_tree_id, :libgit2), Ptr{LibGit2.GitHash}, (Ptr{Cvoid},), tree.ptr)
    oid_ptr == C_NULL && throw("bad tree ID: $tree")
    return unsafe_load(oid_ptr)
end

"""
    get_addable_spec(name, version)

Given a JLL name and registered version, return a `PackageSpec` that, when passed as a
`Dependency`, ensures that exactly that version will be installed.  Example usage:

    dependencies = [
        BuildDependency(get_addable_spec("LLVM_jll", v"9.0.1+0")),
    ]
"""
function get_addable_spec(name::AbstractString, version::VersionNumber)
    Pkg.Registry.update()
    ctx = Pkg.Types.Context()

    # First, resolve the UUID
    uuid = first(Pkg.Types.registry_resolve!(ctx, Pkg.Types.PackageSpec(name=name))).uuid

    # Next, determine the tree hash from the registry
    tree_hashes = Base.SHA1[]
    @static if VERSION >= v"1.4"
        paths = Pkg.Operations.registered_paths(ctx, uuid)
    else
        paths = Pkg.Operations.registered_paths(ctx.env, uuid)
    end
    for path in paths
        vers = Pkg.Operations.load_versions(path; include_yanked = true)
        tree_hash = get(vers, version, nothing)
        tree_hash !== nothing && push!(tree_hashes, tree_hash)
    end

    if isempty(tree_hashes)
        error("Unable to find mapping for $(name) v$(version)")
    end
    tree_hash_sha1 = first(tree_hashes)
    tree_hash_bytes = tree_hash_sha1.bytes

    # Once we have a tree hash, turn that into a git commit sha
    url = "https://github.com/JuliaBinaryWrappers/$(name).jl"
    git_commit_sha = nothing
    mktempdir() do dir
        repo = LibGit2.clone(url, dir)
        LibGit2.with(LibGit2.GitRevWalker(repo)) do walker
            LibGit2.push_head!(walker)
            for oid in walker
                tree = LibGit2.peel(LibGit2.GitTree, LibGit2.GitCommit(repo, oid))
                if all(get_tree_hash(tree).val .== tree_hash_bytes)
                    git_commit_sha = LibGit2.string(oid)
                    break
                end
            end
        end
    end

    if git_commit_sha === nothing
        error("Unable to find corresponding revision in $(name) v$(version) for tree hash $(bytes2hex(tree_hash_bytes))")
    end

    @static if VERSION >= v"1.4"
        repo=Pkg.Types.GitRepo(rev=git_commit_sha, source=url)
    else
        repo=Pkg.Types.GitRepo(rev=git_commit_sha, url=url)
    end
    return Pkg.Types.PackageSpec(
        name=name,
        uuid=uuid,
        version=version,
        tree_hash=tree_hash_sha1,
        repo=repo,
    )
end

ext: jl


text:
# This file holds the definition of the functions pertaining to the Condition
#   type.

# The functions of the Transition type require no other types.
requiredTypes = [ "condition" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# ==============================================================================
# Non-exported methods.
# ==============================================================================

relationFunctions = Dict{String, Function}( "IS" => ==, "IS NOT" => !=,
    "IN" => ∈, "NOT IN" => ∉, ">" => >, ">=" => >=, "<" => <, "<=" => <= )
timeAttrs = [ "age", "tenure", "time in state" ]

function processCondition( attr::String, rel::String, val::T ) where T <: Union{Real, String}

    newCond = Condition{Int}( "", ==, 0 )
    relOp = get( relationFunctions, rel, nothing )

    # If the operator is unknown (shouldn't happen), do not process the
    #   condition.
    if isa( relOp, Void )
        warn( "Unknown operator. This warning should not appear!" )
        return newCond, false
    end  # if isa( relOp, Void )

    attrName = strip( attr )
    attrName = lowercase( attrName ) ∈ timeAttrs ? lowercase( attrName ) :
        attrName

    # If the condition is on the age, only process if the operator makes sense.
    if ( attrName ∈ timeAttrs ) && ( relOp ∉ [ ∈, ∉ ] )
        newCond = Condition{Float64}( attrName, relOp, val * 12.0 )
    else
        attrVal = val
        valType = T

        if relOp ∈ [ ∈, ∉ ]
            attrVal = split( val, "," )
            attrVal = Vector{String}( map( opt -> strip( opt ), attrVal ) )
            valType = Vector{String}
        end  # if relOp ∈ [ ∈, ∉ ]

        newCond = Condition{valType}( attrName, relOp, attrVal )
    end  # if ( attrName ∈ timeAttrs ) && ...

    return newCond, true

end  # processCondition( attr::String, rel::String, val::T )


text:
f = open(readlines, "hightemp.txt")
n = int(ARGS[1])
filelength = length(f)
if n > filelength
  n = filelength
end
s = filelength - n + 1

for l in f[s:filelength]
  print(l)
end


text:
#functions for computing mutual informations (in bits)
#po is the distribution where the KL is averaged over
#the KL measures going from px to pxgy
function mutualinformation(py::Vector, px::Vector, pxgy::Matrix)
    card_y = size(py,1)
    card_x = size(px,1)
    if (size(pxgy,1)!=card_x) || (size(pxgy,2)!=card_y)
        error("Dimensionality of p(x|y) does not match p(x), p(y)!")
    end

    MI = 0
    for i in 1:card_y
        MI += py[i] * kl_divergence_bits(vec(pxgy[:,i]),px)
    end

    return MI
end

#functions for computing the conditional mutual information
#I(X;Y|Z)(in bits)
function conditional_mutualinformation(py::Vector, pzgy::Matrix, pxgz::Matrix, pxgzy)
    #I(X;Y|Z) = ∑_y p(y) DKL( p(x|z,y)||p(x|z) )
    #I(X;Y|Z) = ∑_y p(y) ∑_z p(z|y) ∑_x p(x|z,y) log( p(x|z,y)/p(x|z) )

    card_y = length(py)
    card_z = size(pzgy,1)

    MI = 0
    for j in 1:card_y
        for k in 1:card_z
            MI += py[j]*pzgy[k,j] * kl_divergence_bits( vec(pxgzy[:,k,j]), vec(pxgz[:,k]))
        end
    end

    return MI
end

#function for computing the expected utility
#pagw and U_pre must have the same dimensionality
function expectedutility(pw::Vector, pagw::Matrix, U_pre::Matrix)
    #E[U] = ∑_a,w p(w)p(a|w) U(a,w)
    card_w = size(pw,1)
    card_a = size(pagw,1)
    if (size(pagw,2)!=card_w) || (size(U_pre,1)!=card_a) || (size(U_pre,2)!=card_w)
        error("Dimensionality of p(a|w), U_pre(a,w) and p(w) does not match!")
    end

    EU = 0
    for i in 1:card_w
            EU += pw[i] * sum(pagw[:,i] .* U_pre[:,i])
    end

    return EU
end

#This function assumes that the utility function U(a,w) is not a function of the percept o
function expectedutility(pw::Vector, pogw::Matrix, pagow, U_pre::Matrix)
    #E[U] = ∑_a,o,w p(w)p(o|w)p(a|o,w) U(a,w)

    #reuse the function that computes E[U] for ∑_a,w p(w)p(a|w) U(a,w)
    #to do so, compute p(a|w)
    pagw = marginalizeo(pogw,pagow)

    return expectedutility(pw,pagw,U_pre)
end

#Entropy in bits (using entropy from Distributions.jl)
function entropybits(d::Distribution)
    return Distributions.entropy(d,2) #in bits
end

#Entropy in bits for an discrete distribution represented as a vector
function entropybits(p::Vector)
    return entropybits(Categorical(p))
end

#Explicitly provide a function for the log in bits here.
#The rationale behind this is that it makes it easily possible to spot all the
#places (in the code) to change if one wants to use the nats (natural logarithm) intead of bits.
function log_bits(x)
    return log2(x)
end

#Kullback-Leibler divergence in bits
function kl_divergence_bits(p_x::Vector, p0_x::Vector)
    #D_KL_bits = ∑_x p(x) log2 (p(x)/p0(x))

    if sum(p0_x.==0) > 0
        error("Zeros in denominator before kl_divergence computation!")
    end

    kl_div = kl_divergence(p_x, p0_x)/log(2) #using function from Distances.jl

    return kl_div
end

#compute value of rate-distortion objective (avg ΔF)
function RDobjective(EU,I,β)
    return EU-I/β
end

#Objective-value for the three-variable general case
function ThreeVArRDobjective(EU, I_ow, I_ao, I_awgo, β1, β2, β3)
    if(β3==0)
        #sequential case - I(A;W|O) should be zero and since β3=Inf → ∞ * 0 = NaN
        #handle this by using a tiny β3 - if I(A;W|O) is nonzero, this should be easily detectable through the huge
        #value of the objective
        β3 = eps()
    end
    return EU - (1/β1)*I_ow - (1/β2)*I_ao - (1/β3)*I_awgo
end

#compute I(A;W), H(A), H(A|W), E[U] and E[U]-I(A;W)/β
function analyzeBAsolution(pw::Vector, pa::Vector, pagw::Matrix, U_pre::Matrix, β)
    #compute I(A;W)
    I = mutualinformation(pw,pa,pagw)
    #compute H(A)
    Ha = entropybits(pa)
    #compute H(A|W)
    Hagw = Ha-I
    #compute EU
    EU = expectedutility(pw,pagw,U_pre)
    #compute value of objective
    RDobj = RDobjective(EU,I,β)

    return I, Ha, Hagw, EU, RDobj
end

#compute mutual informations, entropies and value of objective for three-variable general case
function analyze_three_var_BAsolution(pw::Vector, po::Vector, pa::Vector, pogw::Matrix,
    pago::Matrix, pagow, pagw::Matrix, U_pre::Matrix, β1, β2, β3)

    #compute I(O;W)
    I_ow = mutualinformation(pw,po,pogw)

    #compute I(A;O)
    I_ao = mutualinformation(po,pa,pago)

    #compute I(A;W|O)
    I_awgo = conditional_mutualinformation(pw, pogw, pago, pagow)

    #compute I(A;W)
    I_aw = mutualinformation(pw,pa,pagw)

    #compute H(O)
    #println("1: ∑p(o) = $(sum(po))")
    Ho = entropybits(po)

    #compute H(A)
    #println("2: ∑p(a) =  $(sum(pa))")
    Ha = entropybits(pa)

    #compute H(O|W)
    Hogw = Ho - I_ow

    #compute H(A|O)
    Hago = Ha - I_ao

    #compute H(A|O,W)
    Hagow = Hago - I_awgo

    #compute H(A|W)
    Hagw = Ha - I_aw

    #compute EU
    EU = expectedutility(pw,pagw,U_pre)

    #compute value of objective
    ThreeVarRDobj = ThreeVArRDobjective(EU, I_ow, I_ao, I_awgo, β1, β2, β3)

    return I_ow, I_ao, I_awgo, I_aw, Ho, Ha, Hogw, Hago, Hagow, Hagw, EU, ThreeVarRDobj
end


text:
using TimeseriesPrediction
import Statistics: mean
using Test

import Random
Random.seed!(42)

include("system_defs.jl")

@testset "Barkley Const Boundary" begin
    kwargs = ( tskip = 200,
               ssize = (50,50),
               periodic= false,
               a=0.75, b=0.02, ε=0.02, D=1, h=0.75, Δt=0.1)
    Ttrain = 400
    p = 20
    T = Ttrain + p
    U, V = barkley(T+100; kwargs... )
    Vtrain = V[1:Ttrain]
    Utrain = U[1:Ttrain]

    τ = 1
    k = 1
    BC = ConstantBoundary(20.)

    @testset "V, γ=$γ, B=$B" for γ=10, B=1
        Vtest  = V[Ttrain :  T]
        em = cubic_shell_embedding(Vtrain,γ,τ,B,k,BC)
        em = PCAEmbedding(Vtrain, em)
        Vpred = temporalprediction(Vtrain,em,p)
        @test Vpred[1] == Vtrain[end]
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
        end
    end
    @testset "temp. pred. with offset start" for γ=10, B=1
        Vtest  = V[Ttrain+100:T+100]
        Vstart = V[Ttrain+100-11:Ttrain+100]
        em = cubic_shell_embedding(Vtrain,γ,τ,B,k,BC)
        em = PCAEmbedding(Vtrain, em)
        Vpred = temporalprediction(Vtrain,em,p; initial_ts=Vstart)
        @test Vpred[1] == Vstart[end]
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
        end
    end
    @testset "V LightCone" begin
        Vtest  = V[Ttrain :  T]
        em = SpatioTemporalEmbedding(Vtrain, (γ=1,τ=1,r=1,c=0.5,bc=BC))
        Vpred = temporalprediction(Vtrain,em,p)
        @test Vpred[1] == Vtrain[end]
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
        end
    end
    @testset "U, γ=2, B=1" begin
        γ=3; B=1; τ=4
        #high embedding dim, so predict fewer points to save time
        p = 10
        Utest  = U[Ttrain :  T]
        em = cubic_shell_embedding(Utrain,γ,τ,B,k,BC)
        Upred = temporalprediction(Utrain,em,p)
        @test Upred[1] == Utrain[end]
        err = [abs.(Utest[i]-Upred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
        end
    end
    @testset "crosspred V → U" begin
        γ = 4; B = 1
        p=10
        Utest  = U[Ttrain + 1:  T]
        Vtest  = V[Ttrain  - γ*τ+ 1:  T]
        em = cubic_shell_embedding(Vtrain, γ,τ,B,k,BC)
        Upred = crossprediction(Vtrain,Utrain,Vtest, em)
        err = [abs.(Utest[i]-Upred[i]) for i=1:p-1]
        for i in 1:length(err)
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
            @test mean(err[i]) < 0.1
        end
    end
end

@testset "Periodic Barkley" begin
    kwargs = ( tskip = 200,
               ssize = (50,50),
               periodic= true,
               a=0.75, b=0.02, ε=0.02, D=1, h=0.75, Δt=0.1)
    Ttrain = 400
    p = 20
    T = Ttrain + p
    U, V = barkley(T; kwargs... )
    Vtrain = V[1:Ttrain]
    Utrain = U[1:Ttrain]
    τ = 1
    k = 1
    BC = PeriodicBoundary()

    @testset "crosspred V → U" begin
        γ = 3; B = 1
        Utest  = U[Ttrain + 1:  T]
        Vtest  = V[Ttrain  - γ*τ+ 1:  T]
        em = cubic_shell_embedding(Vtrain, γ,τ,B,k,BC)
        Upred = crossprediction(Vtrain,Utrain,Vtest, em)
        err = [abs.(Utest[i]-Upred[i]) for i=1:p-1]
        for i in 1:length(err)
            @test maximum(err[i]) < 0.2
            @test mean(err[i]) < 0.1
            @test minimum(err[i]) < 0.1
        end
    end

    @testset "Periodic, γ=$γ, B=$B" for γ=10, B=1
        Vtest  = V[Ttrain :  T]
        em = cubic_shell_embedding(Vtrain, γ,τ,B,k,BC)
        em = PCAEmbedding(Vtrain, em)
        Vpred = temporalprediction(Vtrain, em, p)
        @test Vpred[1] == Vtrain[end]
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.1
            @test minimum(err[i]) < 0.1
        end
    end

    @testset "Periodic diff. inital, γ=$γ, B=$B" for γ=10, B=1
        U, V = barkley(T; kwargs... )
        Vtrain = V[1:Ttrain]
        Vtest  = V[Ttrain :  T]
        em = cubic_shell_embedding(Vtrain, γ,τ,B,k,BC)
        em = PCAEmbedding(Vtrain, em)
        Vpred = temporalprediction(Vtrain, em, p)
        @test Vpred[1] == Vtrain[end]
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
        end
    end
end

@testset "Half-Precision Tests" begin
    kwargs = ( tskip = 200,
               ssize = (50,50),
               periodic= false,
               a=0.75, b=0.02, ε=0.02, D=1, h=0.75, Δt=0.1)
    Ttrain = 400
    p = 20
    T = Ttrain + p
    U, V = barkley(T; kwargs... )
    U = map(u -> Float32.(u), U)
    V = map(v -> Float32.(v), V)
    Vtrain = V[1:Ttrain]
    Utrain = U[1:Ttrain]

    τ = 1
    k = 1
    BC = ConstantBoundary{Float32}(20.)

    @testset "Temporal Prediction" begin
        γ=10; B=1
        Vtest  = V[Ttrain :  T]
        em = cubic_shell_embedding(Vtrain,γ,τ,B,k,BC)
        em = PCAEmbedding(Vtrain, em)
        Vpred = temporalprediction(Vtrain,em,p)
        @test Vpred[1] == Vtrain[end]
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
        for i in 1:p
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
        end
        @test Float32 == eltype(Vpred[1])
    end
    @testset "Cross Prediction" begin
        γ = 10; B = 0
        p=10
        Utest  = U[Ttrain-γ*τ+1:T]
        Vtest  = V[Ttrain+1:T]
        em = cubic_shell_embedding(Utrain, γ,τ,B,k,BC)
        Vpred = crossprediction(Utrain,Vtrain,Utest, em)
        err = [abs.(Vtest[i]-Vpred[i]) for i=1:p-1]
        for i in 1:length(err)
            @test maximum(err[i]) < 0.2
            @test minimum(err[i]) < 0.1
            @test mean(err[i]) < 0.1
        end
        @test Float32 == eltype(Vpred[1])
    end
end

# Complex Numbers are not compatible with KDTrees
# @testset "Complex Number Test" begin
#     kwargs = ( tskip = 200,
#                ssize = (50,50),
#                periodic= false,
#                a=0.75, b=0.02, ε=0.02, D=1, h=0.75, Δt=0.1)
#     Ttrain = 400
#     p = 20
#     T = Ttrain + p
#     U, V = barkley(T; kwargs... )
#     UV = map(U,V) do u,v
#         u .+ 1im* v
#     end
#     UVtrain = UV[1:Ttrain]
#     UVtest  = UV[Ttrain :  T]
#
#     @testset "Temporal Prediction" begin
#         D=0; τ = 1; r = 1; c = 0
#         BC = ConstantBoundary{ComplexF64}(20.)
#         em = light_cone_embedding(UVtrain,D,τ,r,c,BC)
#         UVpred = temporalprediction(UVtrain,em,p)
#         @test Vpred[1] == Vtrain[end]
#         err = [abs.(Vtest[i]-Vpred[i]) for i=1:p+1]
#         for i in 1:p
#             @test maximum(err[i]) < 0.2
#             @test minimum(err[i]) < 0.1
#         end
#     end
#
#
# end


text:
@testset "Test Joins" begin
    y = rand(10)
    z = rand(10)

    t  = table((x=1:10,   y=y), pkey=:x)
    t2 = table((x=1:2:20, z=z), pkey=:x)

    @testset "how = :inner" begin
        t_inner = table((x = 1:2:9, y = y[1:2:9], z = z[1:5]), pkey = :x)
        @test isequal(join(t, t2; how=:inner), t_inner)
    end
    @testset "how = :left" begin
        # Missing
        z_left = Union{Float64,Missing}[missing for i in 1:10]
        z_left[1:2:9] = z[1:5]
        t_left = table((x = 1:10, y = y, z = z_left))
        @test isequal(join(t, t2; how=:left), t_left)

        # DataValue
        z_left2 = [DataValue{Float64}() for i in 1:10]
        z_left2[1:2:9] = z[1:5]
        t_left2 = table((x=1:10, y = y, z = z_left2))
        @test isequal(join(t, t2, how=:left, missingtype=DataValue), t_left2)
    end
    @testset "how = :outer" begin
        # Missing
        x_outer = union(1:10, 1:2:20)
        y_outer = vcat(y, fill(missing, 5))
        z_left = Union{Float64,Missing}[missing for i in 1:10]
        z_left[1:2:9] = z[1:5]
        z_outer = vcat(z_left, z[6:10])
        t_outer = table((x=x_outer, y=y_outer, z=z_outer); pkey=:x)
        @test isequal(join(t, t2; how=:outer), t_outer)

        # DataValue
        y_outer2 = vcat(y, fill(DataValue{Float64}(), 5))
        z_left2 = [DataValue{Float64}() for i in 1:10]
        z_left2[1:2:9] = z[1:5]
        z_outer2 = vcat(z_left2, z[6:10])
        t_outer2 = table((x=x_outer, y=y_outer2, z=z_outer2); pkey=:x)
        @test isequal(join(t, t2; how=:outer, missingtype=DataValue), t_outer2)
    end
    @testset "how = :anti" begin
        t_anti = table((x=2:2:10, y=y[2:2:10]), pkey=:x)
        @test isequal(join(t, t2; how=:anti), t_anti)
    end
end


text:
module MagicMarshalling

export marshal, unmarshal

using ArgCheck: @check, @argcheck
using ConstructionBase: constructorof, ConstructionBase

KEY_META = "__meta__"
KEY_CONSTRUCTOR = "constructor"
KEY_PROPERTYNAMES = "propertynames"
PKGNAME = "MagicMarshalling"

function __init__()
    # THIS IS HACKY
    # We do this so that ConstructionBase is visible in Main
    Core.eval(Main, :(using ConstructionBase: ConstructionBase))
end

constructor_string(o) = string(constructorof(typeof(o)))
constructor_string(::Tuple) = string(tuple)
constructor_string(o::Array) = string(reshape)

unmarshal(o::Vector) = map(unmarshal, o)
function marshal(o::Array)
    d = Dict{String, Any}()
    meta = Dict(
        KEY_CONSTRUCTOR => constructor_string(o),
        KEY_PROPERTYNAMES => ["data", "size"], # eltype?
    )
    d[KEY_META] = meta
    d["data"] = map(marshal, vec(o))
    d["size"] = marshal(size(o))
    d
end

function marshal_struct(o)
    d = Dict{String, Any}()
    meta = Dict(
        KEY_CONSTRUCTOR => constructor_string(o),
        KEY_PROPERTYNAMES => collect(map(string, propertynames(o))),
    )
    d[KEY_META] = meta
    for prop in propertynames(o)
        key = string(prop)
        d[key] = marshal(getproperty(o, prop))
    end
    d
end

marshal(o) = marshal_struct(o)

SelfMarshalTypes = Union{
    String, Int, Float64,
    }
marshal(o::SelfMarshalTypes) = o
unmarshal(o::SelfMarshalTypes) = o

function unmarshal(d::Dict)
    @argcheck haskey(d, KEY_META)
    meta = d[KEY_META]
    @check haskey(meta, KEY_CONSTRUCTOR)
    @check haskey(meta, KEY_PROPERTYNAMES)

    ctor_str = meta[KEY_CONSTRUCTOR]
    @check ctor_str isa String
    ctor_expr = nothing
    try
        ctor_expr = Meta.parse(ctor_str)
    catch err
        msg = """Error parsing constructor:
        string: $ctor_str
        error: $err
        """
        error(msg)
    end
    ctor = nothing
    try
        ctor = Core.eval(Main, ctor_expr)
    catch err
        msg = """Error evaluating constructor expression:
        string: $ctor_str
        expr: $ctor_expr
        error: $err
        """
        error(msg)
    end
    propnames = meta[KEY_PROPERTYNAMES]
    args = map(propnames) do key
        unmarshal(d[key])
    end
    ctor(args...)
end

end # module


text:
using BlockDiagonals: getblock

#
# Matrix-Matrix product.
#

function LinearAlgebra.mul!(
    C::Matrix{T},
    A::BlockDiagonal{T},
    B::Matrix{T},
    α::T,
    β::T,
) where {T<:Real}
    @assert size(A, 1) == size(A, 2)
    start_row = 1

    @views for n in 1:nblocks(A)

        # Compute end row.
        (p, q) = BlockDiagonals.blocksize(A, n)
        @assert p == q
        end_row = start_row + p - 1

        # Multiply nth block of A by rows of B.
        mul!(C[start_row:end_row, :], getblock(A, n), B[start_row:end_row, :], α, β)

        # Update position.
        start_row += p
    end
    return C
end

function LinearAlgebra.mul!(
    C::Matrix{T},
    A::BlockDiagonal{T},
    B::Adjoint{T, Matrix{T}},
    α::T,
    β::T,
) where {T<:Real}
    @assert size(A, 1) == size(A, 2)
    start_row = 1

    @views for n in 1:nblocks(A)

        # Compute end row.
        (p, q) = BlockDiagonals.blocksize(A, n)
        @assert p == q
        end_row = start_row + p - 1

        # Multiply nth block of A by rows of B.
        mul!(C[start_row:end_row, :], getblock(A, n), B.parent[:, start_row:end_row]', α, β)

        # Update position.
        start_row += p
    end
    return C
end

function LinearAlgebra.mul!(
    C::Matrix{T},
    A::Matrix{T},
    B::BlockDiagonal{T},
    α::T,
    β::T,
) where {T<:Real}
    @assert size(B, 1) == size(B, 2)
    start_col = 1

    @views for n in 1:nblocks(B)

        # Compute end row.
        (p, q) = BlockDiagonals.blocksize(B, n)
        @assert p == q
        end_col = start_col + p - 1

        # Multiply nth col-block A by nth block of B.
        mul!(C[:, start_col:end_col], A[:, start_col:end_col], getblock(B, n), α, β)

        # Update postion.
        start_col += p
    end
    return C
end

function LinearAlgebra.mul!(
    C::Matrix{T},
    A::Adjoint{T, Matrix{T}},
    B::BlockDiagonal{T},
    α::T,
    β::T,
) where {T<:Real}
    @assert size(B, 1) == size(B, 2)
    start_col = 1

    @views for n in 1:nblocks(B)

        # Compute end row.
        (p, q) = BlockDiagonals.blocksize(B, n)
        @assert p == q
        end_col = start_col + p - 1

        # Multiply nth col-block A by nth block of B.
        mul!(C[:, start_col:end_col], A.parent[start_col:end_col, :]', getblock(B, n), α, β)

        # Update postion.
        start_col += p
    end
    return C
end

#
# Matrix-Vector product.
#

function LinearAlgebra.mul!(
    c::Vector{T},
    A::BlockDiagonal{T, Matrix{T}},
    b::Vector{T},
    α::T,
    β::T,
) where {T<:Real}
    @assert size(A, 1) == size(A, 2)
    start_row = 1

    @views for n in 1:nblocks(A)

        # Compute end row.
        (p, q) = BlockDiagonals.blocksize(A, n)
        @assert p == q
        end_row = start_row + p - 1

        # Multiply nth block of A by a block of rows of b.
        mul!(c[start_row:end_row], getblock(A, n), b[start_row:end_row], α, β)

        # Update position.
        start_row += p
    end
    return c
end


text:
# RGI-based runs for validation

using BSON, VAWTools, DataStructures
region = [3, 4, 14, 11, 7][end-1]
runtyp = ["test", "testmid", "prodlow", "prod"][3]
runtyp == :test && println("\n\nTEST RUN !!!!!!!!\n\n")
parallel = true
all_glaciers = parallel
save_geotiff = true
save_sols = false # makes huge files and did not work with NamedTuples...
fit_sigma = false # not implemented yet

today = Dates.today()
region_str = VAWTools.int2str2(region)

if runtyp in ["test", "testmid"]
    dir_output = "output/validation-test/$(lowercase(Dates.monthabbr(today)))-$(VAWTools.int2str2(Dates.day(today)))/"
else
    if fit_sigma
        dir_output = "output/validation-$(region_str)_fit-sigma/$(lowercase(Dates.monthabbr(today)))-$(VAWTools.int2str2(Dates.day(today)))/"
    else
        dir_output = "output/validation-$(region_str)/$(lowercase(Dates.monthabbr(today)))-$(VAWTools.int2str2(Dates.day(today)))/"
    end
end
!isdir(dir_output) && mkpath(dir_output)

@show region, runtyp, parallel, dir_output

include("rgi-setup.jl") # makes sure all is pre-compiles

tmp = BM.get_git_shas(error_on_dirty=false)
description = """
    bugfix with expectation value

    fsl uniform and 0.5

    Landmask for adjacent glacier now activated

    Validation on Svalbard

    Now after some major bug-fix in expectation-prior update calculation.

    BITEModel commit $(tmp["BITEModel"][1])
              branch $(tmp["BITEModel"][2])
    """

if parallel && nprocs()<2
    # if Sys.CPU_CORES>=17

    #     addprocs(17) # /2 to get to physical cores
    # else
    #     addprocs(Sys.CPU_CORES÷2 + 1) # /2 to get to physical cores
    # end
    # addprocs(Sys.CPU_CORES÷5) # five regions
    addprocs(Sys.CPU_CORES)
end

@everywhere begin
    const debugging = true
    include("rgi-setup.jl")
    using NamedTuples
    const progress_bar = true
end

update_cache = true
glnrs_with_radar = get_glaciers(region)[1] #[8:end] # drop the first few, too big
# delete the ones with a bad DEM
if region==7
    for x in [00025, 00028, 00037, 00042, 00228, 00246, 00279, 00409, 00474, 00529, 00531, 00570, 01425, 01466, 01481, 01523, 01559]
        i = findfirst(glnrs_with_radar, x)
        if i>0
            deleteat!(glnrs_with_radar, i)
        end
    end

    # also delete the ones which error (none found currently)

    # Svalbard glaciers with dhdt data as hand-selected in qgis
    svalbard_gls_with_dhdt = [133, 168, 170, 174, 175, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196,
                              197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 210, 213, 225, 266, 267, 268, 269, 270, 271, 272, 273, 274,
                              275, 276, 277, 278, 279, 280, 281, 284, 293, 294, 296, 313, 314, 315, 316, 325, 326, 327, 328, 329, 330, 331, 332, 333,
                              334, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358,
                              359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381,
                              382, 383, 384, 385, 386, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406,
                              407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 426, 427, 428, 429, 430, 431, 432, 433,
                              434, 435, 436, 437, 1403, 1414, 1415, 1461, 1467, 1468, 1469, 1470, 1475, 1476, 1478, 1480, 2, 3, 4, 6, 8, 9, 10, 11,
                              12, 13, 14, 15, 16, 17, 18, 19, 21, 22, 24, 86, 87, 88, 89, 90, 91, 93, 95, 96, 97, 100, 101, 102, 103, 104, 105, 106,
                              107, 108, 109, 112, 113, 114, 115, 116, 118, 119, 120, 121, 124, 125, 126, 128, 129, 130, 131, 132, 134, 135, 136, 137,
                              138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160,
                              161, 162, 163, 164, 165, 166, 167, 169, 172, 173, 176, 177, 179, 208, 209, 226, 227, 228, 229, 230, 231, 232, 233, 234
                              , 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 252, 253, 254, 255, 256, 257, 258, 259, 260,
                              261, 262, 263, 284, 295, 297, 301, 441, 442, 443, 446, 1402, 1406, 1419, 1420, 1421, 1422, 1423, 1438, 1439]
    # 39 have both:
    svalbard_gls_with_dhdt_and_radar = [s for s in svalbard_gls_with_dhdt if s in glnrs_with_radar]
elseif region==11
    # in the Alps not all have IV
    alps_gl_with_iv_n_radar = [1450, 1478, 1328, 2766, 1238, 2739, 1198, 2773, 2630, 1275, 2801, 2709, 1698, 2704, 2715, 2624, 1144, 2596,
                               2793, 2890, 2740, 2507, 2787, 2774, 2746, 1876, 2490, 1791, 1509, 2584, 2006, 2558, 2249, 2755, 1576, 1296, 2884, 2679, 1987, 1986, 1067,
                               2869, 2909, 1786, 1857, 2600, 1928, 2448, 2673, 1492, 2583, 2244, 1894, 2549, 2027, 1199, 2745, 1367, 1344, 1376]
end

glnrs_fit = glnrs_with_radar[1:2:end]
glnrs_test = glnrs_with_radar[2:2:end]
if length(glnrs_fit)>length(glnrs_test)
    # set to same length
    glnrs_fit = glnrs_fit[1:end-1]
end

# To reduce to just a few glacier
if !all_glaciers
    st = [17, 34, 100][2] # 5, 3, 1
    glnrs_fit = glnrs_fit[1:st:end]
    glnrs_test = glnrs_test[1:st:end]
else
    # st = 3
    # glnrs_fit = glnrs_fit[1:st:end][1:20]
    # glnrs_test = glnrs_test[1:st:end][1:20]
end

if region==7
    # 22:
    glnrs_fit_with_dhdt = [s for s in glnrs_fit if s in svalbard_gls_with_dhdt_and_radar]
    # 17:
    glnrs_test_with_dhdt = [s for s in glnrs_test if s in svalbard_gls_with_dhdt_and_radar]
elseif region==11
    # 30 each
    glnrs_fit_with_iv = [s for s in glnrs_fit if s in alps_gl_with_iv_n_radar]
    glnrs_test_with_iv = [s for s in glnrs_test if s in alps_gl_with_iv_n_radar]
end

## MCMC

# parameters to play with
sigma_of_model = @NT(sigma_h_model = 20.0,
                     sigma_iv_model = 20.0)

# number of divisions in elevation band variables
fit_vars = @NT(n_1d_btilde=3,
               n_1d_fsl=3,
               n_1d_temp=1)

@show region, debugging, runtyp, length(glnrs_fit)

# Method errors contain their arguments, which in our case can be huge.
# Remove them in case of error in pmap
@everywhere strip_MethodErrors(err) = err isa MethodError ? MethodError(0, err.f, err.world) : err

# Write description of model run
write(dir_output*"description.txt", description)

# Full fitting Mode 1
#####################

tic = time()

# runtyp=:prodlow
# nr = 2819
# out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(1), dir_output="", retsol=true)
# error()

# out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(1), dir_output="", retsol=true)
# out_fit = map(nr -> fit_it(nr, region, sigma_of_model, fit_vars,
#                         fit_target, runtyp, Val(1)), glnrs_fit)

expect_mean = Dict()
expect_median = Dict()
stdev = Dict()

out_ = 1 # to get "debug" access to below calculation in case of errors
for fit_target = [BM.FitTarget.h, BM.FitTarget.iv, BM.FitTarget.h_iv, BM.FitTarget.length]

    glnrs_fit_ = if region==11 && fit_target==BM.FitTarget.iv
        #  for the Alps only use the ones with IV
        glnrs_fit_with_iv
    else
        glnrs_fit
    end

    out_ = pmap(nr -> fit_it(nr, region, sigma_of_model, fit_vars,
                                fit_target, runtyp, Val(1),
                                dir_output=dir_output, save_geotiff=save_geotiff,
                                save_sol=false, retsol=save_sols, flroot="FIT_$(fit_target)",
                                fit_sigma=fit_sigma,
                                ),
                glnrs_fit_, on_error=strip_MethodErrors, retry_delays = zeros(3));
    out_fit = save_sols ? [ (o isa Exception) ? o : o[1] for o in out_] : out_;

    @assert length(out_fit)>0 "All fitting glaciers errored!"

    toc1 = time() - tic
    println("Iteration over fit glaciers $fit_target completed in $(toc1/60)min")

    # Make the priors of the fitted parameters
    ##########################################
    drop_sigma = true
    _, expect_median_, expect_mean_, stdev_ = calculate_prior_update(out_fit, drop_sigma)
    expect_mean[Symbol(fit_target)] = expect_mean_
    expect_median[Symbol(fit_target)] = expect_median_
    stdev[Symbol(fit_target)] = stdev_

    ## store
    fl = joinpath(dir_output, "reg$(region_str)_FIT_$(fit_target)")
    println("Storing $fl.bson ...")
    @time bson(fl*".bson", Dict(:out_fit => deepcopy(out_fit),
                                :description => description,
                                :commit => BM.get_git_shas(error_on_dirty=false),
                                :glnrs_fit => glnrs_fit_,
                                :glnrs_test => glnrs_test,
                                :runtyp => runtyp,
                                :expect_median => expect_median_,
                                :expect_mean => expect_mean_,
                                :stdev => stdev_,
                                ))

    println("... done.")
    if save_sols
        sols = [o[2] for o in out_ if !(o isa Exception)]

        println("Storing $fl.jls ...")
        sols_fit = OrderedDict(BM.getrgi_nr(o[:rgi])=>s for (o,s) in zip(out_fit,sols))
        open(fl*".jls", "w") do io
            serialize(io, sols_fit)
        end
        println("... done.")
    end

end
#  Base case: no change in prior
ln = length(first(expect_median)[2])
expect_median[:base] = zeros(ln)
expect_mean[:base] = zeros(ln)
stdev[:base] = Inf.*ones(ln)

fl = joinpath(dir_output, "reg$(region_str)_priors")
bson(fl*".bson", Dict(:expect_median => deepcopy(expect_median),
                      :expect_mean => deepcopy(expect_mean),
                      :stdev => deepcopy(stdev),
                      ))

# keys(expect_mean) == [:base, :h, :h_iv, :iv, :length]

# Validation: posterior1d only fitting
############################

tic2 = time()
# validation using updated priors both against test-set and fit-set (the latter is for double-checking)
for (fl_str, glnrs_) in zip(["testgls", "fitgls"], [glnrs_test, glnrs_fit])
    for k in keys(expect_median)
        fit_target = BM.FitTarget.length
        expect = expect_median[k]
        stdev_ = stdev[k]
        if k!=:base
            @assert expect!=nothing
            @assert stdev_!=nothing
        end

        #out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(2), expect, stdev_, dir_output=dir_output, retsol=true)
        out_ = pmap(nr -> fit_it(nr, region, sigma_of_model, fit_vars,
                                 fit_target, runtyp, Val(2),
                                 expect, stdev_,
                                 dir_output=dir_output, save_geotiff=save_geotiff,
                                 flroot="PRIOR_$k-$fl_str",
                                 fit_sigma=false,
                                 save_sol=false, retsol=save_sols),
                    glnrs_, on_error=strip_MethodErrors, retry_delays = zeros(3))
        out_test = save_sols ? [ (o isa Exception) ? o : o[1] for o in out_] : out_
        @assert length(out_test)>0 "All testing glaciers errored!"

        toc2 = time() - tic2

        println("Iteration over validation glaciers $k completed in $(toc2/60)min")

        fl = joinpath(dir_output, "reg$(region_str)_PRIOR_$k-$fl_str")
        println("Storing $fl.bson ...")
        @time bson(fl*".bson", Dict(:out_test => deepcopy(out_test),
                                    :description => description,
                                    :commit => BM.get_git_shas(error_on_dirty=false),
                                    :expect => expect,
                                    :stdev => stdev_,
                                    ))
        println("... done.")
        if save_sols
            println("Storing $fl.jls ...")
            sols = [o[2] for o in out_ if !(o isa Exception)]
            sols_test = OrderedDict(BM.getrgi_nr(o[:rgi])=>s for (o,s) in zip(out_fit,sols))
            open(fl*".jls", "w") do io
                serialize(io, sols_test)
            end
            println("... done.")
        end
    end
end

#####
# Last IV fitting with priors from h_iv and base (all else doesn't make sense)
fit_target = BM.FitTarget.iv
for k in [:base, :h_iv]
    # Here we need to add back the stripped sigma_iv term.  With stdev==Inf it
    # does not modify the prior.
    expect = [expect_median[k]..., 0]
    stdev_ = [stdev[k]..., Inf]

    glnrs_test_ = if region==11 && fit_target==BM.FitTarget.iv
        #  for the Alps only use the ones with IV
        glnrs_test_with_iv
    else
        glnrs_test
    end

    out_ = pmap(nr -> fit_it(nr, region, sigma_of_model, fit_vars,
                             fit_target, runtyp, Val(1),
                             expect, stdev_,
                             dir_output=dir_output, save_geotiff=save_geotiff,
                             flroot="PRIOR_$(k)_FIT_iv",
                             fit_sigma=fit_sigma,
                             save_sol=false, retsol=save_sols),
                glnrs_test_, on_error=strip_MethodErrors, retry_delays = zeros(3))
    out_test = save_sols ? [ (o isa Exception) ? o : o[1] for o in out_] : out_
    @assert length(out_test)>0 "All IV-testing glaciers errored!"

    toc1 = time() - tic

    fl = joinpath(dir_output, "reg$(region_str)_PRIOR_$(k)_FIT_$(fit_target)")
    println("Storing $fl.bson ...")
    @time bson(fl*".bson", Dict(:out_test => deepcopy(out_test),
                                :description => description,
                                :commit => BM.get_git_shas(error_on_dirty=false),
                                :glnrs_fit => glnrs_fit,
                                :glnrs_test => glnrs_test_,
                                :runtyp => runtyp,
                                :expect => expect,
                                :stdev => stdev_,
                                ))

    println("... done.")
    if save_sols
        sols = [o[2] for o in out_ if !(o isa Exception)]
        println("Storing $fl.jls ...")
        sols_test_iv = OrderedDict(BM.getrgi_nr(o[:rgi])=>s for (o,s) in zip(out_fit,sols))
        open(fl*".jls", "w") do io
            serialize(io, sols_test_iv)
        end
        println("... done.")
    end
end

toc = time() - tic
println("Total time $(toc/60/60) hours")


text:
@recipe function f(ta::T) where {T<:TimeArray}
    st = get(plotattributes, :seriestype, :path)
    if in(st, [:candlestick, :heikinashi])
        Candlestick(ta)
    #elseif st == :ohlc #ohlc (meaning sticks with steps on the sides) should be passed on to Plots internal ohlc plot engine
    #    ta, ohlc = extract_ohlc(ta)
    #    collect(zip(ohlc)) # But there are currently issues with that
    else
        label --> get(plotattributes, :label, reshape(String.(colnames(ta)),1,length(colnames(ta))))
        seriestype := st
        timestamp(ta), values(ta)
    end
end

# FIXME: refine Candlestick as a subtype of AbstractTimeSeries
#        it requires the Base.view supports. (#419)
"""
    Candlestick(ta::TimeArray)

# Argument

- There are four required columns from `ta::TimeArray`: `open`, `high`, `low` and
  `close`. The column names is case-insensitive.

# Examples

```julia-repl
julia> using MarketData

julia> TimeSeries.Candlestick(ohlcv)
```
"""
mutable struct Candlestick{D <: TimeType}
    time::Vector{D}
    open::AbstractVector
    high::AbstractVector
    low::AbstractVector
    close::AbstractVector
end

Candlestick(ta::TimeArray) = Candlestick(extract_ohlc(ta)...)

function extract_ohlc(ta::TimeArray)
    C     = ["open", "high", "low", "close"]
    cols  = colnames(ta)
    cols′ = lowercase.(string.(colnames(ta)))
    V     = map(C) do x
        i = findfirst(isequal(x), cols′)
        i ≡ nothing && throw(ArgumentError("the TimeArray did not have column `$x`"))
        values(ta[cols[i]])
    end
    (timestamp(ta), V...)
end

function HeikinAshi!(cs::Candlestick) # some values here are made too high!
    cs.close[1] = (cs.open[1] + cs.low[1] + cs.close[1] + cs.high[1]) / 4
    cs.open[1] = (cs.open[1] + cs.close[1]) / 2
    cs.high[1] = cs.high[1]
    cs.low[1] = cs.low[1]

    for i in 2:length(cs.open)
        cs.close[i] = (cs.open[i] + cs.low[i] + cs.close[i] + cs.high[i]) / 4
        cs.open[i] = (cs.open[i-1] + cs.close[i-1]) / 2
        cs.high[i] = maximum([cs.high[i], cs.open[i], cs.close[i]])
        cs.low[i] = minimum([cs.low[i], cs.open[i], cs.close[i]])
    end
end

@recipe function f(cs::Candlestick)
    st = get(plotattributes, :seriestype, :candlestick)
    st == :heikinashi && HeikinAshi!(cs)

    seriestype := :scatter #ignored
    legend --> false
    linewidth --> 0.7
    grid --> false

    bw = get(plotattributes, :bar_width, nothing)
    bw == nothing && (bw = 0.8)
    bar_width := bw / 2 * minimum(diff(unique(Dates.value.(cs.time))))

    # allow passing alternative colors as a vector
    cols = get(plotattributes, :seriescolor, nothing)
    cols = (isa(cols, Vector{Symbol}) && length(cols) == 2) ? cols : [:red, :blue]

    attributes = [
        Dict(:close_open => <,
             :close_prev => <,
             :bottombox => cs.close,
             :topbox => cs.open,
             :fill => cols[1],
             :line => cols[1],
             :fillalpha => 1),
        Dict(:close_open =>
             <, :close_prev => >=,
             :bottombox => cs.close,
             :topbox => cs.open,
             :fill => cols[2],
             :line => cols[2],
             :fillalpha => 1),
        Dict(:close_open => >=,
             :close_prev => <,
             :bottombox => cs.open,
             :topbox => cs.close,
             :fill => :white,
             :line => cols[1],
             :fillalpha => 0),
        Dict(:close_open => >=,
             :close_prev => >=,
             :bottombox => cs.open,
             :topbox => cs.close,
             :fill => :white,
             :line => cols[2],
             :fillalpha => 0)
    ]

    for att in attributes
        inds = similar(cs.close, Int)
        inds[1] = att[:close_open](cs.close[1], cs.open[1]) & att[:close_prev](cs.close[1], cs.close[1])
        @. inds[2:end] = att[:close_open](cs.close[2:end], cs.open[2:end]) & att[:close_prev]($diff(cs.close), 0)
        inds = findall(Bool.(inds))

        if length(inds) > 0
            @series begin
                linecolor := att[:line]
                fillcolor := att[:fill]
                fillalpha := att[:fillalpha]
                fillto := att[:bottombox][inds]
                seriestype := :bar
                cs.time[inds], att[:topbox][inds]
            end

            for j in 1:2
                @series begin
                    primary := false
                    linecolor := att[:line]
                    seriestype := :sticks
                    fillto := j == 1 ? cs.low[inds] : att[:topbox][inds]
                    cs.time[inds], j == 1 ? att[:bottombox][inds] : cs.high[inds]
                end
            end
        end
    end
end


text: f (1)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using CSV, DataFrames

function rename_csv_col_Order_1_to_Unif(pth)
    for (root, dirs, files) in walkdir(pth)
        for file in files
            fname, extension = splitext(file)
            if ((extension==".csv")&&(!occursin("sim",fname))&&!occursin("sample",fname))
                fpath = joinpath(root,file)
                df = CSV.read(fpath,DataFrame)
                try
                    rename!(df,Dict(:Order_1 => :Unif))
                    CSV.write(fpath,df)
                catch
                    @warn "file at "*fpath*" has col names "*join(names(df),", ")*". Did not overwrite"
                    @show names(df)
                end
            end
        end
    end
end

rename_csv_col_Order_1_to_Unif(".")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libspecfemUtils_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libspecfemUtils")
JLLWrappers.@generate_main_file("libspecfemUtils", UUID("a60b80ba-484f-5667-bf37-cf11dc804be2"))
end  # module libspecfemUtils_jll


text:
## Common Test Helpers. ##

meth(f, args) = first(methods(f, args))

fmeth(f) = first(methods(f))

function rawdocs(entries, m)
    m = entries[m]
    Docile.Interface.data(Docile.Interface.docs(m))
end

function docsmeta(entries, k, m)
    m = entries[m]
    Docile.Interface.metadata(m)[k]
end

macrofunc(mod, s) = getfield(mod, symbol(string("@", s)))

qs(mod, sym) = (mod, Docile.Collector.QualifiedSymbol(mod, sym))


text:
# Eval

const evalerr = :secret_lighttable_eval_err_keyword

# Data gets attached to the result in Julia
# Info gets attached to the result in LT
function eval(editor, mod, code, file, bounds)
  task_local_storage()[:SOURCE_PATH] = file
  file == nothing && (file = "REPL")
  try
    result = include_string(mod, code, file, bounds[1])
    Jewel.isdefinition(code) && (result = nothing)
    return result
  catch e
    showexception(editor, isa(e, LoadError)?e.error:e, catch_backtrace(), bounds)
    return evalerr
  end
end

function evaldisplay(editor, mod, code, file, bounds; data = Dict(), info = Dict())
  try
    result = eval(editor, mod, code, file, bounds)
    result === evalerr && return
    withcurrentresult(register_result(result, data)) do
      displayinline!(result, @d(:editor => editor,
                                :bounds => bounds,
                                :info => merge(info, @d(:id => string(_currentresult_.id)))))
    end
  catch e
    showexception(editor, e, catch_backtrace(), bounds)
  end
end

noselection(data) = data["start"] == data["end"]

LNR.cursor(data::Dict) = cursor(data["line"], data["col"])

handle("eval.selection") do editor, data
  code, bounds =
    noselection(data) ?
      Jewel.getblock(data["code"], data["start"]["line"]) :
      Jewel.getblock(data["code"], cursor(data["start"]), cursor(data["end"]))
  code == "" && return notify_done()
  mod = Jewel.getmodule(data["code"], bounds[1], filemod = data["module"])
  evaldisplay(editor, mod, code, data["path"], bounds)
end

handle("eval.block") do editor, data
  code = data["block"]
  bounds = data["bounds"]
  mod = Jewel.getmodule(data["code"], bounds[1], filemod = data["module"])
  # We need some custom data to enable reevaluation
  evaldisplay(editor, mod, code, data["path"], bounds,
              data = @d(:editor => editor,
                        :mod => mod,
                        :bounds => bounds,
                        :code => code,
                        :path => data["path"]),
              info = @d(:scales => bounds))
end

handle("eval.all") do editor, data
  file = @or data["path"] "REPL"
  code = data["code"]
  mod = Jewel.getthing(data["module"], Main)
  try
    task_local_storage()[:SOURCE_PATH] = data["path"]
    include_string(mod, code, file)
  catch e
    println(STDERR, "Error evaluating $file:")
    showerror(STDERR, e.error, catch_backtrace())
    println(STDERR)
    notify_error(sprint(showerror, e.error))
  end
  notify_done()
end

handle("editor.block") do editor, data
  block, bounds = Jewel.getblock(data["code"], data["line"])
  raise(editor, "return-block", @d("block" => block,
                                   "bounds" => bounds,
                                   "id" => data["id"]))
end

# Reevaluation

handle("result.reval") do _, data
  withcurrentresult(data["id"]) do

    jscall("lt.plugins.reptile.working();")

    code = fillranges(_currentresult_.data[:code],
                      data["vals"],
                      data["locs"],
                      _currentresult_.data[:bounds][1])
    result = eval(_currentresult_.data[:editor],
                  _currentresult_.data[:mod],
                  code,
                  _currentresult_.data[:path],
                  _currentresult_.data[:bounds])

    local html
    try
      html = stringmime(MIME"text/html"(), applydisplayinline(result))
    catch e
      return showexception(_currentresult_.data[:editor], e, catch_backtrace(), _currentresult_.data[:bounds])
    end

    jscall("""
      this.querySelector('.julia.result').innerHTML = '$(jsescapestring(html))';
      lt.plugins.reptile.done_working();
    """)

    queuecmds()
    while length(cmdqueue) > 1 && cmdqueue[1][2] == cmdqueue[2][2] == "result.reval"
      shift!(cmdqueue)
    end
  end
end


text:
using LinearAlgebra
using FermiCG
using Printf
using Test
using LinearMaps
using Arpack
using Random
using Profile

#@testset "fci_matvec" begin
    atoms = []
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,0,1]))
    push!(atoms,Atom(3,"H",[1,0,2]))
    push!(atoms,Atom(4,"H",[1,0,3]))
    push!(atoms,Atom(5,"H",[2,0,4]))
    push!(atoms,Atom(6,"H",[2,0,5]))
    push!(atoms,Atom(7,"H",[3,0,6]))
    push!(atoms,Atom(8,"H",[3,0,7]))
    push!(atoms,Atom(9,"H",[0,0,8]))
    push!(atoms,Atom(10,"H",[0,0,9]))
    push!(atoms,Atom(11,"H",[0,0,10]))
    push!(atoms,Atom(12,"H",[0,0,11]))
    #push!(atoms,Atom(13,"H",[0,0,12]))
    #push!(atoms,Atom(14,"H",[0,0,13]))
    basis = "6-31g"
    basis = "sto-3g"

    mol     = Molecule(0,1,atoms,basis)
    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));

    na = 6
    nb = 6

    e_mf = mf.e_tot - mf.energy_nuc()
    if 1==0
        @printf(" Mean-field energy %12.8f", e_mf)
        @time e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,na,na)
        # @printf(" FCI Energy: %12.8f\n", e_fci)
    end

    norbs = size(ints.h1)[1]

    problem = StringCI.FCIProblem(norbs, 1, 1)
    display(problem)
    nr = 1
    v0 = rand(problem.dim,nr)
    #v0[:,1] .= 0
    v0[1,1] = 1
    v0 = v0 * inv(sqrt(v0'*v0))

    Hmap = StringCI.get_map(ints, problem)
    Random.seed!(3);
    A = Diagonal(rand(20)) + .0001*rand(20,20)
    A = A'+A

    function test_matvec(v,prb,n)
        ket_a = FermiCG.StringCI.DeterminantString(prb.no, prb.na)
        ket_b = FermiCG.StringCI.DeterminantString(prb.no, prb.nb)

        lookup_a = FermiCG.StringCI.fill_ca_lookup2(ket_a)
        lookup_b = FermiCG.StringCI.fill_ca_lookup2(ket_b)

        s = similar(v)
        for i=1:n
            println(i)
            flush(stdout)
            v = reshape(v,ket_a.max,ket_b.max, 1)
            #@time s = FermiCG.StringCI.compute_ab_terms2(v,ints,prb,lookup_a, lookup_b)
            #@profilehtml s = FermiCG.StringCI.compute_ab_terms2(v,ints,prb,lookup_a, lookup_b)
            #@time s = FermiCG.StringCI.compute_ss_terms2(v,ints,prb,lookup_a, lookup_b)
            @time s = FermiCG.StringCI.compute_ss_terms2(v,ints,prb,lookup_a, lookup_b)
            #s = Matrix(H*v)
        end
    end

    test_matvec(v0,problem,1)

    problem = StringCI.FCIProblem(norbs, na, nb)
    display(problem)
    nr = 1
    v0 = rand(problem.dim,nr)
    #v0[:,1] .= 0
    v0[1,1] = 1
    v0 = v0 * inv(sqrt(v0'*v0))

    Random.seed!(3);
    A = Diagonal(rand(20)) + .0001*rand(20,20)
    A = A'+A
    @time test_matvec(v0,problem,3)

#end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# This file is for different constraints

abstract type AbstractLinearProgram end
struct StandardLP            <: AbstractLinearProgram end
struct LinearRelaxedLP       <: AbstractLinearProgram end
struct TriangularRelaxedLP   <: AbstractLinearProgram end
struct BoundedMixedIntegerLP <: AbstractLinearProgram end
struct SlackLP <: AbstractLinearProgram
    slack::Vector{Vector{VariableRef}}
end
SlackLP() = SlackLP([])
struct MixedIntegerLP <: AbstractLinearProgram
    m::Float64
end

# Any encoding passes through here first:
function encode_network!(model::Model,
                         network::Network,
                         zs::Vector{Vector{VariableRef}},
                         δs::Vector,
                         encoding::AbstractLinearProgram)

    for (i, layer) in enumerate(network.layers)
        encode_layer!(encoding, model, layer, zs[i], zs[i+1], δs[i])
    end
    return encoding # only matters for SlackLP
end

# TODO: find a way to eliminate the two methods below.
# i.e. make BoundedMixedIntegerLP(bounds) and TriangularRelaxedLP(bounds) or something
function encode_network!(model::Model,
                         network::Network,
                         zs::Vector{Vector{VariableRef}},
                         δs::Vector,
                         bounds::Vector{Hyperrectangle},
                         encoding::AbstractLinearProgram)

    for (i, layer) in enumerate(network.layers)
        encode_layer!(encoding, model, layer, zs[i], zs[i+1], δs[i], bounds[i])
    end
    return encoding # only matters for SlackLP
end

function encode_network!(model::Model,
                         network::Network,
                         zs::Vector{Vector{VariableRef}},
                         bounds::Vector{Hyperrectangle},
                         encoding::AbstractLinearProgram)

    for (i, layer) in enumerate(network.layers)
        encode_layer!(encoding, model, layer, zs[i], zs[i+1], bounds[i])
    end
    return encoding # only matters for SlackLP
end

# For an Id Layer, any encoding type defaults to this:
function encode_layer!(::AbstractLinearProgram,
                       model::Model,
                       layer::Layer{Id},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       args...)
    @constraint(model, z_next .== layer.weights*z_current + layer.bias)
end

# SlackLP is slightly different, because we need to keep track of the slack variables
function encode_layer!(SLP::SlackLP,
                       model::Model,
                       layer::Layer{Id},
                       z1::Array{VariableRef,1},
                       z2::Array{VariableRef,1},
                       δ...)

    encode_layer!(StandardLP(), model, layer, z1, z2)
    # We need identity layer slack variables so that the algorithm doesn't
    # "get confused", but they are set to 0 because they're not relevant
    slack_vars = @variable(model, [1:n_nodes(layer)])
    @constraint(model, slack_vars .== 0.0)
    push!(SLP.slack, slack_vars)
    return nothing
end

# alternative signature:
# encode_layer!(encoding, model::Model, current_layer::VarLayer{ReLU},  next_layer::VarLayer)
function encode_layer!(::StandardLP,
                       model::Model,
                       layer::Layer{ReLU},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       δ::Vector{Bool})

    # The jth ReLU is forced to be active or inactive,
    # depending on the activation pattern given by δᵢ.
    # δᵢⱼ == true denotes ẑ >=0 (i.e. an *inactive* ReLU)

    ẑ = layer.weights * z_current + layer.bias
    for j in 1:length(layer.bias)
        if δ[j]
            @constraint(model, ẑ[j] >= 0.0)
            @constraint(model, z_next[j] == ẑ[j])
        else
            @constraint(model, ẑ[j] <= 0.0)
            @constraint(model, z_next[j] == 0.0)
        end
    end
end

function encode_layer!(SLP::SlackLP,
                       model::Model,
                       layer::Layer{ReLU},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       δ::Vector{Bool})

    ẑ = layer.weights * z_current + layer.bias
    slack_vars = @variable(model, [1:length(layer.bias)])
    for j in 1:length(layer.bias)
        if δ[j]
            @constraint(model, z_next[j] == ẑ[j] + slack_vars[j])
            @constraint(model, ẑ[j] + slack_vars[j] >= 0.0)
        else
            @constraint(model, z_next[j] == slack_vars[j])
            @constraint(model, 0.0 >= ẑ[j] - slack_vars[j])
        end
    end
    push!(SLP.slack, slack_vars)
    return nothing
end

function encode_layer!(::LinearRelaxedLP,
                       model::Model,
                       layer::Layer{ReLU},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       δ::Vector{Bool})

    ẑ = layer.weights * z_current + layer.bias
    for j in 1:length(layer.bias)
        if δ[j]
            @constraint(model, z_next[j] == ẑ[j])
        else
            @constraint(model, z_next[j] == 0.0)
        end
    end
end

function encode_layer!(::TriangularRelaxedLP,
                       model::Model,
                       layer::Layer{ReLU},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       bounds::Hyperrectangle)

    ẑ = layer.weights * z_current + layer.bias
    ẑ_bound = approximate_affine_map(layer, bounds)
    l̂, û = low(ẑ_bound), high(ẑ_bound)
    for j in 1:length(layer.bias)
        if l̂[j] > 0.0
            @constraint(model, z_next[j] == ẑ[j])
        elseif û[j] < 0.0
            @constraint(model, z_next[j] == 0.0)
        else
            @constraints(model, begin
                                    z_next[j] >= ẑ[j]
                                    z_next[j] <= û[j] / (û[j] - l̂[j]) * (ẑ[j] - l̂[j])
                                    z_next[j] >= 0.0
                                end)
        end
    end
end

function encode_layer!(MIP::MixedIntegerLP,
                       model::Model,
                       layer::Layer{ReLU},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       δ::Vector{VariableRef})
    m = MIP.m

    ẑ = layer.weights * z_current + layer.bias
    for j in 1:length(layer.bias)
        @constraints(model, begin
                                z_next[j] >= ẑ[j]
                                z_next[j] >= 0.0
                                z_next[j] <= ẑ[j] + m * δ[j]
                                z_next[j] <= m - m * δ[j]
                            end)
    end
end

function encode_layer!(::BoundedMixedIntegerLP,
                       model::Model,
                       layer::Layer{ReLU},
                       z_current::Vector{VariableRef},
                       z_next::Vector{VariableRef},
                       δ::Vector,
                       bounds::Hyperrectangle)

    ẑ = layer.weights * z_current + layer.bias
    ẑ_bound = approximate_affine_map(layer, bounds)
    l̂, û = low(ẑ_bound), high(ẑ_bound)

    for j in 1:length(layer.bias) # For evey node
        if l̂[j] >= 0.0
            @constraint(model, z_next[j] == ẑ[j])
        elseif û[j] <= 0.0
            @constraint(model, z_next[j] == 0.0)
        else
            @constraints(model, begin
                                    z_next[j] >= ẑ[j]
                                    z_next[j] >= 0.0
                                    z_next[j] <= û[j] * δ[j]
                                    z_next[j] <= ẑ[j] - l̂[j] * (1 - δ[j])
                                end)
        end
    end
end

#=
Add input/output constraints to model
=#
function add_complementary_set_constraint!(model::Model, output::HPolytope, z::Vector{VariableRef})
    out_A, out_b = tosimplehrep(output)
    # Needs to take the complementary of output constraint
    n = length(constraints_list(output))
    if n == 1
        # Here the output constraint is a half space
        halfspace = first(constraints_list(output))
        add_complementary_set_constraint!(model, halfspace, z)
    else
        error("Non-convex constraints are not supported. Please make sure that the
            output set is a HalfSpace (or an HPolytope with a single constraint) so that the
            complement of the output is convex. Got $n constraints.")
    end
    return nothing
end

function add_complementary_set_constraint!(m::Model, H::HalfSpace, z::Vector{VariableRef})
    a, b = tosimplehrep(H)
    @constraint(m, a * z .>= b)
    return nothing
end
function add_complementary_set_constraint!(m::Model, PC::PolytopeComplement, z::Vector{VariableRef})
    add_set_constraint!(m, PC.P, z)
    return nothing
end

function add_set_constraint!(m::Model, set::Union{HPolytope, HalfSpace}, z::Vector{VariableRef})
    A, b = tosimplehrep(set)
    @constraint(m, A * z .<= b)
    return nothing
end

function add_set_constraint!(m::Model, set::Hyperrectangle, z::Vector{VariableRef})
    @constraint(m, z .<= high(set))
    @constraint(m, z .>= low(set))
    return nothing
end

function add_set_constraint!(m::Model, PC::PolytopeComplement, z::Vector{VariableRef})
    add_complementary_set_constraint!(m, PC.P, z)
    return nothing
end


text:
function make_batch(o, samples, vocab)
    images = make_images_batch(o, map(s->s[1], samples))
    captions, masks = make_captions_batch(o, map(s->s[2], samples), vocab)
    return images, captions, masks
end

function make_images_batch(o, filenames)
    images = h5open(o[:images], "r") do f
        mapreduce(
            x->o[:finetune]?x:reshape(x,1,prod(size(x,1,2)),size(x,3)),
            (x...)->cat(o[:finetune]?4:1, x...),
            map(x->read(f,x), filenames))
    end

    batch = nothing
    if o[:finetune]
        batch = zeros(typeof(images[1]), size(images,1:3...)..., o[:batchsize])
        batch[:,:,:,1:length(filenames)] = images
    else
        # info(size(images))
        batch = zeros(typeof(images[1]), o[:batchsize], size(images)[2:end]...)
        batch[1:length(filenames),:,:] = images
    end

    return batch
end

function make_captions_batch(o, tokens, vocab)
    # captions batch
    vectors = map(t->sen2vec(vocab, t), tokens)
    longest = mapreduce(length, max, vectors)
    pad = word2index(vocab, PAD)
    captions = map(i -> pad * ones(Int, o[:batchsize]), [1:longest...])
    masks = map(i -> falses(o[:batchsize]), [1:(longest-1)...])
    for i = 1:length(tokens)
        map!(t->captions[t][i] = vectors[i][t], [1:length(vectors[i])...])
        map!(t->masks[t][i] = 1, [1:(length(vectors[i])-1)...])
    end
    return captions, masks
end

max_stars_repo_name: ilkerkesen/show-attend-and-tell

max_issues_repo_name: ilkerkesen/show-attend-and-tell

max_forks_repo_name: ilkerkesen/show-attend-and-tell


text:
module Electron

using JSON, URIParser, Sockets, Base64, Pkg.Artifacts, FilePaths, UUIDs

export Application, Window, URI, windows, applications, msgchannel, toggle_devtools, load, ElectronAPI

function conditional_electron_load()
    try
        return artifact"electronjs_app"
    catch error
        return nothing
    end
end

const electronjs_path = conditional_electron_load()

function prep_test_env()
    if haskey(ENV, "GITHUB_ACTIONS") && ENV["GITHUB_ACTIONS"] == "true"
        if Sys.islinux()
            run(Cmd(`Xvfb :99 -screen 0 1024x768x24`), wait=false)
            ENV["DISPLAY"] = ":99"
        end
    end
end

const OptDict = Dict{String, Any}

struct JSError
    msg
end
Base.showerror(io::IO, e::JSError) = print(io, "JSError: ", e.msg)

mutable struct _Application{T} # forward declaration of Application
    connection::IO
    proc
    secure_cookie::Vector{UInt8}
    windows::Vector{T}
    exists::Bool

    global function _Application(::Type{T}, connection::IO, proc, secure_cookie) where {T} # internal constructor
        new_app = new{T}(connection, proc, secure_cookie, T[], true)
        push!(_global_applications, new_app)
        return new_app
    end
end

mutable struct Window
    app::_Application{Window}
    id::Int64
    exists::Bool
    msg_channel::Channel{Any}

    global function _Window(app::_Application{Window}, id::Int64; msg_channel_size=128) # internal constructor
        new_window = new(app, id, true, Channel{Any}(msg_channel_size))
        push!(app.windows, new_window)
        return new_window
    end
end

const Application = _Application{Window}

function Base.show(io::IO, app::Application)
    if app.exists
        if length(app.windows) == 1
            appstate = ", [1 window])"
        else
            appstate = ", [$(length(app.windows)) windows])"
        end
    else
        appstate = ", [dead])"
    end
    print(io, "Application(", app.connection, ", ", app.proc, appstate)
end

const _global_applications = Vector{Application}(undef,0)
const _global_default_application = Ref{Union{Nothing,Application}}(nothing)

function __init__()
    atexit() do # let Electron know we want it to die quietly and sanely
        for app in _global_applications
            if app.exists
                close(app)
            end
        end
    end
    nothing
end

function applications()
    return _global_applications
end

function default_application()
    if _global_default_application[]===nothing || _global_default_application[].exists==false
        _global_default_application[] = Application()
    end

    return _global_default_application[]
end

function windows(app::Application)
    return app.windows
end

function generate_pipe_name(name)
    return if Sys.iswindows()
        "\\\\.\\pipe\\$name"
    elseif Sys.isunix()
        joinpath(tempdir(), name)
    end
end

function get_electron_binary_cmd()
    if electronjs_path===nothing
        return "electron"
    elseif Sys.isapple()
        return joinpath(electronjs_path, "Julia.app", "Contents", "MacOS", "Julia")
    elseif Sys.iswindows()
        return joinpath(electronjs_path, "electron.exe")
    else # assume unix layout
        return joinpath(electronjs_path, "electron")
    end
end

"""
    function Application()

Start a new Electron application. This will start a new process
for that Electron app and return an instance of `Application` that
can be used in the construction of Electron windows.
"""
function Application()
    electron_path = get_electron_binary_cmd()
    mainjs = joinpath(@__DIR__, "main.js")

    id = replace(string(uuid1()), "-"=>"")
    main_pipe_name = generate_pipe_name("jlel-$id")
    server = listen(main_pipe_name)

    id = replace(string(uuid1()), "-"=>"")
    sysnotify_pipe_name = generate_pipe_name("jlel-sn-$id")
    sysnotify_server = listen(sysnotify_pipe_name)

    secure_cookie = rand(UInt8, 128)
    secure_cookie_encoded = base64encode(secure_cookie)
    # proc = open(`$electron_path --inspect-brk=5858 $mainjs $main_pipe_name $sysnotify_pipe_name $secure_cookie_encoded`, "w", stdout)
    proc = open(`$electron_path $mainjs $main_pipe_name $sysnotify_pipe_name $secure_cookie_encoded`, "w", stdout)

    sock = accept(server)
    if read!(sock, zero(secure_cookie)) != secure_cookie
        close(server)
        close(sysnotify_server)
        close(sock)
        error("Electron failed to authenticate with the proper security token")
    end

    let sysnotify_sock = accept(sysnotify_server)
        if read!(sysnotify_sock, zero(secure_cookie)) != secure_cookie
            close(server)
            close(sysnotify_server)
            close(sysnotify_sock)
            close(sock)
            error("Electron failed to authenticate with the proper security token")
        end
        let app = _Application(Window, sock, proc, secure_cookie)
            @async begin
                try
                    try
                        while true
                            try
                                line_json = readline(sysnotify_sock)
                                isempty(line_json) && break # EOF
                                cmd_parsed = JSON.parse(line_json)
                                if cmd_parsed["cmd"] == "windowclosed"
                                    win_index = findfirst(w -> w.id == cmd_parsed["winid"], app.windows)
                                    app.windows[win_index].exists = false
                                    close(app.windows[win_index].msg_channel)
                                    deleteat!(app.windows, win_index)
                                elseif cmd_parsed["cmd"] == "appclosing"
                                    break
                                elseif cmd_parsed["cmd"] == "msg_from_window"
                                    win_index = findfirst(w -> w.id == cmd_parsed["winid"], app.windows)
                                    put!(app.windows[win_index].msg_channel, cmd_parsed["payload"])
                                end
                            catch er
                                bt = catch_backtrace()
                                io = PipeBuffer()
                                print_with_color(Base.error_color(), io, "Electron ERROR: "; bold = true)
                                Base.showerror(IOContext(io, :limit => true), er, bt)
                                println(io)
                                write(stderr, io)
                            end
                        end
                    finally
                        # Cleanup all the windows that are associated with this application
                        for w in app.windows
                            w.exists = false
                        end
                        empty!(app.windows)
                    end
                finally
                    # Cleanup the application instance
                    app.exists = false
                    close(sysnotify_sock)
                    app_index = findfirst(a -> a === app, _global_applications)
                    deleteat!(_global_applications, app_index)
                end
            end
            return app
        end
    end
end

"""
    close(app::Application)

Terminates the Electron application referenced by `app`.
"""
function Base.close(app::Application)
    app.exists || error("Cannot close this application, the application does no longer exist.")
    while length(windows(app))>0
        close(first(windows(app)))
    end
    app.exists = false
    close(app.connection)
end

function req_response(app::Application, cmd)
    connection = app.connection
    json = JSON.json(cmd)
    c = Condition()
    t = @async try
        println(connection, json)
        fetch(c)
    catch ex
        close(connection) # kill Application, since it probably must be in a bad state now
        rethrow(ex)
    end
    retval_json = readline(connection)
    notify(c)
    fetch(t)
    return JSON.parse(retval_json)
end

"""
    run(app::Application, code::AbstractString)

Run the JavaScript code that is passed in `code` in the main
application thread of the `app` Electron process. Returns the
value that the JavaScript expression returns.
"""
Base.run(app::Application, code::AbstractString) = run(app, String(code))
function Base.run(app::Application, code::String)
    app.exists || error("Cannot run code in this application, the application does no longer exist.")
    message = OptDict("cmd" => "runcode", "target" => "app", "code" => code)
    retval = req_response(app, message)
    haskey(retval, "error") && throw(JSError(retval["error"]))
    return retval["data"]
end

"""
    run(win::Window, code::AbstractString)

Run the JavaScript code that is passed in `code` in the render
thread of the `win` Electron windows. Returns the value that
the JavaScript expression returns.
"""
function Base.run(win::Window, code::AbstractString)
    win.exists || error("Cannot run code in this window, the window does no longer exist.")
    message = OptDict("cmd" => "runcode", "target" => "window", "winid" => win.id, "code" => code)
    retval = req_response(win.app, message)
    if haskey(retval, "data")
        return retval["data"]
    else
        @assert haskey(retval, "error")
        error("JSError: $(JSON.json(retval["error"]))")
    end
end

"""
    load(win::Window, uri::URI)

Load `uri` in the Electron window `win`.
"""
function load(win::Window, uri::URI)
    win.exists || error("Cannot load URI in this window, the window does no longer exist.")
    message = OptDict("cmd" => "loadurl", "winid" => win.id, "url" => string(uri))
    req_response(win.app, message)
    return nothing
end

"""
    load(win::Window, path::AbstractPath)

Load `path` in the Electron window `win`.
"""
function load(win::Window, path::AbstractPath)
    win.exists || error("Cannot load path in this window, the window does no longer exist.")
    message = OptDict("cmd" => "loadurl", "winid" => win.id, "url" => string(URI(path)))
    req_response(win.app, message)
    return nothing
end

"""
    load(win::Window, html::AbstractString)

Load `html` in the Electron window `win`.
"""
load(win::Window, html::AbstractString) =
    load(win, URI("data:text/html;charset=utf-8," * escape(html)))

"""
    function Window([app::Application,] options::Dict)

Open a new Window in the application `app`. Pass the content
of `options` to the Electron `electron.BrowserWindow` constructor.

If `app` is not specified, use the default Electron application,
starting one if needed.
"""
function Window(app::Application, options::Dict=OptDict())
    message = OptDict("cmd" => "newwindow", "options" => options)
    retval = req_response(app, message)
    ret_val = retval["data"]
    return _Window(app, ret_val)
end

"""
    function Window([app::Application,] uri::URI)

Open a new Window in the application `app`. Show the content
that `uri` points to in that new window.

If `app` is not specified, use the default Electron application,
starting one if needed.
"""
function Window(app::Application, uri::URI; options::Dict=OptDict())
    internal_options = OptDict()
    merge!(internal_options, options)
    internal_options["url"] = string(uri)
    return Window(app, internal_options)
end

"""
    function Window([app::Application,] path::AbstractPath)

Open a new Window in the application `app`. Show the content
that `path` points to in that new window.

If `app` is not specified, use the default Electron application,
starting one if needed.
"""
function Window(app::Application, path::AbstractPath; options::Dict=OptDict())
    return Window(app, URI(path); options=options)
end

"""
    function Window([app::Application,] content::AbstractString)

Open a new Window in the application `app`. Show the `content`
as a text/html file with utf-8 encoding.

If `app` is not specified, use the default Electron application,
starting one if needed.
"""
function Window(app::Application, content::AbstractString; kwargs...)
    return Window(app, URI("data:text/html;charset=utf-8," * escape(content)); kwargs...)
end

Window(a1::Application, args...; kwargs...) = throw(MethodError(Window, (a1, args...)))
Window(args...; kwargs...) = Window(default_application(), args...; kwargs...)

function toggle_devtools(w::Window)
    run(w.app, "electron.BrowserWindow.fromId($(w.id)).webContents.toggleDevTools()")
end

"""
    close(win::Window)

Close the windows referenced by `win`.
"""
function Base.close(win::Window)
    win.exists || error("Cannot close this window, the window does no longer exist.")
    message = OptDict("cmd" => "closewindow", "winid" => win.id)
    retval = req_response(win.app, message)
    return nothing
end

msgchannel(win::Window) = win.msg_channel

"""
    ElectronAPI

A shim object for calling Electron API functions.

See:
* <https://electronjs.org/docs/api/browser-window>

# Examples
```jldoctest
julia> using Electron

julia> win = Window();

julia> ElectronAPI.setBackgroundColor(win, "#000");

julia> ElectronAPI.show(win);
```
"""
ElectronAPI

struct ElectronAPIType end
const ElectronAPI = ElectronAPIType()

struct ElectronAPIFunction <: Function
    name::Symbol
end

Base.getproperty(::ElectronAPIType, name::Symbol) = ElectronAPIFunction(name)

function (api::ElectronAPIFunction)(w::Window, args...)
    name = api.name
    json_args = JSON.json(collect(args))
    run(w.app, "electron.BrowserWindow.fromId($(w.id)).$name(...$json_args)")
end

end


text:
println("First line of the script")
time_start = time()
using Test
using Logging
using SparseArrays
using ExSup
using Pio3d
using JLD2
@info "Before Distributed"
using Distributed
if nprocs() == 1
    addprocs(3)
end
# using Revise
# using ExSup
using LarSurf
# @everywhere using LarSurf

fn = "exp_surface_extraction_ircad_times.csv"
data = Dict()

@info "before everywhere using"
@info "time from start: $(time()-time_start) [s]"
# using Plasm
# @everywhere using LarSurf
# @everywhere using Distributed

@info "after everywhere using, time from start: $(time()-time_start) [s]"

data_id = 1
stepz = 1
stepxy = 1
# stepxy = 4
block_size = [64, 64, 64]
# block_size = [128, 128, 128]
# block_size = [128, 128, 128]
# block_size = [32, 32, 32]
# data_size1 = 128
# data_size1 = 256
# data_size1 = 512

LarSurf.set_time_data(data)
LarSurf.report_init_row(@__FILE__)

data["block size"] = block_size[1]
data["using done"] = time()-time_start
# segmentation = LarSurf.data234()
@info "Generate data..."
@info "time from start: $(time() - time_start) [s]"
mask_labels=["liver", "portalvein"]
# pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.$data_id/MASKS_DICOM/liver")

for mask_label in mask_labels
	pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.$data_id/MASKS_DICOM/$(mask_label)")
	datap = Pio3d.read3d(pth)
	data3d_full = datap["data3d"]
	data3d_full = data3d_full[1:stepz:end, 1:stepxy:end, 1:stepxy:end]
	    # round(size(data3d_full, 1) / target_size1)
	    # return data3d_full
	segmentation = convert(Array{Int8, 3}, data3d_full .> 0)
	voxelsize_mm = datap["voxelsize_mm"]
	voxelsize_mm[1] = voxelsize_mm[1] * stepz
	voxelsize_mm[2] = voxelsize_mm[2] * stepxy
	voxelsize_mm[3] = voxelsize_mm[3] * stepxy
	LarSurf.Experiments.report_add_data_info(data, segmentation, voxelsize_mm)
	@info "voxelsize mm = $(voxelsize_mm), size=$(size(data3d_full))"
	# segmentation = LarSurf.generate_cube(data_size1; remove_one_pixel=true)
	@info "==== using done, data generated time from start: $(time()-time_start) [s]"
	data["data generated"] = time()-time_start

	@info "Setup..."
	setup_time = @elapsed LarSurf.lsp_setup(block_size;reference_time=time_start)
	println("setup time: $setup_time")
	@info "==== setup done, time from start: $(time()-time_start) [s]"
	data["setup done"] = time()-time_start
	# for wid in workers()
	#     # println("testing on $wid")
	#     ftr = @spawnat wid LarSurf._single_boundary3
	#     @test fetch(ftr) != nothing
	# end

	# @debug "Setup done"
	tmd = @timed larmodel = LarSurf.lsp_get_surface(segmentation; voxelsize=voxelsize_mm)
	val, tm, mem, gc = tmd
	println("Total time per $(mask_label): $tm")
	@info "==== time from start: $(time()-time_start) [s]"
	data["finished"] = time()-time_start
	ExSup.add_to_csv(data, fn)

	V, FV = larmodel
	FVtri = LarSurf.triangulate_quads(FV)

	@JLD2.save "ircad_$(mask_label).jld2" V FV
	using ViewerGL
	println("=== ViewerGL readed")
	ViewerGL.VIEW([
	    ViewerGL.GLGrid(V,FVtri,ViewerGL.Point4d(1,1,1,0.1))
		ViewerGL.GLAxis(ViewerGL.Point3d(-1,-1,-1),ViewerGL.Point3d(1,1,1))
	])
	# FV = FVtri
	@JLD2.save "ircad_$(mask_label)_tri.jld2" V FVtri
	# @JLD2.save "liver01tri.jld2" V FVtri
end
# Plasm.view(val)


text:
module TestDecorate
export @syms, @noenv, @test_decorate
using Compat, Debug.AST, Debug.Meta, Debug.Analysis

macro assert_fails(ex)
    quote
        if (try; $(esc(ex)); true; catch e; false; end)
            error($("@assert_fails: $ex didn't fail!"))
        end
    end
end

type BlockEnv
    defined ::Set{Symbol}
    assigned::Set{Symbol}

    BlockEnv(d, a) = new(Set{Symbol}(d), Set{Symbol}(a))
end

macro syms(args...)
    if length(args) == 0
        Node(BlockEnv([],[]))
    elseif is_expr(args[end], [:hcat, :vcat, :vect])
        Node(BlockEnv(args[1:end-1], args[end].args))
    else
        Node(BlockEnv(args, []))
    end
end
macro noenv()
    Node(NoEnv())
end

rebuild(node::@compat(Union{Node{BlockEnv},Node{NoEnv}})) = node
function rebuild(node::Node)
    if isblocknode(node)
        env = envof(node)
        for arg in argsof(node)
            if isa(arg, Node{BlockEnv})
                benv = valueof(arg)
                if !(env.defined == benv.defined)
                    error("env.defined = $(env.defined) != $(benv.defined)")
                end
                just_assigned = setdiff(env.assigned, env.defined)
                if !(just_assigned == benv.assigned)
                    error("just_assigned = $(just_assigned) != $(benv.assigned)")
                end
            elseif isa(arg, Node{NoEnv})
                @assert isa(env, NoEnv)
            end
        end
        Expr(:block, [rebuild(arg) for arg in argsof(node)]...)
    elseif isa(node, ExNode)
        Expr(headof(node), [rebuild(arg) for arg in argsof(node)]...)
    else
        exof(node)
    end
end

function test_decorate(code)
    ecode = macroexpand(code)
    dcode = analyze(NoEnv(), ecode, false)
    rcode = rebuild(dcode)
    @assert rcode == ecode
end
macro test_decorate(ex)
    quote
        test_decorate($(quot(ex)))
    end
end

@test_decorate let
    @syms [f]
    function f(x::Int)
        @syms x [y]
        y = 0
        while x > 0
            @syms [x y]
            x -= y
            y += 1
        end
        y
    end
end

@test_decorate let
    @syms [f]
    function f(z, x::Int=1; a=3, b=4, c=5)
        @syms z x a b c [y]
        y = 0
    end
end

@assert_fails @test_decorate let
    @syms
    local x = 3
end
@assert_fails @test_decorate let
    @syms
    x = 3
end
@assert_fails @test_decorate let
    @syms x
end
@assert_fails @test_decorate let
    @syms [x]
end

@assert_fails @test_decorate let
    @noenv
end
@assert_fails @test_decorate begin
    @syms
end

# ---- scoping tests ----------------------------------------------------------

# symbol defining/assigning constructs
@test_decorate let
    @syms
    let
        @syms d1 d2 d3 d4 d5 [a1 a2 a3 a4 a5]
        # define
        global d1, d2=3, d3::Int, d4::AbstractString = "foo"
        local d5::Float64 = 3
        # assign
        a1 = 5
        a2, a3::Integer = 6, 7
        a4::Int = 23
        a5 += 2
        # neither
        y
        v[i] = x
    end
end

# while
@test_decorate let
    @syms [i]
    i=1
    while (@syms [i]; i < 3) # condition evaluated in outside scope
        @syms j [i z]
        i=i+1
        local j=i^2
        z = i-j
    end
end

# try
@test_decorate let
    try
        @syms x
        local x
    catch e
        @syms e [y]
        y = 2
    end
    try
        @syms x
        local x
    catch e
        @syms e [y]
        y = 2
    finally
        @syms z
        local z
    end
end

# for
@test_decorate let
    @syms [a]
    for x=(a=11; 1:n)
        @syms x [x2]
        x2 = x^2
        push!(z, x2)
    end
end

# cartesian for
@test_decorate let
    @syms
    for i=1:n, j=(@syms i; 1:i)
        @syms j
    end
end

# let
@test_decorate let
    @syms [a]
    let x, y=3, z::Int, u::Int=11, v=(a=11; 23)
        @syms x y z u v
    end
end

# comprehensions
@test_decorate let
    let
        @syms [a]
        [(@syms x y; x*y+z) for x=(@syms [a]; 1:5), y=(a=5; 1:3)]
    end
    let
        @syms [a]
        Any[(@syms x y; x*y+z) for x=(@syms [a]; 1:5), y=(a=5; 1:3)]
    end
    let
        @syms [a b]
        (b=5;Int)[(@syms x; x+z) for x=(@syms [a b]; a=5; 1:5)]
    end
end

# dict comprehensions
@test_decorate let
    let
        @syms [a]
        [(@syms x y; x*y=>z) for x=(@syms [a]; 1:5), y=(a=5; 1:3)]
    end
    let
        @syms [a]
        Any[(@syms x y; x*y=>z) for x=(@syms [a]; 1:5), y=(a=5; 1:3)]
    end
    let
        @syms [a b]
        (b=5;Int=>Int)[(@syms x; x=>z) for x=(@syms [a b]; a=5; 1:5)]
    end
end

# functions
@test_decorate let
    @syms [f1 f2 f3 f4]
    function f1(x, y::(w=3; Int), args...)
        @syms x y args [z w]
        z = x*y
    end
    f2(x, y::(w=3; Int), args::Int...) = begin
        @syms x y args [z w]
        z = x*y
    end
    f3 = function(x, y::(w=3; Int), args...)
        @syms x y args [z w]
        z = x*y
    end
    f4 = (x, y::(w=3; Int), args...)->begin
        @syms x y args [z w]
        z = x*y
    end
end

# functions with type parameters
@test_decorate let
    @syms [f g]
    function f{S,T<:Int}(x::S, y::T)
        @syms S T x y
        (x, y)
    end
    g{S,T<:Int}(x::S, y::T) = begin
        @syms S T x y
        (x, y)
    end
end

# test that the right hand sides below are really evaluated inside the scope
#@assert_fails eval(:(typealias P{Real} Real)) # currently doesn't fail in Julia 0.4
@assert_fails eval(:(abstract  Q{Real} <: Real))
@assert_fails eval(:(type      R{Real} <: Real; end))

# typealias, abstract
@test_decorate let
    @syms T1 T2 T3
    abstract T1
    abstract T2 <: Q
    typealias T3 T1
end

# type
@test_decorate let
    @syms T [a]
    type T<:(a=3; Integer)
        @syms c T
        x::Int
        y::Int
        c(x,y) = begin
            @syms x y
            new(x,y)
        end
        function T(x)
            @syms x
            c(x,2x)
        end
    end
end

# types with type parameters
@test_decorate let
    @syms P Q
    type P{T}
        @syms T
        x::T
    end
    type Q{S,T} <: Associative{S,T}
        @syms S T
        x::T
    end
end

# abstract/typealias with type parameters
@test_decorate let
    @syms A X
    abstract A{S,T} <: B{(@syms A X; Int)}
    typealias X{Q<:Associative{Int,Int},R<:Real} Dict{(@syms Q R; Int),Int}
end

# global/local scope
@test_decorate begin
    @noenv
    begin
        @noenv
    end
    for x=(@noenv; 1:3)
        @syms x
    end
    let y=(@noenv; 1)
        @syms y
    end
end

end # module