text:
# Compute right eigenvectors of a complex upper triangular matrix TT.
# If Z is nontrivial, multiply by it to get eigenvectors of Z*TT*Z'.
# based on LAPACK::ztrevc
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function _geigvecs!(TT::StridedMatrix{T},
                    Z::StridedMatrix{T}=Matrix{T}(undef,0,0)
                    ) where {T <: Complex}
    n = size(TT,1)
    RT = real(T)
    ulp = eps(RT)
    # Note: LAPACK has
    #   smallnum = safemin(RT) * (n / ulp)
    # but that makes no sense to me and breaks some tests
    smallnum = safemin(RT) * n
    vectors = Matrix{T}(undef,n,n)
    v = zeros(T,n)

    # save diagonal since we modify it to avoid copies
    tdiag = diag(TT)

    # We use the 1-norms of the strictly upper part of TT columns
    # to avoid overflow
    tnorms = zeros(RT,n)
    @inbounds for j=2:n
        for i=1:j-1
            tnorms[j] += abs(TT[i,j])
        end
    end

    for ki=n:-1:1
        smin = max(ulp * abs1(TT[ki,ki]), smallnum)
        #
        # (T[1:k,1:k]-λI) x = b
        # where k=kᵢ-1

        v[1] = one(T) # for ki=1
        @inbounds for k=1:ki-1
            v[k] = -TT[k,ki]
        end
        @inbounds for k=1:ki-1
            TT[k,k] -= TT[ki,ki]
            (abs1(TT[k,k]) < smin) && (TT[k,k] = smin)
        end
        if ki > 1
            vscale = _usolve!(TT,ki-1,v,view(tnorms, 1:ki-1))
            v[ki] = vscale
        else
            vscale = one(RT)
        end
        if size(Z,1) > 0
            # This is done here to avoid allocating a work matrix
            # and to exploit the subspace property to reduce work.
            # Using a work matrix would allow for level-3 ops (cf. ztrevc3).
            @inbounds for j=1:n
                vectors[j,ki] = vscale * Z[j,ki]
                for i=1:ki-1
                    vectors[j,ki] += Z[j,i] * v[i]
                end
            end
        else
            @inbounds for j=1:ki
                vectors[j,ki] = v[j]
            end
            vectors[ki+1:n,ki] .= zero(T)
        end

        # normalize
        t0 = abs1(vectors[1,ki])
        @inbounds for i=2:n; t0 = max(t0, abs1(vectors[i,ki])); end
        remax = one(RT) / t0
        @inbounds for i=1:n; vectors[i,ki] *= remax; end

        # restore diagonal
        @inbounds for k=1:ki-1
            TT[k,k] = tdiag[k]
        end
    end

    vectors
end

# based on LAPACK::ztrevc
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function _gleigvecs!(TT::StridedMatrix{T},
                    Z::StridedMatrix{T}=Matrix{T}(undef,0,0)
                    ) where {T <: Complex}
    n = size(TT,1)
    RT = real(T)
    ulp = eps(RT)
    # replaces LAPACK's smallnum = safemin(RT) * (n / ulp)
    smallnum = safemin(RT) * n
    vectors = Matrix{T}(undef,n,n)
    v = zeros(T,n)

    # save diagonal since we modify it to avoid copies
    tdiag = diag(TT)

    # we use the 1-norm of strictly upper part of TT
    # cols to control overflow
    tnorms = zeros(RT,n)
    @inbounds for j=2:n
        for i=1:j-1
            tnorms[j] += abs1(TT[i,j])
        end
    end

    # logic from ztrevc
    for ki=1:n
        smin = max(ulp * abs1(TT[ki,ki]), smallnum)

        # form RHS
        @inbounds for k=ki+1:n
            v[k] = -conj(TT[ki,k])
        end
        v[1:ki] .= zero(T)
        # solve (T[ki+1:n,ki+1:n] - T[ki,ki])ᴴ x = σ v
        @inbounds for k=ki+1:n
            TT[k,k] -= TT[ki,ki]
            (abs1(TT[k,k]) < smin) && (TT[k,k] = smin)
        end
        if ki < n
            vscale = _cusolve!(view(TT,ki+1:n,ki+1:n),n-ki,view(v,ki+1:n),tnorms)
            v[ki] = vscale
        else
            v[n] = one(T)
            vscale = one(RT)
        end
        if size(Z,1) > 0
            # This is done here to avoid allocating a work matrix
            # and to exploit the subspace property to reduce work.
            # Using a work matrix would allow for level-3 ops (cf. ztrevc3).
            @inbounds for j=1:n
                vectors[j,ki] = vscale * Z[j,ki]
                for i=ki+1:n
                    vectors[j,ki] += Z[j,i] * v[i]
                end
            end
        else
            @inbounds for j=ki:n
                vectors[j,ki] = v[j]
            end
            vectors[1:ki-1,ki] .= zero(T)
        end

        # normalize
        t0 = abs1(vectors[1,ki])
        @inbounds for i=2:n; t0 = max(t0, abs1(vectors[i,ki])); end
        remax = one(RT) / t0
        @inbounds for i=1:n; vectors[i,ki] *= remax; end

        # restore diagonal
        @inbounds for k=ki+1:n
            TT[k,k] = tdiag[k]
        end
    end

    vectors
end


text:
function initializeexasim(version)

# create pde model
pde = Preprocessing.initializepde(version);

# create a mesh structure
mesh = Preprocessing.initializemesh(version);

return pde, mesh

end

ext: jl

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module SPMD

using Distributed

import DistributedArrays: gather, next_did, close
export sendto, recvfrom, recvfrom_any, barrier, bcast, scatter, gather
export context_local_storage, context, spmd

mutable struct WorkerDataChannel
    pid::Int
    rc::Union{RemoteChannel,Nothing}
    lock::ReentrantLock

    WorkerDataChannel(pid) = new(pid, nothing, ReentrantLock())
end

mutable struct SPMDContext
    id::Tuple
    chnl::Channel
    store::Dict{Any,Any}
    pids::Array
    release::Bool

    function SPMDContext(id)
        ctxt = new(id, Channel(typemax(Int)), Dict{Any,Any}(), [], false)
        finalizer(finalize_ctxt, ctxt)
        ctxt
    end
end

function finalize_ctxt(ctxt::SPMDContext)
    ctxt.release && close(ctxt)
end

function context_local_storage()
    ctxt = get_ctxt_from_id(task_local_storage(:SPMD_CTXT))
    ctxt.store
end

function context(pids=procs())
    global map_ctxts
    ctxt = SPMDContext(next_did())
    ctxt.pids = pids
    ctxt.release = true
    ctxt
end

# Every worker is associated with its own RemoteChannel
const map_worker_channels = Dict{Int, WorkerDataChannel}()

# mapping between a context id and context object
const map_ctxts = Dict{Tuple, SPMDContext}()

# Multiple SPMD blocks can be executed concurrently,
# each in its own context. Messages are still sent as part of the
# same remote channels associated with each worker. They are
# read from the remote channel into local channels each associated
# with a different run of `spmd`.

function get_dc(wc::WorkerDataChannel)
    lock(wc.lock)
    try
        if wc.rc === nothing
            if wc.pid == myid()
                myrc = RemoteChannel(()->Channel(typemax(Int)))
                wc.rc = myrc

                # start a task to transfer incoming messages into local
                # channels based on the execution context
                @async begin
                    while true
                        msg = take!(myrc)
                        ctxt_id = msg[1] # First element of the message tuple is the context id.
                        ctxt = get_ctxt_from_id(ctxt_id)
                        put!(ctxt.chnl, msg[2:end]) # stripping the context_id
                    end
                end
            else
                wc.rc = remotecall_fetch(()->get_remote_dc(myid()), wc.pid)
            end
        end
    finally
        unlock(wc.lock)
    end
    return wc.rc
end

function get_ctxt_from_id(ctxt_id)
    global map_ctxts
    ctxt = get(map_ctxts, ctxt_id, nothing)
    if ctxt == nothing
        ctxt = SPMDContext(ctxt_id)
        map_ctxts[ctxt_id] = ctxt
    end
    return ctxt
end

# Since modules may be loaded in any order on the workers,
# and workers may be dynamically added, pull in the remote channel
# handles when accessed for the first time.
function get_remote_dc(pid)
    global map_worker_channels
    if !haskey(map_worker_channels, pid)
        map_worker_channels[pid] = WorkerDataChannel(pid)
    end

    return get_dc(map_worker_channels[pid])
end

function send_msg(to, typ, data, tag)
    ctxt_id = task_local_storage(:SPMD_CTXT)
    @async begin
        dc = get_remote_dc(to)
        put!(dc, (ctxt_id, typ, myid(), data, tag))
#        println("Sent to ", dc)
    end
end

function get_msg(typ_check, from_check=false, tag_check=nothing)
    ctxt_id = task_local_storage(:SPMD_CTXT)
    chnl = get_ctxt_from_id(ctxt_id).chnl

    unexpected_msgs=[]
    while true
        typ, from, data, tag = take!(chnl)

        if (from_check != false && from_check != from) || (typ != typ_check) || (tag != tag_check)
            push!(unexpected_msgs, (typ, from, data, tag))
#            println("Unexpected in get_msg ", unexpected_msgs, " looking for ", typ_check, " ", from_check, " ", tag_check)
        else
            # put all the messages we read (but not expected) back to the local channel
            foreach(x->put!(chnl, x), unexpected_msgs)
            return (from, data)
        end
    end
end

function sendto(pid::Int, data::Any; tag=nothing)
    send_msg(pid, :sendto, data, tag)
end

function recvfrom(pid::Int; tag=nothing)
    _, data = get_msg(:sendto, pid, tag)
    return data
end

function recvfrom_any(; tag=nothing)
    from, data = get_msg(:sendto, false, tag)
    return (from,data)
end

function barrier(;pids=procs(), tag=nothing)
    # send a message to everyone
    for p in sort(pids)
        send_msg(p, :barrier, nothing, tag)
    end
    # make sure we recv a message from everyone
    pending=deepcopy(pids)
    unexpected_msgs=[]

    while length(pending) > 0
        from, _ = get_msg(:barrier, false, tag)
        if from in pending
            filter!(x->x!=from, pending)
        else
            # handle case of 2 (or more) consecutive barrier calls.
            push!(unexpected_msgs, (:barrier, from, nothing, tag))
#            println("Unexpected ", from)
        end
#        length(pending) == 1 && println("Waiting for ", pending)
    end

    ctxt_id = task_local_storage(:SPMD_CTXT)
    chnl = get_ctxt_from_id(ctxt_id).chnl
    foreach(x->put!(chnl, x), unexpected_msgs)
    return nothing
end

function bcast(data::Any, pid::Int; tag=nothing, pids=procs())
    if myid() == pid
        for p in filter(x->x!=pid, sort(pids))
            send_msg(p, :bcast, data, tag)
        end
        return data
    else
        from, data = get_msg(:bcast, pid, tag)
        return data
    end
end

function scatter(x, pid::Int; tag=nothing, pids=procs())
    if myid() == pid
        @assert rem(length(x), length(pids)) == 0
        cnt = div(length(x), length(pids))
        for (i,p) in enumerate(sort(pids))
            p == pid && continue
            send_msg(p, :scatter, x[cnt*(i-1)+1:cnt*i], tag)
        end
        myidx = findfirst(isequal(pid), sort(pids))
        return x[cnt*(myidx-1)+1:cnt*myidx]
    else
        _, data = get_msg(:scatter, pid, tag)
        return data
    end
end

function gather(x, pid::Int; tag=nothing, pids=procs())
    if myid() == pid
        gathered_data = Array{Any}(undef, length(pids))
        myidx = findfirst(isequal(pid), sort(pids))
        gathered_data[myidx] = x
        n = length(pids) - 1
        while n > 0
            from, data_x = get_msg(:gather, false, tag)
            fromidx = findfirst(isequal(from), sort(pids))
            gathered_data[fromidx] = data_x
            n=n-1
        end
        return gathered_data
    else
        send_msg(pid, :gather, x, tag)
        return x
    end
end

function spmd_local(f, ctxt_id, clear_ctxt)
    task_local_storage(:SPMD_CTXT, ctxt_id)
    f()
    clear_ctxt && delete_ctxt_id(ctxt_id)
    return nothing
end

function spmd(f, args...; pids=procs(), context=nothing)
    f_noarg = ()->f(args...)
    clear_ctxt = false
    if context == nothing
        ctxt_id = next_did()
        clear_ctxt = true    # temporary unique context created for this run.
                             # should be cleared at the end of the run.
    else
        ctxt_id = context.id
    end
    @sync for p in pids
        @async remotecall_fetch(spmd_local, p, f_noarg, ctxt_id, clear_ctxt)
    end
    nothing
end

function delete_ctxt_id(ctxt_id)
    global map_ctxts
    haskey(map_ctxts, ctxt_id) && delete!(map_ctxts, ctxt_id)
    nothing
end

function Base.close(ctxt::SPMDContext)
    for p in ctxt.pids
        remote_do(delete_ctxt_id, p, ctxt.id)
    end
    ctxt.release = false
end

end


text:
"""
    AbstractPointProcess{L, M}

The common supertype of all point processes with location type `L` and mark type `M`.
"""
abstract type AbstractPointProcess{L,M} <: AbstractMeasure end

"""
    build_transform(pp)

Return a transformation object from `TransformVariables` that can turn a `Vector{Float64}` into a `NamedTuple` with fields matching those of `pp`.
"""
function build_transform end


text:
"""
Functionality for authenticating Genie users.
"""
module Authentication

using SearchLight, App, Genie, Sessions

export current_user, current_user!!

const USER_ID_KEY = :__auth_user_id

"""
    authenticate(user_id::Union{String,Symbol,Int}, session) :: Sessions.Session
    authenticate(user_id::Union{String,Symbol,Int}, params::Dict{Symbol,Any}) :: Sessions.Session

Stores the user id on the session.
"""
function authenticate(user_id::Union{String,Symbol,Int}, session) :: Sessions.Session
  Sessions.set!(session, USER_ID_KEY, user_id)
end
function authenticate(user_id::Union{String,Symbol,Int}, params::Dict{Symbol,Any}) :: Sessions.Session
  authenticate(user_id, params[:SESSION])
end

"""
    deauthenticate(session) :: Sessions.Session
    deauthenticate(params::Dict{Symbol,Any}) :: Sessions.Session

Removes the user id from the session.
"""
function deauthenticate(session) :: Sessions.Session
  Sessions.unset!(session, USER_ID_KEY)
end
function deauthenticate(params::Dict{Symbol,Any}) :: Sessions.Session
  deauthenticate(params[:SESSION])
end

"""
    is_authenticated(session) :: Bool
    is_authenticated(params::Dict{Symbol,Any}) :: Bool

Returns `true` if a user id is stored on the session.
"""
function is_authenticated(session) :: Bool
  Sessions.is_set(session, USER_ID_KEY)
end
function is_authenticated(params::Dict{Symbol,Any}) :: Bool
  is_authenticated(params[:SESSION])
end

"""
    get_authentication(session) :: Nullable
    get_authentication(params::Dict{Symbol,Any}) :: Nullable

Returns the user id stored on the session, if available.
"""
function get_authentication(session) :: Nullable
  Sessions.get(session, USER_ID_KEY)
end
function get_authentication(params::Dict{Symbol,Any}) :: Nullable
  get_authentication(params[:SESSION])
end

"""
    login(user, session) :: Nullable{Sessions.Session}
    login(user, params::Dict{Symbol,Any}) :: Nullable{Sessions.Session}

Persists on session the id of the user object and returns the session.
"""
function login(user, session) :: Nullable{Sessions.Session}
  authenticate(getfield(user, Symbol(user._id)) |> Base.get, session) |> Nullable{Sessions.Session}
end
function login(user, params::Dict{Symbol,Any}) :: Nullable{Sessions.Session}
  login(user, params[:SESSION])
end

"""
    logout(session) :: Sessions.Session
    logout(params::Dict{Symbol,Any}) :: Sessions.Session

Deletes the id of the user object from the session, effectively logging the user off.
"""
function logout(session) :: Sessions.Session
  deauthenticate(session)
end
function logout(params::Dict{Symbol,Any}) :: Sessions.Session
  logout(params[:SESSION])
end

"""
    current_user(session) :: Nullable{User}
    current_user(params::Dict{Symbol,Any}) :: Nullable{User}

Returns the `User` instance corresponding to the currently authenticated user, wrapped into a Nullable.
"""
function current_user(session) :: Nullable{User}
  auth_state = Authentication.get_authentication(session)
  if isnull(auth_state)
    Nullable{User}()
  else
    SearchLight.find_one(User, Base.get(auth_state))
  end
end
function current_user(params::Dict{Symbol,Any}) :: Nullable{User}
  current_user(params[:SESSION])
end

"""
    current_user!!(session) :: User
    current_user!!(params::Dict{Symbol,Any}) :: User

Attempts to get the `User` instance corresponding to the currently authenticated user - throws error on failure.
"""
function current_user!!(session) :: User
  try
    current_user(session) |> Base.get
  catch ex
    Logger.log("The current user is not authenticated", :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    rethrow(ex)
  end
end
function current_user!!(params::Dict{Symbol,Any}) :: User
  current_user!!(params[:SESSION])
end

"""
    with_authentication(f::Function, fallback::Function, session)
    with_authentication(f::Function, fallback::Function, params::Dict{Symbol,Any})

Invokes `f` only if a user is currently authenticated on the session, `fallback` is invoked otherwise.
"""
function with_authentication(f::Function, fallback::Function, session)
  if ! is_authenticated(session)
    fallback()
  else
    f()
  end
end
function with_authentication(f::Function, fallback::Function, params::Dict{Symbol,Any})
  with_authentication(f, fallback, params[:SESSION])
end

"""
    without_authentication(f::Function, session)
    without_authentication(f::Function, params::Dict{Symbol,Any})

Invokes `f` if there is no user authenticated on the current session.
"""
function without_authentication(f::Function, session)
  ! is_authenticated(session) && f()
end
function without_authentication(f::Function, params::Dict{Symbol,Any})
  without_authentication(f, params[:SESSION])
end

"""
Sets up files and migrations needed by the Authentication module
"""
module FileTemplates

using Genie

function users_migration() :: String
  """
  module CreateTableUsers

  using Genie, SearchLight

  function up()
    SearchLight.query("CREATE SEQUENCE users__seq_id")
    SearchLight.query("
      CREATE TABLE IF NOT EXISTS users (
        id            integer CONSTRAINT users__idx_id PRIMARY KEY DEFAULT NEXTVAL('users__seq_id'),
        name          varchar(100) NOT NULL,
        email         varchar(100) NOT NULL,
        password      varchar(256) NOT NULL,
        role_id       integer DEFAULT NULL,
        CONSTRAINT users__idx_name UNIQUE(email)
      )
    ")
    SearchLight.query("ALTER SEQUENCE users__seq_id OWNED BY users.id")
    SearchLight.query("CREATE INDEX users__idx_role_id ON users (role_id)")
  end

  function down()
    SearchLight.query("DROP INDEX users__idx_role_id")
    SearchLight.query("DROP INDEX users__seq_id")
    SearchLight.query("DROP TABLE users")
  end

  end
  """
end

function roles_migration() :: String
  """
  module CreateTableRoles

  using Genie, SearchLight

  function up()
    SearchLight.query("CREATE SEQUENCE roles__seq_id")
    SearchLight.query("
      CREATE TABLE IF NOT EXISTS roles (
        id            integer CONSTRAINT roles__idx_id PRIMARY KEY DEFAULT NEXTVAL('roles__seq_id'),
        name          varchar(20) NOT NULL,
        CONSTRAINT    roles__idx_id UNIQUE(id),
        CONSTRAINT    roles__idx_name UNIQUE(name)
      )
    ")
    SearchLight.query("ALTER SEQUENCE roles__seq_id OWNED BY roles.id")
  end

  function down()
    SearchLight.query("DROP INDEX roles__seq_id")
    SearchLight.query("DROP TABLE roles")
  end

  end
  """
end

function routes() :: String
  """
  # Authentication
  route("/login", "user_sessions#UserSessionsController.show_login", named = :show_login)
  route("/login", "user_sessions#UserSessionsController.login", method = POST, named = :login)
  route("/logout", "user_sessions#UserSessionsController.logout", named = :logout)
  """
end

end # module FileTemplates

"""
File generation functionality for the Authentication module.
"""
module Generator

using Migration, Authentication, Logger, Genie, Router

function setup() :: Void
  Logger.log("Creating migrations")
  Migration.new("create_table_users", Authentication.FileTemplates.users_migration())
  Migration.new("create_table_roles", Authentication.FileTemplates.roles_migration())

  Logger.log("Creating paths")
  try
    mkpath(joinpath(Genie.APP_PATH, "layouts"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "layouts")) already exists. Skipping.", :warn)
  end
  try
    mkpath(joinpath(Genie.APP_PATH, "resources", "roles"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "resources", "roles")) already exists. Skipping.", :warn)
  end
  try
    mkpath(joinpath(Genie.APP_PATH, "resources", "users"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "resources", "users")) already exists. Skipping.", :warn)
  end
  try
    mkpath(joinpath(Genie.APP_PATH, "resources", "user_sessions"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "resources", "user_sessions")) already exists. Skipping.", :warn)
  end

  Logger.log("Copying files")
  try
    cp(joinpath(Pkg.dir("Genie"), "files", "authentication", "app", "layouts", "login.flax.html"), joinpath(Genie.APP_PATH, "layouts", "login.flax.html"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "layouts", "login.flax.html")) already exists. Skipping.", :warn)
  end
  try
    cp(joinpath(Pkg.dir("Genie"), "files", "authentication", "app", "resources", "roles", "model.jl"), joinpath(Genie.APP_PATH, "resources", "roles", "model.jl"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "resources", "roles", "model.jl")) already exists. Skipping.", :warn)
  end
  try
    cp(joinpath(Pkg.dir("Genie"), "files", "authentication", "app", "resources", "users", "model.jl"), joinpath(Genie.APP_PATH, "resources", "users", "model.jl"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "resources", "users", "model.jl")) already exists. Skipping.", :warn)
  end
  try
    cp(joinpath(Pkg.dir("Genie"), "files", "authentication", "app", "resources", "user_sessions", "controller.jl"), joinpath(Genie.APP_PATH, "resources", "user_sessions", "controller.jl"))
  catch ex
    Logger.log("Destination $(joinpath(Genie.APP_PATH, "resources", "user_sessions", "controller.jl")) already exists. Skipping.", :warn)
  end

  Logger.log("Appending routes")
  Router.append_to_routes_file(Authentication.FileTemplates.routes())

  Logger.log("Success!")

  nothing
end

end # module Generator

end


text:
module Soss

import Base.rand
using Reexport: @reexport

@reexport using Distributions
@reexport using StatsFuns

import MacroTools: prewalk, postwalk, @q, striplines, replace, flatten, @capture
import MLStyle
@reexport using MonteCarloMeasurements

include("statement.jl")
include("model.jl")
# include("weighted.jl")
include("rand.jl")
include("for.jl")
include("dist.jl")
include("flat.jl")
include("iid.jl")
include("utils.jl")
include("examples.jl")
include("graph.jl")
include("nuts.jl")
# include("optim.jl")
include("importance.jl")
include("canonical.jl")
include("symbolic.jl")
# include("sobols.jl")
# include("fromcube.jl")
# include("tocube.jl")
include("particles.jl")
include("xform.jl")
include("toposort.jl")
include("advancedhmc.jl")
include("plots.jl")
include("rejection.jl")
end # module


text:
module Template
    export template

    const cache = Dict{AbstractString, Function}()

    """
    low-level API; template string + argnames -> (args -> rendered string)
    """
    function compile(x::AbstractString, argnames::Symbol...) # TODO: escape """, auto find free variables
        raw = split(x, "\$\$")
        buf = IOBuffer(sizehint = sizeof(x) + length(raw)*32 + 128)

        arglist = join(("$x=throw(ArgumentError(\"missing required argument $x\"))" for x in argnames), ',')

        write(buf, "(_o42CzdTO::IO; ", arglist, ") -> begin\n")

        iscode = false # indicating if i is code
        for i in raw
            if iscode
                write(buf, '\n', i, '\n')
            else
                write(buf, "write(_o42CzdTO, \"\"\"", i, "\"\"\")")
            end

            iscode = !iscode
        end

        write(buf, "\n_o42CzdTO\nend")

        buf |> take! |> String |> Meta.parse |> eval
    end

    """
    render with cache; automatically determin argnames by the arguments of first call
    filepath and args -> rendered string
    """
    function render(file::AbstractString; kwargs...)
        f = Base.@get!(cache, file, compile(read(file, String), keys(kwargs)...))
        Base.invokelatest(f, IOBuffer(); kwargs...)::IOBuffer |> take! |> String
    end

    function template(req, res, route)
        res[:render] = this -> (file; kwargs...) -> begin
            this.status = 200
            this.body = render(file; kwargs...)
            this.content_type = "text/html"
        end

        route.next()
    end
end


text:
import TimeZones: TimeZone, localzone
import Compat: is_linux

# Ensure that the current system's local time zone is supported. If this test fails make
# sure to report it as an issue.
@test isa(localzone(), TimeZone)

if is_linux()
    # Bad TZ environmental variables
    withenv("TZ" => "") do
        @test_throws ErrorException localzone()
    end
    withenv("TZ" => "Europe/Warsaw") do
        @test_throws ErrorException localzone()
    end

    # Absolute filespec
    warsaw_path = joinpath(TZFILE_DIR, "Europe", "Warsaw")
    warsaw_from_file = open(warsaw_path) do f
        TimeZones.read_tzfile(f, "local")
    end
    withenv("TZ" => ":" * abspath(warsaw_path)) do
        @test localzone() == warsaw_from_file
    end

    # Relative filespec
    warsaw = TimeZone("Europe/Warsaw")
    withenv("TZ" => ":Europe/Warsaw") do
        @test localzone() == warsaw
    end

    # Set TZDIR and use time zone unrecognized by TimeZone
    @test_throws ArgumentError TimeZone("Etc/UTC")
    utc = open(joinpath(TZFILE_DIR, "Etc", "UTC")) do f
        TimeZones.read_tzfile(f, "Etc/UTC")
    end
    withenv("TZ" => ":Etc/UTC", "TZDIR" => TZFILE_DIR) do
        @test localzone() == utc
    end

    # Use system installed files
    @test_throws ArgumentError TimeZone("Etc/GMT-9")
    gmt_minus_9 = FixedTimeZone("Etc/GMT-9", 9 * 3600)
    withenv("TZ" => ":Etc/GMT-9") do
        @test localzone() == gmt_minus_9
    end

    # Unable to locate time zone on system
    withenv("TZ" => ":") do
        @test_throws SystemError localzone()
    end
    withenv("TZ" => ":Etc/Foo") do
        @test_throws SystemError localzone()
    end
end


text:
# test/runtests.jl
# Abstract Wallets Testing

using AbstractWallets


text:
using SubModule
using SubModule2

touch(joinpath(@__DIR__, "buildartifact"))

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function assemble_mk(iop::IntegralOperator)
    assert_mk_compatibility(iop)
    X = target_surface(iop)
    Y = source_surface(iop)
    N   = length(X) # number of dof
    out = zero(iop) # output matrix to be filled
    for target_ent in entities(X)
        for source_ent in entities(Y)
            if source_ent == target_ent
                _assemble_mk_self!(out,iop,target_ent)
            else
                # nothing special to do here, just fill up the matrix using the
                # regular quadrature
                target_dof  = dom2dof(X,target_ent)
                source_dof  = dom2dof(Y,source_ent)
                @. out[target_dof,source_dof] = iop[target_dof,source_dof]
            end
        end
    end
    return out
end

# Assemble the self-interaction part for a given entity
function _assemble_mk_self!(out,iop,ent::AbstractEntity)
    X         = target_surface(iop)
    etype2tag = X.ent2elt[ent]
    # loop over individual elements composing ent. Note that because we must
    # handle meshes composed of elements of possiby different types, we must
    # loop over element types, and then over the elements of that type
    # (contained in the entity)
    for (target_E,target_tags) in etype2tag, target_t in target_tags
        for (source_E,source_tags) in etype2tag, source_t in source_tags
            if source_E === target_E && target_t === source_t # same element
                dof = X.elt2dof[source_E][:,source_t] # dof for filling out in global matrix
                el   = X.elements[source_E][source_t]                 # retrieve the element (i.e. access to its parametrization)
                _assemble_mk_self!(out,iop,el,dof)
            else
                target_dof  = X.elt2dof[target_E][:,target_t] # dof for filling out in global matrix
                source_dof  = X.elt2dof[source_E][:,source_t] # dof for filling out in global matrix
                out[target_dof,source_dof] = iop[target_dof,source_dof]
            end
        end
    end
    return out
end

function _assemble_mk_self!(out,iop,el::AbstractElement,idxs)
    X     = target_surface(iop)
    qrule = X.etype2qrule[typeof(el)] # quadrature used on elements of type el
    x̂,ŵ   = qrule()                   # nodes on reference segment [0,1]
    # compute all quantities needed for mk rule, which includes high order
    # derivatives of the parametrization.
    N     = length(idxs)
    dx    =  [derivative(el,u)/(2π)     for u in x̂]
    d2x   =  [derivative2(el,u)/(4*π^2) for u in x̂]
    τ     = [norm(τ) for τ in dx]
    Δs    = 2π / N
    psi   = -0.577215664901532
    x     = qnodes(X)
    ν     = qnormals(X)
    w     = qweights(X)
    R     = nystrom_weights(N)
    k     = iop.kernel.op.k
    K     = kernel(iop)
    φ     = mk_split(K)
    for (iloc,iglob) in enumerate(idxs)
        for (jloc,jglob) in enumerate(idxs)
            lSin = log(4*(sin(Δs*(iloc-jloc)/2))^2) # singular part factored
            if kernel_type(iop) === SingleLayer()
                K1   = φ(x[iglob],x[jglob])*w[jglob]
                if iloc != jloc
                    K2   = K(x[iglob],x[jglob])*w[jglob]-K1*lSin # what is left after factoring out lSin and φ
                else
                    K2   = (im/4+psi/2/pi-1/4/pi*log(k^2/4*τ[jloc]^2))*w[jglob] + 2*log(w[jglob]/(τ[jloc]*Δs))*K1
                end
            elseif kernel_type(iop) == DoubleLayer()
                K1   = φ(x[iglob],x[jglob],ν[jglob])*w[jglob]
                if iloc != jloc
                    K2   = K(x[iglob],x[jglob],ν[jglob])*w[jglob]-K1*lSin
                else
                    K2   = -1/(4*pi)*(dx[jloc][1]*d2x[jloc][2]-dx[jloc][2]*d2x[jloc][1])/τ[jloc]^3*w[jglob] + 2*log(w[jglob]/(τ[jloc]*Δs))*K1
                end
            end
            out[iglob,jglob] = (R[iloc,jloc]*K1 + K2)
        end
    end
    return out
end

# analytical splitting of kernel in 2d Helmholtz
function mk_split(SL::SingleLayerKernel{T,S}) where {T,S<:Helmholtz}
    k = SL.op.k
    ϕ = (x,y) -> begin
        d = norm(x-y)
        (-1/(4*pi))*besselj0(k*d) |> T
    end
    return ϕ
end

function mk_split(DL::DoubleLayerKernel{T,S}) where {T,S<:Helmholtz}
    k = DL.op.k
    ϕ = (x,y,ν) -> begin
        x == y && (return zero(T)) # analitic limit taken by hand to avoid 0/0
        r = x-y
        d = norm(r)
        (-k/(4*pi))*besselj1(k*d)/d*dot(r,ν) |> T
    end
    return ϕ
end

function nystrom_weights(M)
    @assert M%2 == 0  "number of points `M` must be even"
    Md2 = div(M,2)
    R   = zeros(M,M)
    for p=1:M
        tp = pi/Md2*(p-1)
        for j=1:M
            tj = pi/Md2*(j-1)
            R[p,j]=-2*sum((ones(Md2-1)./collect(1:Md2-1)).*cos.(collect(1:Md2-1)*(tp-tj)))-pi/(Md2^2)*cos(Md2*(tp-tj))
        end
    end
    return R
end

function assert_mk_compatibility(iop::IntegralOperator)
    # check 2d
    @assert ambient_dimension(iop) == 2 "Martensen-Kussmaul method only applicable to two-dimensional problems"
    # FIXME: add other checks...
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
function check_grad(g_gpu, g_cpu, atol, rtol)
  @show g_gpu g_cpu
  @test false
end
check_grad(g_gpu::Base.RefValue, g_cpu::Base.RefValue, atol, rtol) =
    check_grad(g_gpu[], g_cpu[], atol, rtol)
check_grad(g_gpu::Nothing, g_cpu::Nothing, atol, rtol) = @test true
check_grad(g_gpu::Float32, g_cpu::Float32, atol, rtol) = @test g_cpu ≈ g_gpu   rtol=rtol atol=atol
check_grad(g_gpu::CuArray{Float32}, g_cpu::Array{Float32}, atol, rtol) =
  @test g_cpu ≈ collect(g_gpu)   rtol=rtol atol=atol

function check_grad(g_gpu::Tuple, g_cpu::Tuple, atol, rtol)
  for (v1, v2) in zip(g_gpu, g_cpu)
    check_grad(v1, v2, atol, rtol)
  end
end

function check_grad(g_gpu::NamedTuple, g_cpu::NamedTuple, atol, rtol)
  for ((k1,v1), (k2,v2)) in zip(pairs(g_gpu), pairs(g_cpu))
    @test k1 == k2
    # @show k2 v2
    check_grad(v1, v2, atol, rtol)
  end
end

function gpu_autodiff_test(f_cpu, xs_cpu::Array{Float32}...;
        test_equal=true, rtol=1e-4, atol=1e-4)

  check_type(x) = false
  check_type(x::Float32) = true
  check_type(x::CuArray{Float32}) = true
  check_type(x::Array{Float32}) = true

  ### GRADIENT WITH RESPECT TO INPUT #####
  # y_cpu, back_cpu = pullback((f, x...) -> f(x...), f_cpu, xs_cpu...)
  y_cpu, back_cpu = pullback((x...) -> f_cpu(x...), xs_cpu...)
  @test check_type(y_cpu)
  Δ_cpu = size(y_cpu) == () ? randn(Float32) : randn(Float32, size(y_cpu))
  gs_cpu = back_cpu(Δ_cpu)

  f_gpu = f_cpu |> gpu
  xs_gpu = gpu.(xs_cpu)
  Δ_gpu = Δ_cpu |> gpu
  # y_gpu, back_gpu = pullback((f, x...) -> f(x...), f_gpu, xs_gpu...)
  y_gpu, back_gpu = pullback((x...) -> f_gpu(x...), xs_gpu...)
  @test check_type(y_gpu)
  gs_gpu = back_gpu(Δ_gpu)

  if test_equal
    @test collect(y_cpu) ≈ collect(y_gpu)   rtol=rtol atol=atol
    for (g_gpu, g_cpu) in zip(gs_gpu, gs_cpu)
      check_grad(g_gpu, g_cpu, atol, rtol)
    end
  end

  ### GRADIENT WITH RESPECT TO f #####
  ps_cpu = Flux.params(f_cpu)
  y_cpu, back_cpu = pullback(() -> f_cpu(xs_cpu...), ps_cpu)
  gs_cpu = back_cpu(Δ_cpu)

  ps_gpu = Flux.params(f_gpu)
  y_gpu, back_gpu = pullback(() -> f_gpu(xs_gpu...), ps_gpu)
  gs_gpu = back_gpu(Δ_gpu)

  if test_equal
    @test collect(y_cpu) ≈ collect(y_gpu)   rtol=rtol atol=atol
    @assert length(ps_gpu) == length(ps_cpu)
    for (p_gpu, p_cpu) in zip(ps_gpu, ps_cpu)
      check_grad(gs_gpu[p_gpu], gs_cpu[p_cpu], atol, rtol)
    end
  end
end


text:
export LBFGS

mutable struct LBFGS{T <: Real} <: AbstractOptimizationMethod{T}
    f_tol::T
    g_tol::T
    iterations::Int
    initial_points::Vector{Vector{T}}
end

"""
    LBFGS(model::Model, args...; kwargs...)

If an `Int` is provided the method will sample that many random initial_points and use them as
initial points for Optim LBFGS method. If a `Vector{Vector{T}}` is provided it will use them as
initial points for Optim LBFGS method.
"""
function LBFGS(model::Model{D, T}, n_initial_points::Int; f_tol::T = T(1e-6), g_tol::T = T(1e-6),
                                               iterations::Int = 10^5, LB::T = 0.0, UB::T = 0.6) where {D, T}

    initial_points = create_initial_points(model, n_initial_points, LB, UB)

    return LBFGS{T}(f_tol, g_tol, iterations, initial_points)
end

function LBFGS(model::Model{D, T}, initial_points::Vector{Vector{T}}; f_tol::T = T(1e-6), g_tol::T = T(1e-6),
                                                           iterations::Int = 10^5) where {D, T}

    ensure_seeds_dimensions(model, initial_points)

    return LBFGS{T}(f_tol, g_tol, iterations, initial_points)
end

function optimize(func::Optim.TwiceDifferentiable, opt_method::LBFGS{T}, verbose::Int, i::Int) where T
    return optimize(func, opt_method, Optim.LBFGS(), verbose, i)
end


text:
using PyPlot, MRIReco

# load fully sampled data
f = ISMRMRDFile(@__DIR__()*"/data/knee_3dFSE_slice170.h5")
acqData = AcquisitionData(f);

# reconstruct
params = Dict{Symbol, Any}()
params[:reco] = "direct"
params[:reconSize] = (320,320) # this size is also contained in acqData.encodingSize

img = reconstruction(acqData, params)
img = sqrt.(sum(img.^2,dims=5))

# To simulate an undersampled reconstruction, we retrospectively undersample the data using a Poisson disk pattern.

redFac = 4.0
acqDataSub = sample_kspace(acqData,redFac,"poisson",calsize=30,profiles=false);

# show sampling pattern
msk = zeros(acqDataSub.encodingSize[1],acqDataSub.encodingSize[2])
msk[acqDataSub.subsampleIndices[1]] .= 1

# Estimate the coil sensitivities using ESPIRiT and reconstruct using SENSE

# coil sensitivities
smaps = espirit(acqData,(6,6),30,eigThresh_1=0.035,eigThresh_2=0.98)

# SENSE reconstruction
params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = (320,320)
params[:senseMaps] = smaps

params[:solver] = "cgnr"
params[:regularization] = "L2"
params[:λ] = 1.e-4
params[:iterations] = 5
params[:normalizeReg] = true

img_cg = reconstruction(acqDataSub, params)

# Using TV regularization recquires us to change some parameters.

params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = (320,320)
params[:senseMaps] = smaps

params[:solver] = "admm"
params[:regularization] = "TV"
params[:λ] = 1.e-1 # 5.e-2
params[:iterations] = 50
params[:ρ] = 0.1
params[:absTol] = 1.e-4
params[:relTol] = 1.e-2
params[:tolInner] = 1.e-2
params[:normalizeReg] = true

img_tv = reconstruction(acqDataSub, params)

# use PyPlot for interactive display
figure(1)
clf()
subplot(2,2,1)
title("Phantom")
imshow(abs.(img[:,:,1,1,1]))
subplot(2,2,2)
title("Mask")
imshow(abs.(msk))
subplot(2,2,3)
title("CG Reconstruction")
imshow(abs.(img_cg[:,:,1,1,1]))
subplot(2,2,4)
title("TV Reconstruction")
imshow(abs.(img_tv[:,:,1,1,1]))

# export images
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/kneeOrig.png")
exportImage(filename, abs.(img[:,:,1,1,1]) )
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/mask.png")
exportImage(filename, abs.(msk[:,:,1,1,1]) )
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/kneeCG.png")
exportImage(filename, abs.(img_cg[:,:,1,1,1]) )
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/kneeTV.png")
exportImage(filename, abs.(img_tv[:,:,1,1,1]) )


text:
@testset "set length and range" begin
    ks = 1:10
    rs = set_range.(Chebyshev(), ks)
    @test rs[1:4] == [1:1, 2:3, 4:5, 6:9]
    @test all(@. last(rs[1:(end - 1)]) + 1 == first(rs[2:end])) # contiguous
    @test length.(rs) == set_length.(Chebyshev(), ks)           # consistent with length
end

@testset "gridpoint sets" begin
    @test gridpoint_set(Chebyshev(), 1) == [0]
    @test gridpoint_set(Chebyshev(), 2) == [-1, 1]
    @test gridpoint_set(Chebyshev(), 3) ≈ [-√2, √2]./2
    @test issorted(gridpoint_set(Chebyshev(), 4))
    x = sort(mapreduce(k -> gridpoint_set(Chebyshev(), k), vcat, 1:5))
    @test allunique(x)
    @test length(x) == 2^4 + 1
    @test x ≈ cospi.((16:-1:0) ./ 16)
end

@testset "Chebyshev polynomials" begin
    J = 19
    x = rand(J) .* 2 .- 1
    N = 7
    Z = @inferred evaluate(N, Chebyshev(), x)
    @test size(Z) == (J, N)
    @test Z ≈ ((t, n) -> cos(n * t)).(acos.(x), (0:(N-1))')
    @test @inferred(evaluate(N, Chebyshev(), first(x))) == Z[1, :]

    @test size(evaluate(0, Chebyshev(), x)) == (J, 0)
    @test evaluate(1, Chebyshev(), x) == ones(J, 1)
    @test evaluate(2, Chebyshev(), x) == hcat(ones(J, 1), x)
end


text:
using StrainRecon
using Base.Test

@testset "StrainRecon" begin
	println("=== testing StrainRecon.jl code ===")
	println("")
	include("testFindw.jl")
	include("testFindM.jl")
	include("testSolveMAP.jl")
	include("testSolveConvexMAP.jl")
	include("testReducedMisfit.jl")
	include("testRandomw.jl")
	include("testRandomM.jl")
	include("testBCD.jl")
	include("testIntegratePosterior.jl")
	include("testASQP.jl")
	include("testIPLP.jl")
	include("testEntropyM.jl")
	include("testSrDistance.jl")
	include("testIntmat2binmat.jl")
end


text:
tree = BitVector([1,1,1])
@test Visualizations.treenodes_matrix(tree) == BitArray([1 1; 1 1])
# plot_tfbdry() test
@test typeof(plot_tfbdry(tree)) == Plots.Plot{Plots.GRBackend}
@test_nowarn plot_tfbdry(tree, start=1)
@test_throws AssertionError plot_tfbdry(tree, start=2)
@test_nowarn plot_tfbdry(tree, nd_col=:red, bg_col=:white, ln_col=:black)
x = randn(16,5)
# wiggle() test
@test typeof(wiggle(x)) == Plots.Plot{Plots.GRBackend}
p = plot()
@test typeof(wiggle(p, x)) == Plots.Plot{Plots.GRBackend}
# wiggle!() test
@test typeof(wiggle!(x)) == Plots.Plot{Plots.GRBackend}
@test typeof(wiggle!(p, x)) == Plots.Plot{Plots.GRBackend}
@test_throws ErrorException("Inconsistent taxis dimension!") wiggle!(p, x, taxis=1:5)
@test_throws ErrorException("Inconsistent zaxis dimension!") wiggle!(p, x, zaxis=1:4)
@test_nowarn wiggle!(p, x, sc=0.3)
@test_nowarn wiggle!(p, x, EdgeColor=:red)
@test_nowarn wiggle!(p, x, FaceColor=:red)
@test_nowarn wiggle!(p, x, Overlap=false)
@test_nowarn wiggle!(p, x, Orient=:down)
@test_throws AssertionError wiggle!(p, x, Orient=:fail)
@test_nowarn wiggle!(p, x, ZDir=:reverse)
@test_throws AssertionError wiggle!(p, x, ZDir=:fail)


text:
const American = OptionExercise("American")
const European = OptionExercise("European")

const USD = Currency("USD")
const GBP = Currency("GBP")
const EUR = Currency("EUR")
const AUD = Currency("AUD")
const NZD = Currency("NZD")
const JPY = Currency("JPY")

const EURUSD = CurrencyPair(EUR,USD)
const EURJPY = CurrencyPair(EUR,JPY)
const AUDUSD = CurrencyPair(AUD,USD)
const NZDUSD = CurrencyPair(NZD,USD)
const GBPUSD = CurrencyPair(GBP,USD)
const USDJPY = CurrencyPair(USD,JPY)

const F = "Jan"
const G = "Feb"
const H = "Mar"
const J = "Apr"
const K = "May"
const M = "Jun"
const N = "Jul"
const Q = "Aug"
const U = "Sep"
const X = "Oct"
const V = "Nov"
const Z = "Dec"


text:
@testset "assemble_flux_matrices.jl" begin

end


text:
using IntervalWavelets
using OffsetArrays
using Test

@testset "Interior scaling functions" begin
    @testset "Dyadic dilation matrix has the expected form" begin
        h1 = interior_filter(3, :min)
        H1 = IntervalWavelets.dyadic_dilation_matrix(h1)

        @test size(H1) == (length(h1) - 2, length(h1) - 2)

        @test H1 == sqrt(2) *
        [ h1[1] h1[0]    0     0  ;
          h1[3] h1[2] h1[1] h1[0] ;
          h1[5] h1[4] h1[3] h1[2] ;
             0     0  h1[5] h1[4] ]

        h2 = interior_filter(3, :symmlet)
        H2 = IntervalWavelets.dyadic_dilation_matrix(h2)

        @test H2 == sqrt(2) *
        [ h2[-1] h2[-2]      0      0  ;
          h2[1]  h2[0]   h2[-1] h2[-2] ;
          h2[3]  h2[2]   h2[1]  h2[0] ;
             0      0    h2[3]  h2[2] ]
    end

    @testset "Construct interior scaling function" begin
        h = interior_filter(2)
        phi = interior_scaling_function(h)

        @test collect(support(h)) == Integer.(support(phi))
        @test vanishing_moments(phi) == 2
        @test resolution(phi) == 0

        phi = interior_scaling_function(h, 1)

        @test vanishing_moments(phi) == 2
        @test resolution(phi) == 1

        @test_throws DomainError interior_scaling_function(h, -1)
    end

    @testset "Specific interior scaling functions" begin
        h = interior_filter(2)
        phi = interior_scaling_function(h, 1)

        @test phi.(support(phi)) ≈ [0.0 ; 0.93301 ; 1.36602 ; 0.0 ; -0.36602 ; 0.06698 ; 0.0] atol = 10.0^-4
    end

    @testset "Haar scaling functions" begin
        haar_filter = interior_filter(1, :min)

        haar = interior_scaling_function(haar_filter)
        @test haar.(support(haar)) == [1.0 ; 0.0]

        haar = interior_scaling_function(haar_filter, 1)
        @test haar.(support(haar)) ≈ [1.0 ; 1.0 ; 0.0]
    end

    @testset "Increase resolution of interior scaling function" begin
        h = interior_filter(2)
        phi0 = interior_scaling_function(h)

        phi0_support = support(phi0)
        @test phi0_support == IntervalWavelets.all_dyadic_rationals(-1, 2, 0)

        phi1 = IntervalWavelets.increase_resolution(phi0)
        @test support(phi1) == IntervalWavelets.all_dyadic_rationals(-1, 2, 1)

        @test phi1.(phi0_support) == phi0.(phi0_support)
    end
end


text:
using Bridge, Distributions, StaticArrays
using ConjugatePriors

PLOT = :winston
include("plot.jl")

Random.seed!(10)

diag2(x, y) = SDiagonal(@SVector [x,y])

struct FitzHughNagumo  <: ContinuousTimeProcess{SVector{2,Float64}}
    α::Float64
    β::Float64
    γ1::Float64
    γ2::Float64
    ϵ::Float64
    s::Float64
    σ::SDiagonal{2,Float64}
    a::SDiagonal{2,Float64}
#    σ::SMatrix{2,2,Float64}
#    a::SMatrix{2,2,Float64}

    FitzHughNagumo(σ, α = 1/3,  β = 0.08*0.7, γ1 = 0.08, γ2 = 0.08*0.8,  ϵ = 1.,  s = 1.) = new(α, β, γ1, γ2, ϵ, s, σ, σ*σ') #'
end
#phi(t, x, P::FitzHughNagumo) = (@SMatrix [-x[1]^3 x[1]-x[2] 1. 0. 0. 0. ; 0. 0. 0. x[1] -x[2] 1.])
phi(t, x, P::FitzHughNagumo) = (@SMatrix [0.  -x[1]^3 + x[1]-x[2] 1. ; x[1]  0.  0. ]) #γ, ϵ, s

Bridge.b(t, x, P::FitzHughNagumo) = (@SVector [-P.α*x[1]^3+ P.ϵ*(x[1]-x[2]) + P.s, P.γ1*x[1]- P.γ2*x[2] + P.β])
Bridge.σ(t, x, P::FitzHughNagumo) = P.σ
Bridge.a(t, x, P::FitzHughNagumo) = P.a
Bridge.Γ(t, x, P::FitzHughNagumo) = inv(P.a)

function conjugateb(YY, th, xi, P)
        n = length(th)
        G = zeros(n,n)
        mu = zeros(th)
        for Y in YY
            for i in 1:length(Y)-1
                phii = phi(Y[i]..., P)
                Gphii = Bridge.Γ(Y[i]..., P)*phii
                zi = phii'*Gphii
                dy = Y.yy[i+1] - Y.yy[i]
                mu = mu + Vector(Gphii'*dy)
                G[:] = G +  Matrix(zi * (Y.tt[i+1]-Y.tt[i]))
            end
        end #for m
        WW = G .+ diagm(xi)
        WL = transpose(chol(WW))
        th° = WL'\(randn(n)+WL\mu)
end

mcstart(yy) = (zeros(yy), zeros(yy), 0)
function mcnext(mc, x)
    m, m2, n = mc
    delta = x - m
    n = n + 1
    m = m + delta*(1/n)
    m2 = m2 + map((x,y)->x.*y, delta, (x - m))
    m, m2, n
end
function mcbandste(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    ste = eltype(m)[sqrt(v) for v in m2 * (1/(k - 1))]*sqrt(1/k)
    m -Q*ste, m +Q*ste
end
function mcband(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    std = eltype(m)[sqrt.(v) for v in m2 * (1/(k - 1))]
    m-Q*std, m+Q*std
end

function param(θ, σ)
    β, γ, ϵ, s = θ
    σ1,σ2 = σ
#   σ,    α, β, γ1,γ2, ϵ,  s
    diag2(σ1,σ2), ϵ, β, γ, 1., ϵ, s
end
simid = 2
simname =["exci", "nonexci"][simid]
mkpath(joinpath("output",simname))
try
    cp(@__FILE__(), joinpath("output",simname,"$simname.jl"); force=true)
catch
end

θ =  [[0.6, 1.4][simid], 1.5, 10., 0.5*10] # [β, γ, ϵ, s]
θtrue=copy(θ)
scaleθ = 0.08*[0., 1., 5., 1.]
σ = [0.3, 0.5] #[σ1]
scaleσ = [0.01]
σtrue = copy(σ)
Ptrue = FitzHughNagumo(param(θ, σ)...)
si = sitrue = 0.1 # observation error

n = 100 # number of segments
m = 100 # number of euler steps per segments
TT = 20.
dt = TT/n/m
tt = range(0., stop=TT, length=n*m+1)
tttrue = range(0., stop=TT, length=n*10*m+1)
ttf = tt[1:m:end]

uu = (@SVector [0., 1.])
Y = euler(uu, sample(tttrue, Wiener{SVector{2,Float64}}()), Ptrue)

Ytrue = Y[1:10m:end] #subsample
assert(endof(Ytrue) == n+1)

L = @SMatrix [ 1. 0. ]
Yobs = SamplePath{SVector{1,Float64}}(Ytrue.tt, [L*Ytrue.yy[i] + si*randn(SVector{1,Float64}) for i in 1:length(Ytrue.yy)])

r = [(:xrange,(-2,2)), (:yrange,(-1,3))]

if false
    r = [(:xrange,(-2,2)), (:yrange,(-1,3))]

    plot(Y, xrange=(-4, 4), yrange=(-4,4))
    plot(Yobs, "o:")
    plot(Yobs[1:2],"o:")

    i = 4
    v = SVector(Matrix(L)\Vector(Yobs.yy[i+1]))
    P° = PBridgeProp(Ptrue, Ytrue[i]..., Yobs.tt[i+1], v, Ytrue[i+2]..., L,si^2*I, Ptrue.a)
    B = bridge(sample(tt[m*(i-1)+1:m*(i+1)+1] , Wiener{SVector{2,Float64}}()),P°)
    # B = eulerb(sample(tt[m*(i-1)+1:m*(i+1)+1] , Wiener{SVector{2,Float64}}()),P°)
    plot(B; r...)
    oplot(Ytrue[i:i+2],"ro:"; r...)
    oplot([Yobs.yy[i+1][1],Yobs.yy[i+1][1]],[-2, 2], ":"; r...)
    oplot(Y[1:2m+1], "b"; r...)
end

#Y2 = SamplePath(tt, copy(Y.yy[1:10:end])) #initialize with truth
Y2 = SamplePath{SVector{2,Float64}}(tt, map(x->x+uu,zeros(SVector{2,Float64}, length(tt)))) #initialize path with constant
Yfil = Y2[1:m:end]

θs = Float64[]
iter = 0
tt2 = [collect(tt[m*(i-1)+1:m*(i+1)+1]) for i in 1:n-1]
BB = SamplePath[Y2[m*(i-1)+1:m*i+1] for i in 1:n]
BBnew = SamplePath[Y2[m*(i-1)+1:m*i+1] for i in 1:n]

# arbitrary starting values
si = 0.3
σ = [0.7, 0.7]
θ = 0.0 + 0.0*rand(length(θ))

for i in 1:length(θ) # start with truth for parameters not to be estimated
    if scaleθ[i] == 0
        θ[i] = θtrue[i]
    end
end

P = FitzHughNagumo(param(θ, σ)...)

if false
    for i in 1:n # reparametrization
        v = SVector(Matrix(L)\Vector(Yobs.yy[i+1]))
        P° = PBridgeProp(P, Ytrue[i]..., Yobs.tt[i+1], v, 1000., (@SVector [0.,0.]), L,  si^2*I, P.a)
        BB[i] = bridge(sample(tt[m*(i-1)+1:m*i+1], Wiener{SVector{2,Float64}}()),P°)
        # BB[i] = eulerb(sample(tt[m*(i-1)+1:m*i+1], Wiener{SVector{2,Float64}}()),P°)
        Ytrue.yy[i+1] = BB[i].yy[m+1]
    end
end
plot2(Y, "r","r" ; yrange=(-3,3),linewidth=0.5)
display(oplot2(Y2, "b", "b"; yrange=(-3,3), linewidth=0.7))

i = 1

# reserve some space
ww = Array{SVector{2,Float64},1}(length(m*(i-1)+1:m*(i+1)+1))
yy = copy(ww)

# Prior
Alpha = 1/100
Beta= 1/100

piσ²(s2) = pdf(InverseGamma(Alpha,Beta), s2[1])*pdf(InverseGamma(Alpha,Beta), s2[2])
lq(x, si) = Bridge.logpdfnormal(x, si^2*I)
PiError = InverseGamma(Alpha,Beta)

xi = 1./[20., 20., 20]

open(joinpath("output",simname,"truth.txt"), "w") do f
    println(f, "beta gamma eps s sigma err")
    println(f, join(round.([θtrue ; σtrue; sitrue],3)," "))
end

open(joinpath("output",simname,"params.txt"), "w") do f
    println(f, "n beta gamma eps s sigma err")
end

open(joinpath("output",simname,"info.txt"), "w") do f
    println(f, "n $n m $m T $TT A $Alpha B $Beta")
    println(f, "Y0 = $uu")
    println(f, "xi = $xi")
end

bacc = 0

mc = mcstart(vcat([BB[i][1:end-1] for i in 1:n]...).yy )
mcparams = mcstart([θ ; σ; si])

Bridge.constdiff(::FitzHughNagumo) = true

while true
    P = FitzHughNagumo(param(θ, σ)...)
    iter += 1

    for j = 1:2

        for i in j:2:n-1
            v = SVector{2,Float64}(Matrix(L)\Vector(Yobs.yy[i+1]))
            cs = Bridge.CSpline(Yobs.tt[i], Yobs.tt[i+2], Bridge.b( Yfil[i]..., P),  Bridge.b( Yfil[i+2]..., P))

            P° = PBridgeProp(P, Yfil[i]..., Yobs.tt[i+1], v, Yfil[i+2]..., L,  si^2*I, P.a, cs)
            B2 = bridge!(SamplePath(tt2[i], yy), sample!(SamplePath(tt2[i],ww), Wiener{SVector{2,Float64}}()),P°)
            # B2 = eulerb!(SamplePath(tt2[i], yy), sample!(SamplePath(tt2[i],ww), Wiener{SVector{2,Float64}}()),P°)

            llold = llikelihood([BB[i]; BB[i+1][2:end]], P°) + lq(Yobs.yy[i+1] -L*Yfil.yy[i+1], si)
            llnew = llikelihood(B2, P°) + lq(Yobs.yy[i+1]-L*B2.yy[m+1], si)
         #   println("$i->$(i+2) ($llnew - $llold) ")
            if rand() < exp(llnew - llold)
                bacc += 1
                Yfil.yy[i+1] = B2.yy[m+1]
                BB[i] = B2[1:m+1]
                BB[i+1] = B2[m+1:end]
            end
        end

    end
    let i = n #do the last bridge
            v = SVector{2,Float64}(Matrix(L)\Vector(Yobs.yy[i+1]))
            P° = PBridgeProp(P, Yfil[i]..., Yobs.tt[i+1], v, 1000., (@SVector [0.,0.]), L,  si^2*I, P.a)
            B2 = euler(P°.v0, sample(tt[m*(i-1)+1:m*i+1], Wiener{SVector{2,Float64}}()),P°)

            llold = llikelihood(BB[i], P°) + lq(Yobs.yy[i+1] -L*Yfil.yy[i+1], si)
            llnew = llikelihood(B2, P°) + lq(Yobs.yy[i+1]-L*B2.yy[m+1], si)
            #println("$i->$(i+1) ($llnew - $llold) ")

            if rand() < exp(llnew - llold)
                Yfil.yy[i+1] = B2.yy[m+1]
                BB[i] = B2
            end
    end

# update error variance
    if size(L,1) == 1 # one dimensional
        residual = Yobs.yy - map(x->L*x,Yfil.yy)
        si = sqrt(rand(ConjugatePriors.posterior_canon(PiError,suffstats(Distributions.NormalKnownMu(0.), Bridge.mat(residual)))))
        si = sqrt(rand(PiError))
    else
        error("todo: multivariate observations")
    end

    # update theta

    if iter % 1 == 1
        BBall = vcat([BB[i][1:end-1] for i in 1:n]...)
        θ° = θ + (2rand(length(θ)) .- 1).*scaleθ
        Pθ = FitzHughNagumo(param(θ, σ)...)
        Pθ° = FitzHughNagumo(param(θ°, σ)...)
        ll = girsanov(BBall, Pθ°, Pθ)
    #    println(ll)

        if rand()<exp(ll)  #iter > 200
            θ = θ°
        end
    end

    # update conjugate theta

    θ[2:end] = conjugateb(BB, θ[2:end],xi, P)

    # update sigma (and theta)
    if iter % 1 == 0
        σ° = σ .* exp.(scaleσ .* randn(length(σ)))
        θ° = θ #+ (2rand(length(θ)) .- 1).*scaleθ/3
        Pσ = FitzHughNagumo(param(θ, σ)...)
        Pσ° = FitzHughNagumo(param(θ°, σ°)...)
        ll = 0.
        for i in 1:n # reparametrization
            cs = Bridge.CSpline(Yobs.tt[i], Yobs.tt[i+1], Bridge.b( Yfil[i]..., P),  Bridge.b( Yfil[i+1]..., P))

            P° = BridgeProp(Pσ, BB[i].tt, (Yfil.yy[i+1], Yfil.yy[i+1]), Pσ.a, cs)
            P°° = BridgeProp(Pσ°, BB[i].tt, (Yfil.yy[i+1], Yfil.yy[i+1]), Pσ°.a, cs)
            Z = innovations(BB[i], P°)

            BBnew[i] = bridge(BridgePre(), Z, P°°)
            # BBnew[i] = eulerb(Z, P°°)
            ll += lptilde(P°°) - lptilde(P°) + llikelihood(BBnew[i], P°°) - llikelihood(BB[i], P°)
    #        println(lptilde(P°°), " - ", ptilde(lP°))
    #        println(llikelihood(BBnew[i], P°°), " - ",  llikelihood(BB[i], P°))
        end
        #print(ll)
        # f(σ²) = log σ²
        # df(σ²) = 1/σ²
        # 2 log σ° = 2 log σ + Z, Z ~ N(0,v)
        # q(σ°²|σ²) = 1/(σ°² sqrt(2piv)) exp(-1/2v |log(σ°²) - log(σ²)|^2)
        # q(σ²|σ°²) / q(σ°²|σ²) = σ°²/σ²

        if rand() < exp(ll) * (piσ²(σ°.^2)/piσ²(σ.^2)) * (prod(σ°.^2)/prod(σ.^2))
            σ = σ°
            θ = θ°
            for i in 1:n
                BB[i] = BBnew[i]
            end
        end
    end
    open(joinpath("output",simname,"params.txt"), "a") do f; println(f, iter, " ", join(round.([θ ; σ; si],8)," ")) end
    println(iter, "\t", join(round.([θ./θtrue; σ./σtrue; si/sitrue; 100bacc/iter/n  ],3),"\t"))

    BBall = vcat([BB[i][1:end-1] for i in 1:n]...)

    mc = mcnext(mc, BBall.yy)
    mcparams = mcnext(mcparams, [θ ; σ; si])

    if iter % 10 == 0
        xr = (5,7)
        @show iter
        plot2(Y, "r","r" ; xrange=xr, yrange=(-3,3),linewidth=0.5)
        oplot2(Yfil, "o","o"; xrange=xr, yrange=(-3,3),linewidth=0.1)
        oplot(Yobs, "ro"; xrange=xr, yrange=(-3,3),linewidth=0.1)
        display(oplot2(BBall, "b", "b"; xrange=xr, yrange=(-3,3), linewidth=0.7))
    end
    if iter >= 100 break end
end

plot2(Y, "r","b" ; yrange=(-3,3),linewidth=0.5)
oplot(Yobs,"+r")
savefig(joinpath("output",simname,"truth.svg"))
savefig(joinpath("output",simname,"truth.pdf"))

plot2(Y, "r","b" ; yrange=(-3,3),linewidth=0.5)
oplot(Yobs,"or";symbolsize=0.3)
mcb = mcband(mc);
oplot2(SamplePath(tt, mcb[1]),"r","b";linewidth=0.5)
oplot2(SamplePath(tt, mcb[2]),"r","b";linewidth=0.5)
hcat(mcbandste(mcparams)..., [θtrue; σtrue; si])

savefig(joinpath("output",simname,"band.svg"))
savefig(joinpath("output",simname,"band.pdf"))

plot2(vcat([BB[i] for i in 1:n]...), "r", "b"; yrange=(-3,3), linewidth=0.5)
oplot(Yobs,"+r")
savefig(joinpath("output",simname,"sample.svg"))
savefig(joinpath("output",simname,"sample.pdf"))


text:
s = BenchmarkGroup()
suite["connectivity"] = s

Random.seed!(1)
for G in GLIST
    for n in [100], d in [2,3]
        g = G(n , 5n, seed=1)
        s["neighborhood d=$d","$g"] = @benchmarkable neighborhood($g, 1, $d)
    end
end


text:
# Ex 2.1
open("ex1.acc") do file
    for ln in eachline(file)
        println(ln)
    end
end

# Ex 2.2
println("Enter a file name:")
fn = chomp(readline())

open(fn) do f
    for ln in eachline(f)
        println(ln)
    end
end

# Ex 2.3
println("Enter a file name:")
fn = chomp(readline())

nl = open(fn) do f
    i = 0
    for ln in eachline(f)
        i += 1
    end
    i
end

println(nl)

#Ex 2.4-9
for i = 1:3
    avg, sum, lc, neg, pos, zer, max_num, min_num = open("ex1.dat.$i") do f
        lc, sum, neg, pos, zer = 0, 0, 0, 0, 0
        max_num, min_num = -Inf, Inf
        for ln in eachline(f)
            ln_num = parse(Float64, chomp(ln))
            (ln_num > 0) ? pos += 1 : (ln_num < 0) ? neg += 1 : zer += 1
            (ln_num > max_num) && (max_num = ln_num)
            (ln_num < min_num) && (min_num = ln_num)
            lc += 1
            sum += ln_num
        end
        sum / lc, sum, lc, neg, pos, zer, max_num, min_num
    end
    println("Column ", i, " - Average: ", avg, ", Sum: ", sum,
    ", Line count: ", lc)
    println("Negative: ", neg, ", Positive: ", pos, ", Zeros: ", zer)
    println("Maximum number: ", max_num, ", Minimum number: ", min_num)
    println()
end

#=
avg = open("1to9") do f
    i, sum = 0, 0
    for ln in eachline(f)
        i += 1
        sum += parse(Float64, chomp(ln))
    end
    sum / i
end
println(avg)
=#

function get_number(num)
    println(string("Is your number ", string(num), "?"))
    chomp(readline())
end

function guess_number(mini, maxi)
    guess = Int(floor( (mini + maxi) / 2 ))
    response = get_number(guess)
    if response == "higher"
        guess_number(guess + 1, maxi)
    elseif response == "lower"
        guess_number(mini, guess)
    elseif response == "yes"
        println("Gotcha!")
        return
    end
end

guess_number(1,10)


text:
using GitHub, GitHub.name
using Base.Test

# The below tests are network-dependent, and actually make calls to GitHub's
# API. They're all read-only, meaning none of them require authentication.

testuser = Owner("julia-github-test-bot")
julweb = Owner("JuliaWeb", true)
ghjl = Repo("JuliaWeb/GitHub.jl")
testcommit = Commit("3a90e7d64d6184b877f800570155c502b1119c15")

hasghobj(obj, items) = any(x -> name(x) == name(obj), items)

# This token has public, read-only access, and is required so that our
# tests don't get rate-limited. The only way a malicious party could do harm
# with this token is if they used it to abuse the rate limit associated with
# the token (not too big of a deal). The token is hard-coded in an obsfucated
# manner in an attempt to thwart token-stealing crawlers.
auth = authenticate(string(circshift(["bcc", "3fc", "03a", "33e",
                                      "c09", "363", "5f1", "bd3",
                                      "fc6", "77b", '5', "9cf",
                                      "868", "033"], 3)...))

@test rate_limit(; auth = auth)["rate"]["limit"] == 5000

##########
# Owners #
##########

# test GitHub.owner
@test name(owner(testuser; auth = auth)) == name(testuser)
@test name(owner(julweb; auth = auth)) == name(julweb)

# test GitHub.orgs
@test hasghobj("JuliaWeb", first(orgs("jrevels"; auth = auth)))

# test GitHub.followers, GitHub.following
@test hasghobj("jrevels", first(followers(testuser; auth = auth)))
@test hasghobj("jrevels", first(following(testuser; auth = auth)))

# test GitHub.repos
@test hasghobj(ghjl, first(repos(julweb; auth = auth)))

################
# Repositories #
################

# test GitHub.repo
@test name(repo(ghjl; auth = auth)) == name(ghjl)

# test GitHub.forks
@test length(first(forks(ghjl; auth = auth))) > 0

# test GitHub.contributors
@test hasghobj("jrevels", map(x->x["contributor"], first(contributors(ghjl; auth = auth))))

# test GitHub.stats
@test stats(ghjl, "contributors"; auth = auth).status < 300

# test GitHub.branch, GitHub.branches
@test name(branch(ghjl, "master"; auth = auth)) == "master"
@test hasghobj("master", first(branches(ghjl; auth = auth)))

# test GitHub.commit, GitHub.commits
@test name(commit(ghjl, testcommit; auth = auth)) == name(testcommit)
@test hasghobj(testcommit, first(commits(ghjl; auth = auth)))

# test GitHub.file, GitHub.directory, GitHub.readme, GitHub.permalink
readme_file = file(ghjl, "README.md"; auth = auth)
src_dir = first(directory(ghjl, "src"; auth = auth))
owners_dir = src_dir[findfirst(c -> get(c.path) == "src/owners", src_dir)]
test_sha = "eab14e1ab7b4de848ef6390101b6d40b489d5d08"
readme_permalink = string(permalink(readme_file, test_sha))
owners_permalink = string(permalink(owners_dir, test_sha))
@test readme_permalink == "https://github.com/JuliaWeb/GitHub.jl/blob/$(test_sha)/README.md"
@test owners_permalink == "https://github.com/JuliaWeb/GitHub.jl/tree/$(test_sha)/src/owners"
@test readme_file == readme(ghjl; auth = auth)
@test hasghobj("src/GitHub.jl", src_dir)

# test GitHub.status, GitHub.statuses
@test get(status(ghjl, testcommit; auth = auth).sha) == name(testcommit)
@test !(isempty(first(statuses(ghjl, testcommit; auth = auth))))

# test GitHub.comment, GitHub.comments
@test name(comment(ghjl, 154431956; auth = auth)) == 154431956
@test !(isempty(first(comments(ghjl, 40; auth = auth))))

# These require `auth` to have push-access (it's currently a read-only token)
# @test hasghobj("jrevels", first(collaborators(ghjl; auth = auth)))
# @test iscollaborator(ghjl, "jrevels"; auth = auth)

##########
# Issues #
##########

state_param = Dict("state" => "all")

# test GitHub.pull_request, GitHub.pull_requests
@test get(pull_request(ghjl, 37; auth = auth).title) == "Fix dep warnings"
@test hasghobj(37, first(pull_requests(ghjl; auth = auth, params = state_param)))

# test GitHub.issue, GitHub.issues
@test get(issue(ghjl, 40; auth = auth).title) == "Needs test"
@test hasghobj(40, first(issues(ghjl; auth = auth, params = state_param)))

############
# Activity #
############

# test GitHub.stargazers, GitHub.starred
@test length(first(stargazers(ghjl; auth = auth))) > 10 # every package should fail tests if it's not popular enough :p
@test hasghobj(ghjl, first(starred(testuser; auth = auth)))

# test GitHub.watched, GitHub.watched
@test hasghobj(testuser, first(watchers(ghjl; auth = auth)))
@test hasghobj(ghjl, first(watched(testuser; auth = auth)))


text:
# Learn the environment!
function learn!(learner, env::Union{EnvMultinomial, EnvGaussian}, callback)
    for i in 1:env.nsteps
        interact!(env)
        update!(learner, env.state)
        callback!(callback, learner, env)
    end
end


text:
"""
An appender object that can be used to append to a table
"""
mutable struct Appender
    handle::duckdb_appender

    function Appender(con::Connection, table::AbstractString)
        handle = Ref{duckdb_appender}()
        if duckdb_appender_create(con.handle, C_NULL, table, handle) != DuckDBSuccess
            error_ptr = duckdb_appender_error(handle)
            if error_ptr == C_NULL
                error_message = string("Opening of Appender for table \"", table, "\" failed: unknown error")
            else
                error_message = unsafe_string(error_ptr)
            end
            duckdb_appender_destroy(handle)
            throw(QueryException(error_message))
        end
        con = new(handle[])
        finalizer(_close_appender, con)
        return con
    end
    function Appender(db::DB, table::AbstractString)
        return Appender(db.main_connection, table)
    end
end

function _close_appender(appender::Appender)
    if appender.handle != C_NULL
        duckdb_appender_destroy(appender.handle)
    end
    appender.handle = C_NULL
    return
end

function close(appender::Appender)
    _close_appender(appender)
    return
end

append(appender::Appender, val::AbstractFloat) = duckdb_append_double(appender.handle, Float64(val));
append(appender::Appender, val::Bool) = duckdb_append_boolean(appender.handle, val);
append(appender::Appender, val::Int8) = duckdb_append_int8(appender.handle, val);
append(appender::Appender, val::Int16) = duckdb_append_int16(appender.handle, val);
append(appender::Appender, val::Int32) = duckdb_append_int32(appender.handle, val);
append(appender::Appender, val::Int64) = duckdb_append_int64(appender.handle, val);
append(appender::Appender, val::UInt8) = duckdb_append_uint8(appender.handle, val);
append(appender::Appender, val::UInt16) = duckdb_append_uint16(appender.handle, val);
append(appender::Appender, val::UInt32) = duckdb_append_uint32(appender.handle, val);
append(appender::Appender, val::UInt64) = duckdb_append_uint64(appender.handle, val);
append(appender::Appender, val::Float32) = duckdb_append_float(appender.handle, val);
append(appender::Appender, val::Float64) = duckdb_append_double(appender.handle, val);
append(appender::Appender, val::Missing) = duckdb_append_null(appender.handle);
append(appender::Appender, val::Nothing) = duckdb_append_null(appender.handle);
append(appender::Appender, val::AbstractString) = duckdb_append_varchar(appender.handle, val);
append(appender::Appender, val::Vector{UInt8}) = duckdb_append_blob(appender.handle, val, sizeof(val));
append(appender::Appender, val::WeakRefString{UInt8}) = duckdb_append_varchar(stmt.handle, i, val.ptr, val.len);

function append(appender::Appender, val::Any)
    println(val)
    throw(NotImplementedException("unsupported type for append"))
end

function end_row(appender::Appender)
    duckdb_appender_end_row(appender.handle)
    return
end

function flush(appender::Appender)
    duckdb_appender_flush(appender.handle)
    return
end

DBInterface.close!(appender::Appender) = _close_appender(db)

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
"""
Experimental module providing several types to store annotated text [`Token`](ref)s,
[`Line`](ref)s.
"""
module Tokens
############################################################
## Tokens
## TODO: move intenring into parsing (creating from a db interning in tokens wastes mem)
export AbstractToken, variable, value
export Token, TokenValue, TokenTuple, TokenString
import Base: ==, hash
using CombinedParsers
using CombinedParsers.Regexp
import CombinedParsers.Regexp: whitespace_maybe, whitespace_horizontal, word
using ..CombinedParserTools
import ..CombinedParserTools: footnote, quotes, delimiter, extension

"""
An abstract type for annotated textual tokens.
"""
abstract type AbstractToken end
## TODO: rename variable to annotation
variable(x::AbstractToken) = error("implement variable(x::$(typeof(x)))")
value(x::AbstractToken) = error("implement value(x::$(typeof(x)))")
label(x::AbstractToken) = error("implement label(x::$(typeof(x)))")

variable_colors=Dict(
    :ext => 36,
    :macro => 36,
    :number => 36,
    :root => :dark_gray,
    :folder => :blue,
    :file => :yellow,
    :syllable => :yellow,
    :operator => :yellow,
    :name => :yellow,
    :footnote => :yellow,
    :type => :red,
    :field => :light_red,
    :ellipsis => :grey,
    :whitespace => :grey,
    :list => :light_black,
    :protocol => :light_blue,
    :domain => :light_blue,
    :Symbol => :yellow,
    :String => :yellow,
    :paren => :light_black,
    :quote => :light_black,
    Symbol("wikt:de") => :light_blue,
    :htmlcomment => :light_black,
    :unknown => :light_red,
    :meaning => :light_black
)

# export value_empty
# #value_empty(x::Pair) = value_empty(x.second) ## needed in tryparsenext
# value_empty(x::Vector) = isempty(x)
# value_empty(x) = false
# value_empty(::Union{Nothing,Missing}) = true
# value_empty(x::String) = x==""
# value_empty(x::AbstractToken) = value(x) === missing || value(x)==""

export isinformative, isvariable
isinformative(i) = true
isvariable(i) = false

function Base.show(io::IO, z::AbstractToken)
    if get(stdout,:color,false)
    color=get(variable_colors,
              Symbol(variable(z)), 36)
    if variable(z) in [ :literal, :capitalized, :delimiter ]
        value(z)!==missing && print(io,value(z))
    elseif !isinformative(z)
        printstyled(io, value(z); bold=true,
                    color=:darkgray)
    elseif value(z)===missing ## || value(z)==""
        printstyled(io, variable(z); bold=true,
                    color=color)
    else
        col=get(variable_colors, variable(z), missing)
        if variable(z) == :macro
            printstyled(
                io, "{{{", value(z),"}}}";
                bold=true, color=color
            )
        elseif col === missing
            if get(io,:compact,false)
                print(io,value(z))
            else
                printstyled(
                    io, "[[",variable(z), "][", value(z),"]]";
                    bold=true, color=color
                )
            end
        else
            printstyled(io, value(z); color=color)
        end
    end
    else
        print(io,value(z))
    end
end

export Token
"""
Represents an annotation of field `value::String` with field `name::Symbol`.
"""
struct Token <: AbstractToken
    name::Symbol ## todo: CategoricalArrays.CategoricalValue
    value::String
    function Token(name::Symbol, value::T) where {T<:AbstractString}
        new(name, value)
    end
end
Token(name::Symbol, value::Union{Missing, Nothing}) = Token(name, "")
Token(name::Symbol) = Token(name, "")
Token(x::Pair) = Token(x.first, x.second)
Token(x::Pair{<:AbstractString,Symbol}) = Token(x.second, x.first)
Token(x::Token) = x
"""
    Token(p::NamedParser)

Parser matching p.parser, transforming the result in a
Token
"""
function Token(p::NamedParser)
    @assert result_type(p)<:AbstractString
    map(v->Token(p.name,v),p)
end
function Token(name::AbstractString, value)
    Token(Symbol(name), value)
end

isinformative(i::Token)  =
    !(variable(i) in [ :delimiter, :indent, :list, :enum, :whitespace ])
isvariable(i::Token)  =
    !(variable(i) in [ :literal, :capitalized ]) && isinformative(i)

import Base: convert
Base.convert(::Type{Token},e::Pair) =
    Token(Symbol(e.first), e.second)

export @l_str, @ws_str, @delim_str, @T_str
"""
    T_str(value, name)

String macro to create a Token.
```jldoctest
julia> T"Jack"name
Jack

julia> T"123"number
123
```
"""
macro T_str(value, name)
    Token(name, string(value))
end
macro l_str(x)
    Token(:literal, x)
end
macro ws_str(x)
    Token(:whitespace, x)
end
macro delim_str(x)
    Token(:delimiter, x)
end
# ws(x) = Token(:whitespace, x)

export TokenPair
"""

"""
struct TokenPair{K,V} <: AbstractToken
    key::K
    value::V
end
==(x::TokenPair,y::TokenPair) =
    x.key==y.key && x.value==y.value
hash(x::TokenPair, h::UInt) = hash(x.key, hash(x.value,h))
parentheses = Dict{Any,Any}(:paren=>("(", ")"),
                            :bracket=>("[", "]"),
                            :curly=>("{", "}"),
                            :angle=>("<", ">"),
                            :quote=>("\"","\""),
                            :squote=> ("'","'"),
                            :german_quote => ("„","“"),
                            :htmlcomment=> ("<!--","-->"),
                            ## :pre=> ("<pre>","</pre>"),
                            ## :nowiki=> ("<nowiki>","</nowiki>"),
                            )
import Base: with_output_color
function Base.show(io::IO, z::TokenPair)
    inner_print(io::IO,x::AbstractVector) =
        for t in x; print(io, t); end
    inner_print(io::IO,x) =
        print(io, x)

    if z.key==:hyphenation
        join(io, z.value, "·")
    elseif z.key==:link
        join(io, z.value)
    elseif z.key==:italics
        with_output_color(inner_print, :underline, io, z.value)
    elseif z.key==:bold
        with_output_color(inner_print, :bold, io, z.value)
    elseif z.key==:bolditalics
        with_output_color(inner_print, :bold, io, z.value)
    else
        open, close = get(parentheses, z.key, (z.key,z.key))
        print(io, open)
        inner_print(io,z.value)
        print(io, close)
    end
end

export open_close

"""
    open_close(name::Symbol, inner, open, close=open, stops::Tuple=tuple())

Parser returning a [`TokenPair`](@ref)`(name,inner_parsing)`
"""
function open_close(name::Symbol,
                    inner,
                    open, close=open,
                    stops::Tuple=tuple())
    with_name(
        name,
        Sequence(open,
                 Repeat_until(
                     inner,
                     sEither(close, stops...))) do v
        TokenPair(name,v[2])
        end)
end

export @annotate
"""
```jldoctest
julia> digit = re"[0-9]";

julia> @annotate [digit, :lletter => re"[a-z]"]
|🗄... Either |> map(NamedString)
├─ [0-9] CharIn |> map(#161) |> with_name(:digit)
└─ [a-z] CharIn |> map(#161) |> with_name(:lletter)
::NamedString

```
"""
macro annotate(x)
    if x.head == :vect
        top = Expr(:vect)
        for e_ in x.args
            e = e_
            while e isa Expr && e.head==:call && e.args[1]==Symbol("!")
                e = e.args[2]
            end
            if e isa Symbol
                push!(top.args,:($(QuoteNode(e)) => $e_))
            elseif e.head==:call && e.args[1] == Symbol("=>")
                push!(top.args,e)
            else
                dump(x)
                error("@annotate supports only variables and `:annotation => parser`, not `$e`.")
            end
        end
        esc(:(Either( tuple(( Token(parser(p)) for p in $top)...))))
    else
        dump(x)
        :()
    end
end

export NamedString
"parametrized Token struct -- dangerously slow!"
struct NamedString{name} <: AbstractToken
    value::String
    function NamedString(name, value)
        new{Symbol(name)}(value)
    end
    function NamedString(x::Token)
        new{variable(x)}(value(x))
    end
end
function NamedString(x::NamedString)
    x
end
Base.convert(::Type{NamedString}, x::Token) =
    NamedString(variable(x),value(x))

Base.isless(x::NamedString, y::NamedString) =
    isless(value(x), value(y))

Base.show(io::IO,x::NamedString{:type}) =
    print(io,Token(:type,x.value),".")

Base.show(io::IO,x::NamedString{:field}) =
    print(io,Token(:field,x.value),"=")

Base.show(io::IO,x::NamedString{:whitespace}) =
    printstyled(io,x.value; color=:underline)

Base.propertynames(x::NamedString) = (:name,:value)
Base.getproperty(x::NamedString, p::Symbol) =
    if p == :name
        variable(x)
    elseif p == :value
        value(x)
    else
        error("no field $p in NamedString")
    end

variable(x::NamedString{name}) where name = name
value(x::NamedString) = getfield(x,1)
==(x::NamedString,y::NamedString) =
    variable(x)==variable(y) && value(x)==value(y)
hash(x::NamedString, h::UInt) = hash(variable(x), hash(value(x),h))

export Node
struct Node{A,T} <: AbstractToken
    name::Symbol
    attributes::Vector{A}
    children::Vector{T}
    function Node(name::Symbol, attrs, value)
        new{eltype(attrs),eltype(value)}(name, attrs,value)
    end
    function Node(name::AbstractString, attrs::Vector{A}, value::Vector{T}) where {A,T}
        new{A,T}(Symbol(name), attrs,value)
    end
    function Node{T}(name::AbstractString, attrs, value) where T
        new{Token,T}(Symbol(name), attrs, _convert(Vector{T},value))
    end
    function Node{A,T}(name, attrs, value) where {A,T}
        new{A,T}(Symbol(name), convert(Vector{A},attrs), convert(Vector{T},value))
    end
end

==(a::Node, b::Node) = a.name==b.name && a.attributes==b.attributes && a.children==b.children
hash(x::Node, h::UInt) = hash(x.name, hash(x.attributes, hash(x.children,h)))
function Base.show(io::IO, x::Node) where {T}
    print(io,"<$(x.name)")
    for a in x.attributes
        if a isa Token
            print(io," ", variable(a), "=\"", value(a), "\"")
        else ## this is for templates and such in wikitext
            print(io,a)
        end
    end
    if !isempty(x.children)
        print(io,">")
        for c in x.children
            print(io,c)
        end
        print(io,"</$(x.name)>")
    else
        print(io,"/>")
    end
end

import InternedStrings: intern
# import ..ParserAlchemy: ParserTypes, instance, map_at, Repeat, Sequence, Either, Optional, alternate, FlatMap, Repeat_until
# import ..ParserAlchemy: result_type, regex_string

export ReferringToken
struct ReferringToken{Tt, Tv, I} <: AbstractToken
    name::Tt
    value::Tv
    reference::I
end
value(x::Union{Token, TokenPair, ReferringToken}) = x.value
variable(x::Union{Token, ReferringToken}) = x.name
variable(x::TokenPair) = x.key

export TokenString
const TokenString = Vector{<:AbstractToken}

function Base.string(x::TokenString)
    b = IOBuffer()
    bc = IOContext(b, :compact => true)
    for e in x
        show(bc, e)
    end
    String(take!(b))
end
# @deprecate TokenString(x...) tokenize(x...)

# const TokenTuple = Tuple{Vararg{Token, N} where N}
# TokenString{Tt, T} = Tuple{Vararg{Token{t,T} where {t <:Tt}, N} where N}
# TokenNest{Tv} = Tuple{Vararg{Union{TokenTuple{t},Token{t,s}} where {t, s <:Tv}, N} where N}
## TokenTuple(x::AbstractToken{Any, Tv}...) where {Tv} = x
# TokenString(x::TokenTuple) = x

export tokens
tokens(x::Vector{<:Union{AbstractToken, AbstractString, Symbol}}) =
    Iterators.repeated(1  => x, 1)
tokens(x::Union{Number,Symbol, AbstractString, AbstractToken}) =
    Iterators.repeated(1 => Iterators.repeated(x, 1), 1)
## tokens(d::Dict{Symbol,<:AbstractString}) =
##     [ 1 => [ Token(x.first, x.second) for x in d ] ]
tokens(d::Dict) =
    [] #Iterators.flatten( tokens(x.second) for x in d )
tokens(x::Vector) =
    Iterators.flatten(tokens(y) for y in x)

"""
    tokens(n::NamedTuple{names,t})

Recursively collect all tokens in a NamedTuple.
"""
function tokens(n::NamedTuple{names,t}) where {names,t}
    val(field) = tokens(getproperty(n,field))
    R = Iterators.flatten( val(field)
                           for field in names
                           # if !isempty(val(field))
                           )
    R
    # Pair[ w => collect(v) for (w,v) in R ]
end

"""
    tokens(n::T)

Recursively collect all tokens in a struct.
"""
function tokens(n::T) where {T}
    val(field) = tokens(getfield(n,field))
    R = Iterators.flatten( val(field)
                           for field in fieldnames(T)
                           # if !isempty(val(field))
                           )
    R
    # Pair[ w => collect(v) for (w,v) in R ]
end

import Base: convert
Base.convert(::Type{TokenString}, x::String) = tokenize(x)

import CombinedParsers: _iterate, MatchState, state_type

export IteratorParser
"""
    IteratorParser{T}(label::String,match::Function,f::Function)

An iterator parser `p` on an AbstractArray Sequence `str`.
filters for `p.match(str[i])`, and return transformed `p.f(str[i],i).

TODO: deprecate for CharIn
"""
struct IteratorParser{T} <: CombinedParser{MatchState,T}
    match::Function
    label::String
    IteratorParser{T}(f::Function, l::String) where {T} =
        new{T}(x->x isa T && f(x),l)
    IteratorParser(f::Function, T::Type, l::String) =
        new{T}(x->x isa T && f(x),l)
    IteratorParser(T::Type, l::String) =
        new{T}(x->x isa T,l)
    IteratorParser(T::Type) =
        new{T}(x->x isa T,"$T")
end
Base.show(io::IO, x::IteratorParser) = print(io, x.label)
## CombinedParsers.state_type(::Type{<:IteratorParser}) = MatchState
# IteratorParser{T}(label::String,match::Function,f::Function) =
#     filter

CombinedParsers._iterate(
    tok::IteratorParser,
    str, till,
    i,after,state::MatchState) = nothing
function CombinedParsers._iterate(tok::IteratorParser,
                  str, till,
                  i,after,state::Nothing)
    ##@show typeof(str[i])
    if i<=lastindex(str) && tok.match(str[i])
        Base.nextind(str,i), MatchState()
    else
        nothing
    end
end

function Base.get(tok::IteratorParser, str, till, after,i,state::MatchState)
    str[i]
end
@inline Base.nextind(str,i::Int,parser::IteratorParser,x) =
    Base.nextind(str,i)
@inline Base.prevind(str,i::Int,parser::IteratorParser,x) =
    Base.prevind(str,i)

include("lines.jl")

include("html.jl")

import CombinedParserTools: word, footnote, quotes, capitalized, delimiter
import CombinedParsers.Regexp: word_char

export default_tokens
default_tokens(mask=AnyChar()) =
    @annotate [ :number     => !!Repeat1(re"[0-9]" .& mask),
                :literal    => !!Repeat1(word_char .& mask),
                :delimiter  => !!(delimiter .& mask),
                :footnote   => !!footnote .& mask,
                :quote      => !!(quotes .& mask)
                ]

export tokenstring,default_token,quotes
tokenstring(unknown=AnyChar()) = Repeat(default_token(unknown) | @annotate [ :unknown => !!unknown])

end


text:
module test_microcoverage_ex
function internalbranching(i)
    x = (i < 1)? 5 : ((i > 2)? 9 : 3)
    x
end

statementfunction(x) = x+1

function andor(k)
    y = 0
    if k < 10 && k*k < 36
        y = 8
    end
    if k > 20 || k * k > 225
        y = 9
    end
end

function hasunreachablestmt(k)
    return k * 19
    x *= 10
end

function runfuncs()
    for j = 0 : 3
        internalbranching(j)
        statementfunction(j)
    end

    for l = 5 : 22
        andor(l)
    end
    hasunreachablestmt(5)
end

end

max_stars_repo_name: StephenVavasis/microcoverage

max_issues_repo_name: StephenVavasis/microcoverage

max_forks_repo_name: StephenVavasis/microcoverage


text:
using Documenter, AWSClusterManagers

makedocs(;
    modules = [AWSClusterManagers],
    authors = "Invenia Technical Computing Corporation",
    repo = "https://github.com/JuliaCloud/AWSClusterManagers.jl/blob/{commit}{path}#L{line}",
    sitename = "AWSClusterManagers.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages = [
        "Home" => "index.md",
        "Docker" => "pages/docker.md",
        "Batch" => "pages/batch.md",
        "ECS" => "pages/ecs.md",
        "Design" => "pages/design.md",
        "API" => "pages/api.md",
    ],
    strict = true,
    checkdocs = :exports,
)

deploydocs(;
    repo = "github.com/JuliaCloud/AWSClusterManagers.jl",
    devbranch = "main",
    push_preview = true,
)


text:
include(joinpath(@__DIR__, "testcommon.jl"))

# Wrap the functions defined in C++
module CppTypes

using CxxWrap

@wrapmodule(Main.libtypes)

const greet = getindex ∘ greet_cref

export enum_to_int, get_enum_b, World
export AConstRef, ReturnConstRef, value, CallOperator, ConstPtrConstruct

end

# Stress test
for i in 1:1000000
  global d = CppTypes.DoubleData()
end

# Default constructor
w = CppTypes.World()
println("Dumping type w...")
dump(w)
@test CppTypes.greet(w) == "default hello"

@show fw = CppTypes.world_factory()
@test CppTypes.greet(fw[]) == "factory hello"

swf = CppTypes.shared_world_factory()
@test CppTypes.greet_shared(swf) == "shared factory hello"
@test CppTypes.greet_shared_const(swf) == "shared factory hello"
@test CppTypes.greet(swf[]) == "shared factory hello" # Explicit dereference
@test CppTypes.greet(swf) == "shared factory hello" # Automatic conversion
swf2 = CppTypes.smart_world_factory()
@test CppTypes.greet_smart(swf2) == "smart factory hello"
@test CppTypes.greet_smart(swf) == "shared factory hello" # auto-convert between pointers
@test CppTypes.greet(swf2[]) == "smart factory hello" # Explicit dereference
@test CppTypes.greet(swf2) == "smart factory hello" # Automatic conversion
@test CppTypes.greet_weak(swf) == "shared factory hello"
@test_throws MethodError CppTypes.greet_weak(swf2) == "shared factory hello"

swfr = CppTypes.shared_world_ref()
@test CppTypes.greet(swfr[]) == "shared factory hello ref"
CppTypes.reset_shared_world!(swfr, "reset shared pointer")
@test CppTypes.greet(swfr[]) == "reset shared pointer"
CppTypes.reset_shared_world!(swfr, "shared factory hello ref")

@test CppTypes.greet(CppTypes.boxed_world_factory()) == "boxed world"
@test CppTypes.greet(CppTypes.boxed_world_pointer_factory()[]) == "boxed world pointer"
@test CppTypes.greet(CppTypes.world_ref_factory()) == "reffed world"

@show uwf = CppTypes.unique_world_factory()
@test CppTypes.greet(uwf) == "unique factory hello"

byval = CppTypes.world_by_value()
@test CppTypes.greet(byval) == "world by value hello"

CppTypes.set(w, "hello")
@show CppTypes.greet(w)
@test CppTypes.greet(w) == "hello"
@test CppTypes.greet_lambda(w) == "hello"

w = CppTypes.World("constructed")
@test CppTypes.greet(w) == "constructed"

w_assigned = w
w_deep = deepcopy(w)

@test w_assigned == w
@test w_deep != w

# Destroy w: w and w_assigned should be dead, w_deep alive
finalize(w)
println("finalized w")
if !(Sys.iswindows() && Sys.WORD_SIZE == 32)
  @test_throws ErrorException CppTypes.greet(w)
  println("throw test 1 passed")
  @test_throws ErrorException CppTypes.greet(w_assigned)
  println("throw test 2 passed")
end
@test CppTypes.greet(w_deep) == "constructed"
println("completed deepcopy test")

wnum = CppTypes.World(1)
@test CppTypes.greet(wnum) == "NumberedWorld"
finalize(wnum)
@test CppTypes.greet(wnum) == "NumberedWorld"

noncopyable = CppTypes.NonCopyable()
other_noncopyable = deepcopy(noncopyable)
@test other_noncopyable.cpp_object == noncopyable.cpp_object

@test CppTypes.value(CppTypes.value(CppTypes.ReturnConstRef())) == 42

wptr = CppTypes.World()
@test CppTypes.greet(CppTypes.ConstPtrConstruct(CxxPtr(wptr))) == "default hello"

call_op = CppTypes.CallOperator()
@test call_op() == 43
@test call_op(42) == 42

mutable struct JuliaTestType
    a::Float64
    b::Float64
end

function julia_test_func(data)
  println("a: ", data.a, ", b: ", data.b)
  @test data.a == 2.0
  @test data.b == 3.0
end

CppTypes.call_testtype_function()

@test CppTypes.enum_to_int(CppTypes.EnumValA) == 0
@test CppTypes.enum_to_int(CppTypes.EnumValB) == 1
@test CppTypes.get_enum_b() == CppTypes.EnumValB
@test CppTypes.EnumValA + CppTypes.EnumValB == CppTypes.EnumValB
@test CppTypes.EnumValA | CppTypes.EnumValB == CppTypes.EnumValB

foovec = Any[CppTypes.Foo(StdWString("a"), [1.0, 2.0, 3.0]), CppTypes.Foo(StdWString("b"), [11.0, 12.0, 13.0])] # Must be Any because of the boxing
@test CppTypes.name(foovec[1]) == "a"
@test CppTypes.data(foovec[1]) == [1.0, 2.0, 3.0]
@test CppTypes.name(foovec[2]) == "b"
@test CppTypes.data(foovec[2]) == [11.0, 12.0, 13.0]
CppTypes.print_foo_array(foovec)

@test !isnull(CppTypes.return_ptr())
@test isnull(CppTypes.return_null())

warr1 = CppTypes.World[CppTypes.World("world$i") for i in 1:5]
warr2 = [CppTypes.World("worldalloc$i") for i in 1:5]
wvec1 = StdVector(warr1)
wvec2 = StdVector(warr2)

for (i,(w1,w2)) in enumerate(zip(wvec1,wvec2))
  @test CppTypes.greet(w1) == "world$i"
  @test CppTypes.greet(w2) == "worldalloc$i"
end

@test CppTypes.greet_vector(wvec1) == string(("world$i " for i in 1:5)...)[1:end-1]

a = [4.0]
return_int() = Int32(3)
return_ptr_double() = pointer(a)
return_world() = CppTypes.World("returned_world")
wptr = CppTypes.World("returned_world_ptr")
wref = CppTypes.World("returned_world_ref")
return_world_ptr() = CxxPtr(wptr)
return_world_ref() = CxxRef(wref)

@test CppTypes.test_unbox() == fill(true,7)


text:
# m4_3t.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
	x̄ = mean(df.weight)
	x = range(minimum(df.weight), maximum(df.weight), length = 100)
end;

@model function ppl4_3(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* (weights .- x̄)
    for i in eachindex(heights)
        heights[i] ~ Normal(μ[i], σ)
    end
end

m4_3t = ppl4_3(df.weight, df.height)
q4_3t = quap(m4_3t, NelderMead())

nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns4_3t = mapreduce(c -> sample(m4_3t, sampler, nsamples), chainscat, 1:nchains)
chns4_3t_df = DataFrame(chns4_3t)

quap4_3t_df = DataFrame(rand(q4_3t.distr, 10_000)', q4_3t.params)

part4_3t =	Particles(chns4_3t[[:a, :b, :σ]])

# End of m4.3t.jl


text:
## NACA Parametrisations
#==========================================================================================#

# NACA 4-digit parameter functions
naca4_thickness(t_by_c, xc, sharp_trailing_edge :: Bool) = 5 * t_by_c * (0.2969 * √xc - 0.1260 * xc - 0.3516 * xc^2 + 0.2843 * xc^3 - (ifelse(sharp_trailing_edge, 0.1036, 0.1015) * xc^4))

naca4_camberline(pos, cam, xc) = ifelse(
                                        xc < pos,
                                        (cam / pos^2) * xc * (2 * pos - xc),
                                        cam / (1 - pos)^2 * ( (1 - 2 * pos) + 2 * pos * xc - xc^2)
                                       )

naca4_gradient(pos, cam, xc) = atan(2 * cam / (ifelse(xc < pos, pos^2, (1 - pos)^2)) * (pos - xc))

"""
    naca4_coordinates(digits :: NTuple{4, <: Real}, n :: Integer, sharp_trailing_edge :: Bool)

Generate the coordinates of a NACA 4-digit series profile with a specified number of points, and a Boolean flag to specify a sharp or blunt trailing edge.
"""
function naca4_coordinates(digits :: NTuple{4, <: Real}, n :: Integer, sharp_trailing_edge :: Bool)
    # Camber
    cam = digits[1] / 100
    # Position
    pos = digits[2] / 10
    # Thickness-to-chord ratio
    t_by_c = (10 * digits[3] + digits[4]) / 100

    # Cosine spacing
    xs = cosine_spacing(0.5, 1.0, n)

    # Thickness distribution
    thickness = naca4_thickness.(Ref(t_by_c), xs, Ref(sharp_trailing_edge))
    if pos == 0 || cam == 0
        x_upper = xs
        y_upper = thickness
        x_lower = xs
        y_lower = -thickness
    else
        # Compute camberline
        camber  = naca4_camberline.(Ref(pos), Ref(cam), xs)
        # Compute gradients
        grads   = naca4_gradient.(Ref(pos), Ref(cam), xs)
        # Upper surface
        x_upper = @. xs - thickness * sin(grads)
        y_upper = @. camber + thickness * cos(grads)
        # Lower surface
        x_lower = @. xs + thickness * sin(grads)
        y_lower = @. camber - thickness * cos(grads)
    end
    coords = [ [x_upper y_upper][end:-1:2,:];
                x_lower y_lower             ]
end

"""
    naca4(digits :: NTuple{4, <: Real}, n :: Integer = 40; sharp_trailing_edge :: Bool)

Generate a `Foil` of a NACA 4-digit series profile with a specified number of points (40 by default), and a named option to specify a sharp or blunt trailing edge.
"""
naca4(digits :: NTuple{4, <: Real}, n = 40; sharp_trailing_edge = true) = Foil(naca4_coordinates(digits, n, sharp_trailing_edge), string("NACA ", digits...))

naca4(a, b, c, d, n = 40; sharp_trailing_edge = true) = naca4((a,b,c,d), n; sharp_trailing_edge = sharp_trailing_edge)


text:
mutable struct LoadingsGibbsOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    if_xml::IntegratedFactorsXMLElement
    tdl_xml::TraitLikelihoodXMLElement
    sparsity::String
    weight::Real

    function LoadingsGibbsOperatorXMLElement(
            if_xml::IntegratedFactorsXMLElement,
            tdl_xml::TraitLikelihoodXMLElement)

        return new(nothing, if_xml, tdl_xml, bn.NONE, 1.0)
    end
end

function make_xml(lgoxml::LoadingsGibbsOperatorXMLElement)
    el = new_element(bn.LOADINGS_GIBBS_OP)
    set_attribute(el, bn.WEIGHT, lgoxml.weight)
    set_attributes(el, [bn.RANDOM_SCAN => bn.FALSE,
                        bn.NEW_MODE => bn.TRUE,
                        bn.CONSTRAINT => bn.NONE,
                        bn.SPARSITY => lgoxml.sparsity])

    make_xml(lgoxml.if_xml)
    make_xml(lgoxml.tdl_xml)

    add_ref_el(el, lgoxml.if_xml.el)
    add_ref_el(el, lgoxml.tdl_xml.el)
    add_ref_el(el, get_normal_prior(lgoxml.if_xml.loadings_prior))

    lgoxml.el = el
end

function sparsity_constraint!(lgo::LoadingsGibbsOperatorXMLElement, constraint::String)
    lgo.sparsity = constraint
end

function get_parameter(op::LoadingsGibbsOperatorXMLElement)
    return get_loadings(op.if_xml)
end


text:
using Test, Lasso, StableRNGs

const testrng = StableRNG(13)

using CSV, DataFrames
readcsvmat(path;header=false, kwargs...) = Matrix{Float64}(DataFrame(CSV.File(path;header=header, kwargs...)))

@testset "Lasso" begin

include("lasso.jl")
include("gammalasso.jl")
include("fusedlasso.jl")
include("trendfiltering.jl")
include("cross_validation.jl")
include("segselect.jl")

end


text:
if !@isdefined(TimeTestExecution)
    const TimeTestExecution = false
end

module BlackBoxOptimTests

using LinearAlgebra, Random
using Printf: @printf, @sprintf
using SpatialIndexing

const SI = SpatialIndexing

TestDir = first(splitdir(@__FILE__()))

# If two arguments the second one if filename of a testset file
# listing the testfiles to use.
if length(ARGS) == 2 && isfile(ARGS[2])
    test_file_list = ARGS[2]
else
    test_file_list = joinpath(TestDir, "testset_normal.txt")
end

TestFiles = filter(fn -> isfile(joinpath(TestDir, fn)),
    readlines(test_file_list))

function latest_changed_file(files, dir = "")
    files[first(sortperm(map(fn -> mtime(joinpath(dir, fn)), files), rev = true))]
end

# If a first argument is given it must be either:
#  all =>           run all test files in the testset
#  latestchanged => run only the latest changed file in the testset
if length(ARGS) > 0
    if ARGS[1] == "all"
        TestFiles = TestFiles # Change nothing so run them all
    elseif ARGS[1] == "latestchanged"
        TestFiles = AbstractString[latest_changed_file(TestFiles, TestDir)]
        println("Testing files: $(TestFiles)")
    end
end

startclocktime = time()
include("helper.jl")

import Compat.String

if Main.TimeTestExecution

function get_git_remote_and_branch()
    lines = split(read(`git remote -v show`, String), "\n")
    remote = match(r"[a-z0-9]+\s+([^\s]+)", lines[1]).captures[1]
    branch = strip(read(`git rev-parse --abbrev-ref HEAD`, String))
    commit = strip(read(`git rev-parse HEAD`, String))
    return remote, branch, commit
end

gitremote, gitbranch, gitcommit = get_git_remote_and_branch()
gitstr = gitremote * "/" * gitbranch * "/" * gitcommit[1:6]
versionstr = string(VERSION)

using CSV, DataFrames

TestTimingFileName = "test/timing_testing.csv"

if isfile("test/timing_testing.csv")
    timing_data = CSV.read("test/timing_testing.csv")
else
    timing_data = DataFrame(TimeStamp = AbstractString[],
                            Julia = AbstractString[],
                            Git = AbstractString[],
                            TestFile = AbstractString[],
                            Elapsed = Float64[])
end

end

using CPUTime

starttime = CPUtime_us()
@testset "BlackBoxOptim test suite" begin

for t in TestFiles
    Main.TimeTestExecution && CPUtic()

    # Including the test file runs the tests in there...
    include(t)

    if Main.TimeTestExecution
        elapsed = CPUtoq()
        datestr = Libc.strftime("%Y%m%d %H:%M.%S", time())
        push!(timing_data, [datestr, versionstr, gitstr, t, elapsed])
    end
    print("."); flush(stdout);
end
println() # So Base.Test summary is correctly aligned...
end
elapsed = float(CPUtime_us() - starttime)/1e6

if Main.TimeTestExecution
    datestr = Libc.strftime("%Y%m%d %H:%M.%S", time())
    using SHA
    hash = bytes2hex(sha512(join(map(fn -> read(open(joinpath("test", fn)), String), TestFiles))))[1:16]
    push!(timing_data, [datestr, versionstr, gitstr, "TOTAL TIME for $(length(TestFiles)) test files, $(hash)", elapsed])
    CSV.write(TestTimingFileName, timing_data)
    println("Wrote $(nrow(timing_data)) rows to file $TestTimingFileName")
end

elapsedclock = time() - startclocktime
println("Tested $(length(TestFiles)) files in $(round(elapsedclock, digits=1)) seconds.")

end # module BlackBoxOptimTests