text:
facts("Timing functions") do

    context("fpswhen") do
        b = Signal(false)
        t = fpswhen(b, 2)
        acc = foldp((x, y) -> x+1, 0, t)
        sleep(0.75)

        @fact queue_size() --> 0
        push!(b, true)

        step() # processing the push to b will start the fpswhen's timer
        # then we wait for two pushes from the timer, which should take ~ 1sec
        dt = @elapsed Reactive.run(2)
        push!(b, false)
        Reactive.run(1) # setting b to false should stop the timer

        sleep(0.75) # no more updates
        @fact queue_size() --> 0

        @show dt
        @fact dt --> roughly(1, atol=0.25) # mac OSX needs a lot of tolerence here
        @fact value(acc) --> 2

    end

    context("every") do
        num_steps = 4
        dt = 0.5
        # gets pushed the `time()` every `dt` seconds
        t = every(dt)

        # append the `time()` to the Float64[] array, once every dt secs
        acc = foldp(push!, Float64[], t)

        # process num_steps pushes (should take num_steps*dt secs)
        Reactive.run(num_steps)
        end_t = time() # should be equal to the last time in acc
        # close(t) to avoid `acc` getting pushed to again
        close(t)

        accval = value(acc)
        @show end_t end_t .- accval
        @fact accval[end] --> roughly(end_t, atol=0.01)

        Reactive.run_till_now()

        @fact [0.5, 0.5, 0.5] --> roughly(diff(accval), atol=0.1)

        sleep(0.75)
        # make sure the `close(t)` above actually also closed the timer
        @fact queue_size() --> 0
    end

    context("throttle") do
        gc()
        x = Signal(0; name="x")
        ydt = 0.5
        y′dt = 1.1
        y = throttle(ydt, x; name="y", leading=false)
        y′ = throttle(y′dt, x, push!, Int[], x->Int[]; name="y′", leading=false) # collect intermediate updates
        z = foldp((acc, x) -> begin
            println("z got ", x)
            acc+1
        end, 0, y)
        z′ = foldp((acc, x) -> begin
            println("z′ got ", x)
            acc+1
        end, 0, y′)
        y′prev = previous(y′)

        i = 0
        sleep_time = 0.15
        t0 = typemax(Float64)
        # push and sleep for a bit, y and y′ should only update every ydt and
        # y′dt seconds respectively
        while time() - t0 <= 2.2
            i += 1
            push!(x, i)

            Reactive.run_till_now()
            t0 == typemax(Float64) && (t0 = time()) # start timer here to match signals
            sleep(sleep_time)
        end
        dt = time() - t0
        sleep(max(ydt,y′dt) + 0.1) # sleep for the trailing-edge pushes of the throttles
        Reactive.run_till_now()

        zcount = ceil(dt / ydt) # throttle should have pushed every ydt seconds
        z′count = ceil(dt / y′dt) # throttle should have pushed every y′dt seconds

        @show i dt ydt y′dt zcount z′count value(y) value(y′) value(z) value(z′)

        @fact value(y) --> i
        @fact value(z) --> roughly(zcount, atol=1)
        @fact value(y′) --> [y′prev.value[end]+1 : i;]
        @fact length(value(y′)) --> less_than(i/(z′count-1))
        @fact value(z′) --> roughly(z′count, atol=1)

        # type safety
        s1 = Signal(3)
        s2 = Signal(rand(2,2))
        m = merge(s1, s2)
        t = throttle(1/5, m; typ=Any)
        r = rand(3,3)
        push!(s2, r)
        Reactive.run(1)
        sleep(0.5)
        # Reactive.run(1)
        Reactive.run_till_now()
        @fact value(t) --> r
    end

    context("debounce") do
        x = Signal(0)
        y = debounce(0.5, x)
        y′ = debounce(1, x, push!, Int[], x->Int[]) # collect intermediate updates
        z = foldp((acc, x) -> acc+1, 0, y)
        z′ = foldp((acc, x) -> acc+1, 0, y′)

        push!(x, 1)
        step()

        push!(x, 2)
        step()

        push!(x, 3)
        t0=time()
        step()

        @fact value(y) --> 0
        @fact value(z) --> 0
        @fact queue_size() --> 0

        sleep(0.55)

        @fact queue_size() --> 1 # y should have been pushed to by now
        step() # run the push to y
        @fact value(y) --> 3
        @fact value(z) --> 1
        @fact value(z′) --> 0
        sleep(0.5)

        @fact queue_size() --> 1 # y′ should have pushed by now
        step() # run the push to y′
        @fact value(z′) --> 1
        @fact value(y′) --> Int[1,2,3]

        push!(x, 3)
        step()

        push!(x, 2)
        step()

        push!(x, 1)
        step()
        sleep(1.1)

        @fact queue_size() --> 2 #both y and y′ should have pushed
        step()
        step()
        @fact value(y) --> 1
        @fact value(z′) --> 2
        @fact value(y′) --> Int[3,2,1]

        # type safety
        s1 = Signal(3)
        s2 = Signal(rand(2,2))
        m = merge(s1, s2)
        t = debounce(1/5, m; typ=Any)
        r = rand(3,3)
        push!(s2, r)
        Reactive.run(1)
        sleep(0.5)
        Reactive.run(1)
        @fact value(t) --> r
    end
end


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

using CausticsEngineering
using Documenter

#
# HTML docs
#
Documenter.makedocs(
    modules = [CausticsEngineering],
    format = Documenter.HTML(),
    build = "build_html",
    sitename = "Amazing Caustics Back-Engineering!",
    pages = ["Table of Contents" => "table_of_contents.md", "Index" => "index.md"],

    # root = "./doc",
    # source = ".",
    # clean = true,
    # doctest = true,
    # repo = "",
    # highlightsig = true,
    # expand = [],
)

#
# pdf docs
#
Documenter.makedocs(
    modules = [CausticsEngineering],
    format = Documenter.LaTeX(),
    build = "build_pdf",
    sitename = "Amazing Caustics Back-Engineering!",
    pages = ["Table of Contents" => "table_of_contents.md", "Index" => "index.md"],
)

ext: jl

max_stars_repo_name: Mehul-Bagaria/causticsEngineering

max_issues_repo_name: Mehul-Bagaria/causticsEngineering

max_forks_repo_name: Mehul-Bagaria/causticsEngineering


text:
struct CharacterTable{Gr,T,O} <: AbstractMatrix{T}
    group::Gr
    conjugacy_classes::Vector{O}
    inv_of::Vector{Int}
    pmap::PowerMap
    values::Matrix{T}
        # rows indexed by irreducible characters,
        # columns indexed by conjugacy classes
end

## Array Interface:
Base.size(chtbl::CharacterTable) = size(chtbl.values)
Base.getindex(chtbl::CharacterTable, i::Integer, j::Integer) = chtbl.values[i, j]

## Accessors
powermap(chtbl::CharacterTable) = chtbl.pmap
Base.parent(chtbl::CharacterTable) = chtbl.group

conjugacy_classes(chtbl::CharacterTable) = chtbl.conjugacy_classes
nconjugacy_classes(chtbl::CharacterTable) = size(chtbl.values, 2)
nirreps(chtbl::CharacterTable) = size(chtbl.values, 1)

## irreps
irreducible_characters(chtbl::CharacterTable) =
    [Character(chtbl, i) for i in 1:size(chtbl, 1)]
irreducible_characters(T::Type, chtbl::CharacterTable) =
    Character{T}[Character{T}(chtbl, i) for i in 1:size(chtbl, 1)]
irreducible_characters(G::Group, cclasses = conjugacy_classes(G)) =
    irreducible_characters(Rational{Int}, G, cclasses)

function irreducible_characters(
    R::Type{<:Rational},
    G::Group,
    cclasses=conjugacy_classes(G)
)
    return irreducible_characters(CharacterTable(R, G, cclasses))
end

function trivial_character(chtbl::CharacterTable)
	# return Character(chtbl, findfirst(r->all(isone, r), eachrow(chtbl)))
	# can't use findfirst(f, eachrow(...)) on julia-1.6
	for i in 1:size(chtbl, 1)
		all(isone, @view(chtbl[i, :])) && return Character(chtbl, i)
	end
	# never hit, to keep compiler happy
	return Character(chtbl, 0)
end

## construcing tables

function CharacterTable(
    Fp::Type{<:FiniteFields.GF},
    G::Group,
    cclasses = conjugacy_classes(G),
)
    Ns = [CMMatrix(cclasses, i) for i in 1:length(cclasses)]
    esd = common_esd(Ns, Fp)
    @assert isdiag(esd)

    tbl = CharacterTable(G, cclasses, _inv_of(cclasses), PowerMap(cclasses), esd.basis)

    tbl = normalize!(tbl)
    return tbl
end

function CharacterTable(R::Type{<:Rational}, G::Group, cclasses=conjugacy_classes(G))
    Fp = FiniteFields.GF{dixon_prime(cclasses)}
    tblFp = CharacterTable(Fp, G, cclasses)
    return complex_character_table(R, tblFp)
end

function complex_character_table(
    R::Type{<:Rational},
    tblFp::CharacterTable{<:Group, <:FiniteFields.GF},
)
    charsFp = irreducible_characters(tblFp)
    mult_c = _multiplicities(charsFp)

    e = size(mult_c, 3) # the exponent

    C = Cyclotomics.Cyclotomic{R,Cyclotomics.SparseVector{R,Int}}
    values = Matrix{C}(undef, size(tblFp))

    Es = [E(e, k) for k in 0:e - 1]
    Threads.@threads for j in 1:size(tblFp, 2) # conjugacy_classes
        for i in 1:size(tblFp, 1) # characters
            # reduced_embedding may prevent overflow sometimes
            values[i, j] = Cyclotomics.reduced_embedding(
                sum(mult_c[i, j, k + 1] * Es[k + 1] for k in 0:e - 1)
            )
        end
    end

    return CharacterTable(
        parent(tblFp),
        conjugacy_classes(tblFp),
        tblFp.inv_of,
        powermap(tblFp),
        values,
    )
end

function _inv_of(cc::AbstractVector{<:AbstractOrbit})
    inv_of = zeros(Int, size(cc))
    for (i, c) in enumerate(cc)
        g = inv(first(c))
        inv_of[i] = something(findfirst(k -> g in k, cc), 0)
    end
    any(iszero, inv_of) && throw(
        ArgumentError(
            "Could not find the conjugacy class for inverse of $(first(cc[findfirst(iszero, inv_of)])).",
        ),
    )
    return inv_of
end

function normalize!(chtbl::CharacterTable{<:Group,<:FiniteFields.GF})
    id = one(parent(chtbl))
    for (i, χ) in enumerate(irreducible_characters(chtbl))
        k = χ(id)
        if !isone(k)
            chtbl.values[i, :] .*= inv(k)
        end
        # ⟨χ, χ⟩ = 1/d²

        deg = sqrt(inv(dot(χ, χ)))
        @debug "normalizing with" dot(χ, χ) χ(id) χ

        # normalizing χ
        chtbl.values[i, :] .*= deg
    end
    return chtbl
end

## "fancy" io

function Base.show(io::IO, ::MIME"text/plain", chtbl::CharacterTable)
    println(io, "Character table of ", parent(chtbl), " over $(eltype(chtbl))")

	if !(get(io, :limit, false)::Bool)
		screenheight = screenwidth = typemax(Int)
	else
		sz = displaysize(io)::Tuple{Int,Int}
		screenheight, screenwidth = sz[1] - 4, sz[2]
	end

    nirr = nirreps(chtbl)
    nccl = nconjugacy_classes(chtbl)
    pre_width = 1 + length(string(nirr))
    sep = " │ "

    if nirr > screenheight - 1
        rows_to_print = [1:screenheight - 2; nirr]
    else
        rows_to_print = [1:nirr;]
    end

    maxpossiblecols = div(screenwidth - pre_width - 3, 3)
    if nccl > maxpossiblecols
        cols_to_print = [1:maxpossiblecols - 1; nccl]
    else
	    cols_to_print = [1:nccl;]
    end

    A = Base.alignment(
        io,
        chtbl.values,
        rows_to_print,
        cols_to_print,
        screenwidth,
        screenwidth,
        2,
    )

	hellipsis = nccl > length(A) ? " …" : ""

    print(io, " "^(pre_width), "   ")
    Base.print_matrix_row(
        io,
        reshape(cols_to_print, (1, length(cols_to_print))),
        A,
        1,
        cols_to_print,
        " ",
    )

    println(io, hellipsis)
    println(io, "─"^pre_width, "─┬─", "─"^(sum(sum, A) + length(A) - 1), hellipsis)

    if nirr > screenheight - 1
        for i in 1:screenheight - 2
            print(io, rpad("χ$(FiniteFields.subscriptify(i))", pre_width), sep)
            Base.print_matrix_row(io, chtbl.values, A, i, cols_to_print, " ")
            println(io, hellipsis)
        end
        print(io, " ⋮", " "^(pre_width - 2), sep)
		Base.print_matrix_vdots(io, "⋮", A, "  ", 2, 1, false)
		println(io)

        print(io, rpad("χ$(FiniteFields.subscriptify(nirr))", pre_width), sep)
        Base.print_matrix_row(io, chtbl.values, A, nirr, cols_to_print, " ")
		print(io, hellipsis)
    else
        for i in 1:nirr
            print(io, rpad("χ$(FiniteFields.subscriptify(i))", pre_width), sep)
            Base.print_matrix_row(io, chtbl.values, A, i, cols_to_print, " ")
            i != nirr && println(io, hellipsis)
        end
    end
end


text:
@testset "Foo test" begin
    v = NewPackage.foo(10,5)
    @test v[1] == 10
    @test v[2] == 5
    @test eltype(v) == Int
    v = NewPackage.foo(10.0, 5)
    @test v[1] == 10
    @test v[2] == 5
    @test eltype(v) == Float64
end


text:
module HullWrapper

using JavaCall

const Collection = @jimport java.util.Collection
const ConvexHull = @jimport org.hipparchus.geometry.hull.ConvexHull
const ConvexHullGenerator = @jimport org.hipparchus.geometry.hull.ConvexHullGenerator
const Point = @jimport org.hipparchus.geometry.Point
const Region = @jimport org.hipparchus.geometry.partitioning.Region

export ConvexHull
export ConvexHullGenerator
export create_region
export generate
export get_vertices

include("convex_hull.jl")
include("convex_hull_generator.jl")

end


text:
export ArxivLinks

struct ArxivLinks
end

function process_source!(A::ArxivLinks,
                         source::AbstractString,
                         fileinfo::FileInfo;
                         args...)
  link_re = r"(arxiv|cond-mat|quant-ph|math|math-ph|physics)[/:]\W*?([\d\.]+)"i
  res = ""
  pos = 1
  match = false
  for m in eachmatch(link_re,source)
    match = true
    res *= source[pos:m.offset-1]
    prefix = m.captures[1]
    number = m.captures[2]
    if lowercase(prefix) == "arxiv"
      res *= "arxiv:[$number](https://arxiv.org/abs/$number)"
    else
      nbprefix = replace(prefix,"-" => "‑") #non-breaking hypen
      res *= "<span>$nbprefix/[$number](https://arxiv.org/abs/$prefix/$number)</span>"
    end
    pos = m.offset+length(m.match)
  end
  if !match
    return source
  else
    res *= source[pos:end]
  end
  return res
end


text:
using BrokenRecord: BrokenRecord, HTTP, playback
using Dates: Day, UTC, now
using RegistryCI: TagBot
using SimpleMock: Mock, called_with, mock
using Test: @test, @testset, @test_logs

const TB = TagBot
const GH = TB.GH

TB.AUTH[] = GH.OAuth2(get(ENV, "GITHUB_TOKEN", "abcdef"))
BrokenRecord.configure!(;
    path=joinpath(@__DIR__, "cassettes"),
    ignore_headers=["Authorization"],
)

@testset "is_merged_pull_request" begin
    @test !TB.is_merged_pull_request(Dict())
    @test !TB.is_merged_pull_request(Dict("pull_request" => Dict("merged" => false)))
    @test TB.is_merged_pull_request(Dict("pull_request" => Dict("merged" => true)))
end

@testset "is_cron" begin
    withenv(() -> @test(!TB.is_cron(())), "GITHUB_EVENT_NAME" => nothing)
    withenv(() -> @test(!TB.is_cron(())), "GITHUB_EVENT_NAME" => "pull_request")
    withenv(() -> @test(TB.is_cron(())), "GITHUB_EVENT_NAME" => "schedule")
end

@testset "repo_and_version_of_pull_request" begin
    body(url) = """
        - Repository: $url
        - Version: v1.2.3
        """
    github = body("https://github.com/Foo/Bar")
    @test TB.repo_and_version_of_pull_request_body(github) == ("Foo/Bar", "v1.2.3")
    ssh = body("git@github.com:Foo/Bar.git")
    @test TB.repo_and_version_of_pull_request_body(ssh) == ("Foo/Bar", "v1.2.3")
    gitlab = body("https://gitlab.com/Foo/Bar")
    @test TB.repo_and_version_of_pull_request_body(gitlab) == (nothing, "v1.2.3")
end

@testset "clone_repo" begin
    mock(run, mktempdir => Mock("a")) do run, _mktempdir
        @test TB.clone_repo("A") == "a"
        @test called_with(run, `git clone --depth=1 https://github.com/A a`)
    end
end

@testset "is_tagbot_enabled" begin
    mock(TB.clone_repo => repo -> joinpath(@__DIR__, "repos", repo)) do _clone
        @test !TB.is_tagbot_enabled("no_actions")
        @test !TB.is_tagbot_enabled("no_tagbot")
        @test TB.is_tagbot_enabled("yes_tagbot")
    end
end

@testset "get_repo_notification_issue" begin
    repo = "christopher-dG/TestRepo"
    playback("get_repo_notification_issue.bson") do
        @test_logs (:info, "Creating new notification issue") begin
            issue = TB.get_repo_notification_issue(repo)
            @test issue.number == 4
        end
        @test_logs (:info, "Found existing notification issue") begin
            issue = TB.get_repo_notification_issue(repo)
            @test issue.number == 4
        end
    end
end

@testset "notification_body" begin
    base = "Triggering TagBot for merged registry pull request"
    @test TB.notification_body(Dict()) == base
    event = Dict("pull_request" => Dict("html_url" => "foo"))
    @test TB.notification_body(event) == "$base: foo"
end

@testset "notify" begin
    playback("notify.bson") do
        comment = TB.notify("christopher-dG/TestRepo", 4, "test notification")
        @test comment.body == "test notification"
    end
end

@testset "collect_pulls" begin
    PR = GH.PullRequest
    prs = [
        [PR(), PR(; merged_at=now(UTC))],
        [PR(; merged_at=now(UTC) - Day(2)), PR(; merged_at=now(UTC) - Day(4))],
    ]
    pages = [Dict("next" => "abc"), Dict()]
    pulls = mock(GH.pull_requests => Mock(collect(zip(prs, pages)))) do _prs
        TB.collect_pulls("JuliaRegistries/General")
    end
    @test pulls == [prs[1][2], prs[2][1]]
end

@testset "tag_exists" begin
    playback("tag_exists.bson") do
        @test TB.tag_exists("JuliaRegistries/RegistryCI.jl", "v0.1.0")
        @test !TB.tag_exists("JuliaRegistries/RegistryCI.jl", "v0.0.0")
    end
end

@testset "maybe_notify" begin
    @test_logs match_mode=:any (:info, r"not enabled") begin
        mock(TB.is_tagbot_enabled => Mock(false)) do _is_tagbot_enabled
            TB.maybe_notify((), "repo", "v")
        end
    end
    @test_logs match_mode=:any (:info, r"already exists") begin
        mock(TB.is_tagbot_enabled => Mock(true), TB.tag_exists => Mock(true)) do ite, te
            TB.maybe_notify((), "repo", "v"; check_tag=true)
        end
    end
    mock(
        TB.is_tagbot_enabled => Mock(true),
        TB.get_repo_notification_issue => Mock(1),
        TB.notify,
    ) do is_enabled, get_issue, notify
        TB.maybe_notify(Dict(), "repo", "v")
        @test called_with(is_enabled, "repo")
        @test called_with(get_issue, "repo")
        msg = "Triggering TagBot for merged registry pull request"
        @test called_with(notify, "repo", 1, msg)
    end
end


text:
using BinaryProvider

const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libPixhawkData"], :PixhawkData),
]

if Sys.iswindows()
    using Pkg
    Pkg.add(PackageSpec(name="Cxx", rev="gn-patch-crash"))
end

# TODO create a repo
bin_prefix = "http://www.github.com/akshayhiregoudar/PixhawkDataBuilder/releases/download/v1.0"

# TODO check version number
download_info = Dict(
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/PixhawkData.v1.0.x86_64-linux-gnu-gcc7.tar.gz", "57f0d9f9e89f05db78ec41c6cb5878b8e2290e9868175f556a34bcfc61df0578"),
    #Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/PixhawkData.v1.0.x86_64-linux-gnu-gcc8.tar.gz", "7caaa648e059c354892a5503b30bf3b2d368469b4295f54564f8a212fa3e2893"),
    #Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/PixhawkData.v1.0.x86_64-w64-mingw32-gcc7.tar.gz", "782f52898e81c51abd8930c34184a6f99d78e4a49e7a6c753b9eb0741e8bcdfc"),
    #Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/PixhawkData.v1.0.x86_64-w64-mingw32-gcc8.tar.gz", "283473593d120f968366f85da52bd4341ac94715c49476adf9bdd1688b445411"),
    #MacOS(:x86_64) => ("$bin_prefix/PixhawkData.v1.0.x86_64-apple-darwin14.tar.gz", "de9236acab1975a06f218d2939f6c319fb0bc67546e89cdc5c85ef90be509cb4"),
)

unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
"""
    log_partition(rbm)

Log-partition of `rbm`, computed by extensive enumeration of visible states
(except for particular cases such as Gaussian-Gaussian RBM).
This is exponentially slow for large machines.

If your RBM has a smaller hidden layer, mirroring the layers of the `rbm` first
(see [`mirror`](@ref)).
"""
function log_partition(rbm::RBM)
    v = ChainRulesCore.ignore_derivatives() do
        collect_states(rbm.visible)
    end
    return logsumexp(-free_energy(rbm, v))
end

"""
    log_likelihood(rbm, v)

Log-likelihood of `v` under `rbm`, with the partition function compued by
extensive enumeration. For discrete layers, this is exponentially slow for large machines.
"""
function log_likelihood(rbm::RBM, v::AbstractArray)
    logZ = log_partition(rbm)
    F = free_energy(rbm, v)
    return -F .- logZ
end

function iterate_states(layer::Binary)
    itr = generate_sequences(length(layer), false:true)
    return map(x -> reshape(x, size(layer)..., 1), itr)
end

function iterate_states(layer::Spin)
    itr = generate_sequences(length(layer), (-Int8(1), Int8(1)))
    return map(x -> reshape(x, size(layer)..., 1), itr)
end

function iterate_states(::Potts)
    error("not implemented")
end

"""
    collect_states(layer)

Returns an array of all states of `layer`.
Only defined for discrete layers.

!!! warning
    Use only for small layers.
    For large layers, the exponential number of states will not fit in memory.
"""
function collect_states(layer::Union{Binary, Spin, Potts})
    reduce(iterate_states(layer)) do x, y
        cat(x, y; dims = ndims(layer) + 1)
    end
end


text:
export RBMSplit

struct RBMSplit{VT,MT} <: MatrixNeuralNetwork
    ar::VT
    ac::VT
    b::VT
    Wr::MT
    Wc::MT
end
@treelike RBMSplit

"""
    RBMSplit([T=Complex{STD_REAL_PREC}], N, α, [initW, initb, inita])

Constructs a Restricted Bolzmann Machine to encode a vectorised density matrix,
with weights of type `T` (Defaults to ComplexF32), `2N` input neurons,
2N⋅α hidden neurons.
This network does not ensure positive-definitness of the density matrix.

`N` must match the size of the lattice.

The initial parameters of the neurons are initialized with a rescaled normal
distribution of width 0.01 for the coupling matrix and 0.05 for the local
biases. The default initializers can be overriden by specifying

initW=(dims...)->rescaled_normal(T, 0.01, dims...),
initb=(dims...)->rescaled_normal(T, 0.05, dims...),
initb=(dims...)->rescaled_normal(T, 0.01, dims...),

Refs:
    https://arxiv.org/abs/1902.07006
"""
RBMSplit(in::Int, α::Number, args...) = RBMSplit(ComplexF32, in, α, args...)
RBMSplit(T::Type, in, α,
         initW=(dims...)->rescaled_normal(T, 0.01, dims...),
         initb=(dims...)->rescaled_normal(T, 0.05, dims...),
         inita=(dims...)->rescaled_normal(T, 0.01, dims...)) =
    RBMSplit(inita(in), inita(in),
             initb(convert(Int, α*in)),
             initW(convert(Int, α*in), in), initW(convert(Int, α*in), in))

input_type(net::RBMSplit)  = real(eltype(net.ar))
weight_type(net::RBMSplit) = out_type(net)
out_type(net::RBMSplit)    = eltype(net.Wr)
input_shape(net::RBMSplit) = (length(net.ar), length(net.ac))
random_input_state(net::RBMSplit) =
    (eltype(net.ar).([rand(0:1) for i=1:length(net.ar)]), eltype(net.ar).([rand(0:1) for i=1:length(net.ar)]))
is_analytic(net::RBMSplit) = true

(net::RBMSplit)(σ::State) = net(config(σ)...)
(net::RBMSplit)(σr, σc)   = transpose(net.ar)*σr .+ transpose(net.ac)*σc .+ sum_autobatch(logℒ.(net.b .+
                                                        net.Wr*σr .+ net.Wc*σc))

function Base.show(io::IO, m::RBMSplit)
    print(io, "RBMSplit($(eltype(m.ar)), n=$(length(m.ar)), α=$(length(m.b)/length(m.ar)))")
end
Base.show(io::IO, ::MIME"text/plain", m::RBMSplit) = print(
"RBMSplit($(eltype(m.ar)), n=$(length(m.ar)), α=$(length(m.b)/length(m.ar)))")

# Cached version
mutable struct RBMSplitCache{VT,VS,VST} <: NNCache{RBMSplit}
    θ::VT
    θ_tmp::VT
    logℒθ::VT
    ∂logℒθ::VT

    # complex sigmas
    res::VS #batch
    res_tmp::VST #batch

    # states
    σr::VT
    σc::VT

    valid::Bool # = false
end

cache(net::RBMSplit) =
    RBMSplitCache(similar(net.b),
                  similar(net.b),
                  similar(net.b),
                  similar(net.b),
                  similar(net.b),
                  similar(net.b),
                  similar(net.b, length(net.ar)),
                  similar(net.b, length(net.ar)),
                  false)

(net::RBMSplit)(c::RBMSplitCache, σ::State) = net(c, config(σ))
(net::RBMSplit)(c::RBMSplitCache, (σr, σc)::Tuple{AbstractArray,AbstractArray}) = net(c, σr, σc)
function (net::RBMSplit)(c::RBMSplitCache, σr_r, σc_r)
    θ = c.θ
    θ_tmp = c.θ_tmp
    logℒθ = c.logℒθ
    T = eltype(θ)

    # copy the states to complex valued states for the computations.
    σr = c.σr; copyto!(σr, σr_r)
    σc = c.σc; copyto!(σc, σc_r)

    #θ .= net.b .+
    #        net.Wr*σr .+
    #            net.Wc*σc
    mul!(θ, net.Wr, σr)
    mul!(θ_tmp, net.Wc, σc)
    θ .+= net.b .+ θ_tmp

    logℒθ .= logℒ.(θ)
    lnψ = dot(σr,net.ar) + dot(σc,net.ac) + sum(logℒθ)
    return lnψ
end

function logψ_and_∇logψ!(∇logψ, net::RBMSplit, c::RBMSplitCache, σr_r, σc_r)
    θ      = c.θ
    θ_tmp  = c.θ_tmp
    logℒθ  = c.logℒθ
    ∂logℒθ = c.∂logℒθ
    T      = eltype(θ)

    # copy the states to complex valued states for the computations.
    σr = c.σr; copyto!(σr, σr_r)
    σc = c.σc; copyto!(σc, σc_r)

    #θ .= net.b .+
    #        net.Wr*σr .+
    #            net.Wc*σc
    mul!(θ, net.Wr, σr)
    mul!(θ_tmp, net.Wc, σc)
    θ .+= net.b .+ θ_tmp

    logℒθ  .= logℒ.(θ)
    ∂logℒθ .= ∂logℒ.(θ)

    ∇logψ.ar .= σr
    ∇logψ.ac .= σc
    ∇logψ.b  .= ∂logℒθ
    ∇logψ.Wr .= ∂logℒθ .* transpose(σr)
    ∇logψ.Wc .= ∂logℒθ .* transpose(σc)

    lnψ = dot(σr,net.ar) + dot(σc,net.ac) + sum(logℒθ)
    return lnψ
end


text:
function updatelungs!(system::CVSystem,n::Int64)
    # 1st arterial compartment
    system.lungs.Va[n+2,1] = (system.lungs.Va[n+1,1] + system.solverparams.h*
        (system.heart.rv.Q[n+1] - system.lungs.Qa[n+1,1]));
    system.lungs.Qa[n+2,1] = ((1 - system.lungs.Ra[1]/system.lungs.La[1]*
        system.solverparams.h)*system.lungs.Qa[n+1,1]+
        system.solverparams.h/system.lungs.La[1]*(system.lungs.Pa[n+1,1]-
        system.lungs.Pa[n+1,2]));
    system.lungs.Pa[n+2,1] = (1/system.lungs.Ca[1]*(system.lungs.Va[n+2,1]-
        system.lungs.V0[1]));

    # 2nd arterial compartment V, P
    system.lungs.Va[n+2,2] = (system.lungs.Va[n+1,2] + system.solverparams.h*
        (system.lungs.Qa[n+1,1] - system.lungs.Qa[n+1,2]));
    system.lungs.Pa[n+2,2] = (1/system.lungs.Ca[2]*(system.lungs.Va[n+2,2]-
        system.lungs.V0[2]));

    # 3rd arterial compartment V, P
    system.lungs.Va[n+2,3] = (system.lungs.Va[n+1,3] + system.solverparams.h*
        (system.lungs.Qa[n+1,2] - system.lungs.Qa[n+1,3]));
    system.lungs.Pa[n+2,3] = (1/system.lungs.Ca[3]*(system.lungs.Va[n+2,3]-
        system.lungs.V0[3]));

    # 1st venous compartment V, P
    system.lungs.Vv[n+2,1] = (system.lungs.Vv[n+1,1] + system.solverparams.h*
        (system.lungs.Qa[n+1,3] - system.lungs.Qv[n+1,1]));
    system.lungs.Pv[n+2,1] = (1/system.lungs.Cv[1]*(system.lungs.Vv[n+2,1]-
        system.lungs.V0[4]));

    # 2nd venous compartment
    system.lungs.Vv[n+2,2] = (system.lungs.Vv[n+1,2] + system.solverparams.h*
        (system.lungs.Qv[n+1,1] - system.lungs.Qv[n+1,2]));
    system.lungs.Qv[n+2,2] = ((1 - system.lungs.Rv[2]/system.lungs.Lv[1]*
        system.solverparams.h)*system.lungs.Qv[n+1,2]+
        system.solverparams.h/system.lungs.Lv[1]*(system.lungs.Pv[n+1,2]-
        system.heart.la.P[n+1]));
    system.lungs.Pv[n+2,2] = (1/system.lungs.Cv[2]*(system.lungs.Vv[n+2,2]-
        system.lungs.V0[5]));

    # remaining Q's
    system.lungs.Qa[n+2,2] = ((system.lungs.Pa[n+2,2] - system.lungs.Pa[n+2,3])/
        system.lungs.Ra[2]);
    system.lungs.Qa[n+2,3] = ((system.lungs.Pa[n+2,3] - system.lungs.Pv[n+2,1])/
        system.lungs.Ra[3]);
    system.lungs.Qv[n+2,1] = ((system.lungs.Pv[n+2,1] - system.lungs.Pv[n+2,2])/
        system.lungs.Rv[1]);

    # root pulmonary artery P
    system.lungs.Pp[n+2] = (system.lungs.Rp*system.heart.rv.Q[n+2] +
        system.lungs.Pa[n+2,1]);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
@testset "Template" begin
    @testset "Template constructor" begin
        @testset "user" begin
            mock(PT.default_user => () -> "") do _du
                @test_throws ArgumentError Template()
                @test isempty(Template(; disable_defaults=[Git]).user)
            end
            mock(PT.default_user => () -> "username") do _du
                @test Template().user == "username"
            end
        end

        @testset "authors" begin
            @test tpl(; authors=["a"]).authors == ["a"]
            @test tpl(; authors="a").authors == ["a"]
            @test tpl(; authors="a,b").authors == ["a", "b"]
            @test tpl(; authors="a, b").authors == ["a", "b"]
        end

        @testset "host" begin
            @test tpl(; host="https://foo.com").host == "foo.com"
        end

        @testset "dir" begin
            @test tpl(; dir="/foo/bar").dir == joinpath(path_separator, "foo", "bar")
            @test tpl(; dir="foo").dir == abspath("foo")
            @test tpl(; dir="~/foo").dir == abspath(expanduser("~/foo"))
        end

        @testset "plugins / disabled_defaults" begin
            function test_plugins(plugins, expected, disabled=DataType[])
                t = tpl(; plugins=plugins, disable_defaults=disabled)
                @test all(map(==, sort(t.plugins; by=string), sort(expected; by=string)))
            end

            defaults = PT.default_plugins()
            test_plugins([], defaults)
            test_plugins([Citation()], union(defaults, [Citation()]))
            # Overriding a default plugin.
            default_g = defaults[findfirst(p -> p isa Git, defaults)]
            g = Git(; ssh=true)
            test_plugins([g], union(setdiff(defaults, [default_g]), [g]))
            # Disabling a default plugin.
            test_plugins([], setdiff(defaults, [default_g]), [Git])
        end

        @testset "Unsupported keywords warning" begin
            @test_logs tpl()
            @test_logs (:warn, r"Unrecognized keywords were supplied") tpl(; x=1, y=2)
        end
    end

    @testset "hasplugin" begin
        t = tpl(; plugins=[TravisCI(), Documenter{TravisCI}()])
        @test PT.hasplugin(t, typeof(first(PT.default_plugins())))
        @test PT.hasplugin(t, Documenter)
        @test PT.hasplugin(t, PT.is_ci)
        @test PT.hasplugin(t, _ -> true)
        @test !PT.hasplugin(t, _ -> false)
        @test !PT.hasplugin(t, Citation)
    end

    @testset "validate" begin
        foreach((GitHubActions, TravisCI, GitLabCI)) do T
            @test_throws ArgumentError tpl(; plugins=[Documenter{T}()])
        end
        mock(LibGit2.getconfig => (_k, _d) -> "") do _gc
            @test_throws ArgumentError tpl(; plugins=[Git()])
        end
    end
end

@testset "Package generation errors" begin
    mktempdir() do dir
        t = tpl(; dir=dirname(dir))
        @test_throws ArgumentError t(basename(dir))
    end

    mktemp() do f, _io
        t = tpl(; dir=dirname(f))
        @test_throws ArgumentError t(basename(f))
    end

    t = tpl()
    pkg = pkgname()
    mock(LibGit2.init => dir -> (@test isdir(dir); error())) do _init
        @test_throws ErrorException @suppress t(pkg)
    end
    @test !isdir(joinpath(t.dir, pkg))
end


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

@doc raw"""ContainerPort represents a network port in a single container.

    IoK8sApiCoreV1ContainerPort(;
        containerPort=nothing,
        hostIP=nothing,
        hostPort=nothing,
        name=nothing,
        protocol=nothing,
    )

    - containerPort::Int32 : Number of port to expose on the pod's IP address. This must be a valid port number, 0 < x < 65536.
    - hostIP::String : What host IP to bind the external port to.
    - hostPort::Int32 : Number of port to expose on the host. If specified, this must be a valid port number, 0 < x < 65536. If HostNetwork is specified, this must match ContainerPort. Most containers do not need this.
    - name::String : If specified, this must be an IANA_SVC_NAME and unique within the pod. Each named port in a pod must have a unique name. Name for the port that can be referred to by services.
    - protocol::String : Protocol for port. Must be UDP, TCP, or SCTP. Defaults to \"TCP\".
"""
mutable struct IoK8sApiCoreV1ContainerPort <: SwaggerModel
    containerPort::Any # spec type: Union{ Nothing, Int32 } # spec name: containerPort
    hostIP::Any # spec type: Union{ Nothing, String } # spec name: hostIP
    hostPort::Any # spec type: Union{ Nothing, Int32 } # spec name: hostPort
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiCoreV1ContainerPort(;containerPort=nothing, hostIP=nothing, hostPort=nothing, name=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("containerPort"), containerPort)
        setfield!(o, Symbol("containerPort"), containerPort)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("hostIP"), hostIP)
        setfield!(o, Symbol("hostIP"), hostIP)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("hostPort"), hostPort)
        setfield!(o, Symbol("hostPort"), hostPort)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1ContainerPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiCoreV1ContainerPort

const _property_map_IoK8sApiCoreV1ContainerPort = Dict{Symbol,Symbol}(Symbol("containerPort")=>Symbol("containerPort"), Symbol("hostIP")=>Symbol("hostIP"), Symbol("hostPort")=>Symbol("hostPort"), Symbol("name")=>Symbol("name"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiCoreV1ContainerPort = Dict{Symbol,String}(Symbol("containerPort")=>"Int32", Symbol("hostIP")=>"String", Symbol("hostPort")=>"Int32", Symbol("name")=>"String", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ContainerPort }) = collect(keys(_property_map_IoK8sApiCoreV1ContainerPort))
Swagger.property_type(::Type{ IoK8sApiCoreV1ContainerPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ContainerPort[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ContainerPort }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ContainerPort[property_name]

function check_required(o::IoK8sApiCoreV1ContainerPort)
    (getproperty(o, Symbol("containerPort")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ContainerPort }, name::Symbol, val)
end


text:
# Export ObjectHandle API
export ObjectHandle,
       readmeta,
       seek, seekstart, skip, start, iostream, position, read, readuntil, eof,
       endianness, is64bit, isrelocatable, isexecutable, islibrary,
       mangle_section_name, mangle_symbol_name, handle, header, format_string,
       section_header_offset, section_header_size, section_header_type,
       segment_header_offset, segment_header_size, segment_header_type,
       symtab_entry_offset, symtab_entry_size, symtab_entry_type,
       path, find_library, find_libraries

export MagicMismatch

# Import methods for extension
import Base: seek, seekstart, skip, start, position, read, readuntil, eof,
             getindex
import StructIO: unpack

"""
    ObjectHandle

The basic type that provides access to object files.  The list of available API
operations is given below, with methods that subclasses must implement marked
in emphasis.  Note tha "must implement" is a bit of a misnomer, if an Object
file does not have need of a certain piece of this API (e.g. `COFF` files have
no concept of `Segment`s), leaving that chunk of the API undefined will simply
cause errors if a user attempts to use methods that use that part of the API
(in the example above, an error will be thrown if the user calls `Segments(oh)`
where `oh <: COFFHandle`).

### Creation
  - *readmeta()*

### IOStream-like operations:
  - seek()
  - seekstart()
  - skip()
  - start()
  - iostream()
  - position()
  - read()
  - readuntil()
  - eof()
  - unpack()

### Format-specific properties
  - *endianness()*
  - *is64bit()*
  - *isrelocatable()*
  - *isexecutable()*
  - *islibrary()*
  - *mangle_section_name()*
  - *mangle_symbol_name()*
  - handle()
  - *header()*
  - *format_string()*

### Section properties
  - *section_header_offset()*
  - *section_header_size()*
  - *section_header_type()*

### Segment properties
  - *segment_header_offset()*
  - *segment_header_size()*
  - *segment_header_type()*

### Symbol properties
  - *symtab_entry_offset()*
  - *symtab_entry_size()*
  - *symtab_entry_type()*

### Misc
  - *path()*
  - show()
  - find_library()
  - find_libraries()
"""
abstract type ObjectHandle end

"""
    MagicMismatch

This is an error type used to denote that `readmeta()` was called on a file
that does not contain the proper magic at the beginning for the type of object
file that was attempting to be loaded.
"""
struct MagicMismatch
    message
end

"""
    ObjTypes::Vector{Type}

`ObjTypes` contains the registry of file formats that will be used to try and
open a object file, (e.g. `ELF`, `MachO`, etc...).  This vector is initialized
at `__init__()` time, and used within `readmeta()`.
"""
const ObjTypes = Type[]

"""
    readmeta(io::IO, ::ObjectHandle)

Read an Object File out from an `IOStream`.  This is the first method you
should call in order to manipulate object files.
"""
function readmeta{T<:ObjectHandle}(io::IO, ::Type{T})
    # Implementing packages such as MachO.jl must explicitly override this
    error("$T must implement readmeta")
end

"""
    readmeta(io::IO)

Read an Object File out from an `IOStream`, guessing at the type of object
within the stream by calling `readmeta(io, T)` for each `T` within `ObjTypes`,
and returning the first that does not throw a `MagicMismatch`.
"""
function readmeta(io::IO)
    pos = position(io)
    for T in ObjTypes
        seek(io,pos)
        try
            return readmeta(io, T)
        catch e
            if !isa(e,MagicMismatch)
                rethrow(e)
            end
        end
    end

    # If the file didn't match anything, error out
    msg = strip("""
    Object file is not any of $(join(ObjTypes, ", "))!
    To force one object file format use readmeta(io, T).
    """)
    error(replace(msg, "\n", " "))
end

"""
    readmeta(path::AbstractStriong)

Read an Object File out from a file `path`, guessing at the type of object
within the stream by calling `readmeta(io, T)` for each `T` within `ObjTypes`,
and returning the first that does not throw a `MagicMismatch`.
"""
function readmeta(file::AbstractString)
    return readmeta(open(file, "r"))
end

## IOStream-like operations
"""
    iostream(oh::ObjectHandle)

Returns the `IOStream` backing the `ObjectHandle`
"""
@mustimplement iostream(oh::ObjectHandle)

"""
    start(oh::ObjectHandle)

Returns the offset within the underlying `IOStream` at which this
`ObjectHandle` is located.
"""
@mustimplement start(oh::ObjectHandle)

unpack(oh::ObjectHandle, T) = unpack(iostream(oh), T, endianness(oh))
position(oh::ObjectHandle) = position(iostream(oh)) - start(oh)
seek(oh::ObjectHandle, pos::Integer) = seek(iostream(oh), start(oh) + pos)
seekstart(oh::ObjectHandle) = seek(oh, 0)

# Generate a bunch of wrappers for ObjectHandles that just call iostream()
# intelligently, so we don't have to continually unwrap
for f in [:skip, :seekstart, :eof, :read, :readuntil, :readbytes, :write]
    @eval $(f)(oh::H, args...) where {H<:ObjectHandle} = $(f)(iostream(oh), args...)
end

"""
    endianness(oh::ObjectHandle)

Returns the endianness of the given `ObjectHandle` (e.g. `:LittleEndian`)
"""
@mustimplement endianness(oh::ObjectHandle)

"""
    is64bit(oh::ObjectHandle)

Returns `true` if the given `ObjectHandle` represents a 64-bit object
"""
@mustimplement is64bit(oh::ObjectHandle)

"""
    isrelocatable(oh::ObjectHandle)

Returns `true` if the given `ObjectHandle` represents a relocatable object
file, e.g. an `.o` file as generated by `gcc -c`
"""
@mustimplement isrelocatable(oh::ObjectHandle)

"""
    isexecutable(oh::ObjectHandle)

Returns `true` if the given `ObjectHandle` represents an executable object
"""
@mustimplement isexecutable(oh::ObjectHandle)

"""
    islibrary(oh::ObjectHandle)

Returns `true` if the given `ObjectHandle` represents a shared library
"""
@mustimplement islibrary(oh::ObjectHandle)

"""
    mangle_section_name(oh::ObjectHandle, name::AbstractString)

Turn a section `name` into the object-format specific naming convention, e.g.
returning `".bss"` for `ELF`/`COFF` files, and `"__bss"` for `MachO` files
"""
@mustimplement mangle_section_name(oh::ObjectHandle, name::AbstractString)

"""
    mangle_symbol_name(oh::ObjectHandle, name::AbstractString)

Mangle a symbol name using the object-format specific naming convention, e.g.
prefixing `"_"` for MachO files.
"""
@mustimplement mangle_symbol_name(oh::ObjectHandle, name::AbstractString)

# Silly handle() fallthrough for an ObjectHandle itself
handle(oh::ObjectHandle) = oh

"""
    header(oh::ObjectHandle)

Return the `ObjectHandle`'s header object, whatever that may be for this
particular object file format.
"""
@mustimplement header(oh::ObjectHandle)

"""
    format_string(::Type{H}) where {H <: ObjectHandle}

Return the string name of the given `ObjectHandle`, examples are "ELF",
"MachO", "COFF", etc...
"""
@mustimplement format_string(oh::Type)
format_string(::H) where {H <: ObjectHandle} = format_string(H)

"""
    section_header_offset(oh::ObjectHandle)

Given an `ObjectHandle`, return the offset (in bytes) at which the sections
start within the containing object file.
"""
@mustimplement section_header_offset(oh::ObjectHandle)

"""
    section_header_size(oh::ObjectHandle)

Given an `ObjectHandle`, return the size of a section header (used for reading
in the sections header when trying to load a `Section` object or iterating over
a `Sections` object)
"""
@mustimplement section_header_size(oh::ObjectHandle)

"""
    section_header_type(oh::ObjectHandle)

Given an `ObjectHandle`, return the type of a section header (used for reading
in the sections header when trying to load a `Section` object or iterating over
a `Sections` object).  For instance, for a 64-bit ELF file, this would return
the type `ELFSection64`
"""
@mustimplement section_header_type(oh::ObjectHandle)

"""
    segment_header_offset(oh::ObjectHandle)

Given an `ObjectHandle`, return the offset (in bytes) at which the segments
start within the containing object file.
"""
@mustimplement segment_header_offset(oh::ObjectHandle)

"""
    segment_header_size(oh::ObjectHandle)

Given an `ObjectHandle`, return the size of a segment header (used for reading
in the segments header when trying to load a `Segment` object or iterating over
a `Segments` object)
"""
@mustimplement segment_header_size(oh::ObjectHandle)

"""
    segment_header_type(oh::ObjectHandle)

Given an `ObjectHandle`, return the type of a segment header (used for reading
in the segments header when trying to load a `Segment` object or iterating over
a `Segments` object).  For instance, for a 64-bit ELF file, this would return
the type `ELFSegment64`
"""
@mustimplement segment_header_type(oh::ObjectHandle)

"""
    symtab_entry_offset(oh::ObjectHandle)

Given an `ObjectHandle`, return the offset (in bytes) at which the symbol
table starts within the containing object file.
"""
@mustimplement symtab_entry_offset(oh::ObjectHandle)

"""
    symtab_entry_size(oh::ObjectHandle)

Given an `ObjectHandle`, return the size of a symbol table entry (used for
reading in the symbol table when trying to load a `SymtabEntry` object or
iterating over a `Symbols` object).
"""
@mustimplement symtab_entry_size(oh::ObjectHandle)

"""
    symtab_entry_type(oh::ObjectHandle)

Given an `ObjectHandle`, return the type of a symbol table entry (used for
reading in the symbol table when trying to load a `SymtabEntry` object or
iterating over a `Symbols` object).  For instance, for a 64-bit ELF file, this
would return the type `ELFSymtabEntry64`
"""
@mustimplement symtab_entry_type(oh::ObjectHandle)

# Misc. stuff
"""
    path(oh::ObjectHandle)

Return the absolute path to the given `ObjectHandle`, if it was a file loaded
from the local disk.  If it was loaded from a general `IOStream` or in some
other way such that the path is unknown or unknowable, return the empty string.
"""
@mustimplement path(oh::ObjectHandle)

"""
    path(io::IO)

Try to guess the path of an `IO` object.  If it cannot be guessed, returns the
empty string.
"""
function path(io::IO)
    if startswith(io.name, "<file ") && endswith(io.name, ">")
        return abspath(io.name[7:end-1])
    end
    return ""
end

function show(io::IO, oh::H) where {H <: ObjectHandle}
    print(io, "$(format_string(H)) Handle ($(is64bit(oh) ? "64" : "32")-bit)")
end

"""
    find_library(oh::ObjectHandle, soname::String)

Return the absolute path to the given `soname`, using the linker search path
that the given `ObjectHandle` would use at runtime.  See the documentation for
`find_library(::RPath, ::String)` for more details.
"""
function find_library(oh::ObjectHandle, soname::AbstractString)
    return find_library(RPath(oh), soname)
end

"""
    find_libraries(oh::ObjectHandle)

Return a mapping from sonames to absolute paths, containing all the sonames
declared as beeing needed by the given `ObjectHandle`.  See the documentation
for `find_library(::RPath, ::String)` and `RPath` for more details.
"""
function find_libraries(oh::ObjectHandle)
    rpath = RPath(oh)
    sonames = [path(dl) for dl in DynamicLinks(oh)]

    # Remote '@rpath/' prefix if it exists
    function strip_rpath(soname)
        if startswith(soname, "@rpath/")
            return soname[8:end]
        end
        return soname
    end

    return Dict(s => find_library(oh, strip_rpath(s)) for s in sonames)
end


text:
using PumasTutorials, NCA, NCAUtilities, CSV, DataFrames, StatsPlots, Statistics, Unitful

data = DataFrame(CSV.File(joinpath(dirname(pathof(PumasTutorials)), "..", "tutorials", "nca", "NCA_tutorial.csv")))
first(data, 6)

single_dose_data  = filter(x -> x.TIME <= 24 , data)

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg";

pop        = read_nca(  single_dose_data,
                        id    = :ID,
                        time  = :TIME,
                        observations  = :CONC,
                        amt   = :AMT,
                        route = :ROUTE,
                        llq   = 0.001concu,
                        timeu = timeu,
                        concu = concu,
                        amtu  = amtu
                      )

NCAUtilities.observations_vs_time(pop[1])

report0     = run_nca(pop, sigdigits=3)

vz        = NCA.vz(pop, sigdigits=3)  # Volume of Distribution/F, in this case since the drug is given orally
cl        = NCA.cl(pop, sigdigits=3)  # Clearance/F, in this case since the drug is given orally
lambdaz   = NCA.lambdaz(pop, threshold=3, sigdigits=3)  # Terminal Elimination Rate Constant, threshold=3 specifies the max no. of time point used for calculation
lambdaz_1 = NCA.lambdaz(pop, slopetimes=[8,12,16].*timeu, sigdigits=3) # slopetimes in this case specifies the exact time point you want for the calculation
thalf     = NCA.thalf(pop[4], sigdigits=3) # Half-life calculation for 4th individual
cmax_d    = NCA.cmax(pop, normalize=true, sigdigits=3) # Dose Normalized Cmax
mrt       = NCA.mrt(pop, sigdigits=3) # Mean residence time
aumc      = NCA.aumc(pop, method=:linlog, sigdigits=3) # AUMC calculation, using :linlog method
rename!(lambdaz_1, Dict(:lambdaz => "lambdaz_specf")) #since we have two lambdaz calculation rename one column to merge in a dataframe
df_1      = innerjoin(vz,cl,lambdaz, lambdaz_1,cmax_d,mrt,aumc, on=[:id], makeunique=true)

auc0_12   = NCA.auc(pop, interval=(0,12).*timeu, method=:linuplogdown, sigdigits=3) #various other methods are :linear, :linlog
auc12_24  = NCA.auc(pop, interval=(12,24).*timeu, method=:linuplogdown, sigdigits=3)
final     = innerjoin(report0.reportdf, auc0_12, auc12_24, on = [:id], makeunique=true)

nca_output_single = select(final, :doseamt, :kel, :half_life, :cmax, :tmax, :auclast, :aucinf_obs, :vz_f_obs, :cl_f_obs, :aumcinf_obs, :auc0_12, :auc12_24)

Stat_report = describe(nca_output_single)

multiple_dose_data           = filter(x -> x.TIME > 24 , data)

pop1       = read_nca( multiple_dose_data,
                       id    = :ID,
                       time  = :TIME,
                       observations  = :CONC,
                       amt   = :AMT,
                       route = :ROUTE,
                       ii    = :II,        # please specify II for Multiple-dose NCA
                       llq   = 0.001concu,
                       timeu = timeu,
                       concu = concu,
                       amtu  = amtu)
report1    = run_nca(pop1, sigdigits=3)

Stat_report1        = summarize(report1.reportdf, parameters = [:doseamt, :kel, :half_life, :tmax, :cmax, :auclast, :auc_tau_obs, :vz_f_obs, :cl_f_obs, :aumcinf_obs, :tau, :cavgss])

CSV.write("NCA_single_dose.csv", final)      # Single-NCA final report
CSV.write("NCA_multiple_dose.csv", report1.reportdf)  # Multiple-NCA final report


text:
"""
    Camera

"""
abstract type Camera end

"""
    IdealCamera <: Camera

A camera with no added noise.

#Fields
- pixelsize
- xpixels
- ypixels
- gain
- offset

```
IdealCamera(;
pixelsize=0.1,
xpixels::Int=256,
ypixels::Int=256,
gain=1.0,
offset=0.0,
)

```
"""
mutable struct IdealCamera <: Camera
    pixelsize
    xpixels
    ypixels
    gain
    offset
end
function IdealCamera(;
    pixelsize=0.1,
    xpixels::Int=256,
    ypixels::Int=256,
    gain=1.0,
    offset=0.0,
    )
    return IdealCamera(pixelsize,xpixels,ypixels,gain,offset)
end


text:
using ShortStrings, BenchmarkTools

# these two should be on par
@benchmark hash($(ShortString{UInt128}("this is a test")), $(zero(UInt)))
@benchmark hash($(rand(UInt128)), $(zero(UInt)))

# faster for some reason
@benchmark hash($"this is a test", $(zero(UInt)))

# these two should be on par
@benchmark hash($(ShortString{UInt128}("this is a test")))
@benchmark hash($(rand(UInt128)))

# faster
@benchmark hash($"this is a test")


text:
module Deps

using DataDeps

function fetch_test_pots()
   register(DataDep(
         "JuLIP_testpots",
         "A few EAM potentials for testing",
         "https://www.dropbox.com/s/leub1c9ft1mm9fg/JuLIP_data.zip?dl=1",
         post_fetch_method = file -> run(`unzip $file`)
         ))
   return joinpath(datadep"JuLIP_testpots", "JuLIP_data") * "/"
end

end


text:
ENV["CUDA_VISIBLE_DEVICES"] = "2,3"
using Revise
using ADSeismic
using ADCME
using MAT
using PyPlot
using DelimitedFiles
using LineSearches
using Random
using Mamba
using ProgressMeter

# matplotlib.use("Agg")
close("all")

data_dir = "../data/acoustic"
if !ispath(data_dir)
  mkpath(data_dir)
end
figure_dir = "figure/"
if !ispath(figure_dir)
  mkpath(figure_dir)
end
result_dir = "result/"
if !ispath(result_dir)
  mkpath(result_dir)
end
# if isfile(joinpath(result_dir, "loss.txt"))
#   rm(joinpath(result_dir, "loss.txt"))
# end

################### Inversion using Automatic Differentiation #####################
reset_default_graph()
model_name = "../models/marmousi2-model-smooth.mat"
vp_true = matread("../models/marmousi2-model-true.mat")["vp"]

## load model setting
params = load_params(model_name)
src = load_acoustic_source(model_name)
rcv = load_acoustic_receiver(model_name)
vp0 = matread(model_name)["vp"]
mean_vp0 = mean(vp0)
std_vp0 = std(vp0)

## construct a neural, no dropout is present
Random.seed!(0);
z = constant(rand(Float32, 1,8))
x = Generator(z, ratio=size(vp0)[1]/size(vp0)[2], base=4)
x = cast(Float64, x)
x = x * std_vp0
vp = constant(vp0[:,:])
if size(x) <= size(vp)
  i = (size(vp0)[1]-size(x)[1])÷2 +1 :(size(vp0)[1]-size(x)[1])÷2 + size(x)[1]
  j = (size(vp0)[2]-size(x)[2])÷2 +1 : (size(vp0)[2]-size(x)[2])÷2 + size(x)[2]
  vp = scatter_add(vp, i, j, x)
elseif  size(x) > size(vp)
  vp = vp + tf.slice(x, (size(x).-size(vp0)).÷2, size(vp0))
else
  error("Size error: ", size(vp), size(x))
end

## assemble acoustic propagator model
x_var = Variable(ones(size(vp0)...))
model = x->AcousticPropagatorSolver(params, x, x_var)
vars = get_collection()

## load data
std_noise = 0
Random.seed!(1234);
Rs = Array{Array{Float64,2}}(undef, length(src))
σv_ = 0.0
for i = 1:length(src)
    Rs[i] = readdlm(joinpath(data_dir, "marmousi-r$i.txt"))
    Rs[i] .+= randn(size(Rs[i])) .* std(Rs[i]) .* std_noise
    global σv_ += std(Rs[i])
end
σv_ /= length(src)

Rs_ = compute_forward_GPU(model, src, rcv)
loss = sum([sum((x-y)^2) for (x,y) in zip(Rs, Rs_)])

sess = Session(); init(sess)
ADCME.load(sess, "NNFWI_05001.mat")
vp_ = run(sess, vp)[:]
@info run(sess, loss, feed_dict = Dict(x_var=>reshape(vp_, size(x_var)...)))

close("all")
visualize_model(run(sess, vp, x_var=>reshape(vp_, size(x_var)...)), params)
savefig("test-vp.png")

σm = std(vp_)
σv = σv_

## apply some scaling
σv *= 30
σm *= 10

function logf(vp_val)
  loss_val = run(sess, loss, feed_dict =
    Dict(x_var=>reshape(vp_val, size(x_var)...))
  )
  p1 = loss_val/(2σv^2)
  p2 = sum((vp_val .- vp_).^2)/(2σm^2)
  -p1-p2
end

n = 1000000
burnin = 1000
sim = Chains(n, 1)
theta = RWMVariate(vp_, 50*ones(length(x_var)), logf, proposal = SymUniform)
results = zeros(n, length(x_var))

v = theta
num_accept = 0
accept = true
for i in 1:n
  t0 = time()
  # sample!(theta)
  if accept
    global prob_v = logf(v.value)
  end
  x = v + v.tune.scale .* rand(v.tune.proposal(0.0, 1.0), length(v))
  prob_x = logf(x)
  if rand() < exp(prob_x - prob_v)
    v[:] = x
    global num_accept += 1
    global accept = true
  else
    global accept = false
  end
  results[i,:] = theta.value
  if mod(i, 100)==1
    writedlm("Data/UQ_VP.txt", results[1:i,:])
  end
  # @info i, "time consumed:", time() - t0
  @info i, "accept rate:", num_accept/i, "prob_v:", exp(prob_v)
  @info "exp(prob_x - prob_v)", exp(prob_x - prob_v)
end


text:
# ---
# title: 958. Check Completeness of a Binary Tree
# id: problem958
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/check-completeness-of-a-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, determine if it is a _complete binary
# tree_.
#
# In a **[complete binary
# tree](http://en.wikipedia.org/wiki/Binary_tree#Types_of_binary_trees)** ,
# every level, except possibly the last, is completely filled, and all nodes in
# the last level are as far left as possible. It can have between `1` and `2h`
# nodes inclusive at the last level `h`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/15/complete-binary-tree-1.png)
#
#
#
#     Input: root = [1,2,3,4,5,6]
#     Output: true
#     Explanation: Every level before the last is full (ie. levels with node-values {1} and {2, 3}), and all nodes in the last level ({4, 5, 6}) are as far left as possible.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/15/complete-binary-tree-2.png)
#
#
#
#     Input: root = [1,2,3,4,5,null,7]
#     Output: false
#     Explanation: The node with value 7 isn't as far left as possible.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 100]`.
#   * `1 <= Node.val <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Tests for Modified DDF DEA Models
@testset "ModifiedDDFDEAModel" begin

    # Data for Modified DDF test
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    # Modified DDF CRS Ones
    deamddfcrs = deamddf(X, Y, Gx = :Ones, Gy = :Ones, rts = :CRS)

    @test typeof(deamddfcrs) == ModifiedDDFDEAModel

    @test nobs(deamddfcrs) == 8
    @test ninputs(deamddfcrs) == 1
    @test noutputs(deamddfcrs) == 1
    @test efficiency(deamddfcrs) ≈ [1.5; 0; 2; 6; 4.5; 10.5; 8.5; 9.412] atol = 1e-5
    @test efficiency(deamddfcrs, :X) ≈ zeros(8, 1) atol = 1e-5
    @test efficiency(deamddfcrs, :Y) ≈ [1.5; 0; 2; 6; 4.5; 10.5; 8.5; 9.412] atol = 1e-5
    @test convert(Matrix, peers(deamddfcrs)) ≈ [
        0 0.5     0    0    0    0    0    0
        0 1.0     0    0    0    0    0    0
        0 2.0     0    0    0    0    0    0
        0 3.0     0    0    0    0    0    0
        0 1.5     0    0    0    0    0    0
        0 3.5     0    0    0    0    0    0
        0 3.5     0    0    0    0    0    0
        0 2.353   0    0    0    0    0    0]  atol = 1e-5

    @test slacks(deamddfcrs, :X) ≈ zeros(8, 1)
    @test slacks(deamddfcrs, :Y) ≈ zeros(8, 1)

    @test efficiency(deamddf(targets(deamddfcrs, :X), targets(deamddfcrs, :Y), Gx = :Ones, Gy = :Ones, rts = :CRS)) ≈ zeros(8,1) atol = 1e-5

    # Test no slacks
    deamddfnoslack = deamddf(X, Y, Gx = :Ones, Gy = :Ones, slack = false)
    @test efficiency(deamddfnoslack) == efficiency(deamddfcrs)
    @test isempty(slacks(deamddfnoslack, :X)) == 1
    @test isempty(slacks(deamddfnoslack, :Y)) == 1

    # Test default is CRS and :Ones
    @test efficiency(deamddf(X, Y, Gx = :Ones, Gy = :Ones)) ≈ efficiency(deamddfcrs)

    # Modified DDF VRS Ones
    deamddfvrs = deamddf(X, Y, Gx = :Ones, Gy = :Ones, rts = :VRS)

    @test efficiency(deamddfvrs) ≈ [0; 0; 0; 0; 4; 22/3; 2; 8.059] atol = 1e-5
    @test efficiency(deamddfvrs, :X) ≈ [0; 0; 0; 0; 2.0; 22/3; 2; 5.412] atol = 1e-5
    @test efficiency(deamddfvrs, :Y) ≈ [0; 0; 0; 0; 2; 0; 0; 2.647] atol = 1e-5
    @test convert(Matrix, peers(deamddfvrs)) ≈ [
        1.0  0     0    0    0    0    0    0
        0    1.0   0    0    0    0    0    0
        0    0     1.0  0    0    0    0    0
        0    0.0   0    1    0    0    0    0
        0    1.0   0    0    0    0    0    0
        0    1/3   2/3  0    0    0    0    0
        0    0     0    1    0    0    0    0
        0    1     0    0    0    0    0    0] atol = 1e-5

    @test slacks(deamddfvrs, :X) ≈ zeros(8, 1)
    @test slacks(deamddfvrs, :Y) ≈  zeros(8, 1)

    @test efficiency(deamddf(targets(deamddfvrs, :X), targets(deamddfvrs, :Y), Gx = :Ones, Gy = :Ones, rts = :VRS)) ≈ zeros(8,1) atol = 1e-5

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deamddf_crs_ref_eff = zeros(size(X, 1))
    deamddf_vrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deamddf_crs_ref_eff[i] = efficiency(deamddf(Xeval, Yeval, Gx = :Ones, Gy = :Ones,  rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deamddf_vrs_ref_eff[i] = efficiency(deamddf(Xeval, Yeval, Gx = :Ones, Gy = :Ones, rts = :VRS, Xref = Xref, Yref = Yref))[1]
    end

    @test deamddf_crs_ref_eff ≈ efficiency(deamddfcrs)
    @test deamddf_vrs_ref_eff ≈ efficiency(deamddfvrs)

    # Modified DDF VRS Observed
    deamddfvrsobs = deamddf(X, Y, Gx = :Observed, Gy = :Observed, rts = :VRS)

    @test efficiency(deamddfvrsobs) ≈ [0; 0; 0; 0; 35/30; 0.571429; 0.142857; 2.549936] atol = 1e-5
    @test efficiency(deamddfvrsobs, :X) ≈ [0; 0; 0; 0; 0; 0.428571; 0.142857; 0.0] atol = 1e-5
    @test efficiency(deamddfvrsobs, :Y) ≈ [0; 0; 0; 0; 35/30; 0.142857; 0.0; 2.549936] atol = 1e-5

    # Modified DDF VRS Mean
    deamddfvrsmean = deamddf(X, Y, Gx = :Mean, Gy = :Mean, rts = :VRS)

    @test efficiency(deamddfvrsmean) ≈ [0; 0; 0; 0; 0.631299; 0.871894; 0.230508; 1.181294] atol = 1e-5
    @test efficiency(deamddfvrsmean, :X) ≈ [0; 0; 0; 0; 0; 0.691523; 0.230508; 0.162738] atol = 1e-5
    @test efficiency(deamddfvrsmean, :Y) ≈ [0; 0; 0; 0; 0.631299; 0.180371; 0.0; 1.018556] atol = 1e-5

    # Test Custom direction
    @test efficiency(deamddf(X, Y, Gx = ones(8,1), Gy = ones(8,1), rts = :VRS)) ≈ efficiency(deamddfvrs)

    # Print
    show(IOBuffer(), deamddfcrs)
    show(IOBuffer(), deamddfnoslack)
    show(IOBuffer(), deamddfvrs)

    # ------------------
    # Test model with zero Directions
    # ------------------
    deamddfzeros = deamddf(X, Y, Gx = zeros(size(X)), Gy = zeros(size(Y)), rts = :VRS)
    @test efficiency(deamddfzeros) == zeros(8)

    # ------------------
    # Test slacks with different data
    # ------------------
    X = [1 1; 2 2; 1.5 1]
    Y = [1 1; 2 2; 1 0.5]

    deamddfslack = deamddf(X, Y, Gx = :Ones, Gy = :Ones, rts = :VRS)
    @test slacks(deamddfslack, :X) ≈ [0 0; 0 0; 0.5 0] atol = 1e-5
    @test slacks(deamddfslack, :Y) ≈ [0 0; 0 0; 0 0.5] atol = 1e-5

    # ------------------
    # Test errors
    # ------------------
    @test_throws DimensionMismatch deamddf([1; 2 ; 3], [4 ; 5], Gx = :Ones, Gy = :Ones) #  Different number of observations
    @test_throws DimensionMismatch deamddf([1; 2], [4 ; 5], Gx = :Ones, Gy = :Ones, Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deamddf([1 1; 2 2], [4 4; 5 5], Gx = :Ones, Gy = :Ones, Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deamddf([1 1; 2 2], [4 4; 5 5], Gx = :Ones, Gy = :Ones, Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deamddf([1; 2; 3], [4; 5; 6], Gx = :Ones, Gy = :Ones, rts = :Error) # Invalid returns to scale
    @test_throws DimensionMismatch deamddf([1 1; 2 2; 3 3], [4; 5; 6], Gx = [1 1 1; 2 2 2; 3 3 3], Gy = [4; 5; 6]) # Different size of inputs direction
    @test_throws DimensionMismatch deamddf([1; 2; 3], [4 4; 5 5; 6 6], Gx = [1; 2; 3], Gy = [4 4 4; 5 5 5; 6 6 6]) # Different size of inputs direction
    @test_throws ArgumentError deamddf([1; 2; 3], [1; 2; 3], Gx = :Error, Gy = :Ones) # Invalid inuts direction
    @test_throws ArgumentError deamddf([1; 2; 3], [1; 2; 3], Gx = :Ones, Gy = :Error) # Invalid outputs direction

    @test_throws ArgumentError efficiency(deamddfvrs, :Error) # Invalid efficiency type

end


text:
"""
Compute a histogram of the angles `θ` using equally sized bins spanning [-π,π]
"""
function angular_histogram(θ::AbstractVector{T}, Δθ::Real) where T <: Real
    nbins = round(Int64,2π/Δθ)
    bins = collect(range(0.0, step=Δθ,length=div(nbins,2)))
    append!(bins, range(-pi, step=Δθ, length=div(nbins,2)))
    nbins = length(bins)
    sort!(bins)
    counts = fill(0, nbins)
    binidx = fill(0, length(θ))
    for i in 1:nbins
        idx = bins[i] - Δθ/2 .<= θ .< bins[i] + Δθ/2
        counts[i] = sum(idx)
        binidx[idx] .= i
    end
    counts, binidx, bins
end


text:
#######################################################################################
#                                   Load Packages
#######################################################################################
# change directory of this file
cd(@__DIR__)
# pakcage manager
using Pkg
# activate this environment of this project
Pkg.activate("../../")
# load dependencies
using Turing, DataFrames, StatsPlots, Revise, ACTRModels
# load stimuli
include("Stimuli.jl")
# load model functions
include("Paired_Model.jl")
# set seed for RNG
Random.seed!(589)
#######################################################################################
#                                   Generate Data
#######################################################################################
# true decay parameter
d = 0.5
# number of block repetitions
n_blocks = 8
# number of trials within a block
n_trials = 20
# fixed parameters
fixed_parms = (
    τ = -2.5,        # retrieval threshold
    noise = true,    # noise "on"
    bll = true,      # base-level learning "on"
    s = 0.2,         # scale parameter for activation noise
    lf = 0.4,        # latency factor
    ter = 0.535      # encoding, motor and conflict resolution time
)
# generate data
data = simulate(all_stimuli, fixed_parms, n_blocks, n_trials; d)
data[1]
#######################################################################################
#                                    Define Model
#######################################################################################
@model model(data, parms) = begin
    d ~ Beta(5, 5)
    data ~ Paired(d, parms)
end
#######################################################################################
#                                 Estimate Parameters
#######################################################################################
# number of samples
n_samples = 1000
delta = 0.85
# number of adaption samples
n_adapt = 1000
# number of independent chains
n_chains = 4
# sampler object
specs = NUTS(n_adapt, delta)
# Start sampling.
chain = sample(model(data, fixed_parms), specs, MCMCThreads(), n_samples, n_chains, progress=true)
#######################################################################################
#                                         Plot
#######################################################################################
pyplot()
post_plot = plot(chain, xaxis=font(5), yaxis=font(5), seriestype=(:density),
  grid=false, titlefont=font(5), size=(200,200))
plot!(title = "", xlabel="d", subplot=1)
#######################################################################################
#                                  Posterior Predictive
#######################################################################################
preds = posterior_predictive(x -> show_learning(all_stimuli, fixed_parms; x...), chain, 100)
preds = vcat(preds..., source=:rep)
groups = groupby(preds, [:rep,:retrieved,:block])
mean_rts = combine(groups, :rt=>mean=>:mean_rt)
sort!(mean_rts, [:rep,:block,:retrieved])
fmean_rts = filter(x->x.retrieved != :truncated, mean_rts)

rt_plot = @df fmean_rts plot(:block, :mean_rt, group=(:rep,:retrieved), xlabel="Block", ylabel="Mean RT", xaxis=font(9),
    yaxis=font(10), color=[:darkred :black], grid=false, size=(300,200), titlefont=font(10), linewidth=.5,
    labeltitle="Retrieved", legendtitlefontsize=8, ylims=(0,5.1), leg=false, alpha=.5)

groups = groupby(preds, :block)
accuracy = combine(groups, :retrieved=>(x->mean(x .== :retrieved))=>:accuracy)
sort!(accuracy)
accuracy_plot = @df accuracy plot(:block, :accuracy, xlabel="Block", ylabel="Mean Accuracy", xaxis=font(10),
    yaxis=font(10), color=:black, grid=false, size=(300,200), titlefont=font(10), linewidth=2,
    leg=false, legendtitlefontsize=8)
savefig(post_plot, "Paired_Posterior.eps")
savefig(rt_plot, "Paired_RTs.eps")
savefig(accuracy_plot, "Paired_Accuracy.eps")


text:
function demoMaxProd()
    #DEMOMAXPROD Max-Product algorithm test :
    # Variable order is arbitary
    variables=1:5
    a, b, c, d, e=variables
    nstates=round(Int64,3*rand(1,5)+2) # random number of states for each variable

    pot=Array(PotArray,5)
    pot[1]=PotArray([a b],rand(nstates[[a b]]...));
    pot[2]=PotArray([b c d],rand(nstates[[b c d]]...));
    pot[3]=PotArray([c],rand(nstates[[c]]...));
    pot[4]=PotArray([e d],rand(nstates[[e d]]...));
    pot[5]=PotArray([d],rand(nstates[[d]]...));

    A = FactorGraph(pot)
    maxstate, maxval, mess=maxprod(pot,A)

    # check if this is correct by brute force max:
    jointpot =prod(pot); V=length(pot)
    maxval2,maxstate2=max(jointpot,variables,ReturnState=true,Ind2Sub=true)
    println("\t\t maxprod \t brute")
    for i=1:V
        println("Variable $i \t $(maxstate[i]...)\t\t $(maxstate2[i])")
    end

end


text:
using PGFPlotsX
import PGFPlotsX: Plot, PlotInc, Options, Axis

for F in (:Expansion,)
    @eval begin
        Axis(F::$F, trailing...; opts...) =
            @pgf Axis({ymode="log"}, Plot(F, trailing...; opts...))
    end
end

for plot in (:Plot, :PlotInc)
    _plot = Meta.parse("_"*string(plot))
    for EXP in (:Expansion,)
        @eval begin
            $(plot)(F::$EXP; opts...) =
                $(plot)(Options(), F; opts...)
            $(plot)(f::Function, F::$EXP; opts...) =
                $(plot)(Options(), f, F; opts...)
            $(plot)(F::$EXP, f::Function; opts...) =
                $(plot)(Options(), f, F)
            $(plot)(F::$EXP, f::$EXP; opts...) =
                $(plot)(Options(), f, F; opts...)
        end
    end
    @eval begin
        function $(plot)(options::Options, F::Expansion; n=200, plot_complex=false, opts...)
            grid = plotgrid(dictionary(F), n)
            if plot_complex
                vals = F(grid)
                options = @pgf {options..., no_markers}
                $(plot)(options, Table([grid, BasisFunctions.postprocess(dictionary(F), grid, real.(vals))])),
                    $(plot)(options, Table([grid, BasisFunctions.postprocess(dictionary(F), grid, imag.(vals))]))
            else
                vals = real.(F(grid))
                options = @pgf {options..., no_markers}
                $(plot)(options, Table([grid, BasisFunctions.postprocess(dictionary(F), grid, vals)]))
            end
        end

        $(plot)(options::Options, F::Expansion, f::Function; opts...) =
            $(_plot)(options, f, F; opts...)

        $(plot)(options::Options, f::Function, F::Expansion; opts...) =
            $(_plot)(options, f, F; opts...)
        $(plot)(options::Options, f::Expansion, F::Expansion; opts...) =
            $(_plot)(options, f, F; opts...)

        function $(_plot)(options::Options, f, F::Expansion; plot_complex=false, n=200, opts...)
            grid = plotgrid(dictionary(F), n)
            vals = abs.(f.(grid) - F.(grid))
            options = @pgf {options..., no_markers}
            $(plot)(options, Table([grid, BasisFunctions.postprocess(dictionary(F), grid, vals)]))
        end
    end
end


text:
include("preamble.jl")

include("hitting_times_model/reflecting_model/model_def.jl")
#model.P_lwr[1], model.P_lwr[2] = [0.0 1.0]
#model.P_upr[1], model.P_upr[2] = [1.0 0.0]
function exp_initial_condition(x,i)
    return (i∈(1:2))*exp(-x)/(1-exp(-model.b))/2
end

t = 2.0
function run_dg(t)
    n_bases_per_cell = 21
    nodes = 0.0:10.0
    mesh = DGMesh(nodes, n_bases_per_cell)

    dq = DiscretisedFluidQueue(model,mesh)
    B = build_full_generator(dq)

    d0 = SFMDistribution(exp_initial_condition, dq)

    t_step = 0.005
    t_integration_scheme = StableRK4(t_step)

    dt = integrate_time( d0, B, t, t_integration_scheme; limiter=NoLimiter)
    return dt
end

n_bases_per_cell = 15
nodes = 0.0:10.0
mesh = DGMesh(nodes, n_bases_per_cell)

dq = DiscretisedFluidQueue(model,mesh)
B = build_full_generator(dq)

d0 = SFMDistribution(exp_initial_condition, dq)# TrapezoidRule, fun_evals=10_001)

dt = run_dg(t)
#dt.coeffs[:] = d0.coeffs'*exp(B*2)
dt_cdf = cdf(dt)

x = range(eps(),10-10*eps(),length=10001)
dt_cdf_eval = dt_cdf.( x, (1:2)')

sum(abs.(dt_cdf_eval[:,1].-sim[1]./sim[1][end]*dt_cdf_eval[end,1]))*10/length(x) + sum(abs.(dt_cdf_eval[:,2].-sim[2]./sim[2][end]*dt_cdf_eval[end,2]))*10/length(x)

##

sim = JSON.parsefile(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/"*
    "reflecting_model/transient_distribution/exp/data/sims/cdf_evaluated.json"
)

plot!( x, sim[2]-dt_cdf_eval[:,2])
plot!( x, sim[1]-dt_cdf_eval[:,1])

tau = fixed_time(t)

rng = StableRNGs.StableRNG(16092021)#StableRNGs.StableRNG(1)
function initial_condition_generator(n, rng)
    i = rand(rng,1:2,n)
    X = Vector{Float64}(undef,n)
    for k in 1:n
        x = Inf
        while x>=model.b
            x = -log(rand(rng))
        end
        X[k]=x
    end
    # X=5.0*ones(n)#10.0*rand(rng,n)
    return (X=X,φ=i)
end
n_sims = 5_000_000
sim_initial_condition = initial_condition_generator(n_sims, rng)

sims = simulate(model, tau, sim_initial_condition)

const idx1 = sims.φ.==1
const idx2 = sims.φ.==2

sim_cdf1(x) = sum(sims.X[idx1].<x)./n_sims
sim_cdf2(x) = sum(sims.X[idx2].<x)./n_sims

cdf1_eval = sim_cdf1.(x)
cdf2_eval = sim_cdf2.(x)

const idx1 = sim_initial_condition.φ.==1
const idx2 = sim_initial_condition.φ.==2

n1 = sum(idx1)
n2 = sum(idx2)

for i in 1:n_sims
    if n1>n2
        if idx1[i]
            idx1[i]=0
            n1 -= 1
        end
    elseif n2>n1
        if idx2[i]
            idx2[i]=0
            n2 -= 1
        end
    else
        break
    end
end

n_new = n1+n2

sim_cdf1(x) = sum(sim_initial_condition.X[idx1].<x)./n_new
sim_cdf2(x) = sum(sim_initial_condition.X[idx2].<x)./n_new

cdf1_eval = sim_cdf1.(x)
cdf2_eval = sim_cdf2.(x)

plot( x, cdf2_eval+cdf1_eval-dt_cdf_eval[:,2]-dt_cdf_eval[:,1],label="2")
plot!( x, cdf1_eval-dt_cdf_eval[:,1],label="1")

sum(abs.(dt_cdf_eval[:,1].-cdf1_eval))*10/length(x) + sum(abs.(dt_cdf_eval[:,2].-cdf2_eval))*10/length(x)

x2 = (x[1:end-1]+x[2:end])./2
plot!( x2, diff(cdf2_eval)*length(x)/10-pdf(d0).(x2,2),label="2")
plot!( x2, diff(cdf1_eval)*length(x)/10-pdf(d0).(x2,1),label="1")

plot( x2, diff(sim_cdf2.(x))*length(x)/10-pdf(dt).(x2,2),label="2")
plot!( x2, diff(sim_cdf1.(x))*length(x)/10-pdf(dt).(x2,1),label="1")
#plot( x, log10.(abs.(cdf(d0).(x,1)-exp_ic_cdf.(x, 1))))
#plot!( x, log10.(abs.(cdf(d0).(x,2)-exp_ic_cdf.(x, 2))))

exp_ic_cdf(x,i) = (1-exp(-x))/(1-exp(-model.b))/2

2*sum(abs.(cdf(d0).(x,1)-exp_ic_cdf.(x, 1)))*10.0/10_001
2*sum(abs.(pdf(d0).(x,1)-exp_initial_condition.(x, 1)))*10.0/10_001

sum(abs.(dt_cdf_eval[:,1]-cdf1_eval))*10.0/10_001+sum(abs.(dt_cdf_eval[:,2]-cdf2_eval))*10.0/10_001
sum(abs.(dt_cdf_eval[:,1]-exp_ic_cdf.(x,1)))*10.0/10_001+
    sum(abs.(dt_cdf_eval[:,2]-exp_ic_cdf.(x,2)))*10.0/10_001

plot(x,abs.(dt_cdf_eval[:,1]-exp_ic_cdf.(x,1)))
plot!(x,abs.(dt_cdf_eval[:,2]-exp_ic_cdf.(x,2)))

plot(x,dt_cdf_eval[:,1])
plot!(x,dt_cdf_eval[:,2])

plot((x[1:end-1]+x[2:end])./2,diff(cdf1_eval)*10_001/10)
plot!(x,pdf(dt).(x,1))
plot!(x,pdf(dt).(x,2))

## qbd
function run1(t,n)
    n_bases_per_cell = n
    nodes = range(0,10,length=11)
    mesh = FRAPMesh(nodes, n_bases_per_cell)
    @show n_bases_per_phase(mesh)

    dq = DiscretisedFluidQueue(model,mesh)
    B = build_full_generator(dq)
    #display(Matrix(B))

    d0 = SFMDistribution(exp_initial_condition, dq, TrapezoidRule; fun_evals = 20_001)
    t_step = diff(zglj(n,0,0))[1]/2/2 # smallest diff between nodes on [-1,1], divided by 2 to get to [0,1], then divided by 2 to be safe
    t_integration_scheme = StableRK4(t_step)
    dt = integrate_time( d0, B, t, t_integration_scheme)
    return dt
end
function run2(t,n)
    mesh2 = DGMesh(range(0,10,length=11), n)
    @show n_bases_per_phase(mesh2)

    dq2 = DiscretisedFluidQueue(model,mesh2)
    B2 = build_full_generator(dq2)

    #display(Matrix(B2))
    d02 = SFMDistribution(exp_initial_condition, dq2)
    t_step = diff(zglj(n,0,0))[1]/2/2 # smallest diff between nodes on [-1,1], divided by 2 to get to [0,1], then divided by 2 to be safe
    t_integration_scheme = StableRK4(t_step)
    dt2 = integrate_time( d02, B2, t, t_integration_scheme; limiter = NoLimiter)
    return dt2
end
function run_all_qbd()
    plot(layout=(2,1))
    dt_pdf_eval_qbd = zeros(10_001,2)
    dt_pdf_eval_dg = zeros(10_001,2)
    for n in 49
        tqbd = 2.0

        dt_dg, dt_qbd = run1(tqbd,n), run2(tqbd,n)

        dt_cdf_qbd = cdf(dt_qbd)

        dt_cdf_dg = cdf(dt_dg)
        dt_pdf_qbd = pdf(dt_qbd)

        dt_pdf_dg = pdf(dt_dg)

        x = range(eps(),10-10*eps(),length=10_001)
        dt_pdf_eval_qbd .= dt_pdf_qbd.( x, (1:2)')
        dt_pdf_eval_dg .= dt_pdf_dg.( x, (1:2)')

        println(n)
        println(sum(abs.(dt_pdf_eval_qbd[:,1].-dt_pdf_eval_dg[:,1]))*10/length(x) +
            sum(abs.(dt_pdf_eval_qbd[:,2].-dt_pdf_eval_dg[:,2]))*10/length(x))
        # display(plot!((abs.(dt_pdf_eval_qbd[:,1].-dt_pdf_eval_dg[:,1])) +
        #     (abs.(dt_pdf_eval_qbd[:,2].-dt_pdf_eval_dg[:,2]))))
        display(plot!(dt_pdf_eval_dg, subplot=1))
        display(plot!(dt_pdf_eval_qbd, subplot=2))
    end
    return dt_pdf_eval_qbd, dt_pdf_eval_dg
end

qbd_d, dg_d = run_all_qbd()

n=25

plot(x,dt_cdf_eval_qbd.-[sim[1][:] sim[2][:]])
plot!(x,dt_cdf_eval.-[sim[1][:] sim[2][:]])

exp_err(p::Float64, n) = sqrt(2/pi)*sqrt(p*(1-p)/n) #sqrt(1/(2*pi*(1-p)*p*n))

n = 5_000_000

sims_pm_transient = JSON.parsefile(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/reflecting_model/"*
    "transient_distribution/point_mass/data/sims/cdf_evaluated.json"
)

p_pm_transient = sims_pm_transient[1][end]

exp_err(p_pm_transient,n)

sims_exp_transient = JSON.parsefile(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/"*
    "reflecting_model/transient_distribution/exp/data/sims/cdf_evaluated.json"
)

p_exp_transient = sims_exp_transient[1][end]

exp_err(p_exp_transient,n)

sims_exp_hitting = CSV.read(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/hitting_times/data/exp/sims.csv",
    DataFrame
)

p_exp_hitting = sum(sims_exp_hitting.φ.==1)/nrow(sims_exp_hitting)

exp_err(p_exp_hitting,n)

sims_pm_hitting = CSV.read(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/hitting_times_model/hitting_times/data/point_mass/sims.csv",
    DataFrame
)

p_pm_hitting = sum(sims_pm_hitting.φ.==1)/nrow(sims_pm_hitting)

exp_err(p_pm_hitting,n)

sims_first_return_discont = CSV.read(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/fluidfluid_discontinuous/data/sims_cdf_evaluated.csv",
    DataFrame
)

p_first_return_discont = sims_first_return_discont[end,end]

exp_err(p_first_return_discont,n)

sims_first_return_cont = CSV.read(
    "/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/fluidfluid_continuous/data/sims_cdf_evaluated.csv",
    DataFrame
)

p_first_return_cont = sims_first_return_cont[end,end]

exp_err(p_first_return_cont,n)

function integrateD(evals,params)
    # evals is an integer specifying how many points to eval the function at
    # params is a CMEParams dictionary entry, i.e. CMEParams[3]
    N = 2*params["n"]+1 # ME order

    α = zeros(N)
    α[1] = params["c"]
    a = params["a"]
    b = params["b"]
    ω =  params["omega"]
    for k in 1:params["n"]
        kω = k*ω
        α[2*k] = (1/2)*( a[k]*(1+kω) - b[k]*(1-kω) )/(1+kω^2)
        α[2*k+1] = (1/2)*( a[k]*(1-kω) + b[k]*(1+kω) )/(1+kω^2)
    end

    period = 2*π/ω # the orbit repeats after this time
    edges = range(0,period,length=evals+1) # points at which to evaluate the fn
    h = period/(evals)

    orbit_LHS = α
    orbit_RHS = zeros(N)
    v_RHS = zeros(N)
    v_RHS[1] = 1
    v_LHS = ones(N)
    D = zeros(N,N)
    for t in edges[2:end]
        orbit_RHS[1] = α[1]
        for k in 1:params["n"]
            kωt = k*ω*t
            idx = 2*k
            idx2 = idx+1
            temp_cos = cos(kωt)
            temp_sin = sin(kωt)
            orbit_RHS[idx] = α[idx]*temp_cos + α[idx2]*temp_sin
            orbit_RHS[idx2] = -α[idx]*temp_sin + α[idx2]*temp_cos
            v_RHS[idx] = temp_cos - temp_sin
            v_RHS[idx2] = temp_sin + temp_cos
        end
        orbit_RHS = orbit_RHS./sum(orbit_RHS)
        orbit = (orbit_LHS+orbit_RHS)./2

        v = exp(-(t-h))*(v_LHS - exp(-h)*v_RHS)

        Dᵢ = v*orbit'
        D += Dᵢ

        orbit_LHS = copy(orbit_RHS)
        v_LHS = copy(v_RHS)
    end
    D = (1/(1-exp(-period)))*D
    return D
end

# numerical approximation of D
k = 100_000_000
T=[]
N=[]
n=184
evals = Int(ceil(k/(n/3)))# number of function evals
CMEParams[n]["D"], t = @timed integrateD(evals,CMEParams[n])
CMEParams[n]["intDevals"] = evals
push!(N,n)
push!(T,t)


text:
using NativeHTML

frag = fragment() do
    div() do
        span("Hello", style="color:red;")
        text(", ")
        b("World!")
    end
end

println(frag)
display(frag)

page = html() do
    head() do
        title("Hello, World!")
    end
    body() do
        div() do
            span("Hello", inline=true, style="color:red;")
            text(", ")
            b("World!", inline=true)
        end
    end
end

println(page)
display(page)

sleep(1)

page2 = html() do
    head() do
        title("Hello again, world!")
    end
    body() do
        div() do
            text("Block element 1")
        end
        div() do
            text("block element 2")
        end
    end
end

println(page2)
display(page2)


text:
using Documenter
using Rimu
using Rimu.ConsistentRNG
using Rimu.BitStringAddresses
using Rimu.StatsTools
using Literate

EXAMPLES_OUTPUT = joinpath(@__DIR__, "src/generated")

Literate.markdown(
    joinpath(@__DIR__, "../scripts/BHM-example.jl"), EXAMPLES_OUTPUT; documenter=true
)
Literate.markdown(
    joinpath(@__DIR__, "../scripts/BHM-example-mpi.jl"), EXAMPLES_OUTPUT; documenter=true
)

makedocs(;
    modules=[Rimu,Rimu.ConsistentRNG],
    format=Documenter.HTML(prettyurls = false),
    pages=[
        "Guide" => "index.md",
        "Examples" => [
            "1D Bose-Hubbard Model" => "generated/BHM-example.md",
            "Using MPI" => "generated/BHM-example-mpi.md",
        ],
        "User documentation" => [
            "StatsTools" => "statstools.md",
        ],
        "Developer documentation" => [
            "Interfaces" => "interfaces.md",
            "Hamiltonians" => "hamiltonians.md",
            "Dict vectors" => "dictvectors.md",
            "BitString addresses" => "addresses.md",
            "Stochastic styles" => "stochasticstyles.md",
            "RMPI" => "RMPI.md",
            "Random Numbers" => "consistentrng.md",
            "Documentation generation" => "documentation.md",
            "Code testing" => "testing.md",
        ],
        "API" => "API.md",
    ],
    sitename="Rimu.jl",
    authors="Joachim Brand <j.brand@massey.ac.nz>",
    doctest=false # Doctests are done while testing.
)

deploydocs(
    repo = "github.com/joachimbrand/Rimu.jl.git",
    push_preview = true,
)