text:
function search!(actr, ex)
    status = :searching
    # search until target found or termination threshold met
    while status == :searching
        # update decay values of objects in iconic memory
        update_decay!(actr)
        # update the finst values for return of inhabition
        update_finst!(actr)
        # update which items are in iconic memory based on what is currently visible
        update_visibility!(actr)
        # update the activations
        compute_activations!(actr)
        ex.visible ? update_window!(actr, ex) : nothing
        # perform search sequence
        status = _search!(actr, ex)
        ex.visible ? update_window!(actr, ex) : nothing
    end
    # add data for trial
    add_data!(ex)
    # add fixations for trial
    push!(ex.fixations, ex.trial_fixations)
    return nothing
end

function _search!(actr, ex)
    ex.trace ? print_trial(ex) : nothing
    data = ex.trial_data
    # production cycle
    cycle_time!(actr, ex)
    # select visual object with highest activation
    status = find_object!(actr, ex)
    if status == :error
        # if not objects are found, add respond time
        cycle_time!(actr, ex)
        motor_time!(actr, ex)
        # add fixation data
        add_no_fixation!(ex, actr)
        ex.trace ? print_response(actr, "absent") : nothing
        # add response to data
        add_response!(actr, data, :absent)
        return status
    end
    # Attending object in abstract-location
    # add cycle time
    cycle_time!(actr, ex)
    # add attend time
    attend_time!(actr, ex)
    # place object in visual buffer and update attend time
    attend_object!(actr, ex)
    # add fixation to fixation data
    add_fixation!(ex, actr)
    # cycle time before response
    cycle_time!(actr, ex)
    # check if current object matches target
    status = check_object(actr, ex)
    if status == :present
        # if matches, respond present and collect data
        cycle_time!(actr, ex)
        motor_time!(actr, ex)
        ex.trace ? print_response(actr, "present") : nothing
        add_response!(actr, data, status)
        return status
    end
    # if does not match, return "searching" and continue
    return status
end

function motor_time!(actr, ex)
    θ = gamma_parms(.065)
    tΔ = rand(Gamma(θ...))
    actr.time += tΔ
    ex.visible ? sleep(tΔ/ex.speed) : nothing
    return nothing
 end

function cycle_time!(actr, ex)
    tΔ = .05
    if actr.parms.rnd_time
        θ = gamma_parms(.05)
        tΔ = rand(Gamma(θ...))
    end
    actr.time += tΔ
    ex.visible ? sleep(tΔ/ex.speed) : nothing
    return nothing
end

function attend_time!(actr, ex)
    # tΔ = saccade_time(model) + visual_encoding(model)
    tΔ = visual_encoding(actr)
    actr.time += tΔ
    ex.visible ? sleep(tΔ/ex.speed) : nothing
    return nothing
end

saccade_time(actr) = saccade_time(actr, actr.visual_location.buffer[1])

function saccade_time(actr, vo)
    @unpack Δexe, β₀exe = actr.parms
    distance = compute_angular_distance(actr, vo)
    μ = β₀exe + distance*Δexe
    θ = gamma_parms(μ)
    return rand(Gamma(θ...))
end

function visual_encoding(actr)
    if actr.parms.rnd_time
        θ = gamma_parms(.085)
        return rand(Gamma(θ...))
    end
    return .085
end

gamma_parms(μ, σ) = (μ/σ)^2,σ^2/μ

gamma_parms(μ) = gamma_parms(μ, μ/3)

function add_fixation!(ex, actr)
    goal = get_buffer(actr, :goal)
    vision = get_buffer(actr, :visual)
    visicon = actr.visicon
    attend_time = vision[1].attend_time
    slots = goal[1].slots
    idx = findfirst(x->x.location == vision[1].location, visicon)
    fixation = Fixation(; attend_time, color=slots.color,
        shape=slots.shape, idx, stop=false)
    push!(ex.trial_fixations, fixation)
    return nothing
end

function add_no_fixation!(ex, actr)
    goal = get_buffer(actr, :goal)
    vision = get_buffer(actr, :visual)
    visicon = actr.visicon
    attend_time = actr.time
    slots = goal[1].slots
    idx = length(visicon) + 1
    fixation = Fixation(; attend_time, color=slots.color,
        shape=slots.shape, idx, stop=true)
    push!(ex.trial_fixations, fixation)
    return nothing
end

function add_data!(ex)
    push!(ex.data, ex.trial_data)
    return nothing
end

function add_response!(actr, data, status)
    data.response = status
    data.rt = actr.time
    set_trial_type!(data)
    return nothing
end

function relevant_object(actr, vo)
    !vo.visible || vo.attended ? (return false) : (return true)
end

function find_object!(actr, ex)
    v = :_
    if actr.parms.noise
        v = find_object2!(actr, ex)
    else
        v = find_object1!(actr, ex)
    end
    return v
end

function find_object1!(actr, ex)
    @unpack iconic_memory = actr.visual_location
    visible_objects = filter(x->relevant_object(actr, x), iconic_memory)
    if isempty(visible_objects)
        ex.trace ? print_abstract_location(actr, "error locating object") : nothing
        return :error
    end
    p = fixation_probs(actr, actr.visicon, visible_objects)
    idx = sample(1:length(p), Weights(p))
    if idx == length(p)
        ex.trace ? print_abstract_location(actr, "termination threshold exceeded") : nothing
        return :error
    end
    max_vo = visible_objects[idx]
    actr.visual_location.buffer = [max_vo]
    ex.trace ? print_abstract_location(actr, "object found") : nothing
    return :searching
end

function find_object2!(actr, ex)
    visible_objects = filter(x->relevant_object(actr, x), get_iconic_memory(actr))
    if isempty(visible_objects)
        ex.trace ? print_abstract_location(actr, "error locating object") : nothing
        return :error
    end
    max_vo = max_activation(visible_objects)
    if terminate(actr, max_vo)
        ex.trace ? print_abstract_location(actr, "termination threshold exceeded") : nothing
        return :error
    end
    actr.visual_location.buffer = max_vo
    ex.trace ? print_abstract_location(actr, "object found") : nothing
    return :searching
end

attend_object!(actr, ex) = attend_object!(actr, ex, actr.visual_location.buffer[1])

function attend_object!(actr, ex, vo)
    ex.trace ? print_vision(actr) : nothing
    actr.visual.buffer = [vo]
    vo.attended = true
    vo.attend_time = actr.time
    actr.visual.focus = vo.location
    ex.trace ? print_visual_buffer(actr) : nothing
    return nothing
end

check_object(actr, ex) = check_object(actr, actr.visual.buffer[1], ex, actr.goal.buffer[1])

function check_object(actr, vo, ex, target)
    for (f,v) in pairs(target.slots)
        if vo.features[f].value ≠ v
            ex.trace ? print_check(actr, "not found") : nothing
            update_threshold!(actr)
            return :searching
        end
    end
    ex.trace ? print_check(actr, "found") : nothing
    return :present
end

function max_activation(vos)
    max_vo = similar(vos, 1)
    mv = -Inf
    for vo in vos
        if vo.activation > mv
            max_vo[1] = vo
            mv = vo.activation
        end
    end
    return max_vo
end

function update_threshold!(actr)
    actr.parms.τₐ += actr.parms.Δτ
end

terminate(actr, max_vo) = terminate(actr, max_vo[1])

function terminate(actr, max_vo::VisualObject)
    τ = actr.parms.τₐ + add_noise(actr)
    τ > max_vo.activation ? (return true) : (return false)
end

function update_decay!(actr)
    map(x->update_decay!(actr, x), get_iconic_memory(actr))
end

function update_decay!(actr, object)
    for f in object.features
        if actr.parms.persistence < (actr.time - f.fixation_time)
            f.visible = false
            f.fixation_time = 0.0
        end
    end
    return nothing
end

function update_finst!(actr)
    attended_objects = filter(x->x.attended, get_iconic_memory(actr))
    isempty(attended_objects) ? (return) : nothing
    map(x->update_finst_span!(actr, x), attended_objects)
    update_n_finst!(actr, attended_objects)
end

function update_finst_span!(actr, vo)
    if actr.parms.finst_span < (actr.time - vo.attend_time)
        vo.attended = false
        vo.attend_time = 0.0
    end
    return nothing
end

function update_n_finst!(actr, vos)
    N = length(vos) - actr.parms.n_finst
    N <= 0 ? (return) : nothing
    sort!(vos, by=x->x.attend_time)
    for i in 1:N
        vos[i].attended = false
        vos[i].attend_time = 0.0
    end
    return nothing
end

function compute_activations!(actr)
    @unpack iconic_memory = actr.visual_location
    iconic_memory = filter(x->x.visible, actr.visicon)
    topdown_activations!(iconic_memory, actr.goal.buffer[1])
    bottomup_activations!(iconic_memory)
    weighted_activations!(actr, iconic_memory)
    actr.parms.noise ? add_noise!(actr, iconic_memory) : nothing
    return nothing
end

function weighted_activations!(model, iconic_memory)
    for object in iconic_memory
        weighted_activation!(model, object)
    end
    return nothing
end

function weighted_activation!(actr, vo)
    vo.activation = weighted_activation(actr, vo)
    return nothing
end

function weighted_activation(actr, vo)
    return actr.parms.bottomup_weight * vo.bottomup_activation +
    actr.parms.topdown_weight * vo.topdown_activation
end

bottomup_activations!(actr::ACTRV) = bottomup_activations!(get_iconic_memory(actr))

function bottomup_activations!(iconic_memory)
    for vo1 in iconic_memory
        activation = 0.0
        for vo2 in iconic_memory
            vo1 == vo2 ? continue : nothing
            activation += bottomup_activation(vo1, vo2)
        end
        vo1.bottomup_activation = activation
    end
    return nothing
end

function bottomup_activation(vo1, vo2)
    distance = compute_distance(vo1, vo2)
    return bottomup_activation(vo1.features, vo2.features, distance)
end

function bottomup_activation(f1, f2, distance)
    activation = 0.0
    for (f,v) in pairs(f1)
        if (f2[f].value ≠ v.value) && v.visible && f2[f].visible
            activation += 1.0/sqrt(distance)
        end
    end
    return activation
end

topdown_activations!(iconic_memory, chunk::Chunk) = topdown_activations!(iconic_memory, chunk.slots)

function topdown_activations!(iconic_memory, target)
    for object in iconic_memory
        topdown_activation!(object, target)
    end
    return nothing
end

function topdown_activation!(object, target)
    activation = 0.0
    for (f,v) in pairs(object.features)
        if v.visible
            if v.value == target[f]
                activation += 1.0
            end
        else
            activation += .5
        end
    end
    object.topdown_activation = activation
    return nothing
end

add_noise(actr) = rand(Normal(0, actr.parms.σ))

add_noise!(actr, vo::VisualObject) = vo.activation += add_noise(actr)

add_noise!(actr, iconic_memory) = map(x->add_noise!(actr, x), iconic_memory)

function update_visibility!(actr)
    map(x->feature_visibility!(actr, x), get_iconic_memory(actr))
    map(x->object_visibility!(x), get_iconic_memory(actr))
end

function feature_visibility!(actr, vo)
    angular_distance = compute_angular_distance(actr, vo)
    for (f,v) in pairs(vo.features)
        parms = actr.parms.acuity[f]
        threshold = compute_acuity_threshold(parms, angular_distance)
        if feature_is_visible(vo, threshold)
            v.visible = true
        end
    end
    return nothing
end

feature_is_visible(vo::VisualObject, threshold) = feature_is_visible(vo.angular_size, threshold)

function feature_is_visible(angular_size, threshold)
    return angular_size > threshold
end

function object_visibility!(object)
    for f in object.features
        if f.visible
            object.visible = true
            return nothing
        end
    end
    object.visible = false
    return nothing
end

function compute_distance(actr::ACTRV, object)
    sqrt(sum((actr.visual.focus .- object.location).^2))
end

function compute_distance(vo1, vo2)
    sqrt(sum((vo1.location .- vo2.location).^2))
end

"""
Angular distance in degress. Also known as eccentricity
"""
function compute_angular_distance(actr, vo)
    distance = compute_distance(actr, vo)
    return pixels_to_degrees(actr, distance)
end

function compute_angular_size!(actr)
    map(x->compute_angular_size!(actr, x), get_visicon(actr))
    return nothing
end

function compute_angular_size!(actr, vo)
    ppi = 72 # pixels per inch
    distance = actr.parms.viewing_distance * ppi
    vo.angular_size = compute_angular_size(distance, vo.width)
end

function compute_angular_size(distance, width)
    radians = 2 * atan(width / (2 * distance))
    return rad2deg(radians)
end

rad2deg(radians) = radians * 180 / pi

function compute_acuity_threshold(parms, angular_distance)
    return parms.a * angular_distance^2 - parms.b * angular_distance
end

function pixels_to_degrees(actr, pixels)
    ppi = 72 # pixels per inch
    radians = atan(pixels/ppi, actr.parms.viewing_distance)
    return rad2deg(radians)
end

function orient!(actr, ex)
    w = ex.array_width / 2
    actr.visual.focus = fill(w, 2)
end

function fixation_prob(actr, visicon, visible_objects, fixation)
    @unpack stop,idx = fixation
    act = map(x->x.activation, visible_objects)
    push!(act, actr.parms.τₐ)
    σ = actr.parms.σ*sqrt(2)
    vo_act = stop ? act[end] : visicon[idx].activation
    return exp(vo_act/σ)/sum(exp.(act/σ))
end

function fixation_probs(actr, visicon, visible_objects)
    act = map(x->x.activation, visible_objects)
    push!(act, actr.parms.τₐ)
    σ = actr.parms.σ * sqrt(2)
    return exp.(act / σ) / sum(exp.(act / σ))
end

# visual_encoding(model) = visual_encoding(model, model.abstract_location[1])
#
# function visual_encoding(model, vo)
#     frequency = model.init_freq
#     distance = compute_angular_distance(model, vo)
#     μ = model.K_encode*-log(frequency)*exp(distance*model.κ_encode)
#     θ = gamma_parms(μ)
#     # vo.encode_start = model.current_time
#     # vo.encode_time = rand(Gamma(a, b))
#     return rand(Gamma(θ...))
# end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

mutable struct IntervalParameters

    precision_type::Type
    precision::Int
    rounding::Symbol
    pi::Interval{BigFloat}

    IntervalParameters() = new(BigFloat, 256, :narrow)  # leave out pi
end

const parameters = IntervalParameters()

## Precision:

"`big53` creates an equivalent `BigFloat` interval to a given `Float64` interval."
function big53(a::Interval{Float64})
    setprecision(Interval, 53) do  # precision of Float64
        atomic(Interval{BigFloat}, a)
    end
end

function big53(x::Float64)
    # BigFloat(x, 53)  # in Julia v0.6

    setprecision(53) do
        BigFloat(x)
    end
end

setprecision(::Type{Interval}, ::Type{Float64}) = parameters.precision_type = Float64
# does not change the BigFloat precision

function setprecision(::Type{Interval}, ::Type{T}, prec::Integer) where T<:AbstractFloat
    #println("SETTING BIGFLOAT PRECISION TO $precision")
    setprecision(BigFloat, prec)

    parameters.precision_type = T
    parameters.precision = prec
    parameters.pi = atomic(Interval{BigFloat}, pi)

    prec
end

setprecision(::Type{Interval{T}}, prec) where T<:AbstractFloat =
    setprecision(Interval, T, prec)

setprecision(::Type{Interval}, prec::Integer) = setprecision(Interval, BigFloat, prec)

function setprecision(f::Function, ::Type{Interval}, prec::Integer)

    old_precision = precision(Interval)
    setprecision(Interval, prec)

    try
        return f()
    finally
        setprecision(Interval, old_precision)
    end
end

# setprecision(::Type{Interval}, precision) = setprecision(Interval, precision)
setprecision(::Type{Interval}, t::Tuple) = setprecision(Interval, t...)

precision(::Type{Interval}) = (parameters.precision_type, parameters.precision)

function Base.setrounding(f::Function, ::Type{Rational{T}},
    rounding_mode::RoundingMode) where T
    setrounding(f, float(Rational{T}), rounding_mode)
end

float(x::Interval{T}) where T = atomic( Interval{float(T)}, x)  # https://github.com/dpsanders/IntervalArithmetic.jl/issues/174

big(x::Interval) = atomic(Interval{BigFloat}, x)

ext: jl


text:
function _short_desc(s::String)
    return length(s) < 80 ? s : strip(s[1:79]*"…")
end

function Base.show(io::IO, dataset::MangalDataset)
    print(io, """🗃  Mangal dataset #$(dataset.id) ($(dataset.name))""")
end

function Base.show(io::IO, network::MangalNetwork)
    print(io, """🕸  Mangal network #$(network.id) ($(network.name))""")
end

function Base.show(io::IO, interaction::MangalInteraction)
    print(io, """🔄  Mangal int° #$(interaction.id): from $(interaction.from.name) to $(interaction.to.name)""")
end

function Base.show(io::IO, user::MangalUser)
    print(io, """🫂  Mangal contributor #$(user.id): $(user.name) - $(user.orcid)""")
end

function Base.show(io::IO, node::MangalNode)
    print(io, """🦝  Mangal node #$(node.id): $(node.name)""")
end

function Base.show(io::IO, taxon::MangalReferenceTaxon)
    print(io, """🌲  Mangal taxon #$(taxon.id): $(taxon.name)""")
end


text:
using BinaryProvider, Libdl

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libmbedcrypto"], :libmbedcrypto),
    LibraryProduct(prefix, String["libmbedx509"], :libmbedx509),
    LibraryProduct(prefix, String["libmbedtls"], :libmbedtls),
]

const juliaprefix = joinpath(Sys.BINDIR, "..")

juliaproducts = Product[
    LibraryProduct(juliaprefix, "libmbedcrypto", :libmbedcrypto),
    LibraryProduct(juliaprefix, "libmbedx509", :libmbedx509),
    LibraryProduct(juliaprefix, "libmbedtls", :libmbedtls)
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaWeb/MbedTLSBuilder/releases/download/v0.17.0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc) => ("$bin_prefix/MbedTLS.v2.16.0.aarch64-linux-gnu.tar.gz", "aeaaabcef59f4587bbd003a3cb891d85e9d723cba641fd004bc3273bf574dad5"),
    Linux(:aarch64, libc=:musl) => ("$bin_prefix/MbedTLS.v2.16.0.aarch64-linux-musl.tar.gz", "dcc4c21affc2e687ae43bac42ff9481a70e9824fc08518987b805b9e318f0068"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf) => ("$bin_prefix/MbedTLS.v2.16.0.arm-linux-gnueabihf.tar.gz", "48cd2b4cbd9b2f4a7ab2f1060a240b1170b5807af236777602cd331c6b29c037"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf) => ("$bin_prefix/MbedTLS.v2.16.0.arm-linux-musleabihf.tar.gz", "25003edb6fb346e4cd299c30b7660f6abb8e493df75fa224951d252cee4b5e28"),
    Linux(:i686, libc=:glibc) => ("$bin_prefix/MbedTLS.v2.16.0.i686-linux-gnu.tar.gz", "113d3e8220ea24d2e075623cebcb2e0a1864bada52d5d9065cddc7b04ed24f1c"),
    Linux(:i686, libc=:musl) => ("$bin_prefix/MbedTLS.v2.16.0.i686-linux-musl.tar.gz", "48531771fcfb69cbc79abde1c06ac74e5e5d9e24367a111eac6cbfe42e5377c8"),
    Windows(:i686) => ("$bin_prefix/MbedTLS.v2.16.0.i686-w64-mingw32.tar.gz", "5be5cae452e513780380f276a19f2cd5eecbecad01728fd180b1c19022369e96"),
    Linux(:powerpc64le, libc=:glibc) => ("$bin_prefix/MbedTLS.v2.16.0.powerpc64le-linux-gnu.tar.gz", "d07a85c710a82bd7ab0d9eabb469e378ce31083161ff864f157962ec018db4e2"),
    MacOS(:x86_64) => ("$bin_prefix/MbedTLS.v2.16.0.x86_64-apple-darwin14.tar.gz", "8939d8b39618d9206feac886298070f5519ab13ac190b44d99da5e62b3d5fef1"),
    Linux(:x86_64, libc=:glibc) => ("$bin_prefix/MbedTLS.v2.16.0.x86_64-linux-gnu.tar.gz", "f8985d86a78e03c7b7cd169c664b93ebe9eb726c652be7bb96707810b0c68bbf"),
    Linux(:x86_64, libc=:musl) => ("$bin_prefix/MbedTLS.v2.16.0.x86_64-linux-musl.tar.gz", "e9a254ff953899f79b79af8f30f9a95f0acc21f44d5a3cf683cfc19b1fb21808"),
    FreeBSD(:x86_64) => ("$bin_prefix/MbedTLS.v2.16.0.x86_64-unknown-freebsd11.1.tar.gz", "8c923544323c6009a46a70cedf7fb714c5d9e0f23650c2ef8b151a84bf74f527"),
    Windows(:x86_64) => ("$bin_prefix/MbedTLS.v2.16.0.x86_64-w64-mingw32.tar.gz", "85b95b43014fc8370c3dd8b41f341c6849fcde12e321c24faa8a86193161252e"),
)

# First, check to see if we're all satisfied
gpl = haskey(ENV, "USE_GPL_MBEDTLS")
forcebuild = parse(Bool, get(ENV, "FORCE_BUILD", "false")) || gpl
done = false
if any(!satisfied(p; verbose=verbose) for p in products) || forcebuild
    if haskey(download_info, platform_key()) && !forcebuild
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
        done = all(satisfied(p; verbose=verbose) for p in products)
        done && @info "using prebuilt binaries"
    end
    if !done && all(satisfied(p; verbose=verbose) for p in juliaproducts) && !forcebuild
        @info "using julia-shippied binaries"
        products = juliaproducts
    elseif !done || forcebuild
        @info "attempting source build"
        VERSION = "2.16.0"
        url, hash = haskey(ENV, "USE_GPL_MBEDTLS") ?
            ("https://tls.mbed.org/download/mbedtls-$VERSION-gpl.tgz", "c8d7a4696287cb9533c455bdd65859106dbdd4472c125194387486e6d4df2799") :
            ("https://tls.mbed.org/download/mbedtls-$VERSION-apache.tgz", "e3dab56e9093c790b7d5e0f7eb19451010fe680649d25cf1dcca9d5441669ae2")
        download_verify(url, hash, joinpath(@__DIR__, "mbedtls.tgz"), force=true, verbose=true)
        unpack(joinpath(@__DIR__, "mbedtls.tgz"), @__DIR__; verbose=true)
        withenv("VERSION"=>VERSION) do
            run(Cmd(`./build.sh`, dir=@__DIR__))
        end
        if any(!satisfied(p; verbose=verbose) for p in products)
            error("attempted to build mbedtls shared libraries, but they couldn't be located (deps/usr/lib)")
        end
    end
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
module HMM
importall BrmlFunctions
#@reexport using BrmlFunctions, HMM, GraphPlot, Winston ## CHECK!!!

export HMMforward
include("HMMforward.jl")

export HMMbackward
include("HMMbackward.jl")

export HMMgamma
include("HMMgamma.jl")

export HMMsmooth
include("HMMsmooth.jl")

export HMMviterbi
include("HMMviterbi.jl")

export HMMem
include("HMMem.jl")

import PyPlot
export PyPlot

end


text:
# import Pkg; Pkg.add("PortAudio")
using Faust, PortAudio

dsp = compile("""
import("stdfaust.lib");

notes = 48, 55, 60, 62, 64, 67, 71;
n = ba.count(notes);

bank = par(i, n, fi.resonbp(ba.take(i+1, notes) : ba.midikey2hz, 1000, 0.007) : sp.spat(2, theta(i), d(i)))
with {
    theta(i) = (i / n) * 2 * ma.PI;
    d(i) = i / n;
};

process = sp.stereoize(_ <: bank :> _);
""")

devices = PortAudio.devices()
dev = filter(x -> x.maxinchans == 2 && x.maxoutchans == 2, devices)[1]

PortAudioStream(dev, dev) do stream
    block_size = 1024
    init!(dsp, block_size=block_size, samplerate=Int(stream.samplerate))
    while true
        dsp.inputs = convert(Matrix{Float32}, read(stream, block_size))
        write(stream, compute!(dsp))
    end
end


text:
## Exercise 6-5
## The Ackermann function, A(m,n), is defined:

##           ⎧ n + 1  if m = 0
## A(m, n) = ⎨ A(m-1, 1) if m>0 &n = 0
##           ⎩ A(m-1, A(m, n-1)) if m>0 and n>0.

## See https://en.wikipedia.org/wiki/Ackermann_function. Write a function named ack that evaluates the Ackermann function. Use your function to evaluate ack(3, 4), which should be 125. What happens for larger values of m and n?
println("Ans: ")

function ack(m, n)
    if m == 0
        return n + 1
    elseif m > 0 && n == 0
        return ack(m -1, 1)
    elseif m > 0 && n > 0
        return ack(m - 1, ack(m, n -1 ))
    end
end

println(ack(3, 4))

println("End.")


text:
bitstring(x::BFloat16) = bitstring(reinterpret(uinttype(BFloat16), x))

show(io::IO, x::BFloat16) = print("BFloat16($(string(Float32(x))))")
print(io::IO, x::BFloat16) = print(string(Float32(x)))


text:
"This module contains all inner parts."
module Internal

export Problem, Options, solve

using Roots
using UnPack

"""
```
Problem(
    α::Function,
    β::Function,
    γ::Function,
    δ::Function,
    p::Function
    q::Function
    a::Float64,
    b::Float64,
    A::Float64,
    y₀::Float64,
    z₀::Float64,
)
```

A Cauchy problem for a linear-stiff system of two ordinary differentials equations.

Expects the problem in the form:

```
[L]  y'(x) = α(A) y(x) + β(A) z(x) + p(x,z(x))
[NL] z'(x) = γ(A) y(x) + δ(A) z(x) + q(x,y(x),z(x))
y(a) = y₀, z(a) = z₀
```

where [L] is a linear equation, and [NL] is a nonlinear equation. The interval is [a, b].

You can call `y(x)` in `q` as `yᵢ₊₁(x, zᵢ₊₁)`.

"""
struct Problem
    α::Function
    β::Function
    γ::Function
    δ::Function

    p::Function
    q::Function

    a::Float64
    b::Float64
    A::Float64

    y₀::Float64
    z₀::Float64
end

"""
`Options(n::Integer)`

Options to provide to the solver. Includes the number of nodes in the grid `n`.
"""
struct Options
    n::Integer
end

"The default set of options. Sets `n` to 1000."
DefaultOptions = Options(1000)

include("Functions/solve.jl")

end


text:
Mcdf(f,fmin,fmax) = (1.0./f - 1.0/fmax) ./ (1.0/fmin - 1.0/fmax)

function selection(λ, f, tend, t1)
    #define the equation for selection as above
    s = (λ .* t1 + log.(f ./ (1 .- f))) ./ (λ .* (tend - t1))
    return s
end

function selection2clone(λ, f1, f2, tend, t1, t2)
    #define the equation for selection as above

    s1 = zeros(Float64, length(f1))
    s2 = zeros(Float64, length(f1))

    for i in 1:length(f1)
      if (f2[i] + f1[i]) < 1.0
        s1[i] = (λ .* t1[i] + log.(f1[i] ./ (1 - f1[i] - f2[i]))) ./ (λ .* (tend[i] - t1[i]))
        s2[i] = (λ .* t2[i] + log.(f2[i] ./ (1 - f1[i] - f2[i]))) ./ (λ .* (tend[i] - t2[i]))
      else
        s1[i] = (λ .* t1[i] + log.((f1[i] - f2[i]) ./ (1 - f1[i]))) ./ (λ .* (tend[i] - t1[i]))
        s2[i] = (λ .* t2[i] + log.(f2[i] ./ (1 - f1[i]))) ./ (λ .* (tend[i] - t2[i]))
      end
    end

    return s1, s2
end

function collectoutput1clone(abcres; Nmax = 10^10)

    scmuts = map(x -> x.other[2], abcres.particles)
    scdivs = map(x -> x.other[3], abcres.particles)
    scfreq = map(x -> x.other[4], abcres.particles)
    mu = abcres.parameters[:, 1]
    weights = abcres.weights

    # eulergamma/log(2) is stochastic correction see Durrett Branching Process Models of Cancer, needed for selection calculation
    t1 = ((shuffle(scmuts) ./ mu) / (2 * log.(2))) .- MathConstants.eulergamma/log.(2) #shuffle so that t1 posterior is not correlated
    tend = (log.(Nmax .* (1 .- scfreq)) / log.(2))
    s = selection(log.(2), scfreq, tend, t1)

    DF = DataFrame(mu = mu,
    clonalmutations = abcres.parameters[:, 2],
    s = s,
    t = t1,
    cellularity = abcres.parameters[:, 5],
    frequency = scfreq,
    scmuts = map(x -> Float64(x), scmuts),
    weight = weights)

    return DF
end

function swapvalues(x1, x2, indeces)

  swap1 = x1[indeces]
  swap2 = x2[indeces]
  x1[indeces] = swap2
  x2[indeces] = swap1

  return x1, x2
end

function clonesize(scfreq1, scfreq2)
  #assume that if scfreq1 + scfreq2 > 1 then clones are nested
  freqfactor = zeros(Float64, length(scfreq1))
  for i in 1:length(scfreq1)
    if (scfreq1[i] + scfreq2[i]) .> 1.0
      freqfactor[i] = 1 - scfreq1[i]
      freqfactor[i] = maximum([0.001, freqfactor[i]])
    else
      freqfactor[i] = 1 - scfreq1[i] - scfreq2[i]
      freqfactor[i] = maximum([0.001, freqfactor[i]])
    end
  end

  return freqfactor
end

function collectoutput2clone(abcres; Nmax = 10^10)

    scmuts1 = map(x -> x.other[2], abcres.particles)
    scmuts2 = map(x -> x.other[3], abcres.particles)
    scdivs1 = map(x -> x.other[4], abcres.particles)
    scdivs2 = map(x -> x.other[5], abcres.particles)
    scfreq1 = map(x -> x.other[6], abcres.particles)
    scfreq2 = map(x -> x.other[7], abcres.particles)
    mu = abcres.parameters[:, 1]

    s1 = abcres.parameters[:, 3]
    t1 = abcres.parameters[:, 4]
    s2 = abcres.parameters[:, 5]
    t2 = abcres.parameters[:, 6]
    weights = abcres.weights

    scmuts1 = map(x -> Float64(x), scmuts1)
    scmuts2 = map(x -> Float64(x), scmuts2)

    indeces = .!(scfreq1 .> scfreq2)

    #identify subclone 1 as the highest frequency subclone
    scfreq1, scfreq2 = swapvalues(scfreq1, scfreq2, indeces)
    scmuts1, scmuts2 = swapvalues(scmuts1, scmuts2, indeces)
    scdivs1, scdivs2 = swapvalues(scdivs1, scdivs2, indeces)
    s1, s2 = swapvalues(s1, s2, indeces)
    t1, t2 = swapvalues(t1, t2, indeces)

    freqfactor = clonesize(scfreq1, scfreq2)

    t1a = ((shuffle(scmuts1) ./ mu) / (2 * log.(2))) .- MathConstants.eulergamma/log.(2)
    t1b = ((shuffle(scmuts2) ./ mu) / (2 * log.(2))) .- MathConstants.eulergamma/log.(2)
    tend = (log.(Nmax .* (freqfactor)) / log.(2))
    s1, s2 = selection2clone(log.(2), scfreq1, scfreq2, tend, t1a, t1b)

    DF = DataFrame(mu = abcres.parameters[:, 1],
    clonalmutations = abcres.parameters[:, 2],
    s1 = s1,
    t1 = t1a,
    s2 = s2,
    t2 = t1b,
    cellularity = abcres.parameters[:, 7],
    frequency1 = scfreq1,
    frequency2 = scfreq2,
    scmuts1 = scmuts1,
    scmuts2 = scmuts2,
    weight = weights)

    return DF
end

function collectoutput0clone(abcres)

  mupost = abcres.parameters[:, 1]
  cmpost = abcres.parameters[:, 2]
  cellularity = abcres.parameters[:, 3]
  weights = abcres.weights

  DFpost = DataFrame(mu = mupost,
  clonalmutations = cmpost,
  cellularity = cellularity,
  weight = weights)
end

function cumulativedistribution(VAF; fmin = 0.1, fmax = 0.3)

    #calculate cumulative sum
    steps = fmax:-0.001:fmin
    cumsum = Array{Int64}(0)
    v = Array{Float64}(0)

    for i in steps
        push!(cumsum, sum(VAF .>= i))
        push!(v, i)
    end
    cumsum = cumsum - cumsum[1]

    DF = DataFrame(cumsum = map(Float64, cumsum), v = v)
    DF[:invf] = 1 ./ DF[:v] - 1 ./ fmax

    DF[:theory] = Mcdf(DF[:v], fmin, fmax)
    DF[:normalized] = DF[:cumsum] ./ maximum(DF[:cumsum])

    lmfit = fit(LinearModel, @formula(cumsum ~ invf + 0), DF)
    DF[:prediction] = predict(lmfit)

    return CancerSeqSim.AnalysedData(DF, VAF)
end

function averagehistogram(particles, model, VAF)

    x = 0.005:0.01:1.005
    y = fit(Histogram, VAF, x, closed=:right)
    DFhist = DataFrame(VAF = x[1:end-1], freq = y.weights)

    particles = particles[map(x -> x.model, particles).==model]
    wts = map(x -> x.weight, particles)
    N = length(particles)

    M = zeros(Int64, 100, N)
    i = 1

    for j in 1:N
        M[:, i] = convert(Array, particles[i].other[1][:freq])
        i = i + 1
    end

    mvalues = Float64[]
    meanvalues = Float64[]
    lquant = Float64[]
    lquart = Float64[]
    uquant = Float64[]
    uquart = Float64[]

    for i in 1:size(M, 1)
      push!(mvalues, median(vec(collect(M[i, :]')), weights(wts)))
      push!(meanvalues, mean(vec(collect(M[i, :]')), weights(wts)))
      push!(lquant, quantile(vec(collect(M[i, :]')), weights(wts), 0.025))
      push!(uquant, quantile(vec(collect(M[i, :]')), weights(wts), 0.975))
      push!(lquart, quantile(vec(collect(M[i, :]')), weights(wts), 0.25))
      push!(uquart, quantile(vec(collect(M[i, :]')), weights(wts), 0.75))
    end

    DFr = DataFrame(median = mvalues,
                    mean = meanvalues,
                    lowerq95 = lquant,
                    upperq95 = uquant,
                    lowerquartile = lquart,
                    upperquartile = uquart,
                    VAF = DFhist[:VAF],
                    truecounts = DFhist[:freq])

    return DFr
end

function saveresults(res::Results; resultsdirectory = "output")
  makedirectories(joinpath(resultsdirectory, res.SampleName, "finalpopulation"))
  getresults(res.ABCresults, resultsdirectory, res.SampleName, res.VAF, save = true)
  return
end

function getresults(abcres, resultsdirectory, sname, VAF; save = false, Nmaxinf = 10^10, savepopulations = false, popnum = 1)

  if savepopulations
      resultsdirectory = joinpath(resultsdirectory, sname, "populations", "population_$popnum")
  else
      resultsdirectory = joinpath(resultsdirectory, sname, "finalpopulation")
  end
  posteriors = Posterior[]
  #save model posterior
  DFmp = DataFrame(Model = map(x -> string(x),0:length(abcres.modelprob) - 1), Probability = abcres.modelprob)

  if abcres.modelprob[1] > 0.0
    DFpost0 = collectoutput0clone(getmodel(abcres, 1))
    DFr = averagehistogram(abcres.particles, 1, VAF)
    if save == true
      CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-parameters-clone0.csv"), DFpost0)
      CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-histogram-clone0.csv"), DFr)
    end
    push!(posteriors, Posterior(DFr, DFpost0, abcres.modelprob[1]))
  else
    push!(posteriors, Posterior(DataFrame(), DataFrame(), abcres.modelprob[1]))
  end

  if abcres.modelprob[2] > 0.0
    DFpost1 = collectoutput1clone(getmodel(abcres, 2), Nmax = Nmaxinf)
    DFr = averagehistogram(abcres.particles, 2, VAF)
    if save == true
      CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-histogram-clone1.csv"), DFr)
      CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-parameters-clone1.csv"), DFpost1)
    end
    push!(posteriors, Posterior(DFr, DFpost1, abcres.modelprob[2]))
  else
    push!(posteriors, Posterior(DataFrame(), DataFrame(), abcres.modelprob[2]))
  end

  if (length(abcres.modelprob) > 2) && (abcres.modelprob[3] > 0.0)
    DFpost2 = collectoutput2clone(getmodel(abcres, 3), Nmax = Nmaxinf)
    DFr = averagehistogram(abcres.particles, 3, VAF)
    if save == true
      CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-parameters-clone2.csv"), DFpost2)
      CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-histogram-clone2.csv"), DFr)
    end
    push!(posteriors, Posterior(DFr, DFpost2, abcres.modelprob[3]))
  else
    if (length(abcres.modelprob) > 2)
      push!(posteriors, Posterior(DataFrame(), DataFrame(), abcres.modelprob[3]))
    end
  end

  DF = DataFrame(Model = map(x -> string(x),0:length(abcres.modelprob) - 1), Probability = abcres.modelprob)
  if save == true
    CSV.write(joinpath(resultsdirectory, "posterior", "$(sname)-modelprobabilities.csv"), DF)
  end

  return posteriors, DFmp
end

function makedirectories(resultsdirectory)

  if isdir(resultsdirectory) == false
    mkpath(resultsdirectory)
  end

  if isdir(joinpath(resultsdirectory, "plots")) == false
    mkpath(joinpath(resultsdirectory, "plots"))
  end
  if isdir(joinpath(resultsdirectory, "processed")) == false
    mkpath(joinpath(resultsdirectory, "processed"))
  end

  if isdir(joinpath(resultsdirectory, "posterior")) == false
    mkpath(joinpath(resultsdirectory, "posterior"))
  end

  if isdir(joinpath(resultsdirectory, "data")) == false
    mkpath(joinpath(resultsdirectory, "data"))
  end

end

function makeplotsdirectories(resultsdirectory)
  if isdir(joinpath(resultsdirectory)) == false
    mkpath(joinpath(resultsdirectory))
  end
  if isdir(joinpath(resultsdirectory, "plots")) == false
    mkpath(joinpath(resultsdirectory, "plots"))
  end
end

function makedirectory(resultsdirectory)
  if isdir(joinpath(resultsdirectory)) == false
    mkpath(joinpath(resultsdirectory))
  end
end

function show(io::IO, res::Results)
  show(res.ABCresults)
end

function show(io::IO, ABCresults::ApproxBayes.ABCSMCmodelresults)

  posteriors, DFmp = getresults(ABCresults, "nothing", "nothing", [1.0, 2.0]; save = false, Nmaxinf = ABCresults.setup.other)

  @printf("Total number of simulations: %.2e\n", sum(ABCresults.numsims))
  println("Cumulative number of simulations = $(cumsum(ABCresults.numsims))")
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)
  println("Tolerance schedule = $(round.(ABCresults.ϵ, digits = 2))\n")

  print("Model probabilities:\n")
  for j in 1:length(ABCresults.modelprob)
    @printf("\tModel %d (%d subclones): %.3f\n", j, j-1, ABCresults.modelprob[j])
  end

  print("\nParameters:\n\n")

  for j in 1:length(ABCresults.parameters)
    if ABCresults.modelprob[j] > 0.0
      print("Model $j ($(j-1) subclones)\n")

      upperci = zeros(Float64, size(posteriors[j].Parameters, 2) - 1)
      lowerci = zeros(Float64, size(posteriors[j].Parameters, 2) - 1)
      parametermeans = zeros(Float64, size(posteriors[j].Parameters, 2) - 1)
      parametermedians = zeros(Float64, size(posteriors[j].Parameters, 2) - 1)

      for i in 1:(size(posteriors[j].Parameters, 2) - 1)
        parametermeans[i] = mean(posteriors[j].Parameters[:, i],
        weights(posteriors[j].Parameters[:weight]))
        parametermedians[i] = median(posteriors[j].Parameters[:, i],
        weights(posteriors[j].Parameters[:weight]))
        (lowerci[i], upperci[i]) = quantile(posteriors[j].Parameters[:, i],
        weights(posteriors[j].Parameters[:weight]), [0.025,0.975])
      end

      print("\tMedian (95% intervals):\n")
      for i in 1:length(parametermeans)
          if j == 1
            parameternames = ["μ/β", "Clonal Mutations", "Cellularity"]
          elseif j == 2
            parameternames = ["μ/β", "Clonal Mutations", "Fitness", "Time (tumour doublings)", "Cellularity", "Subclone Frequency", "Subclone Mutations"]
          elseif j == 3
            parameternames = ["μ/β", "Clonal Mutations", "Fitness - Subclone 1", "Time (tumour doublings) - Subclone 1", "Fitness - Subclone 2", "Time (tumour doublings) - Subclone 2", "Cellularity", "Subclone 1 Frequency", "Subclone 2 Frequency", "Subclone 1 Mutations", "Subclone 2 Mutations"]
          end
          @printf("\tParameter %d - %s: %.2f (%.2f,%.2f)\n", i, parameternames[i], parametermedians[i], lowerci[i], upperci[i])
      end
    end
  end
end


text:
using Highlights.Tokens, Highlights.Themes

import MyterialColors:
    green,
    green_light,
    purple_light,
    blue,
    red,
    green_light,
    green,
    blue_light,
    orange_light,
    pink_light,
    orange,
    yellow,
    yellow_light,
    white,
    green,
    amber,
    green_dark,
    orange

import Parameters: @with_kw

"""
    Theme

Stores colors for different semantically relevant items, used to
style outputs to terminal.
"""
@with_kw mutable struct Theme
    name::String = "default"
    docstring::String = green_dark
    string::String = "#64b565"
    type::String = purple_light
    emphasis::String = "$blue  bold"
    emphasis_light::String = yellow_light
    code::String = "$(yellow) italic"
    multiline_code::String = "$(yellow) italic"
    symbol::String = orange
    expression::String = amber
    number::String = blue_light
    operator::String = red
    func::String = yellow

    # loggin levels
    info::String = "#7cb0cf"
    debug::String = "#197fbd"
    warn::String = orange
    error::String = "bold #d13f3f"

    # Tree objects
    tree_title_style::String = "$orange italic"
    tree_node_style::String = "$yellow italic"
    tree_leaf_style::String = yellow_light
    tree_guide_style::String = blue
    tree_max_width::Int = 44

    # Repr
    repr_accent_style = "bold #e0db79"
    repr_name_style = "#e3ac8d"
    repr_type_style = "#bb86db"
    repr_values_style = "#b3d4ff"
    repr_line_style = "dim #7e9dd9"
    repr_panel_style = "#9bb3e0"
end

function Base.show(io::IO, ::MIME"text/plain", theme::Theme)
    fields = fieldnames(Theme)
    N = length(fields)
    values = map(f -> getfield(theme, f), fields)

    fields = map(v -> if v[2] isa String
        RenderableText(string(v[1]); style = v[2])
    else
        RenderableText(string(v[1]))
    end, zip(fields, values))

    values = map(
        v -> v isa String ? RenderableText("■■"; style = v) : RenderableText(string(v)),
        values,
    )

    content =
        hLine(30, "Base"; style = "#9bb3e0") /
        (rvstack(values[2:13]) * Spacer(2, 1) * lvstack(fields[2:13]))
    content /= "" / hLine(30, "Logging"; style = "#9bb3e0")
    content /= (rvstack(values[14:17]) * Spacer(2, 1) * lvstack(fields[14:17]))
    content /= "" / hLine(30, "Tree"; style = "#9bb3e0")
    content /= (rvstack(values[18:22]) * Spacer(2, 1) * lvstack(fields[18:22]))
    content /= "" / hLine(30, "REPL"; style = "#9bb3e0")
    content /= (rvstack(values[23:end]) * Spacer(2, 1) * lvstack(fields[23:end]))

    return print(
        io,
        Panel(
            content;
            width = 44,
            justify = :center,
            title = "Theme: {bold}$(theme.name){/bold}",
            padding = (4, 4, 1, 1),
            style = "#9bb3e0",
            title_style = "white",
        ),
    )
end

set_theme(theme::Theme) = (term_theme[] = theme)

# ------------------------------ Highlighters.jl ----------------------------- #

"""
Custom hilighting theme for Highlighters.jl
https://juliadocs.github.io/Highlights.jl/stable/man/theme/
"""
abstract type CodeTheme <: AbstractTheme end

@theme CodeTheme Dict(
    :style => S"",
    :tokens => Dict(
        TEXT => S"fg: dedede;",

        # yellow
        NAME_FUNCTION => S"fg: e8d472;",
        NAME_OTHER => S"fg: e8d472;",

        # red
        KEYWORD => S"fg: 7a93f5;",
        OPERATOR => S"fg: de6d59;",
        PUNCTUATION => S"fg: e38864",

        # green
        STRING => S"fg: 50ad5f",
        COMMENT => S"fg: 287a36; italic",
        STRING_DOC => S"fg: 50ad5f",

        # blue
        NUMBER => S"fg: 90CAF9",
    ),
)


text:
# ============================================================
# Neural & Universal Closures
# ============================================================

export neural_closure, create_neural_closure

"""
    neural_closure(X, model)

Neural closure model computation

@args X: model input
@args model: neural network model
"""
function neural_closure(model, X::T) where {T<:AbstractArray}
    return Zygote.gradient(x -> first(model(x)), X)[1]
end

function neural_closure(model::T, p, X::T1) where {T<:DiffEqFlux.FastLayer,T1<:AbstractArray}
    return Zygote.gradient(x -> first(model(x, p)), X)[1]
end

function neural_closure(model::PyObject, X::T) where {T<:AbstractArray}
    py"""
    import tensorflow as tf

    def dnn(model, input):
        '''
        Input: input.shape = (nCells,nMaxMoment), nMaxMoment = 4 in case of MK7
        Output: Gradient of the network wrt input
        '''
        x_model = tf.Variable(input)
        with tf.GradientTape() as tape:
            predictions = model(x_model) # model.predict(x_model) doesn't work
        gradients = tape.gradient(predictions, x_model)
        return gradients.numpy()
    """

    return py"dnn"(model, X)
end

"""
    create_neural_closure(imputDim, outputDim; acfun = relu)

Create neural closure model
"""
function create_neural_closure(
    Din::T,
    Dout::T,
    Dhid = 10::T,
    Nhid = 1::T;
    acfun = relu,
    mode = :icnn,
) where {T<:Integer}
    return ICNN(Din, Dout, ones(Int, Nhid) * Dhid, acfun)
end


text:
#=
util:
- Julia version:
- Author: bramb
- Date: 2019-01-24
=#

function _test_normalized_strings_are_equal(string1::String, string2::String)
    n1 = _normalize(string1)
    n2 = _normalize(string2)
    @test n1 == n2
end

function _normalize(string::String)
    trimmed_lines = map(l->strip(l), split(string,"\n"))
    trim = join(trimmed_lines,"\n")
    return replace(trim, r"\n\n+" => "\n")
end

function _print_dot(dot::String)
   out = """
   ---------- ✂ ------------------------------------------------------------------
   $(replace(dot, r"\n([\t ]*\n)+" => "\n"))
   -------------------------------------------------- ✂ --------------------------
   """
   println(out)
end


text:
# Autogenerated wrapper script for Couenne_jll for i686-linux-gnu-cxx11
export amplexe, libcouenne

using ASL_jll
using Bonmin_jll
using Cbc_jll
using Ipopt_jll
JLLWrappers.@generate_wrapper_header("Couenne")
JLLWrappers.@declare_library_product(libcouenne, "libCouenne.so.1")
JLLWrappers.@declare_executable_product(amplexe)
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, Bonmin_jll, Cbc_jll, Ipopt_jll)
    JLLWrappers.@init_library_product(
        libcouenne,
        "lib/libCouenne.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        amplexe,
        "bin/couenne",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DynamicHMCModels

Random.seed!(1)

cd(@__DIR__)
include("simulateLR.jl")

 struct RegressionProb
    x::Array{Float64,2}
    y::Array{Float64,1}
    Nd::Int64
    Nc::Int64
end

function (problem::RegressionProb)(θ)
    @unpack x,y,Nd,Nc = problem   # extract the data
    @unpack B0,B,sigma = θ
    μ = B0 .+x*B
    sum(logpdf.(Normal.(μ,sigma),y)) + logpdf(Normal(0,10),B0) +
    loglikelihood(Normal(0,10),B) + logpdf(Truncated(Cauchy(0,5),0,Inf),sigma)
end

# Define problem with data and inits.

x, y, Nd, Nc = simulateLR()
p = RegressionProb(x, y, Nd, Nc)
p((B0 = 0.0, B = fill(0.0, Nc), sigma = 1.0))

# Write a function to return properly dimensioned transformation.

problem_transformation(p::RegressionProb) =
    as((B0=asℝ, B=as(Array, asℝ, Nc), sigma = asℝ₊))

# Use Flux for the gradient.

P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P)

#import Zygote
#∇P = ADgradient(:Zygote, P)

# Sample from the posterior.

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 4000);

# Undo the transformation to obtain the posterior from the chain.

posterior = TransformVariables.transform.(Ref(problem_transformation(p)), get_position.(chain));

chns = nptochain(posterior,NUTS_tuned)


text:
export save_network, load_network

############################################################
# General IO utilities
############################################################
function mkdir_p(path::AbstractString)
  # I bet $5 this bad recursion is not a problem for our usages here
  if isempty(path) || isdir(path)
    return
  end

  if !isdir(dirname(path))
    mkdir_p(dirname(path))
  end
  mkdir(path)
end

# add this function since the built-in tempname() function is not
# working properly on Windows. See https://github.com/JuliaLang/julia/issues/9053
function temp_filename()
  tmp_dir = tempdir()
  joinpath(tmp_dir, "Mocha-$(getpid())-$(randstring(32))")
end

function glob(path::AbstractString, pattern::Regex; sort_by :: Symbol = :none)
  list = filter(x -> ismatch(pattern, x), readdir(path))
  if sort_by == :none
    return list
  elseif sort_by == :name
    return sort(list)
  elseif sort_by == :mtime
    return sort(list, by = fn -> begin
      stat(joinpath(path, fn)).mtime
    end)
  else
    error("Unknown sort_by: $sort_by")
  end
end

############################################################
# JLD IO procedures to save and load network parameters
############################################################
using HDF5, JLD
const NETWORK_SAVE_NAME = "params_all"

function save_network(file::JLD.JldFile, net)
  params_all = Dict{AbstractString, Vector{Array}}()
  for i = 1:length(net.layers)
    if has_param(net.layers[i])
      key = net.layers[i].name
      if haskey(params_all, key)
        error("Duplicated names ($key) for multiple layers, cannot save parameters")
      end

      @debug("Saving parameters for layer $(key)")
      params = map(param -> to_array(param.blob), net.states[i].parameters)
      params_all[key] = params
    end
  end
  write(file, NETWORK_SAVE_NAME, params_all)
end

function is_similar_shape(x, y)
  function compact_shape(obj)
    shape = [size(obj)...]
    mask  = ones(Bool, length(shape))
    # strip singleton dimension at the two end
    for i = 1:length(shape)
      if shape[i] == 1
        mask[i] = false
      else
        break
      end
    end
    for i = length(shape):-1:1
      if shape[i] == 1
        mask[i] = false
      else
        break
      end
    end
    return shape[mask]
  end

  compact_shape(x) == compact_shape(y)
end

function load_network(file::JLD.JldFile, net)
  params_all = read(file, NETWORK_SAVE_NAME)
  for i = 1:length(net.layers)
    if has_param(net.layers[i])
      key = net.layers[i].name
      if !haskey(params_all, key)
        error("Cannot find saved parameters for layer $key")
      end

      @debug("Loading parameters for layer $(key)")
      params = params_all[key]
      @assert length(params) == length(net.states[i].parameters)

      for j = 1:length(params)
        @assert is_similar_shape(params[j], net.states[i].parameters[j].blob)
        copy!(net.states[i].parameters[j].blob, params[j])
        net.states[i].parameters[j].initializer = NullInitializer()
      end
    end
  end
end

############################################################
# Load network parameters from a HDF5 file
# When mapped correctly, this could be used to import
# trained network from other tools (e.g. Caffe)
############################################################
function load_network(file::HDF5File, net, die_if_not_found=true)
  for i = 1:length(net.layers)
    if has_param(net.layers[i])
      layer_name = net.layers[i].name
      @debug("Loading parameters from HDF5 for layer $layer_name")
      for j = 1:length(net.states[i].parameters)
        param_obj  = net.states[i].parameters[j]
        param_name = param_obj.name
        key = "$(layer_name)___$(param_name)"
        if !has(file, key)
          if param_name == "bias"
            @warn("No bias found for $layer_name, use default initialization")
            init(param_obj.initializer, param_obj.blob)
          else
            if die_if_not_found
              error("No saved parameter $param_name not found for layer $layer_name")
            else
              # show error message but do not die
              @error("No saved parameter $param_name not found for layer $layer_name, use default initialization")
              init(param_obj.initializer, param_obj.blob)
            end
          end
        else
          param = read(file, key)
          if !is_similar_shape(param, param_obj.blob)
            error("Dimension for $param_name not match: got $(size(param)), expect $(size(param_obj.blob))")
          end

          if eltype(param) != eltype(param_obj.blob)
            @warn("Automatic converting saved $param_name from $(eltype(param)) to $(eltype(param_obj.blob))")
            param = convert(Array{eltype(param_obj.blob)}, param)
          end
          copy!(param_obj.blob, param)
        end
        param_obj.initializer = NullInitializer()
      end
    end
  end
end


text:
"""
    FEMMDeforSurfaceDampingModule

Module for operations on the damping associated with absorbing boundary
conditions (ABC) representation of the effect of infinite extent
of inviscid fluid next to the surface.
"""
module FEMMDeforSurfaceDampingModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, nodesperelem, manifdim
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobiansurface
import FinEtools.FieldModule: ndofs, gatherdofnums!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.FEMMBaseModule: FEMMAbstractBase
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!
import FinEtools.MatrixUtilityModule: add_nnt_ut_only!, complete_lt!, locjac!
import FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
import LinearAlgebra: norm, cross

"""
    FEMMDeforSurfaceDamping{S<:FESet, F<:Function}

Type for surface damping model.
"""
mutable struct FEMMDeforSurfaceDamping{S<:FESet, F<:Function} <: FEMMAbstractBase
    integdomain::IntegDomain{S, F} # geometry data
end

"""
    dampingABC(self::FEMMDeforSurfaceDamping, assembler::A,
                  geom::NodalField{FFlt}, u::NodalField{T1},
                  impedance::T2, surfacenormal::SurfaceNormal) where {A<:SysmatAssemblerBase, T1<:Number, T2<:Number}

Compute the damping matrix associated with absorbing boundary conditions (ABC) representation of the effect of infinite extent of inviscid fluid next to the surface.
"""
function dampingABC(self::FEMMDeforSurfaceDamping, assembler::A,
    geom::NodalField{FFlt}, u::NodalField{T1},
    impedance::T2, surfacenormal::SurfaceNormal) where {A<:SysmatAssemblerBase, T1<:Number, T2<:Number}
    fes = self.integdomain.fes
    # Constants
    nfes = count(fes); # number of finite elements
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes per finite element
    sdim = ndofs(geom); # spatial dimension
    mdim = manifdim(fes); # manifold dimension of the finite elements
    Cedim = ndn * nne; # size of damping element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc  =  integrationdata(self.integdomain);
    loc = zeros(FFlt, 1, sdim); # quadrature point coordinate -- used as a buffer
    J = zeros(FFlt, sdim, mdim); # Jacobian matrix -- used as a buffer
    Ce = zeros(T2, Cedim, Cedim); # element damping matrix -- used as a buffer
    Nn = zeros(FFlt, Cedim); # column vector
    dofnums = zeros(FInt, Cedim); # degree of freedom array -- used as a buffer
    # Prepare assembler and temporaries
    startassembly!(assembler, Cedim, Cedim, nfes, u.nfreedofs, u.nfreedofs);
    for i = 1:nfes # loop over finite elements
        fill!(Ce, 0.0); # Initialize element damping matrix
        for j = 1:npts # loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            n = updatenormal!(surfacenormal, loc, J, self.integdomain.fes.label[i]);
            for k = 1:nne
                Nn[(k-1)*ndn+1:k*ndn] = n * Ns[j][k]
            end
            add_nnt_ut_only!(Ce, Nn, (-1.0)*impedance*Jac*w[j]);
        end # end loop over quadrature points
        complete_lt!(Ce);
        gatherdofnums!(u, dofnums, self.integdomain.fes.conn[i]); # retrieve degrees of freedom
        assemble!(assembler, Ce, dofnums, dofnums); # assemble the element damping matrix
    end # end loop over finite elements
    return makematrix!(assembler);
end

function dampingABC(self::FEMMDeforSurfaceDamping, geom::NodalField{FFlt},
    u::NodalField{T1}, impedance::T2, surfacenormal::SurfaceNormal) where {T1<:Number, T2<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return dampingABC(self, assembler, geom, u, impedance, surfacenormal);
end

end


text:
#=
    Figure S2
=#

using Model # To get colours

using CSV
using DataFrames
using DataFramesMeta
using Plots
using StatsPlots

##############################################################
## Data
##############################################################

data = CSV.read("Data/CSV/DualMarker.csv",DataFrame)

# Get autofluorescence information
E = [@subset(data,:Time .== 0.0, :Quenched .== false, :Temp .== 37.0).Signal_Cy5,
     @subset(data,:Time .== 0.0, :Quenched .== false, :Temp .== 37.0).Signal_BDP]

##############################################################
## Figure S2
##############################################################

figS2 = density2d(E[1],E[2],
    colorbar=false,
    xlim=(0.0,20.0),
    ylim=(0.0,20.0),
    box=:on,axis=:all,aspect_ratio=:equal,lw=0.0,c=:RdPu_9,levels=8,fill=true,
    xlabel="Cy5 autofluoresence",
    ylabel="BDP autofluoresence"
)

##############################################################
## Save
##############################################################

savefig("$(@__DIR__)/figS2.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
using PyPlot, JUDI, LinearAlgebra, FFTW, JOLI, DSP
using IterativeSolvers

using Shampoo

n = (401, 401)   # (x,y,z) or (x,z)
d = (5f0, 5f0)
o = (0f0, 0f0)

extentx = n[1]*d[1]
extentz = n[2]*d[2]

v = 2f0*ones(Float32,n)

# Slowness squared [s^2/km^2]
m = convert(Array{Float32,2},(1f0 ./ v).^2)	# true model

model = Model(n, d, o, m; nb = 200)

# Model structure

dtS = 1f0
timeS = 3000f0
nt = Int64(timeS/dtS)+1

fmin = 10f0
fmax = 60f0
cfreqs = (fmin, fmin+10f0, fmax-10f0, fmax) # corner frequencies
wavelet = reshape(cfreq_wavelet(500, nt, dtS/1f3, cfreqs; edge=hamming),nt,1)
#wavelet = ricker_wavelet(timeS, dtS, 0.04f0)[:,1]

################################### Source/receiver geometry ######################################

# Set up receiver geometry

domain_x = (model.n[1] - 1)*model.d[1]
domain_z = (model.n[2] - 1)*model.d[2]

nxrec = 512

function circleShape(h, k ,r)
	theta = LinRange(0,2*pi, nxrec)
	h .+ r*sin.(theta), k .+ r*cos.(theta)
end

xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-5*d[2])
yrec = 0f0 #2d so always 0

# receiver sampling and recording time
timeR = 3000f0   # receiver recording time [ms]
dtR   = 1f0
# Set up receiver structure
nsrc = 1
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

mute_matrix = zeros(Float32,n)

for i = 1:n[1]
	for j = 1:n[2]
		mute_matrix[i,j] = ((i-201)^2+(j-201)^2<190^2)
	end
end

M = MaskOp(n,mute_matrix)

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

# Return shots as julia vector (needed for tracker)
opt = Options(free_surface=false,isic=false)

Pr  = judiProjection(info, recGeometry)
F  = judiModeling(info, model; options=opt)	# modeling operator initial model
Pw = judiLRWF(info, wavelet)
# Combined operators
F = Pr*F*adjoint(Pw)

#P = FractionalIntegrationOp(nt,nsrc,nxrec,0.95)
P = DiffOp(nt,nsrc,nxrec)

# Extended modeling CNN layers
q = zeros(Float32,n)
q[201,201] = 1f0
q = judiWeights(q)

PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
figure(); imshow(reshape(sqrt.(1f0./m), n[1], n[2])',extent=(0,extentx,extentz,0));
PyPlot.scatter(xrec,zrec,marker=".",color="red",label="receivers",s=1);
PyPlot.scatter(201*d[1],201*d[2],label="sources",marker="x",color="white");
legend(loc=3,fontsize=8);
xlabel("Lateral position [m]", fontsize=8);ylabel("Horizontal position [m]", fontsize=8);
fig = PyPlot.gcf();
title("Experimental setup")

######################################### Set up JUDI operators ###################################

d_obs = F*q
q_new = F'*d_obs

q_precond = F'*P'*P*d_obs

figure();
subplot(2,1,1);
PyPlot.plot(wavelet[1:1000]);title("zoomed-in wavelet in time domain");
subplot(2,1,2);
PyPlot.plot(abs.(fftshift(fft(wavelet))));title("wavelet in frequency domain");

figure();
subplot(2,2,1)
imshow(reshape(q_new.weights[1],n)[130:270,130:270]',vmin=-0.8*norm(q_new.weights[1],Inf),vmax=0.8*norm(q_new.weights[1],Inf));title("F'*F*q")
subplot(2,2,2)
imshow(reshape(q_precond.weights[1],n)[130:270,130:270]',vmin=-0.8*norm(q_precond.weights[1],Inf),vmax=0.8*norm(q_precond.weights[1],Inf));title("F'*P'*P*F*q")
subplot(2,2,3)
imshow(abs.(fftshift(fft(reshape(q_new.weights[1],n)')))/norm(abs.(fftshift(fft(reshape(q_new.weights[1],n)'))),Inf),cmap="jet",vmin=0,vmax=1)
subplot(2,2,4)
imshow(abs.(fftshift(fft(reshape(q_precond.weights[1],n)')))/norm(abs.(fftshift(fft(reshape(q_precond.weights[1],n)'))),Inf),cmap="jet",vmin=0,vmax=1)

maxit = 10

q1 = 0f0 .* q
q1,his1 = lsqr!(q1,F*M,d_obs,atol=0f0,btol=0f0,conlim=0f0,maxiter=maxit,log=true,verbose=true)
q2 = 0f0 .* q
q2,his2 = lsqr!(q2,P*F*M,P*d_obs,atol=0f0,btol=0f0,conlim=0f0,maxiter=maxit,log=true,verbose=true)
res1 = his1[:resnorm]
res2 = his2[:resnorm]

u1 = 20*log.(res1./norm(d_obs))
u2 = 20*log.(res2./norm(P*d_obs))

figure();
xlabel("iterations")
ylabel("ratio")
PyPlot.plot(res1/norm(d_obs),label="LSQR")
PyPlot.plot(res2/norm(P*d_obs),label="P-LSQR")
title("normalized least squares residual")
legend()

figure();
xlabel("iterations")
ylabel("dB")
PyPlot.plot(u1,label="LSQR")
PyPlot.plot(u2,label="P-LSQR")
title("normalized log-based least squares residual")
legend()

figure();
subplot(2,2,1)
imshow(reshape((M*q1).weights[1],n)[130:270,130:270]',vmin=-0.8*norm((M*q1).weights[1],Inf),vmax=0.8*norm((M*q1).weights[1],Inf));title("LSQR")
subplot(2,2,2)
imshow(reshape((M*q2).weights[1],n)[130:270,130:270]',vmin=-0.8*norm((M*q2).weights[1],Inf),vmax=0.8*norm((M*q2).weights[1],Inf));title("P-LSQR")
subplot(2,2,3)
imshow(abs.(fftshift(fft(reshape((M*q1).weights[1],n)'))),cmap="jet")
subplot(2,2,4)
imshow(abs.(fftshift(fft(reshape((M*q2).weights[1],n)'))),cmap="jet")


text:
using GtkSourceWidget, Test
using Gtk

@testset "GtkSourceView" begin

    w = GtkWindow()
    b = GtkSourceBuffer()
    b.text[String] = "test"
    v = GtkSourceView(b)
    push!(w,v)
    showall(w)
    @test b.text[String] == "test"

    search_settings = GtkSourceSearchSettings()
    search_context = GtkSourceSearchContext(b, search_settings)
    set_search_text(search_settings, "test")

    it = Gtk.GtkTextIter(b, 1)
    found, its, ite = search_context_forward(search_context, it)
    @test found == true
    @test (its:ite).text[String] == "test"

    search_context_replace(search_context, its, ite, "it worked!")
    set_search_text(search_settings, "it worked!")

    it = Gtk.GtkTextIter(b, 1)
    found, its, ite = search_context_forward(search_context, it)
    @test found == true
    @test (its:ite).text[String] == "it worked!"

    mark = create_mark(b, it)
    scroll_to(v, mark, 0, true, 0.5, 0.5)

end


text:
export main

# import Base:
# REPL,
# Terminals,
# Vector,
# println,
# STDIN,
# STDOUT,
# STDERR,
# pushdisplay,
# VersionNumber,
# print,
# source_path,
# Dict

using ..conf
using ..util

using ArgParse
# # test if the program works
# # load all necessary files
# load("maths.jl")
# load("maths.jl")
# load("config.jl")
# load("wordstate.jl")
# load("helper.jl")
# load("lexicon.jl")
# load("gibbs_sampler.jl")
# load("run_all.jl")
# load("blade_runner.jl")

# data_dir = "data-sample/"
# out_prefix = "model-"
# out_dir = "models/"

# book = "new_wordlist.eng"
# filename = "$(data_dir)$(book)"
# outdict = "$(out_dir)$(out_prefix)$(book)"
# gold_cols = read_data(filename)
# cluster_file = "cluster-embedding.eng"

# gold_tags = nothing
# hyperparameters = Dict()
# seq = false
# suffix = false
# prefix = false
# num_processes = 1
# numit = 30
# num_tags = 1
# fix_tags = false
# fix_segs = false
# fix_stem = false
# init_tag = "uniform"
# init_seg = "whole-word"
# init_stem = "first"
# state0 = nothing
# sep_lex_size = true
# flags = nothing
# phrase_counts = gold_cols[1]
# phrases = gold_cols[2]
# gold = gold_cols[2]
# cluster = nothing
# outdicts_state = "$(outdict)-state"
# params = Params(phrase_counts, cluster, phrases, gold, gold_tags, seq, suffix, prefix, numit, fix_tags, fix_segs, fix_stem, init_tag, init_seg, init_stem, num_tags, state0, outdicts_state, sep_lex_size)
# run_all(params)

const arg_parse_settings = ArgParseSettings(
    prog=XSH_BIN,
    description=XSH_ARG_DESCRIPTION,
    epilog=XSH_ARG_EPILOGUE,
    usage = XSH_ARG_USAGE,
    version = "$(VERSION)",
    add_help = false,
    add_version = false,
    commands_are_required = false,
    autofix_names = true,
    suppress_warnings = true,
    allow_ambiguous_opts = true
    )

arg_parse_data = Dict{Symbol, UTF8String}()

function init_arg_table()

    @add_arg_table arg_parse_settings begin
        "file"
        help = "Evaluate xsh script"
        "cmd"
        help = "Run command"
        action = :command
        "--help", "-h"
        help = "Print this message"
        action = :show_help
        "--version", "-v"
        help = "Display version information"
        action = :show_version
        "--quiet", "-q"
        help = "Quiet operation, only errors"
        action = :store_true
        "--verbose", "-V"
        help = "Specify verbosity level, 0≦V≦5 Default level [1]"
        action = :store_true
        "--debug", "-d"
        help = "Enable debug output"
        action = :store_true
        "--home", "-H"
        help = "Set the location of xsh home"
        action = :store_true
        "--interactive", "-i"
        help = "Run in interactive mode"
        action = :store_true
        "--login", "-l"
        help = "Run as a login shell"
        action = :store_true
        "--daemon", "-m"
        help = "Run and detach xsh daemon, then exit"
        action = :store_true
        "--ncpu", "-n"
        help = "Run [n] local processes"
        action = :store_true
        "--print", "-E"
        help = "Evaluate and show <cmd>"
        action = :store_true
        "--color", "-c"
        help = "Enable or disable color text"
        action = :store_true
        "--eval", "-e"
        help = "Evaluate <expr>"
        action = :store_true
        "--no-exec"
        help = "Do not execute, only perform syntax checking"
        action = :store_true
        "--no-history-file"
        help = "Do not archive command history"
        action = :store_true
        "--no-startup"
        help = "Do not load user start-up script [~/.xshrc]"
        action = :store_true
    end
end

function active(val::ANY)
    if val == false
        return false
    end
    if val == nothing
        return false
    end
    return true
end

function parse_argument(key::ASCIIString, value::ANY)
    @switch key begin
        @case "home"
        # xsh_config.home_path = value
        break
        @case "login"
        # xsh_config.login = value
        break
        @case "no_exec"
        # xsh_config.no_exec = value
        break
        @case "no_history_file"
        # xsh_config.no_history = value
        break
        @case "quiet"
        # xsh_config.quiet = value
        break
        @case "interactive"
        # xsh_config.interactive = value
        break
        @case "eval"
        # println("eval")
        break
        @case "debug"
        # xsh_config.debug_mode = value
        break
        @case "verbose"
        # xsh_config.verbose = value
        break
        @case "no_startup"
        # xsh_config.no_startup = value
        break
        @case "daemon"
        # xsh_config.daemon_only = value
        break
        @case "print"
        # println("print")
        break
        @case "file"
        print("file : $(value)")
        # conf.file_name = value
        # print("$(conf.file_name) : $(value)")
        # print("$(conf.file_name) : $(value)")
        break
        @case "color"
        # xsh_config.color_mode = value
        break
        @case "ncpu"
        # xsh_config.ncpu = value
        break
        @default
        println("Not Found [$(key)][$(value)]")
    end
end

function parse_args(args::Vector{UTF8String}=UTF8String[])
    parsed_args =ArgParse.parse_args(args, arg_parse_settings)
    for (key,val) in parsed_args
        if active(val)
            parse_argument(key, val)
        end
    end
end

function init()
    init_arg_table()

end

function main(args::Vector{UTF8String}=UTF8String[])
    init()
    println(XSH_DESCRIPTION)
    parsed_args = parse_args(args)
    # println(conf.file_name)
    # println(file_name)
    # dump(config)
end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
"""
    kfolds(n::Integer, k = 5) -> Tuple

Compute the train/validation assignments for `k` repartitions of
`n` observations, and return them in the form of two vectors. The
first vector contains the index-vectors for the training subsets,
and the second vector the index-vectors for the validation subsets
respectively. A general rule of thumb is to use either `k = 5` or
`k = 10`. The following code snippet generates the indices
assignments for `k = 5`

```julia
julia> train_idx, val_idx = kfolds(10, 5);
```

Each observation is assigned to the validation subset once (and
only once). Thus, a union over all validation index-vectors
reproduces the full range `1:n`. Note that there is no random
assignment of observations to subsets, which means that adjacent
observations are likely to be part of the same validation subset.

```julia
julia> train_idx
5-element Array{Array{Int64,1},1}:
 [3,4,5,6,7,8,9,10]
 [1,2,5,6,7,8,9,10]
 [1,2,3,4,7,8,9,10]
 [1,2,3,4,5,6,9,10]
 [1,2,3,4,5,6,7,8]

julia> val_idx
5-element Array{UnitRange{Int64},1}:
 1:2
 3:4
 5:6
 7:8
 9:10
```
"""
function kfolds(n::Integer, k::Integer = 5)
    2 <= k <= n || throw(ArgumentError("n must be positive and k must to be within 2:$(max(2,n))"))
    # Compute the size of each fold. This is important because
    # in general the number of total observations might not be
    # divideable by k. In such cases it is custom that the remaining
    # observations are divided among the folds. Thus some folds
    # have one more observation than others.
    sizes = fill(floor(Int, n/k), k)
    for i = 1:(n % k)
        sizes[i] = sizes[i] + 1
    end
    # Compute start offset for each fold
    offsets = cumsum(sizes) .- sizes .+ 1
    # Compute the validation indices using the offsets and sizes
    val_indices = map((o,s)->(o:o+s-1), offsets, sizes)
    # The train indices are then the indicies not in validation
    train_indices = map(idx->setdiff(1:n,idx), val_indices)
    # We return a tuple of arrays
    train_indices, val_indices
end

"""
    kfolds(data, [k = 5])

Repartition a `data` container `k` times using a `k` folds
strategy and return the sequence of folds as a lazy iterator.
Only data subsets are created, which means that no actual data is copied until
[`getobs`](@ref) is invoked.

Conceptually, a k-folds repartitioning strategy divides the given
`data` into `k` roughly equal-sized parts. Each part will serve
as validation set once, while the remaining parts are used for
training. This results in `k` different partitions of `data`.

In the case that the size of the dataset is not dividable by the
specified `k`, the remaining observations will be evenly
distributed among the parts.

```julia
for (x_train, x_val) in kfolds(X, k=10)
    # code called 10 times
    # nobs(x_val) may differ up to ±1 over iterations
end
```

Multiple variables are supported (e.g. for labeled data)

```julia
for ((x_train, y_train), val) in kfolds((X, Y), k=10)
    # ...
end
```

By default the folds are created using static splits. Use
[`shuffleobs`](@ref) to randomly assign observations to the
folds.

```julia
for (x_train, x_val) in kfolds(shuffleobs(X), k = 10)
    # ...
end
```

See [`leavepout`](@ref) for a related function.
"""
function kfolds(data, k::Integer)
    n = numobs(data)
    train_indices, val_indices = kfolds(n, k)

    ((obsview(data, itrain), obsview(data, ival))
        for (itrain, ival) in zip(train_indices, val_indices))
end

kfolds(data; k) = kfolds(data, k)

"""
    leavepout(n::Integer, [size = 1]) -> Tuple

Compute the train/validation assignments for `k ≈ n/size`
repartitions of `n` observations, and return them in the form of
two vectors. The first vector contains the index-vectors for the
training subsets, and the second vector the index-vectors for the
validation subsets respectively. Each validation subset will have
either `size` or `size+1` observations assigned to it. The
following code snippet generates the index-vectors for `size = 2`.

```julia
julia> train_idx, val_idx = leavepout(10, 2);
```

Each observation is assigned to the validation subset once (and
only once). Thus, a union over all validation index-vectors
reproduces the full range `1:n`. Note that there is no random
assignment of observations to subsets, which means that adjacent
observations are likely to be part of the same validation subset.

```julia
julia> train_idx
5-element Array{Array{Int64,1},1}:
 [3,4,5,6,7,8,9,10]
 [1,2,5,6,7,8,9,10]
 [1,2,3,4,7,8,9,10]
 [1,2,3,4,5,6,9,10]
 [1,2,3,4,5,6,7,8]

julia> val_idx
5-element Array{UnitRange{Int64},1}:
 1:2
 3:4
 5:6
 7:8
 9:10
```
"""
function leavepout(n::Integer, p::Integer = 1)
    1 <= p <= floor(n/2) || throw(ArgumentError("p must to be within 1:$(floor(Int,n/2))"))
    k = floor(Int, n / p)
    kfolds(n, k)
end

"""
    leavepout(data, p = 1)

Repartition a `data` container using a k-fold strategy, where `k`
is chosen in such a way, that each validation subset of the
resulting folds contains roughly `p` observations. Defaults to
`p = 1`, which is also known as "leave-one-out" partitioning.

The resulting sequence of folds is returned as a lazy
iterator. Only data subsets are created. That means no actual
data is copied until [`getobs`](@ref) is invoked.

```julia
for (train, val) in leavepout(X, p=2)
    # if nobs(X) is dividable by 2,
    # then numobs(val) will be 2 for each iteraton,
    # otherwise it may be 3 for the first few iterations.
end
```

See[`kfolds`](@ref) for a related function.
"""
function leavepout(data, p::Integer)
    n = numobs(data)
    1 <= p <= floor(n/2) || throw(ArgumentError("p must to be within 1:$(floor(Int,n/2))"))
    k = floor(Int, n / p)
    kfolds(data, k)
end

leavepout(data; p::Integer=1) = leavepout(data, p)


text:
@testset "EndgameTracker" begin
    @testset "Tracking and PathResult" begin
        @var x y
        f = [2.3 * x^2 + 1.2 * y^2 + 3x - 2y + 3, 2.3 * x^2 + 1.2 * y^2 + 5x + 2y - 5]
        tracker, starts = total_degree(System(f); gamma = 1.3im + 0.4)
        S = collect(starts)
        @test !isempty(sprint(show, tracker.options))

        res = map(i -> track(tracker, S[i]; path_number = i), 1:4)

        @test is_success(res[1])
        @test path_number.(res) == 1:4
        @test start_solution.(res) == S
        @test all(isfinite, solution(res[1]))
        @test is_finite(res[1])
        @test isfinite(res[1])
        @test multiplicity(res[1]) == nothing
        @test winding_number(res[1]) == nothing
        @test !is_failed(res[1])
        @test accuracy(res[1]) < 1e-12
        @test residual(res[1]) < 1e-12
        @test steps(res[1]) < 10
        @test accepted_steps(res[1]) < 10
        @test rejected_steps(res[1]) == 0
        @test is_success(res[2])
        @test is_nonsingular(res[2])
        @test cond(res[2]) < 1e3
        @test is_at_infinity(res[3])
        @test is_at_infinity(res[4])
        @test !isempty(sprint(show, res[1]))
        @test valuation(res[3]) ≈ [-1, -1] rtol = 1e-3
        # singular stuff
        @var x
        f = [(x - 10)^2]
        tracker, starts = total_degree(System(f))
        res = track.(tracker, starts)
        @test winding_number(res[1]) == 2
        @test last_path_point(res[1]) isa Tuple{Vector{ComplexF64},Float64}
        @test 0 < last(last_path_point(res[1])) < 0.1
        @test all(is_real, res)
        @test all(isreal, res)
        @test all(is_singular, res)
    end

    @testset "Overdetermined tracking" begin
        @var x y z

        p₁ = (x^2 + y^2 + z^2 - 1) * (x - 0.5)
        p₂ = (x^2 + y^2 + z^2 - 1) * (y - 0.5)
        p₃ = (z - x^2 - 2) * (x^2 + y^2 + z^2 - 1) * (z - 0.5)
        F = System([p₁, p₂, p₃])

        L₁ = rand_affine_subspace(3; codim = 2)
        L₂ = rand_affine_subspace(3; codim = 2)
        F_L₁ = System([x^2 + y^2 + z^2 - 1]) ∩ extrinsic(L₁)([x, y, z])
        res = track.(total_degree(F_L₁)...)
        @test all(is_success, res)

        H = AffineSubspaceHomotopy(F, L₁, L₂)
        res2 = track.(EndgameTracker(H), solution.(res))
        @test all(is_success, res2)
    end
end


text:
@testset "QC" begin
    let t = read_test_data()
        e, n = t
        s = splitting(e, n, e.meta.SAC_user0, e.meta.SAC_user2)
        @test snr_restivo_helffrich(s) ≈ 19.695 atol=2
        s′ = splitting(Seis.rotate_through(n, e, rand(1:360))...,
                       e.meta.SAC_user0, e.meta.SAC_user2)
        @test snr_restivo_helffrich(s) ≈ snr_restivo_helffrich(s′)
        @test quality(s) ≈ 0.902 atol=0.02
        @test quality(s) ≈ quality(s′)
        s″ = splitting(n, e, xcorr=false)
        @test_throws ArgumentError quality(s″)
    end
end

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
# Current version has flaws since update to methods will fixed in future
#=
"""
--------------------------------------------------------------------------------
Function: Generalized_Param_Bisection
--------------------------------------------------------------------------------
Description:
Uses the interval arithmetic presented in ValidatedNumerics to bound unique
implicit functions x:P->X defined by a system of equations h(z,p)=0 via
Gauss-Siedel Newton.
--------------------------------------------------------------------------------
Inputs:
* Xin:       Vector{Interval{Float64}} - Bounds for dependent variables
* Pin:       Vector{Interval{Float64}} - Bounds for independent variables
* h:         function - Equations defining potential implicit function
* hj:        function - Jacobian of h(x,p) with respect to x
* opt:       Param_Bisect_Opts - Parametric Bisection Options
--------------------------------------------------------------------------------
Returns:
A tuple (sol_set, ind_set, NIT, hist) where
* sol_set:     Array - Interval boxes known to contain unique implicit function
* ind_set:     Array - Interval boxes known to contain unique implicit function
* NIT:         Int64 - Final iteration number
* hist:        Array - Storage object with info on problem solution history
--------------------------------------------------------------------------------
"""
#=
function Generalized_Param_Bisection(Xin,Pin,h,hj,opt)

  # unpacks options for general bisection
  DAGflag::Bool = opt.DAGflag
  LPflag::Bool = opt.LPflag
  max_iter = opt.kmax_main
  max_iter_cntr = opt.kmax_cntr
  style = opt.style
  display = opt.display
  ptol = opt.ptol
  etol = opt.etol
  rtol = opt.rtol
  pbisect = opt.p_rel_bisect
  Pstart = copy(Pin)
  bp = 1

  if (DAGflag)
    println("started generating DAG parameters")
    DAGr = opt.DAGpass
    DAGh = opt.DAGh
    DAGg = opt.DAGg
    DAGsym = opt.DAGsym
    DAGpack = vcat(Xin,Pin)
    DAGparam = Contractor_Params(DAGh,DAGg,DAGpack,DAGsym)
    println("generated DAG parameters")
  end

  inclusionHigh = [false for i=1:length(Xin)]
  inclusionLow = [false for i=1:length(Xin)]
  opt = [max_iter_cntr,etol,rtol,style]

  # sets dimensions
  np = length(Pin)
  nx = length(Xin)

  # initializes count variables
  NIT = 0
  NMEX = 0
  NPINMEX = 0
  NS = 0

  # creates working stack & storage objects, pushes initial box onto stack
  sol_set = []
  ind_set = []
  stack = []
  hist = []
  push!(stack,[Xin,Pin])
  disp = 0
  # calculates DAG contractor parameters

  # Check Termination
  while (~isempty(stack) && (NIT < max_iter))
    #if (display)
      println("Iteration  ||  Remaining  ||  Solution Size  ||  Indeterminat Size")
      println(NIT, "  ||  ",length(stack),"  ||  ",length(sol_set),"  ||  ",length(ind_set))
    #end
    #### Resets the tests variables and unpacks the stack ####
    Iflag = false
    Eflag = false
    PIflag = false
    PIcert = false
    eDflag = false
    cnode = pop!(stack)
    Xw = copy(cnode[1])
    X0 = copy(cnode[1])
    Pw = copy(cnode[2])
    println("Current Node (X,P): ", Xw, " ",Pw)

    #### Check for exclusion via Miranda's Test ####
    Eflag = MirandaExc(h,Xw,Pw,Eflag)
    if (Eflag)
      println("Miranda Exclusion:  ",Eflag,"  h(Z,P):  ",h(Xw,Pw))
    end
    Eflag && (NMEX += 1)
    disp = 1
    histout = []

    if (~Eflag)
      #=
      if (DAGflag)
        DAGpack = vcat(Xw,Pw)
        DAGpack = Contractor_Eval(DAGh,DAGg,DAGpack,DAGsym,DAGr,DAGparam)
        #println("typeof(DAGpack): ",typeof(DAGpack))
        #println("typeof(Xw): ",typeof(Xw))
        for i=1:length(DAGpack)
          if (isempty(DAGpack[i]))
            println("Fathom by DAG")
            Eflag = true
            break
          end
        end
        for i=1:length(Xw)
          Xw[i] = DAGpack[i]
        end
        count = 1
        for i=(length(Xw)+1):(length(Xw)+length(Pw))
          Pw[count] = DAGpack[i]
          count += 1
        end
      end

      println("ran DAG contractor")
    end

    if (~Eflag)
      #println("parametric iterations start")
      if (style == "NewtonGS")
        Xw,Xw2,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh = PI_NewtonGS(Xw,Pw,hj,h,opt,Eflag,Iflag,eDflag)
        #println("Xw:   ", Xw)
        #println("Xw2:   ", Xw2)
        #println("Eflag:   ", Eflag)
        #println("Iflag:   ", Iflag)
        #println("eDflag:   ", eDflag)
        #println("inclusionLow:   ", inclusionLow)
        #println("inclusionHigh:   ", inclusionHigh)
        if eDflag
          push!(stack,[Xw2,Pw])
        end
      end
      if (style == "KrawczykCW")
      #  println("Xw Initial: ",Xw)
        Xw,Eflag,Iflag,inclusionLow,inclusionHigh =  PI_KrawczykCW(Xw,Pw,hj,h,opt,Eflag,Iflag)
      #  println("Xw Final: ",Xw)
      end
      #println("eDflag: ", eDflag)
    #  println("Eflag: ", Eflag)
      #println("Iflag: ", Iflag)
    #  println("PIflag: ", PIflag)
    #  println("PIcert: ", PIcert)

      #println("~Iflag: ", ~Iflag)
    #  println("~PIflag: ", ~PIflag)
      #println("(~PIflag || PIcert): ", (~PIflag || PIcert))
    #  println("(~Iflag && (~PIflag || PIcert)): ", (~Iflag && (~PIflag || PIcert)))

      println("start boundary test:")
      if (Eflag==true)
        disp = 2
        NPINMEX += 1
      elseif (~Iflag && (~PIflag || PIcert))
        PIcert,PIflag,Iflag,Eflag = BoundaryTest(h,hj,X0,Xw,Pw,opt,PIcert,PIflag,
                                                 Iflag,Eflag,inclusionLow,inclusionHigh)
        println("performed boundary test")
        println("boundary test results")
        println("Eflag: ", Eflag)
        println("Iflag: ", Iflag)
        println("PIflag: ", PIflag)
        println("PIcert: ", PIcert)
      end

      if (Iflag)
        println("start re-refine")
        if (style == "NewtonGS")
          #println("trace 4")
          Xw,Xw2,Eflag,Iflag,eDflag,inclusionLow,inclusionHigh =  MC_NewtonGS(Xw,Pw,hj,h,opt,Eflag,Iflag,eDflag)
          #println("trace 4a")
          #println("Xw:   ", Xw)
          #println("Xw2:   ", Xw2)
          #println("Eflag:   ", Eflag)
          #println("Iflag:   ", Iflag)
          #println("eDflag:   ", eDflag)
          #println("inclusionLow:   ", inclusionLow)
          #println("inclusionHigh:   ", inclusionHigh)
        end
        if (style == "KrawczykCW")
          #println("trace 4")
          Xw,Eflag,Iflag,inclusionLow,inclusionHigh =  MC_KrawczykCW(Xw,Pw,hj,h,opt,Eflag,Iflag)
          #println("trace 4a")
        end
        disp = 3
        push!(sol_set,[Xw,Pw])
        NS += 1
      elseif (~Eflag)
        println("start XP bisection")
        NextNode = copy(Xw)
        NextNodeP =  copy(Pw)
        Ptemp =  copy(Pw)
        PIflag = false
        bp,CNodeX,CNodeP,NNodeX,NNodeP,Ptemp = XP_Bisection(h,hj,Xw,Pw,NextNode,NextNodeP,Ptemp,
                                                            PIflag,Iflag,Eflag,ptol,pbisect,Pstart)
        println("trace 5 - Bisection Results")
        println("bp:  ", bp)
        println("CNodeX:  ", CNodeX)
        println("CNodeP:  ", CNodeP)
        println("NNodeX:  ", NNodeX)
        println("NNodeP:  ", NNodeP)
        println("Ptemp:  ", Ptemp)
        if (bp == 2)
          disp = 4
          push!(stack,[CNodeX,CNodeP],[CNodeX,NNodeP])
        elseif (bp == 3)
          disp = 5
          push!(stack,[CNodeX,CNodeP],[CNodeX,NNodeP],[CNodeX,Ptemp])
        elseif (bp == 1)
          disp = 6
          push!(stack,[CNodeX,CNodeP],[NNodeX,CNodeP])
        else
          disp = 7
          push!(ind_set,[CNodeX,CNodeP])
        end
      end
    end
    #println("trace end loop")
    push!(hist,[NIT,cnode,disp,histout])
    NIT += 1
  end
  #println("trace end")
  return sol_set, ind_set, NIT,hist
end
=#


text:
using Plots
pyplot() # Use PyPlot as a backend (may already be the default)

#x = collect(1:7)
#y = []

#for item in x
    #push!(y, 2 - 2*item + item^2/4)
#end

#plot(x,y)
#plot!(x, y, marker = :diamond, linewidth=2)
#plot!(title = "Sample plot", leg=false)

x = [1 2 3 4 5 6]'
#y = (x-3).^2/4
y = []

for item in x
    push!(y, (item-3).^2/4)
end

plot(x,y, marker = :hex, leg=false, linewidth = 2, linecolor=:black)
plot!(title="Plot for graded quiz")
gui()    # To show plot on atom

max_stars_repo_name: aliyyousuff/MOOC

max_issues_repo_name: aliyyousuff/MOOC

max_forks_repo_name: aliyyousuff/MOOC


text:
@testset "216.combination-sum-iii.jl" begin
    @test compute_area([-3, 0, 3, 4, 0, -1, 9, 2]) == 45
    @test compute_area([0, -1, 9, 2, -3, 0, 3, 4]) == 45
    @test compute_area([-2, -2, 2, 2, -2, -2, 2, 2]) == 16
    @test compute_area([-2, -2, 2, 2, 2, 2, 4, 4]) == 0
end