text:
module Testset

using Test: AbstractTestSet, Broken, DefaultTestSet, Error, Fail, Pass, Test,
            TestSetException, get_testset, get_testset_depth,
            parse_testset_args, pop_testset, push_testset

import Test: finish, record

import Random

using Printf: @sprintf

using Distributed: myid, nworkers

import InlineTest: @testset

using ..ReTest: Pattern, matches

# mostly copied from Test stdlib
# changed from Test: pass nothing as file in ip_has_file_and_func
#-----------------------------------------------------------------------

# (this has to be copied from Test, because `@__FILE__` is hardcoded)

# Backtrace utility functions
function ip_has_file_and_func(ip, file, funcs)
    return any(fr -> ((file === nothing || string(fr.file) == file) &&
                      fr.func in funcs),
               StackTraces.lookup(ip))
end

function scrub_backtrace(bt)
    do_test_ind = findfirst(ip -> ip_has_file_and_func(ip, nothing,
                                                       (:do_test, :do_test_throws)), bt)
    if do_test_ind !== nothing && length(bt) > do_test_ind
        bt = bt[do_test_ind + 1:end]
    end
    name_ind = findfirst(ip -> ip_has_file_and_func(ip, @__FILE__, (Symbol("macro expansion"),)), bt)
    if name_ind !== nothing && length(bt) != 0
        bt = bt[1:name_ind]
    end
    return bt
end

function scrub_exc_stack(stack)
    return Any[ (x[1], scrub_backtrace(x[2])) for x in stack ]
end

mutable struct Format
    stats::Bool
    desc_align::Int
    pass_width::Int
    fail_width::Int
    error_width::Int
    broken_width::Int
    total_width::Int
end

Format(stats, desc_align) = Format(stats, desc_align, 0, 0, 0,0 ,0)

mutable struct ReTestSet <: AbstractTestSet
    mod::String # enclosing module
    description::AbstractString
    id::Int64
    overall::Bool # TODO: could be conveyed by having self.mod == ""
    results::Vector
    n_passed::Int
    anynonpass::Bool
    verbose::Bool
    timed::NamedTuple
    exception::Union{TestSetException,Nothing}
end

ReTestSet(mod, desc, id::Integer=0; overall=false, verbose = true) =
    ReTestSet(mod, desc, id, overall, [], 0, false, verbose, NamedTuple(), nothing)

# For a non-passed result, simply store the result
record(ts::ReTestSet, t::Union{Broken,Fail,Error}) = (push!(ts.results, t); t)
# For a passed result, do not store the result since it uses a lot of memory
record(ts::ReTestSet, t::Pass) = (ts.n_passed += 1; t)

# When a ReTestSet finishes, it records itself to its parent
# testset, if there is one. This allows for recursive printing of
# the results at the end of the tests
record(ts::ReTestSet, t::AbstractTestSet) = push!(ts.results, t)

function print_test_errors(ts::ReTestSet)
    for t in ts.results
        if isa(t, Error) || isa(t, Fail)
            printstyled(ts.description, ": ", color=:white)

            # don't print for interrupted tests
            if t isa Fail || t.test_type !== :test_interrupted
                show(t)
            end
            if t isa Fail # if not gets printed in the show method
                # Base.show_backtrace(stdout, scrub_backtrace(backtrace()))
            end
            println()
        elseif isa(t, ReTestSet)
            print_test_errors(t)
        end
    end
end

function print_test_results(ts::ReTestSet, fmt::Format;
                            depth::Int=0, bold::Bool=false, hasbroken::Bool=false,
                            maxidw::Int)
    # Calculate the overall number for each type so each of
    # the test result types are aligned
    upd = false

    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    total_pass   = passes + c_passes
    total_fail   = fails  + c_fails
    total_error  = errors + c_errors
    total_broken = broken + c_broken
    total = total_pass + total_fail + total_error + total_broken

    # if no tests, count that as a pass:
    dig_pass   = total_pass   > 0 || total == 0 ? ndigits(total_pass)   : 0
    dig_fail   = total_fail   > 0               ? ndigits(total_fail)   : 0
    dig_error  = total_error  > 0               ? ndigits(total_error)  : 0
    dig_broken = total_broken > 0 || hasbroken  ? ndigits(total_broken) : 0

    # max(1, ...) : we always print something at least in "Pass", even when no tests
    nprinted = max(1, (total_pass > 0) + (total_fail > 0) +
                      (total_error > 0) + (total_broken > 0))
    if nprinted == 1 && total_pass > 0 && !hasbroken
        # do not print "Total" when only "Pass" column is printed
        total = 0
    end
    dig_total = total > 0 || hasbroken ? ndigits(total) : 0

    # For each category, take max of digits and header width if there are
    # tests of that type
    fail_width   = dig_fail   > 0 ? max(6,   dig_fail) : 0
    error_width  = dig_error  > 0 ? max(6,  dig_error) : 0
    broken_width = dig_broken > 0 ? max(6, dig_broken) : 0
    total_width  = dig_total  > 0 ? max(6,  dig_total) : 0
    # heuristic: if we print total because of a broken test, still print "Pass",
    # as it's likely that later on a passing test will be printed
    pass_width   = dig_pass   > 0 || total_width > 0 ? max(6,   dig_pass) : 0

    if pass_width > fmt.pass_width
        upd = true
        fmt.pass_width = pass_width
    else
        pass_width = fmt.pass_width
    end
    if fail_width > fmt.fail_width
        upd = true
        fmt.fail_width = fail_width
    else
        fail_width = fmt.fail_width
    end
    if error_width > fmt.error_width
        upd = true
        fmt.error_width = error_width
    else
        error_width = fmt.error_width
    end
    if broken_width > fmt.broken_width
        upd = true
        fmt.broken_width = broken_width
    else
        broken_width = fmt.broken_width
    end
    if total_width > fmt.total_width
        upd = true
        fmt.total_width = total_width
    else
        total_width = fmt.total_width
    end

    # Calculate the alignment of the test result counts by
    # recursively walking the tree of test sets
   if !ts.overall # otherwise, we don't print recursively
       align = get_alignment(ts, 0)
   end

    if !ts.overall && align > fmt.desc_align
        upd = true
        fmt.desc_align = align
    else
        align = fmt.desc_align
    end

    # Print the outer test set header once
    if upd
        if maxidw > 0
            print(' '^(maxidw + 2)) # +2 for "| " after number
        end
        printstyled(rpad("", align, " "), "  ", " "; bold=true, color=:white)
        if pass_width > 0
            printstyled(lpad("Pass", pass_width, " "), "  "; bold=true, color=:green)
        end
        if fail_width > 0
            printstyled(lpad("Fail", fail_width, " "), "  "; bold=true, color=Base.error_color())
        end
        if error_width > 0
            printstyled(lpad("Error", error_width, " "), "  "; bold=true, color=Base.error_color())
        end
        if broken_width > 0
            printstyled(lpad("Broken", broken_width, " "), "  "; bold=true, color=Base.warn_color())
        end
        if total_width > 0
            printstyled(lpad("Total", total_width, " "), fmt.stats ? "  " : ""; bold=true, color=Base.info_color())
        end
        if fmt.stats
            # copied from Julia/test/runtests.jl
            compile_header = VERSION >= v"1.6-" ? " Compile /" : ""
            printstyled("|  Time /$compile_header GC |   Alloc   ΔRSS |", color=:white)
        end
        println()
    end

    # Recursively print a summary at every level
    print_counts(ts, fmt, depth, align,
                 pass_width, fail_width, error_width, broken_width, total_width;
                 bold=bold, maxidw=maxidw)
end

# Called at the end of a @testset, behaviour depends on whether
# this is a child of another testset, or the "root" testset
function finish(ts::ReTestSet, chan)
    # If we are a nested test set, do not print a full summary
    # now - let the parent test set do the printing
    if get_testset_depth() != 0
        # Attach this test set to the parent test set
        parent_ts = get_testset()
        record(parent_ts, ts)
        return ts
    end
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    total_pass   = passes + c_passes
    total_fail   = fails  + c_fails
    total_error  = errors + c_errors
    total_broken = broken + c_broken
    total = total_pass + total_fail + total_error + total_broken

    # Finally throw an error as we are the outermost test set
    if total != total_pass + total_broken
        # Get all the error/failures and bring them along for the ride
        efs = filter_errors(ts)
        ts.exception = TestSetException(total_pass, total_fail, total_error,
                                        total_broken, efs)
    end

    put!(chan.out, ts)
    if myid() == 1
        take!(chan.compute)
    end

    # return the testset so it is returned from the @testset macro
    ts
end

# Recursive function that finds the column that the result counts
# can begin at by taking into account the width of the descriptions
# and the amount of indentation. If a test set had no failures, and
# no failures in child test sets, there is no need to include those
# in calculating the alignment
function get_alignment(ts::ReTestSet, depth::Int)
    # The minimum width at this depth is
    ts_width = 2*depth + length(ts.description)
    # If all passing, no need to look at children
    !ts.anynonpass && return ts_width
    # Return the maximum of this width and the minimum width
    # for all children (if they exist)
    isempty(ts.results) && return ts_width
    child_widths = map(t->get_alignment(t, depth+1), ts.results)
    return max(ts_width, maximum(child_widths))
end
get_alignment(ts, depth::Int) = 0

# Recursive function that fetches backtraces for any and all errors
# or failures the testset and its children encountered
function filter_errors(ts::ReTestSet)
    efs = []
    for t in ts.results
        if isa(t, ReTestSet)
            append!(efs, filter_errors(t))
        elseif isa(t, Union{Fail, Error})
            append!(efs, [t])
        end
    end
    efs
end

# Recursive function that counts the number of test results of each
# type directly in the testset, and totals across the child testsets
function get_test_counts(ts::ReTestSet)
    passes, fails, errors, broken = ts.n_passed, 0, 0, 0
    c_passes, c_fails, c_errors, c_broken = 0, 0, 0, 0
    for t in ts.results
        isa(t, Fail)   && (fails  += 1)
        isa(t, Error)  && (errors += 1)
        isa(t, Broken) && (broken += 1)
        if isa(t, ReTestSet)
            np, nf, ne, nb, ncp, ncf, nce , ncb = get_test_counts(t)
            c_passes += np + ncp
            c_fails  += nf + ncf
            c_errors += ne + nce
            c_broken += nb + ncb
        end
    end
    ts.anynonpass = (fails + errors + c_fails + c_errors > 0)
    return passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken
end

# Recursive function that prints out the results at each level of
# the tree of test sets
function print_counts(ts::ReTestSet, fmt::Format, depth, align,
                      pass_width, fail_width, error_width, broken_width, total_width;
                      bold, maxidw)
    # Count results by each type at this level, and recursively
    # through any child test sets
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    subtotal = passes + fails + errors + broken + c_passes + c_fails + c_errors + c_broken
    # Print test set header, with an alignment that ensures all
    # the test results appear above each other

    style = bold ? (bold=bold, color=:white) : NamedTuple()
    if maxidw > 0
        if ts.id != 0
            printstyled(lpad(ts.id, maxidw), "| ", color = :light_black, bold=true)
        else
            print(' '^(maxidw+2))
        end
    end

    printstyled(rpad(string("  "^depth, ts.description), align, " "); style...)

    np = passes + c_passes
    nf = fails + c_fails
    ne = errors + c_errors
    nb = broken + c_broken

    print_total = true

    if ts.overall && endswith(ts.description, ':')
        # header, we don't print the '|' nor anything afterwards
        @assert np == 0 && nf == 0 && ne == 0 && nb == 0
        print_total = false
    elseif np > 0 || np == 0 && nf == 0 && ne == 0 && nb == 0
        # print `0` in warn color instead of "No tests" like in Test module,
        # which messes up alignments (and am too lazy to fix)
        printstyled(" | ", bold=false)
        printstyled(lpad(string(np), pass_width, " "), "  ",
                    color = np > 0 ? :green : Base.warn_color(), bold=bold)
    elseif pass_width > 0 # TODO: isn't this condition always true here?
        # No passes at this level, but some at another level
        printstyled(" | ", bold=false)
        print(lpad(" ", pass_width), "  ")
    end

    if nf > 0
        printstyled(lpad(string(nf), fail_width, " "), "  ", color=Base.error_color(), bold=bold)
    elseif fail_width > 0
        # No fails at this level, but some at another level
        print(lpad(" ", fail_width), "  ")
    end

    if ne > 0
        printstyled(lpad(string(ne), error_width, " "), "  ", color=Base.error_color(), bold=bold)
    elseif error_width > 0
        # No errors at this level, but some at another level
        print(lpad(" ", error_width), "  ")
    end

    if nb > 0
        printstyled(lpad(string(nb), broken_width, " "), "  ", color=Base.warn_color(), bold=bold)
    elseif broken_width > 0
        # None broken at this level, but some at another level
        print(lpad(" ", broken_width), "  ")
    end

    if total_width > 0 && print_total
        printstyled(lpad(string(subtotal), total_width, " "), fmt.stats ? "  " : "", color=Base.info_color(), bold=bold)
    end

    if fmt.stats && print_total # copied from Julia/test/runtests.jl
        ts.overall && set_timed!(ts)
        timed = ts.timed

        # we don't want to report zeros, which makes it hard to spot non-zeros
        function hide_zero(str, unit)
            strip(str, ' ') in ("0.0", "0.00") ?
                ' '^(1+length(str)) : # +1 for unit
                str * unit
        end

        time_str = hide_zero(@sprintf("%6.2f", timed.time), "s")
        printstyled("| ", time_str, " ", color=:white)
        if VERSION >= v"1.6-"
            compile_str = all(==(' '), time_str) ?
                ' '^6 : # print percentages only if time itself is shown!
                        # (also, this can result in weird things, like "30663.3%",
                        # if e.g. there was no @test in the @testset)
                hide_zero(@sprintf("%5.1f", timed.compile_time / 10^7 / timed.time), "%")
            # can be >= 100% !?
            printstyled(compile_str, " ", color=:white)
        end
        gc_str = all(==(' '), time_str) ?
            ' '^5 :
            hide_zero(@sprintf("%4.1f", 100 * timed.gctime / timed.time), "%")
        printstyled(gc_str, " | ", color=:white)

        alloc_str = hide_zero(@sprintf("%6.1f", timed.bytes / 2^20), "M")
        printstyled(alloc_str, " ", color=:white)

        rss_str = hide_zero(@sprintf("%5.1f", timed.rss / 2^20), "M")
        printstyled(rss_str, " |", color=:white)
    end
    println()

    # Only print results at lower levels if we had failures or if the user
    # wants.
    if !ts.overall && ((np + nb != subtotal) || (ts.verbose))
        for t in ts.results
            if isa(t, ReTestSet)
                print_counts(t, fmt, depth + 1, align,
                             pass_width, fail_width, error_width, broken_width, total_width,
                             bold=false, maxidw=maxidw)
            end
        end
    end
end

#-----------------------------------------------------------------------

default_rng() = isdefined(Random, :default_rng) ?
    Random.default_rng() :
    Random.GLOBAL_RNG

function get_testset_string(remove_last=false)
    testsets = get(task_local_storage(), :__BASETESTNEXT__, Test.AbstractTestSet[])
    join('/' * ts.description for ts in (remove_last ? testsets[1:end-1] : testsets))
end

# HACK: we re-use the same macro name `@testset` for actual execution (like in `Test`)
# as for the one documented in ReTest (deferred execution), because otherwise
# packages wouldn't know about the former, say named @retestset (which has to be evaluated
# within packages' modules). We could export @retestset but this is fragile;
# is there a simple alternative?
# In the meantime, we create this semi-catch-all method, called from make_ts,
# so that when arguments are messed up (forgot to update call in make_ts), we don't end-up
# calling the real-catch-all method from InlineTest (which is again deferred), which
# regularly leads to confusing behavior

# try to not mess up first three arguments!
macro testset(mod::String, isfinal::Bool, pat::Pattern, x...)
    error("invalid arguments")
end

# non-inline testset with regex filtering support
macro testset(mod::String, isfinal::Bool, pat::Pattern, id::Int64, desc::Union{String,Expr}, options,
              stats::Bool, chan, body)
    Testset.testset_beginend(mod, isfinal, pat, id, desc, options, stats, chan,  body, __source__)
end

macro testset(mod::String, isfinal::Bool, pat::Pattern, id::Int64, desc::Union{String,Expr}, options,
              stats::Bool, chan, loops, body)
    Testset.testset_forloop(mod, isfinal, pat, id, desc, options,
                            stats, chan, loops, body, __source__)
end

"""
Generate the code for a `@testset` with a `begin`/`end` argument
"""
function testset_beginend(mod::String, isfinal::Bool, pat::Pattern, id::Int64, desc, options,
                          stats::Bool, chan, tests, source)
    # Generate a block of code that initializes a new testset, adds
    # it to the task local storage, evaluates the test(s), before
    # finally removing the testset and giving it a chance to take
    # action (such as reporting the results)
    desc = esc(desc)
    ex = quote
        local current_str
        if $isfinal
            current_str = string(get_testset_string(), '/', $desc)
        end
        if !$isfinal || matches($pat, current_str, $id)
            local ret
            local ts = ReTestSet($mod, $desc, $id; verbose=$(options.transient_verbose))
            if nworkers() == 1 && get_testset_depth() == 0 && $(chan.preview) !== nothing
                put!($(chan.preview), $desc)
            end
            push_testset(ts)
            # we reproduce the logic of guardseed, but this function
            # cannot be used as it changes slightly the semantic of @testset,
            # by wrapping the body in a function
            local RNG = default_rng()
            local oldrng = copy(RNG)
            try
                # RNG is re-seeded with its own seed to ease reproduce a failed test
                Random.seed!(RNG.seed)
                let
                    ts.timed = @stats $stats $(esc(tests))
                end
            catch err
                err isa InterruptException && rethrow()
                # something in the test block threw an error. Count that as an
                # error in this test set
                record(ts, Error(:nontest_error, Expr(:tuple), err,
                                 Base.catch_stack(), $(QuoteNode(source))))
            finally
                copy!(RNG, oldrng)
                pop_testset()
                ret = finish(ts, $chan)
            end
            ret
        end
    end
    # preserve outer location if possible
    if tests isa Expr && tests.head === :block &&
        !isempty(tests.args) && tests.args[1] isa LineNumberNode

        ex = Expr(:block, tests.args[1], ex)
    end
    return ex
end

"""
Generate the code for a `@testset` with a `for` loop argument
"""
function testset_forloop(mod::String, isfinal::Bool, pat::Pattern, id::Int64,
                         desc::Union{String,Expr}, options, stats, chan, loops, tests, source)

    desc = esc(desc)
    blk = quote
        local current_str
        if $isfinal
            current_str = string(get_testset_string(!first_iteration), '/', $desc)
        end
        if !$isfinal || matches($pat, current_str, $id)
            # Trick to handle `break` and `continue` in the test code before
            # they can be handled properly by `finally` lowering.
            if !first_iteration
                pop_testset()
                push!(arr, finish(ts, $chan))
                # it's 1000 times faster to copy from tmprng rather than calling Random.seed!
                copy!(RNG, tmprng)
            end
            ts = ReTestSet($mod, $desc, $id; verbose=$(options.transient_verbose))
            if nworkers() == 1 && get_testset_depth() == 0 && $(chan.preview) !== nothing
                put!($(chan.preview), ts.description)
            end
            push_testset(ts)
            first_iteration = false
            try
                let
                    ts.timed = @stats $stats $(esc(tests))
                end
            catch err
                err isa InterruptException && rethrow()
                # Something in the test block threw an error. Count that as an
                # error in this test set
                record(ts, Error(:nontest_error, Expr(:tuple), err, Base.catch_stack(), $(QuoteNode(source))))
            end
        end
    end
    quote
        local arr = Vector{Any}()
        local first_iteration = true
        local ts
        local RNG = default_rng()
        local oldrng = copy(RNG)
        Random.seed!(RNG.seed)
        local tmprng = copy(RNG)
        try
            let
                $(Expr(:for, Expr(:block, [esc(v) for v in loops]...), blk))
            end
        finally
            # Handle `return` in test body
            if !first_iteration
                pop_testset()
                push!(arr, finish(ts, $chan))
            end
            copy!(RNG, oldrng)
        end
        arr
    end
end

function set_timed!(ts)
    if ts.overall
        foreach(get_timed!, ts.results)
    end
    ts.timed = NamedTuple{(:time, :bytes, :gctime, :rss, :compile_time)}(
        ntuple(Val(5)) do i
            # init=0.0 kwarg not available on old Julia
            Float64(sum(tsc.timed[i] for tsc in ts.results))
        end)
end

get_timed!(ts) = isempty(ts.timed) ? set_timed!(ts) : ts

# adapted from @timed in Julia/base/timing.jl
# also, @timed inserts a `while false; end` compiler heuristic, which destroys perfs here
macro stats(yes, ex)
    quote
        if $yes
            local stats = Base.gc_num()
            local elapsedtime = time_ns()
            local rss = Sys.maxrss()
            local compile_time = cumulative_compile_time_ns()
        end
        local val = $(esc(ex))
        if $yes
            elapsedtime = time_ns() - elapsedtime
            local diff = Base.GC_Diff(Base.gc_num(), stats)
            rss = Sys.maxrss() - rss
            compile_time = cumulative_compile_time_ns() - compile_time
            # COMPAT: on Julia 1.1, the form `(time=..., bytes=..., ...)` doesn't work
            # (macro name mangling with `#`, e.g. (#115#time = ..., ))
            NamedTuple{(:time, :bytes, :gctime, :rss, :compile_time)}(
                (elapsedtime/1e9, diff.allocd, diff.total_time/1e9, rss, compile_time))
        else
            NamedTuple()
        end
    end
end

cumulative_compile_time_ns() =
    isdefined(Base, :cumulative_compile_time_ns) ?
        Base.cumulative_compile_time_ns() :
        UInt(0)

end # module


text:
struct ElementaryPathInstance{T, O <: CombinatorialObjective} <:
       CombinatorialInstance
    graph::AbstractGraph{T}
    rewards::Dict{Edge{T}, Float64}
    src::T
    dst::T
    objective::O

    function ElementaryPathInstance(
        graph::AbstractGraph{T},
        rewards::Dict{Edge{T}, Float64},
        src::T,
        dst::T,
        objective::O=Maximise(),
    ) where {T <: Real, O <: CombinatorialObjective}
        # Source and destination vertices are allowable.
        if src < 0 || src > nv(graph)
            error("The requested source is not a vertex of the given graph.")
        end

        if dst < 0 || dst > nv(graph)
            error("The requested destination is not a vertex of the given " *
                  "graph.")
        end

        # Source is different from destination.
        if dst == src
            error("The requested source and destination are the same vertex.")
        end

        # Each reward is associated to an existing edge.
        for e in keys(rewards)
            if e ∉ edges(graph)
                error("The edge $(e) has a reward, but is not present in the " *
                      "graph.")
            end
        end

        # Return a new instance.
        return new{T, O}(graph, rewards, src, dst, objective)
    end
end

dimension(i::ElementaryPathInstance{T, O}) where {T, O} = ne(i.graph)

function copy(
    i::ElementaryPathInstance{T, O};
    graph::AbstractGraph{T}=i.graph,
    rewards::Dict{Edge{T}, Float64}=i.rewards,
    src::T=i.src,
    dst::T=i.dst,
    objective::CombinatorialObjective=i.objective,
) where {T, O <: CombinatorialObjective}
    return ElementaryPathInstance(graph, rewards, src, dst, objective)
end

# Solution.

struct ElementaryPathSolution{T, O <: CombinatorialObjective} <:
       CombinatorialSolution
    instance::ElementaryPathInstance{T, O}
    variables::Vector{Edge{T}}
    solutions::Dict{T, Vector{Edge{T}}}
end

function ElementaryPathSolution(
    instance::ElementaryPathInstance{T, O},
    variables::Vector{Edge{T}},
) where {T, O <: CombinatorialObjective}
    return ElementaryPathSolution(
        instance,
        variables,
        Dict{T, Vector{Edge{T}}}(),
    )
end

function make_solution(
    i::ElementaryPathInstance{T, O},
    path::Dict{Edge{T}, Float64},
) where {T, O}
    path_edges = Edge{T}[]
    for (k, v) in path
        if v >= 0.5
            push!(path_edges, k)
        end
    end

    return ElementaryPathSolution(i, path_edges)
end

# Budgeted solution.

struct BudgetedElementaryPathSolution{T, O} <: MultipleMinBudgetedSolution
    instance::MinimumBudget{ElementaryPathInstance{T, O}, T}
    variables::Vector{Edge{T}}
    solutions::Dict{Int, Vector{Edge{T}}}
end

function BudgetedElementaryPathSolution(
    instance::MinimumBudget{ElementaryPathInstance{T, O}, T},
    variables::Vector{Edge{T}},
) where {T, O <: CombinatorialObjective}
    return BudgetedElementaryPathSolution(
        instance,
        variables,
        Dict{Int, Vector{Edge{T}}}(),
    )
end

function make_solution(
    i::MinimumBudget{ElementaryPathInstance{T, O}, T},
    path::Dict{Edge{T}, Float64},
) where {T, O}
    path_edges = Edge{T}[]
    for (k, v) in path
        if v >= 0.5
            push!(path_edges, k)
        end
    end

    return BudgetedElementaryPathSolution(i, path_edges)
end

function _check_and_warn_budget_too_high(
    s::BudgetedElementaryPathSolution{T, O},
    max_budget::Int,
) where {T, O}
    if max_budget > s.instance.min_budget
        @warn "The requested maximum budget $max_budget is higher than the instance's minimum budget $(s.instance.min_budget). Therefore, some values have not been computed and are unavailable."
    end
end

function paths_all_budgets(
    s::BudgetedElementaryPathSolution{T, O},
    max_budget::Int,
) where {T, O}
    _check_and_warn_budget_too_high(s, max_budget)
    mb = min(max_budget, s.instance.min_budget)
    return Dict{Int, Vector{Edge{T}}}(
        budget => s.solutions[budget] for budget in 0:mb
    )
end

function paths_all_budgets_as_tuples(
    s::BudgetedElementaryPathSolution{T, O},
    max_budget::Int,
) where {T, O}
    _check_and_warn_budget_too_high(s, max_budget)
    mb = min(max_budget, s.instance.min_budget)
    return Dict{Int, Vector{Tuple{T, T}}}(
        budget => [(src(e), dst(e)) for e in s.solutions[budget]] for
        budget in 0:mb
    )
end

ext: jl


text:
# DEFINE FUNCTIONS THAT ARE PART OF THE MODEL

#*====== SMOOTH SATURATION FUNCTION =======

# true nonsmooth saturation function
function saturation(x::Float64,x_min::Float64,x_max::Float64)

	x_out::Float64 = median([x, x_min, x_max])
	return x_out
end

# smooth saturation function
function smoothSaturation(x::Float64,x_min::Float64,x_max::Float64,tau_sat::Float64)

# scale smoothing parameter relative to x_max and x_min
	tau::Float64 = tau_sat*(x_max - x_min)

	x_out::Float64 = x_max

	if(x < x_min - tau)
		x_out = x_min
	elseif(x <= x_min + tau)
		x_out = x_min + (1/(4*tau))*(x - (x_min - tau))^2
	elseif(x < x_max - tau)
		x_out = x
	elseif(x <= x_max + tau)
		x_out = x_max - (1/(4*tau))*(x - (x_max + tau))^2
	end

	return x_out
end

# return derivative of smooth saturation function
function smoothSaturationDeriv(x::Float64,x_min::Float64,x_max::Float64,tau_sat::Float64)

# scale smoothing parameter relative to x_max and x_min
	tau::Float64 = tau_sat*(x_max - x_min)

	deriv::Float64 = 0

	if(x < x_min - tau)
		# do nothing
	elseif(x <= x_min + tau)
		deriv = (1/(2*tau))*(x - (x_min - tau))
	elseif(x < x_max - tau)
		deriv = 1
	elseif(x <= x_max + tau)
		deriv = - (1/(2*tau))*(x - (x_max + tau))
	end

	return deriv
end

#*====== LINE PENALTY FUNCTION =======

# line penalty function
function linePen(x::Float64,x_max::Float64,gamma_line::Float64,delta_line::Float64)

	line_viol::Float64 = maximum([abs(x) - delta_line*x_max, 0])
	pen_val::Float64 = gamma_line*(line_viol)^2
	return pen_val
end

# derivative of line penalty function
function linePenDeriv(x::Float64,x_max::Float64,gamma_line::Float64,delta_line::Float64)

	line_viol::Float64 = maximum([abs(x) - delta_line*x_max, 0])
	deriv::Float64 = 2*gamma_line*line_viol*sign(x)
	return deriv
end

#*====== RESERVE PENALTY FUNCTION =======

# returns smooth version of max{x, 0}
function smoothPos(x::Float64,tau_pos::Float64)

	if(x <= 0)
		return (tau_pos*log(1+exp(x/tau_pos)))
	else
	# rewrite for numerical stability
		return (x + tau_pos*log(1 + exp(-x/tau_pos)))
	end

end

# derivative of smooth positive function
function smoothPosDeriv(x::Float64,tau_pos::Float64)

	deriv::Float64 = 1/(1+exp(-x/tau_pos))
	return deriv
end

# reserves penalty function
function reserveSmoothPen(x::Float64,costCoeff::Float64,gamma_res::Float64,tau_pos::Float64)

	res_pen::Float64 = gamma_res*costCoeff*smoothPos(x,tau_pos)
	return res_pen
end

# derivative of reserves penalty function
function reserveSmoothPenDeriv(x::Float64,costCoeff::Float64,gamma_res::Float64,tau_pos::Float64)

	deriv::Float64 = gamma_res*costCoeff*smoothPosDeriv(x,tau_pos)
	return deriv
end

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
#=
LVE.jl
    Written by Matthew White
    5/21/2019
=#
function LVE(surf::TwoDSurf,curfield::TwoDFlowField,nsteps::Int64 = 500,dtstar::Float64 = 0.015,frameCnt = 20,view = "square",LEV_Stop = 1000,anim = false,animStep = 30,wflag = false,writeInterval = 50)
    ## Initialize variables
    mat = zeros(8 ,0) # loads output matrix
    prevCirc = zeros(surf.ndiv-1)
    circChange = zeros(surf.ndiv-1)
    vor_loc = zeros(surf.ndiv-1,2) # (x1,z1 ; x2,z2 ...) Inertial frame
    coll_loc = zeros(surf.ndiv-1,2) # (x1,z1 ; x2,z2 ...)
    global RHS = zeros(surf.ndiv)
    global relVel = zeros(surf.ndiv-1,2)
    IFR_vor = 0 .* vor_loc # Global frame
    IFR_field = TwoDFlowField()
    IFR_surf = TwoDSurf(surf.coord_file,surf.pvt,surf.kindef,surf.lespcrit; ndiv = surf.ndiv, camberType = "linear")
    global dp = zeros(surf.ndiv-1,1) # change in panel pressure
    frameStep = round((nsteps-1) / frameCnt) # Frames to capture animation on
    frameStep < 1 ? frameStep = 1 : frameStep
    frames = 0
    prevNow = 0
    LEVhist = hcat([],[])
    TEVhist = hcat([],[])
    animMat = hcat([],[],[],[],[],[],[],[])

    #   length of each panel
    ds = sqrt.(( surf.x[1:end-1]-surf.x[2:end]).^2 + (surf.cam[1:end-1]-surf.cam[2:end]).^2)
    #   Surf cam_slope does not give the correct panel locations
    cam_slope = asind.((surf.cam[2:end]-surf.cam[1:end-1])./ds) # [deg]

    # Normal vectors for each panel
    n = hcat(-sind.(cam_slope),cosd.(cam_slope))
    # Tangential vectors for each panel
    tau = hcat(cosd.(cam_slope),sind.(cam_slope))

    ## Vortex Locations
    # Located at 25% of each panel
    #   Vortex loactions (at 25% of panel length)
    vor_loc[:,1] = surf.x[1:end-1] + .25 .* ds .* cosd.(cam_slope)
    vor_loc[:,2] = surf.cam[1:end-1] + .25 .* ds .* sind.(cam_slope)
    #   Collocation points (at 75% panel chordwise)
    coll_loc[:,1] = surf.x[1:end-1] + .75 .* ds .* cosd.(cam_slope)
    coll_loc[:,2] = surf.cam[1:end-1] + .75 .* ds .* sind.(cam_slope)

    ## Critical gamma (LEV)
    x = 1 - 2*ds[1] / surf.c
    gamma_crit = surf.lespcrit / 1.13 * ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) )))
    gamma_crit = gamma_crit[1]

    ## Refresh vortex values
    refresh_vortex(surf,vor_loc)

    ## Time Stepping Loop
    t = -dtstar

    for i = 0:nsteps
        t += dtstar

        # Update Kinematics
        update_kinem(surf, t)

        # Inertial reference frame
        IFR_vor[:,1],IFR_vor[:,2] = IFR(surf,vor_loc[:,1],vor_loc[:,2],t)
        IFR_surf = refresh_vortex(IFR_surf,IFR_vor)

        ## TEV setup
        place_tev2(IFR_surf,IFR_field,dtstar,t)

        ## Influence Coeffcients
        x_w, z_w = BFR(surf, IFR_field.tev[end].x, IFR_field.tev[end].z, t)
        global a = influence_coeff(surf,curfield,coll_loc,n,dtstar,x_w,z_w)

        ## RHS Vector
        RHS[end] = sum(prevCirc) # previous total circulation
        # u_w, w_w
        if length(curfield.tev) > 0 # Trailing edge vortexs exist
            surf.uind[1:end-1], surf.wind[1:end-1] = ind_vel([curfield.tev; curfield.lev],coll_loc[:,1],coll_loc[:,2])
        end
        for j = 1:surf.ndiv-1
            alpha = surf.kinem.alpha
            # relVel = [U(t) ; W(t)]
            global relVel[j,:] = [cos(alpha) -sin(alpha) ; sin(alpha) cos(alpha)] * [surf.kinem.u ; -surf.kinem.hdot] + [-surf.kinem.alphadot*coll_loc[j,2] ; surf.kinem.alphadot * (coll_loc[j,1] - surf.pvt)]

            # RHS = -[U_t + u_w , W_t + w_w] dot n
            global RHS[j] = -((relVel[j,1] + surf.uind[j])*n[j,1] + (relVel[j,2] + surf.wind[j])*n[j,2])
        end

        ## Vortex Strenghths a*circ = RHS
        global circ = a\RHS

        ## Circulation changes before each panel
        for j = 1:surf.ndiv-1
            gamma1 = sum(prevCirc[1:j-1])
            gamma2 = sum(circ[1:j-1])
            circChange[j] = (gamma2-gamma1) ./ dtstar
        end

        ## LEV Ejection
        #   Test for LEV criteria
        if abs(circ[1]) > gamma_crit
            if surf.levflag[1] == 0 && t < LEV_Stop # if lev is first in batch
                gamma_crit_use = abs(circ[1])/circ[1] * gamma_crit # = gamma_crit w/ sign of circ[1]
                global t_start = t
                #println("LEV start batch:")
            else # if lev sheds but is not first in batch
                m_slp = -surf.lespcrit / (LEV_Stop - t_start)
                c_slp = surf.lespcrit - (m_slp*t_start)
                LESP_crit_use = (m_slp*t) + c_slp
                x = 1 - 2*ds[1] / surf.c
                gamma_crit_use = abs(circ[1])/circ[1] * LESP_crit_use / 1.13 * ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) ))) # abs(circ[1])/circ[1] keeps sign of original LEV
            end
            #println("LEV Ejection: Step $i, time $t")
            gamma_crit_use = gamma_crit_use[1]

            # Place LEV
            place_lev2(surf,IFR_field,dtstar,t)

            surf.levflag[1] = 1 # Set LEV flag to track constant LEV streams

            # Convert coords to BFR
            x_lev, z_lev = BFR(surf, IFR_field.lev[end].x, IFR_field.lev[end].z, t)
            global a,a1 = mod_influence_coeff(surf,curfield,coll_loc,n,dtstar,x_w,z_w,x_lev,z_lev)

            # Recalculate RHS
            # u_w, w_w
            if length(curfield.tev) > 0 # Trailing edge vortexs exist
                surf.uind[1:end-1], surf.wind[1:end-1] = ind_vel([curfield.tev; curfield.lev],coll_loc[:,1],coll_loc[:,2])
            end
            for j = 1:surf.ndiv-1
                # RHS = -[U_t + u_w , W_t + w_w] dot n
                global RHS[j] = -((relVel[j,1] + surf.uind[j])*n[j,1] + (relVel[j,2] + surf.wind[j])*n[j,2]) - gamma_crit_use*a1[j]
            end
            RHS[end] -= gamma_crit_use # first RHS - gamma_crit_use*a1[end] (which is always 1)

            global circ = a\RHS # Solve new circluations

            IFR_field.lev[end].s = circ[end] #store LEV strength

            circ[2:end] = circ[1:end-1] # shift circs
            circ[1] = gamma_crit_use

            # LEV history data [t,s]
            stamp = hcat(t,IFR_field.lev[end].s)
            LEVhist = vcat(LEVhist,stamp)
        else
            surf.levflag[1] = 0
        end

        prevCirc = circ[1:end-1]

        for j = 1:surf.ndiv-1 # Set bv circs
            surf.bv[j].s = circ[j]
            IFR_surf.bv[j].s = circ[j]
        end
        if length(IFR_field.tev) > 1 # Set TEV circ
            IFR_field.tev[end].s = circ[end]
        end

        # LEV history data [t,s]
        stamp = hcat(t,IFR_field.tev[end].s)
        TEVhist = vcat(TEVhist,stamp)

        #= Position mesh
        if i % frameStep == 0. && Int(i / frameStep) < frameCnt

            temp = Int(i / frameStep) + 1
            now = Dates.format(Dates.now(), "HH:MM")
            if prevNow == now
                println("$temp / $frameCnt")
            else
                println("$temp / $frameCnt,  $now")
            end
            prevNow = now
            # Velocities at mesh points
            global mesh = meshgrid(surf,IFR_field.tev,IFR_field.lev,.25,t,100,view)
            vorts = [IFR_field.tev ; IFR_field.lev ; IFR_surf.bv]
            for j = 1:size(vorts,1)
                u,w = ind_vel(vorts[j],mesh.x,mesh.z)
                global mesh.uMat = mesh.uMat .+ u
                global mesh.wMat = mesh.wMat .+ w
                global mesh.velMag = sqrt.(mesh.uMat.^2 + mesh.wMat.^2)
                global mesh.t = t
            end
            if frames == 0
                frames = [mesh;]
            else
                push!(frames,mesh)
            end
        end
        =#
        if i > 0 # Force Calculation
            ## Pressures (Change in pressure for each panel)
            for j = 1:length(dp)
                global dp[j] = surf.rho * ( ( (relVel[j,1] + surf.uind[j])*tau[j,1] + (relVel[j,2] + surf.wind[j])*tau[j,2] ) * circ[j]/ds[j] + circChange[j])
            end
            ## Loads
            # Cn normal
            ndem = .5 * surf.rho * surf.kinem.u^2 * surf.c # non-dimensionalization constant
            cn = sum( dp[j]*ds[j] for j = 1:length(dp) ) / ndem
            # LESP
            x = 1 - 2*ds[1] / surf.c
            surf.a0[1] = 1.13 * circ[1] / ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) )))
            # Cs suction
            cs = 2*pi*surf.a0[1]^2
            # Cl lift
            alpha = surf.kinem.alpha
            cl = cn*cos(alpha) + cs*sin(alpha)
            # Cd drag
            cd = cn*sin(alpha) - cs*cos(alpha)
            # Cm moment
            ndem = ndem * surf.c # change ndem to moment ndem
            # Cm = Cmz + Cmx
            cm = -sum( dp[j]*ds[j]*cosd(cam_slope[j])*(vor_loc[j,1]-surf.pvt) for j = 1:length(dp) ) / ndem + sum( dp[j]*ds[j]*sind(cam_slope[j])*vor_loc[j,2] for j = 1:length(dp) ) / ndem

            mat = hcat(mat,[t, surf.kinem.alpha*180/pi, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm]) # Pressures and loads

        end
        ## Wake Rollup
        IFR_field = wakeroll(IFR_surf,IFR_field,dtstar)

        # Convert IFR_field values to curfield (IFR -> BFR)
        push!(curfield.tev,TwoDVort(0,0,IFR_field.tev[end].s,0.02*surf.c,0.,0.))
        if size(IFR_field.lev,1) > 0 && surf.levflag[1] == 1
            push!(curfield.lev,TwoDVort(0,0,IFR_field.lev[end].s,0.02*surf.c,0.,0.))
        end
        vor_BFR(surf,t+dtstar,IFR_field,curfield) # Calculate IFR positions for next iteration

        # Draw animation forces and vorticies alonglside calculation
        if anim
            if i % animStep == 0. && i > 0
                clf()
                # Force Data
                animMat = [animMat;mat[:,end]']
                cl = animMat[:,6]
                cd = animMat[:,7]
                cm = animMat[:,8]
                tvec = animMat[:,1]
                tmax = nsteps*dtstar

                # Vorticies
                animMesh = meshgrid(surf,IFR_field.tev,IFR_field.lev,.25,t,100,view)

                #   Calculate velocities

                # Plotting
                subplot2grid((3,3),(0,0)) # Top left (Cl)
                plot(tvec,cl,color = "black")
                xlim(0,tmax)
                ylabel("Cl")

                subplot2grid((3,3),(1,0)) # Middle left (Cd)
                plot(tvec,cd,color = "black")
                xlim(0,tmax)
                ylabel("Cd")

                subplot2grid((3,3),(2,0)) # Bottom left (Cm)
                plot(tvec,cm,color = "black")
                xlim(0,tmax)
                ylabel("Cm")

                subplot2grid((3,3),(0,1), colspan=2,rowspan=3,aspect=1) # Vortex plot

                if length(animMesh.tev) > 0
                    scatter(animMesh.tev[:,1], animMesh.tev[:,2],c = "red", marker = "o", s = 1)
                end
                if length(animMesh.lev) > 0
                    scatter(animMesh.lev[:,1], animMesh.lev[:,2],c = "blue", marker = "o", s = 1)
                end
                plot(animMesh.camX,animMesh.camZ,color = "black", linewidth = 3)
                xmin = animMesh.x[1,1]
                xmax = animMesh.x[1,end]
                zmin = animMesh.z[end,1]
                zmax = animMesh.z[1,1]
                xlim(xmin,xmax)
                ylim(zmin,zmax)
                axis("off")

                tight_layout()
                pause(0.0001)
                show(block = false)
                print("") # needed becasue for some reason the plot wont update without printing to terminal
            end
        end

        # Write iteration step variables out if wflag == true
        if wflag
            if i % writeInterval == 0 && i > 0
                dirname = "$(round(t, sigdigits=6))"
                writeStamp(dirname, t, IFR_surf, IFR_field)
            end
        end

    end

    mat = mat'
    test1 = TEVhist
    test2 = LEVhist
    return frames,IFR_field,mat,test1,test2
end


text:
# Background image subtraction

function BG_correct(img, filename, GRIN)
    (x, y, z) = size(img)
    b_mean = zeros(z, 1)
    @showprogress @sprintf("Background calculating for Record %s...", filename) for i = 1:z
        b_mean[i] = mean(img[:, :, i])
    end
    base = maximum(b_mean)

    filtered = zeros(x, y, z)
    @showprogress @sprintf("Light sheet filtering of Record %s...", filename) for t =
        1:size(img, 3)
        img2 = img[:, :, t]

        for i = 1:x
            for j = 1:y
                filtered[i, j, t] = img2[i, j] - b_mean[t] + base
                if (filtered[i, j, t] >= 1)
                    filtered[i, j, t] = 1
                end
            end
        end
    end
    if (GRIN == 1)
        b_mean = zeros(x, y)
        for i = 1:x
            for j = 1:y
                b_mean[i, j] = mean(filtered[i, j, :])
            end
        end
        base = maximum(b_mean)

        @showprogress @sprintf("GRIN lens correction of Record %s...", filename) for t =
            1:size(img, 3)

            img2 = img1[:, :, t]
            for i = 1:x
                for j = 1:y
                    filtered[i, j, t] = filtered[i, j, t] - b_mean[i, j] + base
                    if (filtered[i, j, t] >= 1)
                        filtered[i, j, t] = 1
                    end
                end
            end
        end
    end
    img_save(filtered, "/home/jchang/image/result/", @sprintf("%s-b.tif", filename))
end


text:
import LinearAlgebra

@testset "$(format("Misc: check[:rank]"))" begin
    let
        tc = ones(5,5)
        @test check(Val{:rank},tc,1)
        @test check(:rank,tc,1)
        @test check(:rank,tc) == 1
        @test check(:rank,tc,2) == false
    end
end

@testset "$(format("Misc: check[:dimension]"))" begin
    let
        tc = rand(5,5)
        @test check(:dimension,tc) == (5,5)
        @test check(:dimension,tc,(5,5)) == true
        @test check(:dimension,tc,(4,5)) == false
    end
end

@testset "$(format("Misc: check[:l2difference]"))" begin
    # number case
    let
        for i = 1:5
            tc_a = rand() * 10
            tc_b = rand() * 10
            @test abs(check(:l2difference,tc_a,tc_b) - LinearAlgebra.norm(tc_a-tc_b,2)) < 1e-5
            @test check(:l2difference,tc_a,tc_b,LinearAlgebra.norm(tc_a-tc_b,2))
            @test abs(check(:l2diff,tc_a,tc_b) - LinearAlgebra.norm(tc_a-tc_b,2)) < 1e-5
            @test check(:l2diff,tc_a,tc_b,LinearAlgebra.norm(tc_a-tc_b,2))
        end
    end
    # vector case
    let
        for i = 1:5
            tc_a = rand(10) .* 10
            tc_b = rand(10) .* 10
            @test abs(check(:l2difference,tc_a,tc_b) - LinearAlgebra.norm(tc_a-tc_b,2)) < 1e-5
            @test check(:l2difference,tc_a,tc_b,LinearAlgebra.norm(tc_a-tc_b,2))
            @test abs(check(:l2diff,tc_a,tc_b) - LinearAlgebra.norm(tc_a-tc_b,2)) < 1e-5
            @test check(:l2diff,tc_a,tc_b,LinearAlgebra.norm(tc_a-tc_b,2))

        end
    end
    # matrix case
    let
        for i = 1:5
            tc_a = rand(10,10) .* 10
            tc_b = rand(10,10) .* 10
            @test abs(check(:l2difference,tc_a,tc_b) - LinearAlgebra.norm(tc_a-tc_b,2)) < 1e-5
            @test check(:l2difference,tc_a,tc_b,LinearAlgebra.norm(tc_a-tc_b,2))
            @test abs(check(:l2diff,tc_a,tc_b) - LinearAlgebra.norm(tc_a-tc_b,2)) < 1e-5
            @test check(:l2diff,tc_a,tc_b,LinearAlgebra.norm(tc_a-tc_b,2))
        end
    end
end

@testset "$(format("Misc: zeros"))" begin
    let
        # matrix case
        tc = rand(100,100)
        @test check(:l2diff,zeros(tc),zeros(100,100)) < 1e-3

    end
end


text:
function whitespacenumber(n)
    result = Char[]
    append!(result, n<0 ? '\t' : ' ' )
    if n ≠ 0
      for b in Iterators.reverse(digits(n,base=2))
          append!( result, b==0 ? ' ' : '\t' )
      end
    end
    append!(result,'\n')
    String(result)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using YAActL, .Threads, Distributed

function comp(i)  # a "heavy" computation
    sleep(2)
    return (id=i, thrd=threadid())
end

println("sequential execution:")
@time map(comp, 1:5)

println("with Julia's Threads:")
@time begin
    t = map(x->(Threads.@spawn comp(x)), 1:5)
    map(fetch, t)
end

println("with YAActL:")
A = map(i->Actor(parallel(), comp, i), 1:5); # start actors
@time begin
    foreach(a->call!(a, USR), A)
    map(x->receive!(USR).y, 1:5)
end

println("with Distributed:")
nworkers() < 5 && addprocs(4);  # add processes
@everywhere function comp(i)    # a "heavy" computation
    sleep(2)
    return (id=i, prc=myid())
end
@time begin
    f = map(i->(@spawnat i comp(i)), 1:5)
    map(fetch, f)
end

rmprocs(workers());


text:
using StanBase
using Test

if haskey(ENV, "JULIA_CMDSTAN_HOME")

  stan_prog = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
    real empty[0];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  ";

  println("\nRunning StanBase.jl tests")

  @testset "Basic HelpModel" begin

    sm = HelpModel( "help", stan_prog)

    return_codes = stan_sample(sm; n_chains=4)

    if success(return_codes)
      println()
      @test sm.method == StanBase.Help(:sample)
      @test StanBase.get_n_chains(sm) == 4
    end

  end
else
  println("\nJULIA_CMDSTAN_HOME not set. Skipping tests")
end


text:
@testset "Triangular" begin

    @test ndims(Triangular((10, 5))) == 2

    @testset "Constructors" begin
        # test default BC
        @test Triangular((10, 5)) == Triangular((10, 5), Periodic())

        for B1 in [Periodic, Open, Helical]
            @test (
                Triangular((10, 5), B1())
                == Triangular((10, 5), (B1(), B1()))
                == Triangular((10, 5), B1())
                == Triangular((10, 5), (B1(), B1()))
            )
            if B1 != Helical
                for B2 in [Periodic, Open]
                    @test (
                        Triangular((10, 5), (B1(), B2()))
                        == Triangular((10, 5), (B1(), B2()))
                    )
                end
                @test_throws ArgumentError Triangular((10, 5), (B1(), Helical()))
            end
        end
    end

    @testset "Translation Vectors" begin
        l = Triangular((10, 10))

        odd_vecs = Set([
            Coordinate(1, 0),
            Coordinate(-1, 0),
            Coordinate(0, 1),
            Coordinate(0, -1),
            Coordinate(1, -1),
            Coordinate(-1, 1)
        ])

        even_vecs = Set([
            Coordinate(2, -1),
            Coordinate(-2, 1),
            Coordinate(1, 1),
            Coordinate(-1, 2),
            Coordinate(1, -2),
            Coordinate(-1, -1),
        ])

        # run k=4 first to fill cache
        @test Set(translation_vectors(l, 4)) == Set([
            Coordinate(2, 1),
            Coordinate(-2, -1),
            Coordinate(3, -1),
            Coordinate(-3, 1),

            Coordinate(2, -3),
            Coordinate(-2, 3),
            Coordinate(3, -2),
            Coordinate(-3, 2),

            Coordinate(1, -3),
            Coordinate(-1, 3),
            Coordinate(-1, -2),
            Coordinate(1, 2),
        ])

        @test Set(translation_vectors(l, 1)) == odd_vecs
        @test Set(translation_vectors(l, 2)) == even_vecs
        @test Set(translation_vectors(l, 3)) == Set(2 .* odd_vecs)

    end

    @testset "Nearest Neighbors" begin
        l = Triangular((10, 10), Periodic())
        @test Set(neighbors(l, Coordinate(10, 10))) == Set([Coordinate(10, 9), Coordinate(10, 1), Coordinate(1, 10), Coordinate(9, 10), Coordinate(1, 9), Coordinate(9, 1)])
        @test Set(neighbors(l, Coordinate(5, 10))) == Set([Coordinate(5, 9), Coordinate(5, 1), Coordinate(4, 10), Coordinate(6, 10), Coordinate(6, 9), Coordinate(4, 1)])
        @test Set(neighbors(l, Coordinate(10, 5))) == Set([Coordinate(10, 4), Coordinate(10, 6), Coordinate(9, 5), Coordinate(1, 5), Coordinate(1, 4), Coordinate(9, 6)])
        @test Set(neighbors(l, Coordinate(5, 5))) == Set([Coordinate(5, 4), Coordinate(5, 6), Coordinate(4, 5), Coordinate(6, 5), Coordinate(4, 6), Coordinate(6, 4)])

        l = Triangular((10, 10), Open())
        @test Set(neighbors(l, Coordinate(10, 10))) == Set([Coordinate(10, 9), Coordinate(9, 10)])
        @test Set(neighbors(l, Coordinate(5, 10))) == Set([Coordinate(5, 9), Coordinate(4, 10), Coordinate(6, 10), Coordinate(6, 9)])
        @test Set(neighbors(l, Coordinate(10, 5))) == Set([Coordinate(10, 4), Coordinate(10, 6), Coordinate(9, 5), Coordinate(9, 6)])
        @test Set(neighbors(l, Coordinate(5, 5))) == Set([Coordinate(5, 4), Coordinate(5, 6), Coordinate(4, 5), Coordinate(6, 5), Coordinate(4, 6), Coordinate(6, 4)])
    end

    @testset "Next Nearest Neighbors" begin
        l = Triangular((10, 10), Periodic())
        @test Set(neighbors(l, Coordinate(10, 10), 2)) == Set([Coordinate(2, 9), Coordinate(8, 1), Coordinate(1, 1), Coordinate(9, 2), Coordinate(1, 8), Coordinate(9, 9)])
        @test Set(neighbors(l, Coordinate(5, 10), 2)) == Set([Coordinate(7, 9), Coordinate(3, 1), Coordinate(6, 1), Coordinate(4, 2), Coordinate(6, 8), Coordinate(4, 9)])
        @test Set(neighbors(l, Coordinate(10, 5), 2)) == Set([Coordinate(2, 4), Coordinate(8, 6), Coordinate(1, 6), Coordinate(9, 7), Coordinate(1, 3), Coordinate(9, 4)])
        @test Set(neighbors(l, Coordinate(5, 5), 2)) == Set([Coordinate(7, 4), Coordinate(3, 6), Coordinate(6, 6), Coordinate(4, 7), Coordinate(6, 3), Coordinate(4, 4)])

        l = Triangular((10, 10), Open())
        @test Set(neighbors(l, Coordinate(10, 10), 2)) == Set([Coordinate(9, 9)])
        @test Set(neighbors(l, Coordinate(5, 10), 2)) == Set([Coordinate(7, 9), Coordinate(6, 8), Coordinate(4, 9)])
        @test Set(neighbors(l, Coordinate(10, 5), 2)) == Set([Coordinate(8, 6), Coordinate(9, 7), Coordinate(9, 4)])
        @test Set(neighbors(l, Coordinate(5, 5), 2)) == Set([Coordinate(7, 4), Coordinate(3, 6), Coordinate(6, 6), Coordinate(4, 7), Coordinate(6, 3), Coordinate(4, 4)])
    end
end


text:
module VisitingICF2021
    include("Krav.jl")
    include("Wignerd.jl")
    include("LagKrav.jl")
    include("ImageAcq.jl")
    include("ImgRots.jl")
end


text:
__precompile__(true)

module Strategems
    using Dates
    using Temporal
    using Indicators
    using Random
    using ProgressMeter

    export
        # universe definitions
        Universe, gather!, get_overall_index,
        # parameter sets
        ParameterSet, count_runs, generate_combinations, generate_dict,
        # indicators
        Indicator, calculate,
        # signals
        Signal, prep_signal, ↑, ↓, @signal,
        # rules
        Rule, @rule, →,
        # portfolios
        Portfolio,#, update_portfolio!,
        # order
        AbstractOrder, MarketOrder, LimitOrder, StopOrder, liquidate, long, buy, short, sell,
        # strategy results
        Backtest,
        # summary statistic calculations
        cum_pnl,
        # strategies
        Strategy, generate_trades, generate_trades!, backtest, backtest!, optimize, optimize!, summarize_results

    include("model/universe.jl")
    include("model/paramset.jl")
    include("model/indicator.jl")
    include("model/signal.jl")
    include("model/rule.jl")
    include("model/portfolio.jl")
    include("model/order.jl")
    include("model/backtest.jl")
    include("model/strategy.jl")
    include("compute/backtest.jl")
    include("compute/optimize.jl")

end


text:
mutable struct Tuner
    n_blk_max::Int  ## maximum number of blocks
    blk_size ::Int  ## block size
    cur_blk  ::Int  ## current block
    cur_iter ::Int  ## number of iterations within current block
    ac_cnt   ::Int  ## number of acceptances within current block
    blk_acr  ::Vector{Float64} ## acceptance rates
    theta    ::Vector{Float64} ## tuning parameter (standard deviation)
end

```
       Tuner(warmup::Int)

Create Tuner object. Tuning aims to adjust the tuning parameters during the `warmup` iterations such that an acceptance rate of 25% is achieved.
It does this by fitting a linear regression to blockwise acceptance rates and the tuning parameters. No tuning is performed when there are less
than 100 warmup iterations.
```
function Tuner(warmup::Int)
    if warmup < 100
        ## no tuning if fewer than 100 warmup iterations
        Tuner(1, -1, 1, 1, 0, [0.0], [2.38])
    else
        blk_size = Int(floor(warmup/25))
        n_blk_max = Int(floor(warmup/blk_size))
        theta = zeros(Float64, n_blk_max+1)
        theta[1:2] = [2.38, 0.1]

        Tuner(n_blk_max, blk_size, 1, 1, 0, zeros(Int, n_blk_max), theta)
    end
end

get_tune_param(x::Tuner) = x.theta[x.cur_blk]

```
    update!(x::Tuner, ac::Bool)

Update the tuner `x` with acceptance status `ac`. It simply counts the number of acceptances in the current block. If the end of the
current block is reached the acceptance rate is calculated and stored and a new block is started. After at least two blocks have been recorded the
tuning parameter is updated via linear regression for every new block.
```
function update!(x::Tuner, ac::Bool)
    x.ac_cnt += ac

    b = x.blk_size
    cb = x.cur_blk

    if x.cur_blk <= x.n_blk_max
        if x.cur_iter == b ## update tuning parameter every b iterations
            x.blk_acr[cb] = x.ac_cnt/b ## acceptance rate of last b iteration
            x.ac_cnt = 0 ## reset acceptance counter
            if cb >= 2 ## need at least two blocks for linear regression
                x.theta[cb+1] = linreg(x.blk_acr[1:cb], x.theta[1:cb])

                ## gamma_theta is standard deviation
                if x.theta[cb+1] <= 0
                    x.theta[cb+1] = x.theta[cb]
                end
            end
            x.cur_iter = 1
            x.cur_blk += 1
        else
            x.cur_iter += 1
        end
    end
end

function linreg(accept::Vector{Float64}, gamma::Vector{Float64})
    @assert all(accept .<= 1.0)

    X = hcat(ones(Float64, length(accept)), accept)

    ## U = cholesky!(X'X).U
    ## w = U' \ X'Y
    ## beta = U \ X'w

    Z = X'X

    if rank(Z) == 2
        beta = cholesky!(Z) \ X'gamma
        beta[1] + 0.25*beta[2]
    else
        mean(gamma)
    end
end


text:
module PolymakeOscar

using Polymake, Hecke, Markdown

Hecke.nrows(A::Polymake.MatrixAllocated) = Int(size(A)[1])
Hecke.ncols(A::Polymake.MatrixAllocated) = Int(size(A)[2])

function _polytope(; A::fmpz_mat=zero_matrix(FlintZZ, 1, 1), b::fmpz_mat=zero_matrix(FlintZZ, ncols(A), 1), C::fmpz_mat=zero_matrix(FlintZZ, 1, 1))
  if !iszero(A)
    bA = Matrix{BigInt}(hcat(-b, A))
    z = findall(i->!iszero_row(bA, i), 1:nrows(bA))
    zbA = Matrix{BigInt}(bA[z, :])
  else
    zbA = Matrix{BigInt}(undef, 0, 0)
  end
  if !iszero(C)
    z = findall(i->!iszero_row(C, i), 1:nrows(C))
    zI = Matrix{BigInt}(hcat(zero_matrix(FlintZZ, nrows(C), 1), C))[z, :]
  else
    zI = Matrix{BigInt}(undef, 0, 0)
  end
  if length(zbA) == 0
    p =  polytope.Polytope(INEQUALITIES = zI)
  else
    if nrows(zI) == 0
      p = polytope.Polytope(EQUATIONS = zbA)
    else
      p = polytope.Polytope(EQUATIONS = zbA, INEQUALITIES = zI)
    end
  end
  return p
end

@doc Markdown.doc"""
    solve_ineq(A::fmpz_mat, b::fmpz_mat)

Solves $Ax<=b$, assumes finite set of solutions.
"""
function solve_ineq(A::fmpz_mat, b::fmpz_mat)
  p = _polytope(C = hcat(b, -A))
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

@doc Markdown.doc"""
    solve_non_negative(A::fmpz_mat, b::fmpz_mat)

Finds all solutions to $Ax = b$, $x>=0$. Assumes a finite set of solutions.
"""
function solve_non_negative(A::fmpz_mat, b::fmpz_mat)
  p = _polytope(A = A, b = b, C = identity_matrix(FlintZZ, ncols(A)))
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

@doc Markdown.doc"""
    solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)

Solves $Ax = b$ under $Cx >= 0$, assumes a finite solution set.
"""
function solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)  # Ax == b && Cx >= 0
  p = _polytope(A = A, b = b, C = C)
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    if out[i,1] != 1
      println("unbounded polytope!!")
      global last_in = (A, b, C)
    end
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    if inner[i,1] != 1
      println("unbounded polytope!!")
      global last_in = (A, b, C)
    end
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

@doc Markdown.doc"""
    solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat, d::fmpz_mat)

Solves $Ax = b$ under $Cx >= d$, assumes a finite solution set.
"""
function solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat, d::fmpz_mat)
  n = ncols(A)
  A = cat(A, identity_matrix(FlintZZ, ncols(d)), dims=(1,2))
  b = vcat(b, identity_matrix(FlintZZ, ncols(d)))
  C = [C -d; zero_matrix(FlintZZ, ncols(d), ncols(C)) identity_matrix(FlintZZ, ncols(d))]
  s = solve_mixed(A, b, C)
  return s[:, 1:n]
end

function Hecke.valuation(a::FacElem{fmpz, FlintIntegerRing}, p::fmpz)
  return sum(k*valuation(b, p) for (b, k) = a.fac)
end

function norm_equation2_fac_elem(R::NfAbsOrd, k::fmpz; abs::Bool = false)
  @assert Hecke.ismaximal(R)
  lp = factor(k*R)
  s, ms = Hecke.sunit_mod_units_group_fac_elem(collect(keys(lp)))
  C = vcat([matrix(FlintZZ, 1, ngens(s), [valuation(ms(s[i]), p) for i=1:ngens(s)]) for p = keys(lp)])

  lp = factor(k)
  A = vcat([matrix(FlintZZ, 1, ngens(s), [valuation(Hecke.factored_norm(ms(s[i])), p) for i=1:ngens(s)]) for p = keys(lp.fac)])
  b = matrix(FlintZZ, length(lp.fac), 1, [valuation(k, p) for p = keys(lp.fac)])

  so = solve_mixed(A, b, C)
  sol = [ms(s(sub(so, i:i, 1:ncols(so)))) for i=1:nrows(so)]

  if !abs
    u, mu = unit_group_fac_elem(R)
    i = findfirst(x -> norm(mu(x)) == -1, gens(u))
    ns = [norm(x) for x = sol]
    if i === nothing
      return [sol[i] for i in 1:length(sol) if ns[i] == k]
    end
    U = mu(u[i])
    return [ ns[i] == k ? sol[i] : sol[i] * U for i = 1:length(sol)]
  end
  return sol
end

function norm_equation_fac_elem(R::NfAbsOrd, k::fmpz; abs::Bool = false)
  @assert Hecke.ismaximal(R)
  lp = factor(k)
  S = []
  for (p, k) = lp.fac
    P = prime_decomposition(R, p)
    s = solve_non_negative(matrix(FlintZZ, 1, length(P), [degree(x[1]) for x = P]), matrix(FlintZZ, 1, 1, [k]))
    push!(S, (P, [view(s, i:i, 1:ncols(s)) for i=1:nrows(s)]))
  end
  sol = []
  for x in Base.Iterators.ProductIterator(Tuple(t[2] for t = S))
    I = ideal(R, 1)
    for i = 1:length(S)
      I *= prod(S[i][1][j][1]^Int(x[i][j]) for j=1:length(S[i][1]))
    end
    fl, g = Hecke.isprincipal_fac_elem(I)
    if fl
      push!(sol, g)
    end
  end
  if !abs
    u, mu = unit_group_fac_elem(R)
    i = findfirst(x -> norm(mu(x)) == -1, gens(u))
    ns = [norm(x) for x = sol]
    if i === nothing
      return [sol[i] for i in 1:length(sol) if ns[i] == k]
    end
    U = mu(u[i])
    return [ ns[i] == k ? sol[i] : sol[i] * U for i = 1:length(sol)]
  end
  return sol
end

norm_equation_fac_elem(R::NfAbsOrd, k::Integer; abs::Bool = false) =
                            norm_equation_fac_elem(R, fmpz(k), abs = abs)

function norm_equation(R::NfAbsOrd, k::fmpz; abs::Bool = false)
  s = norm_equation_fac_elem(R, k, abs = abs)
  return [R(evaluate(x)) for x = s]
end

norm_equation(R::NfAbsOrd, k::Integer; abs::Bool = false) = norm_equation(R, fmpz(k), abs = abs)

function norm_equation_fac_elem(R::Hecke.NfRelOrd{nf_elem,Hecke.NfOrdFracIdl}, a::NfAbsOrdElem{AnticNumberField,nf_elem})

  @assert Hecke.ismaximal(R)
  Ka, mKa, mkK = collapse_top_layer(nf(R))
  Ra = maximal_order(Ka)
  class_group(Ra)
  k = nf(parent(a))
  class_group(parent(a))

  lp = factor(Ra(mkK(k(a)))*Ra)
  la = factor(a*parent(a))
  S, mS = Hecke.sunit_mod_units_group_fac_elem(collect(keys(lp)))
  s, ms = Hecke.sunit_mod_units_group_fac_elem(collect(keys(la)))
  No = hom(S, s, elem_type(s)[preimage(ms, norm(mkK, mS(g))) for g = gens(S)])
  q, mq = quo(S, kernel(No)[1])
  q, mms = snf(q)
  mq = mq*inv(mms)

  C = vcat([matrix(FlintZZ, 1, ngens(q), [valuation(mS(preimage(mq, q[i])), p) for i=1:ngens(q)]) for p = keys(lp)])

  A = vcat([matrix(FlintZZ, 1, ngens(q), [valuation(norm(mkK, mS(preimage(mq, g))), p) for g in gens(q)]) for p = keys(la)])
  b = matrix(FlintZZ, length(la), 1, [valuation(a, p) for p = keys(la)])

  so = solve_mixed(A, b, C)
  u, mu = Hecke.unit_group_fac_elem(parent(a))
  U, mU = Hecke.unit_group_fac_elem(Ra)
  No = hom(U, u, elem_type(u)[preimage(mu, norm(mkK, mU(g))) for g = gens(U)])
  sol = []
  for i = 1:nrows(so)
    aa = mS(preimage(mq, q(sub(so, i:i, 1:ncols(so)))))
    b = norm(mkK, aa)
    c = b*inv(FacElem(k(a)))
    d = preimage(mu, c)
    fl, p = haspreimage(No, d)
    if fl
      push!(sol, FacElem(Dict(mKa(x) => v for (x, v) = (aa*inv(mU(p))).fac)))
    end
  end

  return sol
end

function Hecke.isirreducible(a::NfAbsOrdElem{AnticNumberField,nf_elem})
  if iszero(a)
    return false
  end
  O = parent(a)
  S = collect(keys(factor(a*O)))
  if length(S) == 0
    return false
  end
  s, ms = Hecke.sunit_mod_units_group_fac_elem(S)
  V = matrix([fmpz[valuation(ms(x), y) for y = S] for x = gens(s)])
  b = matrix([fmpz[valuation(a, y) for y = S]])
  sol = solve(V, b)

  #want to write sol = x+y where
  # Cx, Cy > 0
  # if this is possible, then a is not irreducible as a
  # is then ms(Ax) * ms(Ay) and neither is trivial.

  I = identity_matrix(FlintZZ, length(S))
  A = hcat(I, I)
  #so A*(x|y) = x+y = sol is the  1. condition
  C = cat(V, V, dims=(1,2))
  # C(x|y) >=0 iff Cx >=0 and Cy >=0
  #Cx <> 0 iff (1,..1)*Cx >= 1
  one = matrix(FlintZZ, 1, length(S), [1 for p = S])
  zer = matrix(FlintZZ, 1, length(S), [0 for p = S])
  C = vcat(C, hcat(one, zer), hcat(zer, one))
  d = matrix(FlintZZ, 2*length(S)+2, 1, [0 for i = 1:2*length(S) + 2])
  d[end-1, 1] = 1
  d[end, 1] = 1
  global last_irr = a
  pt = solve_mixed(A, sol, C, d)
  return nrows(pt) == 0
end

@doc Markdown.doc"""
    irreducibles(S::Vector{NfAbsOrdIdl{AnticNumberField,nf_elem}}) -> Vector{NfAbsOrdElem}

Computes all irreducibles whose support is contained in $S$.
"""
function irreducibles(S::Vector{NfAbsOrdIdl{AnticNumberField,nf_elem}})
  if length(S) == 0
    return []
  end
  @assert all(isprime, S)
  #TODO: try to get a better bound - in general if S is too large
  #      it cannot work

  O = order(S[1])
  @assert all(x-> order(x) == O, S)

  s, ms = Hecke.sunit_mod_units_group_fac_elem(S)
  if length(S) == 1
    return [O(evaluate(ms(s[1])))]
  end
  c, mc = class_group(O)

  V = matrix([[valuation(ms(x), y) for y = S] for x = gens(s)])

  p = _polytope(C = V)
  z = p.HILBERT_BASIS_GENERATORS
  @assert nrows(z[2]) == 0 #no idea....
  res = [O(evaluate(ms(s(map(fmpz, Array(z[1][i, 2:end])))))) for i=1:nrows(z[1]) if z[1][i,1] == 0]
  return res
end

@doc Markdown.doc"""
    factorisations(a::NfAbsOrdElem{AnticNumberField,nf_elem}) -> Vector{Fac{OrdElem}}

Computes all factorisations of $a$ into irreducibles.
"""
function factorisations(a::NfAbsOrdElem{AnticNumberField,nf_elem})
  O = parent(a)
  S = collect(keys(factor(a*O)))
  if length(S) == 0
    return []
  end
  irr = irreducibles(S)
  A = matrix([fmpz[valuation(x, y) for y = S] for x = irr])
  b = matrix([fmpz[valuation(a, y) for y = S]])
  sol = solve_non_negative(A, b)
  res = Fac{NfAbsOrdElem{AnticNumberField,nf_elem}}[]
  for j=1:nrows(sol)
    x = Dict{typeof(a), Int}()
    y = a
    for i=1:length(irr)
      if sol[j,i] != 0
        x[irr[i]] = sol[j,i]
        y = divexact(y, irr[i]^sol[j,i])
      end
    end
    push!(res, Fac(y, x))
  end
  return res
end

function Base.lastindex(a::fmpz_mat, i::Int)
  i==1 && return nrows(a)
  i==2 && return ncols(a)
  error("illegal dimension")
end

export irreducibles, factorisations

end

using Main.PolymakeOscar


text:
# adding using TimeModels to this file is odd, but it's a trick to get travis to pass

using Compat.Test
using TimeModels

include("kalman_filter.jl")
include("parameter_estimation.jl")


text:
using Alistair

# PERFORMANCE COMPARISON OPTIONS
const TEST_SIZE_1 = 1000

Profile.clear()

#=
println("---- LINEAR REGRESSION ----")

beta = [0.7; 0.45]
X = hcat(ones(TEST_SIZE_1,1), collect(1.0:TEST_SIZE_1))
Y = X * beta + 0.1*randn(TEST_SIZE_1, 1)

olsreg = OLS()

olsreg(X, Y)
@time res = olsreg(X, Y)

linregress(X, Y, olsreg)
@time res = linregress(X, Y, olsreg)

linregress(X, Y, OLS(intercept=true, robust=HCEVariance()))
@time res = linregress(X, Y, OLS(intercept=true, robust=HCEVariance()))

linregress(X, Y, FGLS(intercept=true, robust=HCEVariance()))
@time res = linregress(X, Y, FGLS(intercept=true, robust=HCEVariance()))

linregress(X, Y, IteratedFGLS(10, intercept=true, robust=HCEVariance()))
@time res = linregress(X, Y, IteratedFGLS(10, intercept=true, robust=HCEVariance()))

println("---- NONLINEAR REGRESSION ----")

f(X_, beta_) = exp.(beta_[1]*X_[:,1] + beta_[2]*log.(X_[:,2]))

Ynlin = f(X, beta) + 0.1*randn(TEST_SIZE_1, 1)

OptimizeRegression = Optimize(f, [1.0, 1.0])

res = nlinfit(X, Ynlin, OptimizeRegression)
@time res = nlinfit(X, Ynlin, OptimizeRegression)
=#

println("---- GENERAL LINEAR MODEL REGRESSION ----")

beta_cat = [0.5; 2.0]
X_cat = hcat(ones(TEST_SIZE_1,1), collect(linspace(-2.2,2.0,TEST_SIZE_1)))
Y_cat = round.(Alistair.logistic(X_cat * beta_cat))

res_logit = glmregress(X_cat, Y_cat, Logit())
@time res_logit = glmregress(X_cat, Y_cat, Logit())
println(res_logit)


text:
@testset "Annotations" begin

    @testset "Empty annotations" begin

        annot = Annotations()

        @test length(annot.file) == 0
        @test length(annot.sequences) == 0
        @test length(annot.columns) == 0
        @test length(annot.residues) == 0

        @test length(annot) == 0

        @test ncolumns(annot) == -1

        @test isempty(annot)
    end

    @testset "Getters & Setters" begin

        annot = Annotations()
        example_str = "CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH"

        setannotresidue!(annot, "O31698/18-71", "SS", example_str)

        @test_throws AssertionError setannotresidue!(annot, "O31698/18-71",
                                                    String(rand('A':'Z',51)), example_str)

        @test_throws AssertionError setannotresidue!(annot, "O31698/18-71", "Feature Name",
                                                    example_str)

        @test ncolumns(annot) == 37

        setannotfile!(annot, "AC", "PF00571")
        setannotcolumn!(annot, "SS_cons", example_str)
        setannotsequence!(annot, "O31698/88-139", "OS", "Bacillus subtilis")

        @test getannotfile(annot, "AC") == "PF00571"
        @test getannotcolumn(annot, "SS_cons") == example_str
        @test getannotsequence(annot, "O31698/88-139", "OS") == "Bacillus subtilis"
        @test getannotresidue(annot, "O31698/18-71", "SS") == example_str

        @test getannotfile(annot, "An", "Default") == "Default"
        @test getannotcolumn(annot, "Other", "Default") == "Default"
        @test getannotsequence(annot, "O31698/1-88", "OS", "Default") == "Default"
        @test getannotresidue(annot, "O31698/1-88", "SS", "Default") == "Default"

        @test ncolumns(annot) == 37

        @test_throws DimensionMismatch  setannotresidue!(annot,
                                        "O31698/18-71", "AS", "__*__")

        @test_throws DimensionMismatch  setannotcolumn!(annot, "SS_cons",
                                        "---CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH---")
    end

    @testset "Copy, deepcopy and empty!" begin

        annot = Annotations()
        setannotresidue!(annot,"O31698/18-71","SS","CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH")
        setannotfile!(annot, "AC", "PF00571")
        setannotcolumn!(annot, "SS_cons", "CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH")
        setannotsequence!(annot, "O31698/88-139", "OS", "Bacillus subtilis")

        copy_annot = copy(annot)
        @test copy_annot == annot
        empty!(copy_annot)
        @test ncolumns(annot) == 37
        @test ncolumns(copy_annot) == -1

        @test !isempty(annot)
        @test isempty(copy_annot)

        deepcopy_annot = deepcopy(annot)
        @test deepcopy_annot == annot
        empty!(deepcopy_annot)
        @test ncolumns(annot) == 37
        @test ncolumns(deepcopy_annot) == -1

        @test !isempty(annot)
        @test isempty(deepcopy_annot)
    end

    @testset "Filter" begin

        annot = Annotations()
        setannotresidue!(annot,"O31698/18-71","SS","CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH")
        setannotfile!(annot, "AC", "PF00571")
        setannotcolumn!(annot, "SS_cons", "CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH")
        setannotsequence!(annot, "O31698/88-139", "OS", "Bacillus subtilis")

        @test_throws AssertionError filtercolumns!(annot, [true, false, true])

        #filtersequences!(annot, IndexedArray(["O31698/88-139", "O31698/18-71"]), [false, true])
        #@test length( getannotsequence(annot) ) == 0
        #filtersequences!(annot, IndexedArray(["O31698/88-139", "O31698/18-71"]), [true, false])
        #@test length( getannotresidue(annot) ) == 0

        mask = collect("CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH") .!= Ref('E')
        filtercolumns!(annot, mask)
        @test ncolumns(annot) == 19
        @test getannotcolumn(annot, "SS_cons") == "CCCCCHHHHHHHHHHHHHH"

        filtercolumns!(annot, [1,2,19])
        @test ncolumns(annot) == 3
        @test getannotcolumn(annot, "SS_cons") == "CCH"
    end

    @testset "Print" begin

        output_string = """
                        #=GF AC	PF00571
                        #=GS O31698/88-139	OS	Bacillus subtilis
                        #=GR O31698/18-71	SS	CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH
                        #=GC SS_cons			CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH
                        """

        io = IOBuffer()

        annot = Annotations()
        setannotresidue!(annot,"O31698/18-71","SS","CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH")
        setannotfile!(annot, "AC", "PF00571")
        setannotcolumn!(annot, "SS_cons", "CCCCCHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEH")
        setannotsequence!(annot, "O31698/88-139", "OS", "Bacillus subtilis")

        print(io, annot)
        @test String(take!(io)) == output_string

        show(io, annot)
        @test String(take!(io)) == output_string
    end
end


text:
# I am still not sure LRUStore should actually a store or live
# somewhere else. We are actually not extending the storage interface
# because we cache uncompressed chunks of data, which is already ones
# step further down the processing pipeline. Maybe there should rather
# be a mini-interface for the step of copying uncompressed data into
# the resulting array. However, let's try how this implementation
# works out and modify later if necessary
# TODO: add writing of data

using LRUCache: LRU
struct LRUStore{S<:AbstractStore} <: AbstractStore
    parent::S
    lru::LRU{String,Any}
end
function LRUStore(s::AbstractStore; maxsize=5)
    lru = LRU{String,Any}(maxsize=maxsize)
    LRUStore(s,lru)
end
storagesize(d::LRUStore,p) = storagesize(d.parent,p)
Base.getindex(d::LRUStore,i::AbstractString) = d.parent[i]
Base.setindex!(d::LRUStore,v,i::AbstractString) = error("Writing to LRU stores is currently not supported")
Base.delete!(d::LRUStore, i::AbstractString) = error("Writing to LRU stores is currently not supported")
subdirs(d::LRUStore,p) = subdirs(d.parent,p)
#subkeys(d::LRUStore,p) = subkeys(d.parent,p)

function readchunk!(a::DenseArray,z::ZArray{<:Any,N,<:Zarr.Compressor,<:LRUStore},i::CartesianIndex{N}) where N
    length(a) == prod(z.metadata.chunks) || throw(DimensionMismatch("Array size does not equal chunk size"))
    k = _concatpath(z.path,citostring(i))
    if haskey(z.storage.lru,k)
        a .= z.storage.lru[k]
        return a
    end
    curchunk = z.storage[z.path,i]
    if curchunk === nothing
        fill!(a, z.metadata.fill_value)
    else
        zuncompress!(a, curchunk, z.metadata.compressor, z.metadata.filters)
    end
    z.storage.lru[k] = a
    a
end


text:
const BLOCKQUOTE_ACC = union((:BLOCKQUOTE_LINE,), INLINE_BLOCKS)
const LIST_ACC       = union((:ITEM_U_CAND, :ITEM_O_CAND), INLINE_BLOCKS)
const TABLE_ACC      = union((:TABLE_ROW_CAND,), [b for b in INLINE_BLOCKS if b != :TEXT])

"""
    aggregate!(blocks, items, acc, case)

Merge a bunch of blocks into a parent block. For instance at this point there
may be multiple BLOCKQUOTE_LINE, this function will aggregate them into one
BLOCKQUOTE block.

## Arguments

    * blocks: the current vector of blocks we're working with
    * items:  list of names of blocks that would trigger the aggregation
    * acc:    list of names of blocks that would be taken in the aggregation
    * case:   name of the block resulting from the aggregation
"""
function aggregate!(blocks::Vector{Block}, items::Vector{Symbol},
                    acc::Vector{Symbol}, case::Symbol)
    i = 1
    j = 1
    ps = parent_string(first(blocks))
    @inbounds while i ≤ length(blocks)
        bi = blocks[i]
        if bi.name in items
            # look ahead
            j = i + 1
            while j ≤ length(blocks)
                if blocks[j].name ∉ acc
                    break
                end
                j += 1
            end
            blocks[i] = Block(
                case,
                subs(ps, from(bi), to(blocks[j-1]))
            )
        end
        i = max(j, i + 1)
    end
end

form_blockquotes!(blocks::Vector{Block}) =
    aggregate!(blocks, [:BLOCKQUOTE_LINE], BLOCKQUOTE_ACC, :BLOCKQUOTE)

form_lists!(blocks::Vector{Block}) =
    aggregate!(blocks, [:ITEM_U_CAND, :ITEM_O_CAND], LIST_ACC, :LIST)

form_tables!(blocks::Vector{Block}) =
    aggregate!(blocks, [:TABLE_ROW_CAND], TABLE_ACC, :TABLE)

form_refs!(blocks::Vector{Block}) =
    aggregate!(blocks, [:REF], INLINE_BLOCKS, :REF)


text:
__precompile__()

module Potree

    using Common
    using LasIO
    using LazIO
    using JSON
    using PlyIO
    using DataStructures
    using Dates
    using Printf
    using Distributed
    #using PointClouds

    # util
    include("struct.jl")
    include("utilities.jl")

    # load
    include("Load/hierarchy.jl")
    include("Load/json.jl")
    include("Load/las.jl")
    include("Load/bin.jl")

    export potree2trie, las2pointcloud, DataStructures, LasIO, LazIO, JSON, getmodel, CloudMetadata,
    HEADER_SIZE, DATA_OFFSET, SIZE_DATARECORD, bin2pointcloud, bin2pointcloudNoMultithreading, las2larpoints, las2color, bin2pointcloudNEW

end


text:
export quadraticFieldmap, stairsFieldmap, polarFieldmap

"""
    quadraticFieldmap(Nx::Int64, Ny::Int64, maxOffresonance::Float64=125.0)

Computes a parabolic fieldmap.
"""
function quadraticFieldmap(Nx::Int64, Ny::Int64, maxOffresonance::Float64=125.0)
  xx = range(-1, stop=1, length=Nx)
  yy = range(-1, stop=1, length=Ny)
  fieldmap = zeros(Nx,Ny, 1)
  for nx=1:Nx
    for ny=1:Ny
      fieldmap[nx,ny,1] = maxOffresonance*(xx[nx]^2 + yy[ny]^2) - maxOffresonance
    end
  end
  return fieldmap
end

"""
Computes a polaroid fieldmap.
"""
function polarFieldmap(Nx::Int64, Ny::Int64,
                        maxOffresonance::Float64=125.0;
                        xCenter::Int64=0,
                        yCenter::Int64=0
                       )
  fieldmap = zeros(Nx,Ny,1)
  xx = range(-1, stop=1, length=Nx)
  yy = range(-1, stop=1, length=Ny)
  for nx=1:Nx
    for ny=1:Ny
      fieldmap[nx,ny,1] = maxOffresonance*sqrt(xx[nx]^2 + yy[ny]^2)- maxOffresonance
    end
  end
  return fieldmap
end

"""
Computes a fieldmap consisting of stairs, which has discrete linear
ascending/descending values
"""
function stairsFieldmap(Nx::Int64, Ny::Int64;
                        numStairs::Int64=8,
                        minOffresonance::Float64=-125.0,
                        maxOffresonance::Float64=125.0,
                        ascending=true)
  # Input check for the number of stairs
  if numStairs < 2 || numStairs > (Nx > Ny ? Ny : Nx)
    error("Number of stairs not used correctly")
  end
  # Allocating output
  fieldmap = zeros(Nx,Ny)
  # Determine the height and width of the stairs (e.g for NxN -> deltaX = deltaY)
  deltaX = Int64(floor((Nx / (numStairs*2) )))
  deltaY = Int64(floor((Ny / (numStairs*2) )))

  if ascending == true
    scaling = range(minOffresonance,stop=maxOffresonance,length=numStairs)
  else
    scaling = -range(minOffresonance,stop=maxOffresonance,length=numStairs)
  end

  # Filling (Quarter) Matrix for constructing the stairs
  partialFieldMap = zeros(Int64(floor((Nx/2))),Int64(floor((Ny/2))))

  for i=1:numStairs
    partialFieldMap[(i-1)*deltaX+1 : end     , (i-1)*deltaY+1 : i*deltaY  ] .= -1 * scaling[i]
    partialFieldMap[(i-1)*deltaX+1 : i*deltaX, (i-1)*deltaY+1 : end ]       .= -1 * scaling[i]
  end

  # Correction if numOfStairs is not divisible by 4 fill the lower right corner of
  # the upper left (quarter) partial matrix
  if numStairs % 4 != 0
    partialFieldMap[numStairs*deltaX:end , numStairs*deltaY:end] .= -1 * scaling[numStairs]
  end

  # If width or height is not divisible by two the middle row or column(or both)
  # must be filled manually
  partSizeX , partSizeY = size(partialFieldMap)
  xRemainder, yRemainder = Nx%2, Ny%2

  # Rotating and mirroring the filled quarter of the matrix to complete the stairs
  fieldmap[1:partSizeX                , 1: partSizeY]                 = partialFieldMap
  fieldmap[1:partSizeX                , partSizeY+1+yRemainder : end] = partialFieldMap[:,end:-1:1]
  fieldmap[partSizeX+1+xRemainder:end , 1: partSizeY]                 = partialFieldMap[end:-1:1,:]
  fieldmap[partSizeX+1+xRemainder:end , partSizeY+1+yRemainder : end] = rot180(partialFieldMap)

  # Filling the row or column,needed when number of rows/columns is odd
  if xRemainder != 0 && yRemainder != 0
    midXIndx, midYIndx  = Int64(ceil(Nx/2)), Int64(ceil(Ny/2))
    fieldmap[midXIndx,: ] = fieldmap[midXIndx-1,:]
    fieldmap[:,midYIndx ] = fieldmap[:,midYIndx-1]
  elseif xRemainder != 0
    midXIndx = Int64(ceil(Nx/2))
    fieldmap[midXIndx,: ] = fieldmap[midXIndx-1,:]
  elseif yRemainder != 0
    midYIndx = Int64(ceil(Ny/2))
    fieldmap[:,midYIndx ] = fieldmap[:,midYIndx-1]
  end

  return reshape(fieldmap,Nx,Ny,1)
end


text:
"""
    regularized_fit_data

Returns a dictionary `D` with data for the regularized data fitting
example on page 329 of VMLS.  The items in the dictionary are:
* `D["xtrain"]`: vector of length 10
* `D["ytrain"]`: vector of length 10
* `D["xtest"]`: vector of length 20
* `D["ytest"]`: vector of length 20.

"""
function regularized_fit_data()
    X = [
    1.06830e-01  1.82422e+00;
    2.55780e-01  8.13549e-01;
    2.71120e-01  9.24060e-01;
    3.02990e-01  4.20670e-01;
    5.02040e-01  4.46178e-01;
    6.58450e-01 -3.73407e-02;
    6.84230e-01 -2.10935e-01;
    7.02590e-01 -1.03327e-01;
    7.34060e-01  3.32097e-01;
    9.40350e-01  2.29278e+00
    ];
    xtrain = X[:,1];
    ytrain = X[:,2];
    Xtest = [
    3.96674e-01 -1.05772e+00;
    7.77517e-01  8.79117e-01;
    1.18400e-01  1.98136e+00;
    2.23266e-01  8.67012e-01;
    9.01463e-01  2.13650e+00;
    3.58033e-01 -7.01948e-01;
    2.60402e-01  9.41469e-01;
    8.04281e-01  1.49755e+00;
    6.31664e-01  5.50205e-01;
    1.49704e-01  1.34245e+00;
    5.51333e-01  1.21445e+00;
    6.63999e-01  4.49111e-03;
    1.64948e-01  9.57535e-01;
    6.51698e-01  7.78640e-02;
    1.23026e-01  1.73558e+00;
    3.37066e-01 -3.25244e-01;
    8.32080e-02  2.56555e+00;
    2.04422e-01  1.10081e+00;
    9.78000e-01  2.70644e+00;
    4.03676e-01 -1.10049e+00
    ];
    xtest = Xtest[:,1];
    ytest = Xtest[:,2];
    return Dict([ ("xtrain", xtrain), ("ytrain", ytrain),
        ("xtest", xtest), ("ytest", ytest) ])
end

max_issues_repo_name: lynnlangit/julia-linear-algebra


text:
using Test, Logging, Random, JuMP
using DecisionProgramming

function influence_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int}, n_inactive::Int)
    diagram = InfluenceDiagram()
    random_diagram!(rng, diagram, n_C, n_D, n_V, m_C, m_D, states)
    for c in diagram.C
        random_probabilities!(rng, diagram, c; n_inactive=n_inactive)
    end
    for v in diagram.V
        random_utilities!(rng, diagram, v; low=-1.0, high=1.0)
    end

    # Names needed for printing functions only
    diagram.Names = ["node$j" for j in 1:n_C+n_D+n_V]
    diagram.States = [["s$s" for s in 1:n_s] for n_s in diagram.S]

    diagram.P = DefaultPathProbability(diagram.C, diagram.I_j[diagram.C], diagram.X)
    diagram.U = DefaultPathUtility(diagram.I_j[diagram.V], diagram.Y)

    return diagram
end

function test_decision_model(diagram, n_inactive, probability_scale_factor, probability_cut)
    model = Model()

    @info "Testing DecisionVariables"
    z = DecisionVariables(model, diagram)

    @info "Testing PathCompatibilityVariables"
    if probability_scale_factor > 0
        x_s = PathCompatibilityVariables(model, diagram, z; probability_cut = probability_cut, probability_scale_factor = probability_scale_factor)
    else
        @test_throws DomainError x_s = PathCompatibilityVariables(model, diagram, z; probability_cut = probability_cut, probability_scale_factor = probability_scale_factor)
    end

    x_s = PathCompatibilityVariables(model, diagram, z; probability_cut = probability_cut, probability_scale_factor = 1.0)

    @info "Testing probability_cut"
    lazy_probability_cut(model, diagram, x_s)

    @info "Testing expected_value"
    EV = expected_value(model, diagram, x_s)

    @info "Testing conditional_value_at_risk"
    if probability_scale_factor > 0
        CVaR = conditional_value_at_risk(model, diagram, x_s, 0.2; probability_scale_factor = probability_scale_factor)
    else
        @test_throws DomainError conditional_value_at_risk(model, diagram, x_s, 0.2; probability_scale_factor = probability_scale_factor)
    end

    @test true
end

function test_analysis_and_printing(diagram)
    @info("Creating random decision strategy")
    Z_j = [LocalDecisionStrategy(rng, diagram, d) for d in diagram.D]
    Z = DecisionStrategy(diagram.D, diagram.I_j[diagram.D], Z_j)

    @info "Testing CompatiblePaths"
    @test all(true for s in CompatiblePaths(diagram, Z))
    @test_throws DomainError CompatiblePaths(diagram, Z, Dict(diagram.D[1] => State(1)))
    node, state = (diagram.C[1], State(1))
    @test all(s[node] == state for s in CompatiblePaths(diagram, Z, Dict(node => state)))

    @info "Testing UtilityDistribution"
    U_distribution = UtilityDistribution(diagram, Z)

    @info "Testing StateProbabilities"
    S_probabilities = StateProbabilities(diagram, Z)

    @info "Testing conditional StateProbabilities"
    S_probabilities2 = StateProbabilities(diagram, Z, node, state, S_probabilities)

    @info "Testing printing functions"
    print_decision_strategy(diagram, Z, S_probabilities)
    print_decision_strategy(diagram, Z, S_probabilities, show_incompatible_states=true)
    print_utility_distribution(U_distribution)
    print_state_probabilities(diagram, S_probabilities, [diagram.Names[c] for c in diagram.C])
    print_state_probabilities(diagram, S_probabilities, [diagram.Names[d] for d in diagram.D])
    print_state_probabilities(diagram, S_probabilities2, [diagram.Names[c] for c in diagram.C])
    print_state_probabilities(diagram, S_probabilities2, [diagram.Names[d] for d in diagram.D])
    print_statistics(U_distribution)
    print_risk_measures(U_distribution, [0.0, 0.05, 0.1, 0.2, 1.0])

    @test true
end

@info "Testing model construction"
rng = MersenneTwister(4)
for (n_C, n_D, states, n_inactive, probability_scale_factor, probability_cut) in [
        (3, 2, [2, 3, 4], 0, 1.0, true),
        (3, 2, [2, 3], 0, -1.0, true),
        (3, 2, [3], 1, 100.0, true),
        (3, 2, [2, 3], 0, -1.0, false),
        (3, 2, [4], 1, 10.0, false)
    ]
    diagram = influence_diagram(rng, n_C, n_D, 2, 2, 2, states, n_inactive)
    test_decision_model(diagram, n_inactive, probability_scale_factor, probability_cut)
    test_analysis_and_printing(diagram)
end


text:
#===============================================================================

    The main object defined in this file is: SamplingUnit.
    It is the smallest `composite unit` (i.e. a unit consisting of at least
    two `GuidProp`/containers for paths etc.) and it can be used for sampling
    conditioned diffusions based on conditioning on multiple points.

===============================================================================#
"""
    struct SamplingUnit{TGP,TGPb,TW,TWn,TX}
        PP::Vector{TGP}
        PPb::Vector{TGPb}
        WW::Vector{TW}
        Wnr::TWn
        XX::Vector{TX}
    end

Smallest composite unit with containers needed for sampling of conditioned
diffusions via Guided Proposals.

# Fields
---
- `PP`: a vector of `GuidProp`
- `PPb`: a vector of `GuidProp` that can be used in a blocking schemes as
         `GuidProp`s on terminal subintervals
- `WW`: a vector of containers for sampled Wiener process
- `Wnr`: a flag for sampling Wiener processes
- `XX`: a vector of containers for a sampled process

    SamplingUnit(
        aux_laws, recording, tts, args=tuple();
        aux_laws_blocking=aux_laws, artificial_noise=1e-11,
        solver_choice_blocking=args
    )

Base constructor.

# Arguments
---
- `aux_laws`:
- `recording`:
- `tts`:
- `args`:
- `aux_laws_blocking`:
- `artificial_noise`:
- `solver_choice_blocking`:
"""
struct SamplingUnit{TGP,TGPb,TW,TWn,TX}
    PP::Vector{TGP}
    PPb::Vector{TGPb}
    WW::Vector{TW}
    Wnr::TWn
    XX::Vector{TX}

    function SamplingUnit(
            aux_laws, recording, tts, args=tuple();
            aux_laws_blocking=aux_laws, artificial_noise=1e-11,
            solver_choice_blocking=args
        )
        PP = build_guid_prop(aux_laws, recording, tts, args...)
        PPb = guid_prop_for_blocking(
            PP,
            aux_laws_blocking,
            artificial_noise,
            (args == tuple() ? solver_choice_blocking : args )...
        )
        XX, WW = trajectory(PP)
        Wnr = Wiener(PP)

        init_paths!(PP, WW, Wnr, XX, recording.x0_prior)
        new{eltype(PP),eltype(PPb),eltype(WW),typeof(Wnr),eltype(XX)}(
            PP, PPb, WW, Wnr, XX
        )
    end
end

"""
    init_paths!(P, WW, Wnr, XX, x0_prior)

Sample paths of guided proposals without using the preconditioned Crank–Nicolson
scheme.
"""
function init_paths!(PP, WW, Wnr, XX, x0_prior)
    while true
        forward_guide!(PP, XX, WW, rand(x0_prior); Wnr=Wnr) && return
    end
end

function all_eltypes(
        ::SamplingUnit{TGP,TGPb,TW,TWn,TX}
    ) where {TGP,TGPb,TW,TWn,TX}
    TGP,TGPb,TW,TWn,TX
end

"""
    GP.recompute_guiding_term!(u::SamplingUnit)

Recompute the guiding term assuming `u.PP` is a vector with guided proposal laws
"""
function GP.recompute_guiding_term!(u::SamplingUnit)
    recompute_guiding_term!(u.PP)
end

"""
    GP.loglikhd(u::SamplingUnit)

Return log-likelihood evaluated at a sampled path
"""
GP.loglikhd(u::SamplingUnit) = loglikhd(u.PP, u.XX)

"""
    draw_proposal_path!(u::SamplingUnit)

Sample a proposal path, compute log-likelihood along the way. Assumes
`u.XX[1].x[1]` is a starting point. No preconditioned Crank-Nicolson scheme is
used.
"""
function draw_proposal_path!(u::SamplingUnit)
    rand!(u.PP, u.XX, u.WW, _LL, u.XX[1].x[1]; Wnr=u.Wnr)
end


text:
# Latin hypercube sampling of values from given mins and maxs, per column.
function latin_hypercube_sampling(mins, maxs, numSamples)
  dims = length(mins)
  result = zeros(numSamples, dims)
  for(i in 1:dims)
    interval_len = (maxs[i] - mins[i]) / numSamples
    interval_starts = linspace(mins[i], maxs[i], numSamples+1)[1:numSamples]
    samples_for_dim = interval_starts + interval_len * rand(numSamples, 1)
    result[:,i] = samples_for_dim[shuffle([1:numSamples])]
  end
  result
end


text:
@doc doc"""
    latexalign()
Generate a ``LaTeX`` align environment from an input.

# Examples
## use with arrays

```julia
lhs = [:(dx/dt), :(dy/dt), :(dz/dt)]
rhs = [:(y-x), :(x*z-y), :(-z)]
latexalign(lhs, rhs)
```

```LaTeX
\begin{align}
\frac{dx}{dt} =& y - x \\\\
\frac{dy}{dt} =& x \cdot z - y \\\\
\frac{dz}{dt} =& - z \\\\
\end{align}
```

## use with ParameterizedFunction

```julia-repl
julia> using DifferentialEquations
julia> ode = @ode_def foldChangeDetection begin
    dm = r_m * (i - m)
    dy = r_y * (p_y * i/m - y)
end i r_m r_y p_y

julia> latexalign(ode)
```
```LaTeX
\begin{align}
\frac{dm}{dt} =& r_{m} \cdot \left( i - m \right) \\\\
\frac{dy}{dt} =& r_{y} \cdot \left( \frac{p_{y} \cdot i}{m} - y \right) \\\\
\end{align}
```

"""
function latexalign end

function latexalign(arr::AbstractMatrix; separator=" =& ", double_linebreak=false, starred=false, kwargs...)
    (rows, columns) = size(arr)
    eol = double_linebreak ? " \\\\\\\\\n" : " \\\\\n"
    arr = latexraw(arr; kwargs...)

    str = "\\begin{align$(starred ? "*" : "")}\n"
    for i in 1:(rows-1)
        str *= join(arr[i,:], separator) * eol
    end
    str *= join(arr[end,:], separator) * "\n"  # No EOL on the last line.
    str *= "\\end{align$(starred ? "*" : "")}\n"
    latexstr = LaTeXString(str)
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

function latexalign(lhs::AbstractArray, rhs::AbstractArray; kwargs...)
    return latexalign(hcat(lhs, rhs); kwargs...)
end

function latexalign(lhs::Tuple, rhs::Tuple; kwargs...)
    return latexalign(hcat(collect(lhs), collect(rhs)); kwargs...)
end

latexalign(args::Tuple...; kwargs...) = latexalign(hcat([collect(i) for i in args]...); kwargs...)

latexalign(arg::Tuple; kwargs...) = latexalign(hcat([collect(i) for i in arg]...); kwargs...)

function latexalign(nested::AbstractVector{AbstractVector}; kwargs...)
    return latexalign(hcat(nested...); kwargs...)
end

function latexalign(d::AbstractDict; kwargs...)
    latexalign(collect(keys(d)), collect(values(d)); kwargs...)
end

"""
    latexalign(vec::AbstractVector)

Go through the elements, split at any = sign, pass on as a matrix.
"""
function latexalign(vec::AbstractVector; kwargs...)
    lvec = latexraw(vec; kwargs...)
    ## turn the array into a matrix
    lmat = hcat(split.(lvec, " = ")...)
    ## turn the matrix ito arrays of left-hand-side, right-hand-side.
    larr = [lmat[i,:] for i in 1:size(lmat, 1)]
    length(larr) < 2 && error("Invalid intput to latexalign().")
    return latexalign( hcat(larr...) ; kwargs...)
end


text:
function test_include()
println("version 1")
end


text:
using JSON

# Setup GHA logger in CI
using Logging: global_logger
using GitHubActions: GitHubActionsLogger, set_output
get(ENV, "GITHUB_ACTIONS", "false") == "true" && global_logger(GitHubActionsLogger())

# Dir to store generated test results in
results_dir = mktempdir(; prefix="jl_test-runner_results_", cleanup=false)

# In CI, set results directory as output
set_output("results-path", results_dir)

for exercise_type in ("concept", "practice")
    for exercise in readdir(joinpath("exercises", exercise_type))
        # Allow only testing specified exercises
        if !isempty(ARGS) && !(exercise in ARGS)
            continue
        end

        exercise_path = joinpath("exercises", exercise_type, exercise)
        isdir(exercise_path) || continue

        # Concept exercises have an exemplar, practice exercises merely an example for CI purposes
        example_or_exemplar = exercise_type == "concept" ? "exemplar" : "example"

        @info "[$(uppercase(exercise_type))] Testing $exercise"

        # Determine file name of the solution
        meta_cfg = JSON.parsefile(joinpath(exercise_path, ".meta", "config.json"))
        solution_file = meta_cfg["files"]["solution"][1]

        # Copy solution to temporary directory
        tmp = mktempdir(; prefix="jl_$(exercise)_")
        cp(joinpath(exercise_path, ".meta", example_or_exemplar * ".jl"), joinpath(tmp, solution_file))
        cp(joinpath(exercise_path, "runtests.jl"), joinpath(tmp, "runtests.jl"))

        # Run test-runner
        run(`docker run
                --network none
                --read-only
                --mount type=bind,src=$tmp,dst=/solution
                --mount type=bind,src=$tmp,dst=/output
                --mount type=tmpfs,dst=/tmp
                exercism/julia-test-runner $exercise /solution /output`
        )

        cp(joinpath(tmp, "results.json"), joinpath(results_dir, "$(exercise_type)_$exercise.json"))
    end
end

max_stars_repo_name: stepinski/julia

max_issues_repo_name: stepinski/julia

max_forks_repo_name: stepinski/julia


text:
module Reactive

include("core.jl")
include("operators.jl")
include("async.jl")
include("time.jl")

end # module


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Random
Random.seed!(145)
x, y = 1:2:20, 5*rand(10)
yerr, xerr = 0.4*abs.(randn(10)), abs.(randn(10))
@gp x y yerr 1.5xerr "w boxerrorbars t 'test' lc 'black'"
@gp :- xrange= (0,22) yrange = (0,6)
save(term="pngcairo size 600,400", output="./code/plt_errs2.png") #hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using PDESolver
using AdvectionEquationMod

mesh, sbp, eqn, opts = run_advection(ARGS[1])


text:
function isdivisiblebyall(n, tab)
	for v in tab
		if n % v != 0
			return false
		end
	end
	return true
end

import Base.factorize
function factorize(n,F)
	f = zeros(eltype(n), F)
	f[1] = 1
	i = n
	while i > 1
		if isprime(i) && n % i == 0
			f[i] += 1
			n = div(n,i)
			i = n
		else
			i -= 1
		end
	end
	return f
end

function main()
	N = 20
	f = zeros(eltype(N), N, N)
	for i = N:-1:1
		f[:,i] = factorize(i, N)
	end
	println(f)
	resf = maximum(f, 2)
	println(resf')
	n = 1
	for i = 1:N
		for j = 1:resf[i]
			n *= i
		end
	end
	println(n)
	assert(isdivisiblebyall(n, [1:N]))
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

push!(LOAD_PATH, "../src/")

using Documenter, PiecewiseQuadratics

DocMeta.setdocmeta!(PiecewiseQuadratics, :DocTestSetup, :(using PiecewiseQuadratics); recursive=true)

makedocs(
    sitename = "PiecewiseQuadratics.jl",
    authors = "Nick Moehle, Ellis Brown, Mykel Kochenderfer",
    repo = "github.com/JuliaFirstOrder/PiecewiseQuadratics.jl.git",
    modules = [PiecewiseQuadratics],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "index.md",
        "api.md"
    ]
)

deploydocs(;
    repo = "github.com/JuliaFirstOrder/PiecewiseQuadratics.jl.git",
    devbranch = "main",
    push_preview = true,
)