text:
#Gadfly.set_default_plot_size(24cm, 12cm)
# Loss curves
function plot_loss_history(net, last_train_every, full_train_every, val_every)

    # Number of batches run
    num_batches = length(net.last_train_loss) * last_train_every

    # Last batch training data
    xids = last_train_every:last_train_every:num_batches
    layer3 = layer(x=xids, y=net.last_train_loss, Geom.point, Geom.line,
    Theme(default_color=colorant"lightblue"))

    # Full training data
    xids = full_train_every:full_train_every:(num_batches - 1)
    xids = [xids; num_batches]
    layer2 = layer(x=xids, y=net.full_train_loss, Geom.point, Geom.line,
    Theme(default_color=colorant"blue"))

    # Validation data
    xids = val_every:val_every:(num_batches - 1)
    xids = [xids; num_batches]
    layer1 = layer(x=xids, y=net.val_loss, Geom.point, Geom.line,
    Theme(default_color=colorant"red"))

    # Plot
    plot(layer1, layer2, layer3,
         Guide.xlabel("epoch"), Guide.ylabel("loss"),
         Coord.cartesian(ymin=0.0),
         Guide.manual_color_key("", ["last batch training loss    ",
                                     "full training loss    ",
                                     "validation loss"],
                                    ["lightblue", "blue", "red"]),
         Theme(key_position = :top)
         )
end

# Display colour images that are stored in columns
function display_rgb_cols(A; extrema = (0.0, 1.0), scale::Int = 1)

    # Size and number of images to display
    num_input, num_imgs = size(A)
    @assert num_input % 3 == 0
    channel_size = div(num_input, 3)

    # Error if too big - max 450
    num_imgs > 450 && throw("Maximum number of images is 450. Input a smaller subset.")

    # Dimensions of individual images - must be squares
    dim = Int(sqrt(channel_size))
    channel_size == dim^2 || throw("Individual images must be square")
    dimp1 = dim + 1    # Useful for padding i.e. adding a pixel strip between individual images

    # Shift to range [0, 1] using provided extrema
    minval, maxval = extrema
    diffval = maxval - minval
    A = (A .- minval) ./ diffval

    # Split the input into 3 arrays - one for each RGB channel
    R = A[1:channel_size, :]                             # 1st channel_size rows
    G = A[(channel_size + 1):(channel_size * 2), :]      # 2nd channel_size rows
    B = A[(channel_size * 2 + 1):(channel_size * 3), :]  # 3rd channel_size rows

    # Initialise image canvas rectangle
    cols = round(Int, sqrt(num_imgs * 2))   # Tried also φ
    rows = ceil(Int, num_imgs / cols)
    I = ones(RGB, dim * rows + rows - 1, dim * cols + cols - 1)

    # Transfer features to the canvas
    for x = 0:(rows - 1)
        for y = 0:(cols - 1)

            # Number of the image - count through the previous rows and cols
            img = x * cols + y + 1

            # Stop after through all images - will leave rest of canvas blank
            img > num_imgs && break

            # Indices of this image on the main canvas
            x_start, x_end = x * dimp1 + 1, x * dimp1 + dim
            y_start, y_end = y * dimp1 + 1, y * dimp1 + dim

            # Fill the RGB array
            I[x_start:x_end, y_start:y_end] = RGB.(R[:, img], G[:, img], B[:, img])

        end
    end

    # If not rescaling then we are done
    scale == 1 && return I

    # Optionally rescale the image up for bigger display
    Irows, Icols = size(I)
    ISrows, IScols = Irows * scale, Icols * scale
    IS = zeros(RGB, ISrows, IScols)
    for j = 1:IScols
        for i = 1:ISrows
            IS[i, j] = I[ceil(Int, i / scale), ceil(Int, j / scale)]
        end
    end

    # Return image
    return IS
end

# Function to display network
function display_rgb_weights(A; scale::Int = 1)

    # Size and number of images to display
    num_input, num_imgs = size(A)
    @assert num_input % 3 == 0
    channel_size = div(num_input, 3)

    # Error if too big - max 450
    num_imgs > 450 && throw("Maximum number of images is 450. Input a smaller subset.")

    # Dimensions of individual images - must be squares
    dim = Int(sqrt(channel_size))
    channel_size == dim^2 || throw("Individual images must be square")
    dimp1 = dim + 1    # Useful for padding i.e. adding a pixel strip between individual images

    # Split the input into 3 arrays - one for each RGB channel
    R = A[1:channel_size, :]                             # 1st channel_size rows
    G = A[(channel_size + 1):(channel_size * 2), :]      # 2nd channel_size rows
    B = A[(channel_size * 2 + 1):(channel_size * 3), :]  # 3rd channel_size rows

    # Max neuron firing of each channel and normalise
    R = (R ./ maximum(abs(R), 1) .+ 1) ./ 2
    G = (G ./ maximum(abs(G), 1) .+ 1) ./ 2
    B = (B ./ maximum(abs(B), 1) .+ 1) ./ 2

    # Initialise image canvas rectangle
    cols = round(Int, sqrt(num_imgs * 2))   # Tried also φ
    rows = ceil(Int, num_imgs / cols)
    I = ones(RGB, dim * rows + rows - 1, dim * cols + cols - 1)

    # Transfer features to the canvas
    for x = 0:(rows - 1)
        for y = 0:(cols - 1)

            # Number of the image - count through the previous rows and cols
            img = x * cols + y + 1

            # Stop after through all images - will leave rest of canvas blank
            img > num_imgs && break

            # Indices of this image on the main canvas
            x_start, x_end = x * dimp1 + 1, x * dimp1 + dim
            y_start, y_end = y * dimp1 + 1, y * dimp1 + dim

            # Fill the RGB array
            I[x_start:x_end, y_start:y_end] = RGB.(R[:, img], G[:, img], B[:, img])

        end
    end

    # If not rescaling then we are done
    scale == 1 && return I

    # Optionally rescale the image up for bigger display
    Irows, Icols = size(I)
    ISrows, IScols = Irows * scale, Icols * scale
    IS = zeros(RGB, ISrows, IScols)
    for j = 1:IScols
        for i = 1:ISrows
            IS[i, j] = I[ceil(Int, i / scale), ceil(Int, j / scale)]
        end
    end

    # Return image
    return IS
end


text:
# SPDX-License-Identifier: X11
# 2020-11-30
# Hands (300pt)

a,b,x,y = parse.(Int, (split ∘ readline)())
d = abs(a - b) - (b < a ? 1 : 0)
println(min(y, 2x) * d + x)

ext: jl

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
for (func, elty) in [(:CLBlastCherk, ComplexF32), (:CLBlastZherk, ComplexF64)]

    @eval function $func(layout::CLBlastLayout,
                         triangle::CLBlastTriangle, a_transpose::CLBlastTranspose,
                         n::Integer, k::Integer,
                         alpha::$elty,
                         a_buffer::cl.CL_mem, a_offset::Integer, a_ld::Integer,
                         beta::$elty,
                         c_buffer::cl.CL_mem, c_offset::Integer, c_ld::Integer,
                         queue::cl.CmdQueue, event::cl.Event)
        err = ccall(
            ($(string(func)), libCLBlast),
            cl.CL_int,
            (Cint, Cint, Cint, Csize_t, Csize_t, $elty, Ptr{Cvoid}, Csize_t, Csize_t,
              $elty, Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}),
            Cint(layout), Cint(triangle), Cint(a_transpose), n, k, alpha,
              a_buffer, a_offset, a_ld, beta, c_buffer, c_offset, c_ld, Ref(queue), Ref(event)
        )
        if err != cl.CL_SUCCESS
            println(stderr, "Calling function $(string($func)) failed!")
            throw(CLBlastError(err))
        end
        return err
    end

    @eval function herk!(uplo::Char, trans::Char, α::Number,
                         A::cl.CLArray{$elty,2}, β::Number, C::cl.CLArray{$elty,2};
                         queue::cl.CmdQueue=cl.queue(C))
        # check and convert arguments
        if trans == 'N'
            a_transpose = CLBlastTransposeNo
        elseif trans == 'C'
            a_transpose = CLBlastTransposeConjugate
        else
            throw(ArgumentError("Transpose marker `trans` is $(trans) but only 'N' and 'C' are allowed."))
        end
        if uplo == 'U'
            triangle = CLBlastTriangleUpper
        elseif uplo == 'L'
            triangle = CLBlastTriangleLower
        else
            throw(ArgumentError("Upper/lower marker `uplo` is $(uplo) but only 'U' and 'L' are allowed."))
        end
        n = size(C,1)
        if n != size(C,2)
            throw(DimensionMismatch("`C` has dimensions $(size(C)) but must be square."))
        end
        nn = size(A, trans == 'N' ? 1 : 2)
        if nn != n
            throw(DimensionMismatch("`C` has size ($n,$n), corresponding dimension of `A` is $nn."))
        end
        k  = size(A, trans == 'N' ? 2 : 1)
        alpha = convert($elty, α)
        beta  = convert($elty, β)
        layout = CLBlastLayoutColMajor

        # output event
        event::cl.Event = cl.Event(C_NULL)

        $func(layout, triangle, a_transpose,
              n, k,
              alpha,
              pointer(A), 0, size(A,1),
              beta,
              pointer(C), 0, size(C,1),
              queue, event)

        # wait for kernel
        cl.wait(event)

        C
    end

end


text:
using DrWatson
using Plots
@quickactivate :Swalbe

sys = Swalbe.SysConst(Lx=100, Ly=100, g=-0.001, γ=0.0005)
h = Swalbe.run_rayleightaylor(sys, "CPU"; h₀=1.0, ϵ=0.01, verbos=true)

heatmap(h[1], aspect_ratio=1, c=:viridis)
gui()
readline()


text:
module Client # client

using Sockets

using ..Common: Error, SSConfig, Cipher, SSConnection, close, gensubkey, read, write, ioCopy, init_write, init_read, eof

@inline function handShake(conn::TCPSocket, buff::Array{UInt8})
    nbytes, err = read(conn, buff)
    if err != nothing
        return err
    end

    if buff[1] != 0x05
        return "Not a Socks5 Client"
    end

    if 0x00 in buff[3:nbytes]
        err = write(conn, [0x05; 0x00], 2)
        if err != nothing
            return err
        end
    else
        write(conn, [0x05; 0xFF], 2) != nothing && return "Write Error"
        return "Not a Valid Authentication"
    end

    nbytes, err = read(conn, buff, 3)
    if err != nothing
        return err
    end

    if buff[1:3] != [0x05; 0x01; 0x00]
        return "Not a Supported CMD"
    end

    nbytes, err = read(conn, buff)
    if err != nothing
        return err
    end

    if getipaddr() isa IPv4
        err = write(conn, [0x05; 0x00; 0x00; 0x01; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00], 10)
        if err != nothing
            return err
        end
    else
        err = write(conn, [0x05; 0x00; 0x00; 0x04; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00; 0x00], 22)
        if err != nothing
            return err
        end
    end

    return nbytes
end

@inline function handleConnection(conn::TCPSocket, ssConn::SSConnection) # Client
    while true
        buff = Array{UInt8}(undef, 272)
        nbytes = handShake(conn, buff)
        if !(nbytes isa Integer)
            break
        end

        err = init_write(ssConn)
        if err != nothing
            break
        end

        err = write(ssConn, buff, nbytes)
        if err != nothing
            break
        end

        @async ioCopy(conn, ssConn)

        err = init_read(ssConn)
        if err != nothing
            break
        end

        buff = nothing
        ioCopy(ssConn, conn) # First

        close(conn)
        while !eof(conn)
            sleep(1)
        end

        break
    end

    close(conn)
    close(ssConn)
end

@inline function tcpClient(config::SSConfig, cipher::Cipher)
    server = try
        listen(getipaddr() isa IPv4 ? IPv4(0) : IPv6(0), config.lisPort)
    catch err
        return
    end

    while isopen(server)
        conn = try
            accept(server)
        catch err
            break
        end

        client = try
            connect(config.host, config.port)
        catch err
            close(conn)
            continue
        end

        @async handleConnection(
            conn,
            SSConnection(client, cipher, zeros(UInt8, cipher.ivlen), zeros(UInt8, cipher.ivlen), Array{UInt8}(undef, cipher.taglen + 2), nothing, nothing)
        )
    end

    close(server)
end

@inline function tcpClient(configs::Array{SSConfig}, ciphers::Array{Cipher})
    server = try
        listen(getipaddr() isa IPv4 ? IPv4(0) : IPv6(0), configs[1].lisPort)
    catch err
        return
    end

    nServers = [1:length(configs)...]

    while isopen(server)
        conn = try
            accept(server)
        catch err
            break
        end

        n = rand(nServers, 1)
        cipher = ciphers[n]
        config = configs[n]

        client = try
            connect(config.host, config.port)
        catch err
            close(conn)
            continue
        end

        @async handleConnection(
            conn,
            SSConnection(client, cipher, zeros(UInt8, cipher.ivlen), zeros(UInt8, cipher.ivlen), Array{UInt8}(undef, cipher.taglen + 2), nothing, nothing)
        )
    end

    close(server)
end

@inline function runClient(config::SSConfig)
    cipher = Cipher(config)

    # @async udpClient()
    tcpClient(config, cipher)
end

@inline function runClient(configs::Array{SSConfig})
    ciphers = Cipher.(configs)

    # @async udpclient()
    tcpClient(configs, ciphers)
end

end # module


text:
import .Utils

abstract type AbstractSum{OpT, StringT, CoeffT} end

const ZERO_ATOL = 1e-12

####
#### Constructors
####

## Should we use this ?
op_strings(as::AbstractSum) = as.strings

function Base.similar(ps::AbstractSum{OpT, T, V}, n=0) where {OpT, W, C, V <:Vector{C}, T <: Vector{Vector{W}}}
    m = size(ps, 2)
    strings = [Vector{W}(undef, m) for i in 1:n]
    coeffs = Vector{C}(undef, n)
    return strip_typeof(ps)(strings, coeffs; already_sorted=true)
end

function Base.copy(as::AbstractSum)
    (new_strings, new_coeffs) = (copy.(as.strings), copy(as.coeffs))
    return strip_typeof(as)(new_strings, new_coeffs; already_sorted=true)
end

####
#### IO
####

function Base.show(io::IO, asum::AbstractSum)
    (m, n) = size(asum)
    print(io, m, "x", n, " ", typeof(asum), ":\n")
    for i::Integer in eachindex(asum)
        _show_abstract_term(io, asum[i])
        if i != lastindex(asum)
            print(io, "\n")
        end
    end
end

Base.show(m::MIME{Symbol("text/input")}, asum::AbstractSum) = show(stdout, m, asum)
function Base.show(io::IO, mime::MIME{Symbol("text/input")}, asum::AbstractSum)
    print(io, strip_typeof(asum), "(")
    print(io, "[")
    for term in asum
        show(io, mime, term)
        print(io, ",")
    end
    print(io, "])")
end

####
#### Canonicalization / sorting
####

function sort_and_sum_duplicates!(asum::AbstractSum)
    sort_and_sum_duplicates!(asum.strings, asum.coeffs)
    return asum
end

sum_duplicates!(asum::AbstractSum) = sum_duplicates!(asum.strings, asum.coeffs)

function sort_and_sum_duplicates!(terms, coeffs)
    sort_sums!(terms, coeffs)
    sum_duplicates!(terms, coeffs)
    remove_zeros!(terms, coeffs)
    return nothing
end

## This is expensive. Most time is spent in sortperm.
## There is no ThreadsX.sortperm, only sort.
Base.sort!(asum::AbstractSum; alg=MergeSort) = (sort_sums!(asum.strings, asum.coeffs; alg=alg); asum)

function sort_sums!(strings, coeffs; alg=MergeSort)
    p = sortperm(strings; alg=alg)
    permute!(strings, p)
    permute!(coeffs, p)
    return nothing
end

# TODO: may want to make this a method of zchop
"""
    remove_zeros!(asum::AbstractSum)
    remove_zeros!(terms, coeffs)

Remove terms from `asum` with coefficient (approximately) equal to zero.
If `terms` and `coeffs` are supplied, then elements are deleted from both `terms`
and `coeffs` at indices corresponding to vanishing elements of `coeff`.
"""
function remove_zeros!(asum::AbstractSum, atol=ZERO_ATOL)
    remove_zeros!(asum.strings, asum.coeffs, atol)
    return asum
end

function remove_zeros!(terms::Vector, coeffs::Vector, atol=ZERO_ATOL)
    # ThreadsX is very slow for small arrays
    # The findall(iszero, coeffs) is 500ns for two non-zero floats. What is wrong?
    # Appears to be this: iszero.(array) is taking almost all the time.
    # The following is what Base does, but writing it out is faster. A bug?
    # if length(coeffs) > 10^10
    #     inds = ThreadsX.findall(Utils.isapprox_zero.(coeffs))
    # else
    inds::Vector{Int} = findall(Utils.isapprox_zero.(coeffs, atol))
    #    end
    if ! isempty(inds)
        deleteat!(coeffs::Vector, inds)
        deleteat!(terms::Vector, inds)
    end
    return nothing
end

## Modeled on code in unique! for sorted input
"""
    sum_duplicates!(_op_strings, coeffs)

Find groups of terms whose members differ only in the coefficient.
Replace each group by one term carrying the sum of the coefficients
in that group. This routine assumes that `_op_strings` are sorted.
"""
function sum_duplicates!(_op_strings, coeffs)
    last_pauli::eltype(_op_strings) = first(_op_strings)
    coeff = first(coeffs)
    k = 2
    @inbounds for j in 2:length(_op_strings)
        if _op_strings[j] != last_pauli
            last_pauli = _op_strings[k] = _op_strings[j]
            coeffs[k] = coeffs[j]
            k += 1
        else
            coeffs[k-1] += coeffs[j]
        end
    end
    resize!(_op_strings, k-1)
    resize!(coeffs, k-1)
    return nothing
end

####
#### Container interface
####

## Fails for empty asum
function Base.size(asum::AbstractSum)
    n = isempty(asum) ? 0 : length(first(asum))
    return (length(asum), n)
end
Base.size(asum::AbstractSum, i::Integer) = size(asum)[i]

# Enables using `findall`, for instance.
# Fallback methods for `values` and `pairs` are OK.
Base.keys(asum::AbstractSum) = eachindex(asum)

# Iterate uses getindex to return `AbstractTerm`s.
function Base.iterate(asum::AbstractSum, state=1)
    state > lastindex(asum) && return nothing
    return (asum[state], state + 1)
end

for func in (:length, :eachindex, :lastindex, :firstindex)
    @eval begin
        Base.$func(as::AbstractSum, args...) = return $func(as.coeffs, args...)
    end
end

Base.getindex(asum::AbstractSum, j::Integer) =
    term_type(typeof(asum))(asum.strings[j], asum.coeffs[j])

Base.getindex(asum::AbstractSum, j::Integer, k::Integer) = asum.strings[j][k]
Base.getindex(asum::AbstractSum, i1::Integer, inds2) = asum.strings[i1][inds2]

## TODO: Use already_sorted flag ?
## I think no, because asum[5:-1:1] should sort the reversed result.
## Unless, someone wants to use indexing to get something other than canonical order.
Base.getindex(asum::AbstractSum, inds) = strip_typeof(asum)(asum.strings[inds], asum.coeffs[inds]) # already_sorted=true)

## Should this return a sum or an array of arrays ?
## Currently it returns a sum. Note that selected substrings may be identical, and coeffs combined
Base.getindex(asum::AbstractSum, inds1, inds2) = strip_typeof(asum)([x[inds2] for x in asum.strings[inds1]], asum.coeffs[inds1])

####
#### Compare / predicates
####

# This will fail for an empty `psum`. Use type info instead.
# There no well-defined `one` for `PauliSum`. It depends on the
# width of the string.
function Base.one(asum::AbstractSum)
    t = one(first(asum))
    strip_typeof(asum)([op_string(t)], [t.coeff]; already_sorted=true)
end

function Base.isone(asum::AbstractSum)
    return length(asum) == 1 && isone(first(asum))
end

function Base.:(==)(asum1::AbstractSum, asum2::AbstractSum)
    if length(asum1) != length(asum2)
        return false
    end
    if length(asum1) > 10^3
        # This is 8x-12x faster for large arrays, 10^4 or 5.
        return ThreadsX.all(i -> asum1[i] == asum2[i], eachindex(asum1))
    end
    return all(i -> asum1[i] == asum2[i], eachindex(asum1))
end

"""
    reverse(ps::AbstractSum)

Reverse qubit order in `ps` and sort terms.
"""
Base.reverse(as::AbstractSum) = reverse!(copy(as))

"""
    reverse!(ps::AbstractSum)

Reverse qubit order in `ps` in place and sort terms.
"""
function Base.reverse!(as::AbstractSum)
    strings = as.strings
    @inbounds for i in eachindex(strings)
        reverse!(strings[i])
    end
    return Base.sort!(as)
end

"""
    add!(psum::AbstractSum, pt::AbstractTerm...)

Add `AbstractTerm`s to `psum` in place, assuming `psum` is sorted and has no repeated
strings. Either a new term is inserted, or the coefficient is added to an existing
term. After adding the `pt`, `psum` will be left sorted, with no duplicates, and
no zero coefficients. Use `push!` to insert a term at the end of `psum` with no
simplification performed.
"""
add!(asum::AbstractSum, term::AbstractTerm) = add!(asum, op_string(term), term.coeff)

## TODO: pass higher order function to allow + and - instead of hard coded +=
function add!(asum::AbstractSum, op_string, coeff)
    inds = searchsorted(asum.strings, op_string)
    if length(inds) == 0 # op_string not found, add a new term
        insert_tuple!(asum, first(inds), (op_string, coeff))
    elseif length(inds) == 1 # one element equal to op_string
        i = first(inds) # get the (single) index
        @inbounds asum.coeffs[i] += coeff # add p to existing term
        @inbounds if Utils.isapprox_zero(asum.coeffs[i])
            @inbounds deleteat!(asum, i)
        end
    else
        throw(ErrorException("Duplicate terms found in operator sum."))
    end
    return asum
end

"""
    add!(to::AbstractSum, from::AbstractSum)

Adds the terms in `from` to `to` in place. `to` is mutated. `from` is not.
"""
function add!(to::AbstractSum, from::AbstractSum)
    for p in from
        add!(to, p)
    end
    return to
end

####
#### Updating / adding elements
####

"""
    insert!(ps::AbstractSum, ind, p::AbstractTerm)

Insert `p` into `ps` without sorting resulting `ps`.
"""
Base.insert!(ps::AbstractSum, ind, p::AbstractTerm) = insert_tuple!(ps, ind, (op_string(p), p.coeff))

@inline function insert_tuple!(ps::AbstractSum, ind, (paulis, coeff))
    insert!(ps.strings, ind, paulis)
    insert!(ps.coeffs, ind, coeff)
    return ps
end

function Base.deleteat!(ps::AbstractSum, args...)
    deleteat!(ps.coeffs, args...)
    deleteat!(ps.strings, args...)
    return ps
end

function Base.deleteat!(ps::AbstractSum, ind)
    deleteat!(ps.coeffs, ind)
    deleteat!(ps.strings, ind)
    return ps
end

"""
    push!(psum::AbstractSum, ps::AbstractTerm...)

Push `ps` to the end of `psum` without regard to order or possible duplication.

See `sort_and_sum_duplicates!`.
"""
function Base.push!(psum::AbstractSum, ps::AbstractTerm...)
    for p in ps
        push!(psum.strings, op_string(p))
        push!(psum.coeffs, p.coeff)
    end
    return psum
end

## TODO: This should probably not be a method of push!
function Base.push!(psum::AbstractSum, (string, coeff))
    push!(psum.strings, string)
    push!(psum.coeffs, coeff)
end

function Base.append!(to::AbstractSum, from::AbstractSum)
    for t in from
        push!(to, t)
    end
    return to
end

####
#### Algebra / mathematical operations
####

Base.:+(term::T, terms::T...) where {T <: AbstractTerm} = sum_type(T)([term, terms...])

function Base.:+(ps0::AbstractSum, pss::AbstractSum...)
    ps_out = copy(ps0)
    for ps in pss
        add!(ps_out, ps)
    end
    return ps_out
end

function Base.sum(v::Vector{<:AbstractSum})
    isempty(v) && throw(ArgumentError("Can't sum empty vector of `AbstractSum`."))
    v0, vrest = Iterators.peel(v)
    ps_out = copy(v0)
    for vi in vrest
        add!(ps_out, vi)
    end
    return ps_out
end

## TODO: Do something more efficient here.
function Base.:-(ps1::T, ps2::T) where T <: AbstractSum
    return ps1 + -ps2
end

## TODO: Do something more efficient here.
function Base.:-(pt1::T, pt2::T) where T <: AbstractTerm
    return sum_type(T)([pt1, -one(pt2.coeff) * pt2])
end

function Base.:-(psum::AbstractSum)
    strip_typeof(psum)(psum.strings, -one(eltype(psum.coeffs)) .* psum.coeffs; already_sorted=true)
end

function Base.:*(n::Number, asum::AbstractSum)
    strip_typeof(asum)(asum.strings, n .* asum.coeffs; already_sorted=true)
end

Base.:*(asum::AbstractSum, n::Number) = n * asum

function Base.:/(psum::AbstractSum, n::Number)
    strip_typeof(psum)(psum.strings, psum.coeffs ./ n; already_sorted=true)
end

function mul!(asum_out::T, term::AbstractTerm, asum::T) where T <: AbstractSum
    @inbounds for j in eachindex(asum)
        new_term = term * asum[j]
        asum_out.coeffs[j] = new_term.coeff
        asum_out.strings[j] = op_string(new_term)
    end
    return asum_out
end

function Base.:*(term::AbstractTerm, asum::AbstractSum)
    asum_out = strip_typeof(asum)(similar(asum.strings), similar(asum.coeffs); already_sorted=true)
    return sort!(mul!(asum_out, term, asum))
end

"""
    *(as1::AbstractSum, as2::AbstractSum)

Multiply `as1` and `as2` returning another `AbstractSum`

# Examples
```jldoctest
julia> a = rand_op_sum(Pauli, 5, 3); ma = Matrix(a);

julia> b = rand_op_sum(Pauli, 5, 3); mb = Matrix(b);

julia> Matrix(a * b * a) == ma * mb * ma
true
```
"""
function Base.:*(as1::AbstractSum, as2::AbstractSum)
    ## Using  a buffer and mul! ought to be faster, but it is not.
    ## because mul! takes no time here compared to add!
    asum_out = similar(as1)
    asum_temp = strip_typeof(as2)(similar(as2.strings), similar(as2.coeffs); already_sorted=true)
    @inbounds for i in 1:length(as1)
        mul!(asum_temp, as1[i], as2)
        append!(asum_out, asum_temp)
    end
    return sort_and_sum_duplicates!(asum_out)
end

"""
    filter(f, as::AbstractSum)

Return a copy of `as` keeping only terms for which `f` is true.

# Examples
```julia
julia> filter(x -> count_bodies(x) == 1, fermi_sum)
IIIN * -0.47594871522096355
IINI * -0.47594871522096355
INII * -1.2524635735648981
NIII * -1.2524635735648981
```
"""
function Base.filter(f, as::AbstractSum)
    inds = findall(f, as)
    return strip_typeof(as)(as.strings[inds], as.coeffs[inds])
end


text:
## Define the derivative of a qnm function

#Define the derivative of a HeunConfluentRadial function
function ∂r(ψᵣ::HeunConfluentRadial)
    η = ψᵣ.η;
    α = ψᵣ.α;
    ξ = ψᵣ.ξ;
    ζ = ψᵣ.ζ;
    r₊ = ψᵣ.r₊
    r₋ = ψᵣ.r₋
    aₙ = ψᵣ.coeffs
    """Add a sum over different copies of qnm with some
    change"""
    Ψη = HeunConfluentRadial(η-1,α,ξ,ζ,r₊,r₋,aₙ)
    Ψξ = HeunConfluentRadial(η,α,ξ-1,ζ,r₊,r₋,aₙ)
    aₙshift = convert(Vector{eltype(aₙ)},circshift(aₙ,-1))
    aₙshift[end] = zero(eltype(aₙshift))
    nn = 1:length(aₙshift)
    aₙshift = aₙshift .*nn
    aₙshift_static = similar_type(aₙ)(aₙshift)
    Ψaₙ = HeunConfluentRadial(η-1,α+1,ξ,ζ,r₊,r₋,aₙshift_static)
    (im*(η-α))*Ψη + (im*ξ)*Ψξ + ζ*ψᵣ - Ψaₙ
end

#Define the derivation of the qnm function
function ∂r(Ψ::QuasinormalModeFunction)
    s = Ψ.s; l = Ψ.l; m = Ψ.m; n = Ψ.n; a = Ψ.a; ω = Ψ.ω; Alm = Ψ.Alm;
    ψᵣ = Ψ.R
    η = ψᵣ.η;
    α = ψᵣ.α;
    ξ = ψᵣ.ξ;
    ζ = ψᵣ.ζ;
    r₊ = ψᵣ.r₊
    r₋ = ψᵣ.r₋
    aₙ = ψᵣ.coeffs
    """Add a sum over different copies of qnm with some
    change"""
    Ψη = HeunConfluentRadial(η-1,α,ξ,ζ,r₊,r₋,aₙ)
    Ψηf = QuasinormalModeFunction(s,l,m,n,a,ω,Alm,Ψη,Ψ.S)
    Ψξ = HeunConfluentRadial(η,α,ξ-1,ζ,r₊,r₋,aₙ)
    Ψξf = QuasinormalModeFunction(s,l,m,n,a,ω,Alm,Ψξ,Ψ.S)
    aₙshift = convert(Vector{eltype(aₙ)},circshift(aₙ,-1))
    aₙshift[end] = zero(eltype(aₙ))
    nn = 1:length(aₙshift)
    #print("here")
    aₙshift = aₙshift .*nn
    aₙshift_static = similar_type(aₙ)(aₙshift)
    Ψaₙ = HeunConfluentRadial(η-1,α+1,ξ,ζ,r₊,r₋,aₙshift_static)
    Ψaₙf = QuasinormalModeFunction(s,l,m,n,a,ω,Alm,Ψaₙ,Ψ.S)
    (im*(η-α))*Ψηf + (im*ξ)*Ψξf + ζ*Ψ - Ψaₙf
end

function ∂θ(S::SpinWeightedSpheroidal)
    s = S.s; m = S.m; l=S.l; Cllʼ = S.Cllʼ;
    lmin = S.lmin; lmax = S.lmax;
    ψm1 = SpinWeightedSpheroidal(s-1,l,m,Cllʼ,lmin,lmax)
    ψp1 = SpinWeightedSpheroidal(s+1,l,m,Cllʼ,lmin, lmax)
    Ap1 = -0.5*sqrt((l-s)*(l+s+1))
    Am1 = +0.5*sqrt((l+s)*(l-s+1))
    Ap1*ψp1+Am1*ψm1
end

#Define the derivation of the qnm function
function ∂θ(Ψ::QuasinormalModeFunction)
    s = Ψ.s; l = Ψ.l; m = Ψ.m; n = Ψ.n; a = Ψ.a; ω = Ψ.ω; Alm = Ψ.Alm;
    Cllʼ = Ψ.S.Cllʼ;
    lmin = Ψ.S.lmin; lmax = Ψ.S.lmax;
    ψm1 = SpinWeightedSpheroidal(s-1,l,m,Cllʼ,lmin,lmax)
    ψp1 = SpinWeightedSpheroidal(s+1,l,m,Cllʼ,lmin, lmax)
    Ap1 = -0.5*sqrt((l-s)*(l+s+1))
    Am1 = +0.5*sqrt((l+s)*(l-s+1))
    """Add a sum over different copies of qnm with some
    change"""
    Ψm1 = QuasinormalModeFunction(s-1,l,m,n,a,ω,Alm,Ψ.R,ψm1)
    Ψp1 = QuasinormalModeFunction(s+1,l,m,n,a,ω,Alm,Ψ.R,ψp1)
    Ap1*Ψp1+Am1*Ψm1
end

function ∂t(Ψ::QuasinormalModeFunction)
    ω = Ψ.ω;
    -im*ω*Ψ
end

function ∂ϕ(Ψ::QuasinormalModeFunction)
    m = Ψ.m;
    im*m*Ψ
end

function ∂r(Ψ::LinearCombinationOf{T}) where T
    sum(v*∂r(k) for (k,v) in Ψ.dict)
end

function ∂θ(Ψ::LinearCombinationOf{T}) where T
    sum(v ≈ zero(v) ? v : v*∂θ(k) for (k,v) in Ψ.dict)
end

function ∂t(Ψ::LinearCombinationOf{T}) where T
    sum(v ≈ zero(v) ? v : v*∂t(k) for (k,v) in Ψ.dict)
end

function ∂ϕ(Ψ::LinearCombinationOf{T}) where T
    sum(v ≈ zero(v) ? v : v*∂ϕ(k) for (k,v) in Ψ.dict)
end

export ∂r, ∂θ, ∂t, ∂ϕ


text:
# ---
# title: 1559. Detect Cycles in 2D Grid
# id: problem1559
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/detect-cycles-in-2d-grid/description/>
# hidden: true
# ---
#
# Given a 2D array of characters `grid` of size `m x n`, you need to find if
# there exists any cycle consisting of the **same value**  in `grid`.
#
# A cycle is a path of **length 4  or more** in the grid that starts and ends at
# the same cell. From a given cell, you can move to one of the cells adjacent to
# it - in one of the four directions (up, down, left, or right), if it has the
# **same value** of the current cell.
#
# Also, you cannot move to the cell that you visited in your last move. For
# example, the cycle `(1, 1) -> (1, 2) -> (1, 1)` is invalid because from `(1,
# 2)` we visited `(1, 1)` which was the last visited cell.
#
# Return `true` if any cycle of the same value exists in `grid`, otherwise,
# return `false`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/15/1.png)**
#
#
#
#     Input: grid = [["a","a","a","a"],["a","b","b","a"],["a","b","b","a"],["a","a","a","a"]]
#     Output: true
#     Explanation: There are two valid cycles shown in different colors in the image below:
#     ![](https://assets.leetcode.com/uploads/2020/07/15/11.png)
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/15/22.png)**
#
#
#
#     Input: grid = [["c","c","c","a"],["c","d","c","c"],["c","c","e","c"],["f","c","c","c"]]
#     Output: true
#     Explanation: There is only one valid cycle highlighted in the image below:
#     ![](https://assets.leetcode.com/uploads/2020/07/15/2.png)
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/15/3.png)**
#
#
#
#     Input: grid = [["a","b","b"],["b","z","b"],["b","b","a"]]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m <= 500`
#   * `1 <= n <= 500`
#   * `grid` consists only of lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
struct ParallelSimulation{ET <: Equation, IT <: Tuple, ST <: Solution}
    equation::ET
    integrators::IT
    solution::ST
    ncycle::Int
    run_id::String
    filename::String
end

@inline equation(sim::ParallelSimulation) = sim.equation
@inline integrators(sim::ParallelSimulation) = sim.integrators
@inline integrator(sim::ParallelSimulation, id) = sim.integrators[id]
@inline solution(sim::ParallelSimulation) = sim.solution
@inline cycles(sim::ParallelSimulation) = 1:sim.ncycle
@inline CommonFunctions.eachsample(sim::ParallelSimulation) = eachsample(solution(sim))
@inline CommonFunctions.eachsample(sim::ParallelSimulation, id::Int) = begin
    nthreads = Threads.nthreads()
    ns = nsamples(solution(sim))
    ns_thread = div(ns, nthreads)
    i1 = (id-1)*ns_thread + 1
    i2 =  id   *ns_thread
    id == nthreads ? (i1:ns) : (i1:i2)
end

function ParallelSimulation(equ::ET, ints::IT, sol::ST, run_id::String, filename::String) where {ET,IT,ST}
    @assert mod(sol.ntime, sol.nwrite) == 0
    ncycle = div(sol.ntime, sol.nwrite)
    ParallelSimulation{ET,IT,ST}(equ, ints, sol, ncycle, run_id, filename)
end

function ParallelSimulation(equ::Equation, ints::Tuple, Δt, run_id, filename, ntime; nsave=DEFAULT_NSAVE, nwrite=DEFAULT_NWRITE)
    ParallelSimulation(equ, ints, ParallelSolution(equ, Δt, ntime; nsave=nsave, nwrite=nwrite), run_id, filename)
end

function ParallelSimulation(equ::Equation, tableau::AbstractTableau, Δt, run_id, filename, ntime; kwargs...)
    ints = Tuple(Integrator(equ, tableau, Δt) for i in 1:Threads.nthreads())
    ParallelSimulation(equ, ints, Δt, run_id, filename, ntime; kwargs...)
end

function ParallelSimulation(equ::Equation, integrator, tableau::AbstractTableau, Δt, run_id, filename, ntime; kwargs...)
    ints = Tuple(integrator(equ, tableau, Δt) for i in 1:Threads.nthreads())
    ParallelSimulation(equ, ints, Δt, run_id, filename, ntime; kwargs...)
end

function run!(sim::ParallelSimulation)

    println("Running ", sim.run_id, "...")

    create_hdf5!(solution(sim), sim.filename)

    try
        # loop over integration cycles showing progress bar
        @showprogress 5 for c in cycles(sim)
            # loop over samples
            Threads.@threads for m in eachsample(sim)
                id = Threads.threadid()

                # create atomic solution
                asol = AtomicSolution(equation(sim))

                # get cache from solution
                get_initial_conditions!(solution(sim), asol, m)

                # initilize integrator
                initialize!(integrator(sim, id), asol)

                # loop over time steps
                for n in eachtimestep(solution(sim))
                    integrate!(integrator(sim, id), solution(sim), asol, m, n)
                end
            end

            write_to_hdf5(solution(sim))
            c == sim.ncycle || reset!(solution(sim))
        end
    catch ex
        if isa(ex, DomainError)
            @warn("DOMAIN ERROR")
        elseif isa(ex, ErrorException)
            @warn("Simulation exited early.")
            @warn(ex.msg)
        else
            close(solution(sim))
            throw(ex)
        end
    end

    close(solution(sim))

    return solution(sim)
end


text:
using Plots, Printf, Statistics
# pyplot()
do_viz = false

function compute_flux!(qx, qy, T, λ, dx, dy, nx, ny)
    # Threads.@threads for iy=2:ny-1
    for iy=2:ny-1
        for ix=1:nx-1
            qx[ix,iy-1] = -λ*(T[ix+1,iy]-T[ix,iy])/dx
        end
    end
    # Threads.@threads for iy=1:ny-1
    for iy=1:ny-1
        for ix=2:nx-1
            qy[ix-1,iy] = -λ*(T[ix,iy+1]-T[ix,iy])/dy
        end
    end
    return
end

function update_T!(T, qx, qy, dt, ρCp, dx, dy, nx, ny)
    # Threads.@threads for iy=2:ny-1
    for iy=2:ny-1
        for ix=2:nx-1
            T[ix,iy] = T[ix,iy] - dt/ρCp*((qx[ix,iy-1]-qx[ix-1,iy-1])/dx + (qy[ix-1,iy]-qy[ix-1,iy-1])/dy)
        end
    end
    return
end

@views function heat_2D()
    @show Threads.nthreads()
    # Physics
    lx, ly = 10.0, 10.0
    λ      = 1.0
    ρCp    = 1.0
    nt     = 200
    # Numerics
    nx, ny = 128, 128
    nout   = 10
    # Derived numerics
    dx, dy = lx/nx, ly/ny
    dt     = min(dx^2,dy^2)/ρCp/λ/4.1
    xc     = LinRange(dx/2, lx-dx/2, nx)
    yc     = LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qx     = zeros(nx-1,ny-2)
    qy     = zeros(nx-2,ny-1)
    # Initial condition
    T      = exp.(.-(xc.-lx./2.0).^2 .-(yc.-ly./2.0)'.^2)
    # Time loop
    for it = 1:nt
        if (it==11) global t0 = Base.time() end
        compute_flux!(qx, qy, T, λ, dx, dy, nx, ny)
        update_T!(T, qx, qy, dt, ρCp, dx, dy, nx, ny)
        # Visualise
        if mod(it,nout)==0 && do_viz
            display(heatmap(xc, yc, T', xlabel="lx", ylabel="ly", title="heat diffusion, it=$it", clims=(0.,1.)))
            # sleep(.01)
        end
    end
    time_s = (Base.time()-t0)
    @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", time_s, round((2/1e9*nx*ny*sizeof(lx))/(time_s/(nt-10)), sigdigits=2))
    return
end

@time heat_2D()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
module Cadmium

using HermiteDualContouring
using RegionTrees
using ConstructiveSolidGeometry
using StaticArrays
using Plots

import Plots.Shape

export Plots

function closest_key(d::Dict, new_key)
    if isempty(d)
        return new_key
    else
        for old_key in keys(d)
            if isapprox(old_key, new_key)
               return old_key
            end
        end
        return new_key
    end
end

function next_key(d::Dict, half_key, spent_half_key)
    for key in keys(d)
        if any(i->i==half_key,key) && !any(i->i==spent_half_key,key)
            println("This is reached")
            return key
        end
    end
    return []
end

function assemble_shape(hdc::HermiteDualContour)
#     p_dict = Dict{Array{Float64,1}, Array{Float64,1}}()
    p_dict=Dict()
    index_dict = Dict()
    empty_array = []#Array{Array{Float64,1},1}()
    for leaf in allleaves(hdc.root)
        if !isempty(leaf.data.residual)
            clean_keys=copy(empty_array)
            for p in leaf.data.p
                dict_key = closest_key(p_dict, p)
                push!(clean_keys, dict_key)
                p_dict[dict_key]=push!(get(p_dict, dict_key, copy(empty_array)), getindex(leaf))
            end
            index_dict[getindex(leaf)] = clean_keys
        end
    end

    previous_p = first(keys(p_dict))
    previous_idcs =pop!(p_dict, previous_p)
    previous_idx = previous_idcs[1] #hardcoding for 2D, update with more sophisticated code for 3D

    if isfinite(previous_idcs[1].data.residual)
        verts=[previous_idcs[1].data.qef_min]
    else
        verts=[]
    end

    is_closed = false

    while !isempty(p_dict) && !is_closed
        next_idx = filter(i->i!=previous_idx, previous_idcs)
        if length(next_idx)>1
            min_idx=find(i->isfinite(i.data.residual), next_idx)
            if isempty(min_idx)
                min_idx=indmax([maximum(i.boundary.widths) for i in next_idx])
            else
                min_idx = first(min_idx)
            end
        else
            min_idx = 1
        end
        #allows faces to be traversed even if they don't contain qef points, like planes
        if isfinite(next_idx[min_idx].data.residual)
            push!(verts,next_idx[min_idx].data.qef_min)
        end
        next_ps = index_dict[next_idx[min_idx]] #may be replaceable by a cell reference
        next_p = filter(i->!isapprox(previous_p,i), next_ps)
        try
            previous_idcs = pop!(p_dict, next_p...)
        catch
            is_closed = true #replace with more rigorous logic; perhaps tracking the original p and seeing when it's called again
        end
        previous_idx = next_idx[min_idx]
        previous_p=first(next_p)
    end
    return hcat(verts...)
end

function Plots.Shape(hdc::HermiteDualContour)
    verts = assemble_shape(hdc)
    Plots.Shape(verts[1,:], verts[2,:])
end

function Plots.Shape(
        csg::ConstructiveSolidGeometry.Surface,
        origin::AbstractArray,
        widths::AbstractArray,
        rtol=1e-3,
        atol=1e-3,
        surfcellmax=1e-1)

    hdc=HermiteDualContour(csg, origin, widths, rtol, atol, surfcellmax)
    Plots.Shape(hdc)
end

end # module


text:
mutable struct DescriptorSetLayout
    mVkDevice::vk.VkDevice
    mCreateInfo::DescriptorSetLayoutCreateInfo
    mPool::DescriptorPool

    mHandleRef::Ref{vk.VkDescriptorSetLayout}

    function DescriptorSetLayout(device::vk.VkDevice,
                                   info::DescriptorSetLayoutCreateInfo,
                               poolSize::UInt32)
        this = new()
        this.mVkDevice = device
        this.mCreateInfo = info

        this.mHandleRef = Ref{vk.VkDescriptorSetLayout}()
        if vk.vkCreateDescriptorSetLayout(this.mVkDevice, handleRef(this.mCreateInfo), C_NULL, this.mHandleRef) != vk.VK_SUCCESS
            error("Failed to create descriptor set layout!")
        end

        if poolSize != 0
            sizes = Dict{vk.VkDescriptorType, UInt32}()
            for b in info.mBindings
                count = get!(sizes, b.descriptorType, 0)
                count += b.descriptorCount
                sizes[b.descriptorType] = count
            end
            pinfo = DescriptorPoolCreateInfo()
            allowFreeing(pinfo)
            for s in sizes
                addSize(pinfo, s.first, s.second * poolSize);
            end
            setMaxSets(pinfo, poolSize);

            this.mPool = createDescriptorPool(this.mVkDevice, pinfo);
        end
        return this
    end

end
@class DescriptorSetLayout

function destroy(this::DescriptorSetLayout)
    vk.vkDestroyDescriptorSetLayout(this.mVkDevice, this.mHandleRef[], C_NULL)
    this.mPool.destroy()
end

function handleRef(this::DescriptorSetLayout)::Ref{vk.VkDescriptorSetLayout}
    return this.mHandleRef
end

function createDescriptorSet(this::DescriptorSetLayout)::DescriptorSet
    return DescriptorSet(this.mVkDevice, this.mPool, this)
end


text:
module IntSets

import Base: similar, copy, eltype, sizehint!, push!, union, union!, pop!,
             delete!, empty!, setdiff, setdiff!, symdiff!, symdiff, copy!, in,
             start, next, done, first, last, isempty, shift!, length, show,
             intersect, intersect!, hash, findnextnot,
             findprevnot, unsafe_getindex, unsafe_setindex!
importall Base.Operators
# export IntSet

include("intset.jl")

end # module


text:
"""
    bds(session::Session, security::AbstractString, field::AbstractString;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)

Runs a query for reference data of a security. Returns a `Vector` of named tuples.

Internally, it issues a `ReferenceDataRequest` in `//blp/refdata` service.

See also [`bdh`](@ref).

# Example

```julia
using BLPData, DataFrames
session = BLPData.Session()
result = BLPData.bds(session, "PETR4 BS Equity", "COMPANY_ADDRESS")
df = DataFrame(result)
```
"""
function bds(session::Session, security::AbstractString, field::AbstractString;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        )

    queue, corr_id = send_request(session, "//blp/refdata", "ReferenceDataRequest") do req
        push!(req["securities"], security)
        push!(req["fields"], field)

        if options != nothing
            for (k, v) in options
                req[k] = v
            end
        end
    end

    result = Vector()

    for_each_response_message_element(queue, corr_id, timeout_milliseconds=timeout_milliseconds, verbose=verbose) do element
        @assert has_name(element, "ReferenceDataResponse")
        response_element = get_choice(element)

        if has_name(response_element, "responseError")
            error("Got responseError. \n$response_element")
        end

        @assert has_name(response_element, "securityData")
        for security_data_element in get_element_value(response_element)
            @assert get_element_value(security_data_element["security"]) == security
            push_named_tuples!(result, security_data_element["fieldData"][field])
        end
    end

    return result
end

"""
    bds(session::Session, securities::Vector{T}, field::AbstractString;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        ) where {T<:AbstractString}

Runs a query for reference data of a security.
Returns a `Dict` where the key is the security name and value is a `Vector` of named tuples.

Internally, BLPData will process a `ReferenceDataRequest` request for each security in parallel.
"""
function bds(session::Session, securities::Vector{T}, field::AbstractString;
            options=nothing, # expects key->value pairs or Dict
            verbose::Bool=false,
            timeout_milliseconds::Integer=UInt32(0)
        ) where {T<:AbstractString}

    result = Dict()

    @sync for security in securities
        @async result[security] = bds($session, $security, $field, options=$options, verbose=$verbose, timeout_milliseconds=$timeout_milliseconds)
    end

    return result
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
@testset "module validation" begin
    using Atom: toplevelitems

    path = subjunkspath
    text = read(path, String)

    # basic -- finds every toplevel items with default arguments
    @test filter(toplevelitems(text)) do item
        item isa Atom.ToplevelBinding &&
        item.bind.name == "imwithdoc"
    end |> length === 3

    # don't enter non-target modules, e.g.: submodules
    @test filter(toplevelitems(text; mod = "Junk", inmod = true)) do item
        item isa Atom.ToplevelBinding &&
        item.bind.name == "imwithdoc"
    end |> length === 1 # should only find the `imwithdoc` in Junk module

    # don't include items outside of a module
    @test filter(toplevelitems(text; mod = "SubJunk", inmod = false)) do item
        item isa Atom.ToplevelBinding &&
        item.bind.name == "imwithdoc"
    end |> length === 1 # should only find the `imwithdoc` in SubJunk module
end

@testset "outline edge cases" begin
    # don't leak something in a call
    for t in ["f(call())", "f(call(), @macro)", "f(@macro)", "f(@macro, arg; kwarg = kwarg)"]
        items = toplevelitems(t)
        @test length(items) === 1
        @test items[1] isa Atom.ToplevelCall
    end

    # https://github.com/JunoLab/Juno.jl/issues/502
    @test isempty(toplevelitems("ary[ind] = some"))

    # don't leak rhs of the assignemt
    let items = toplevelitems("lhs = rhscall(somearg)")
        @test length(items) === 1
        @test items[1].bind.name == "lhs"
    end
    let items = toplevelitems("lhs = @rhsmacrocall somearg")
        @test length(items) === 1
        @test items[1].bind.name == "lhs"
    end
end


text:
using TestExtras
using Test

@timedtestset "constinferred tests" begin

    mysqrt(x; complex = true) = x >= 0 ? sqrt(x) : (complex ? im*sqrt(-x) : throw(DomainError(x, "Enable complex return values to take square roots of negative numbers")))

    @constinferred mysqrt(+3)
    @constinferred mysqrt(-3)
    for x = -1.5:0.5:+1.5
        @constinferred mysqrt($x)
        @constinferred mysqrt($(rand() < 0 ? x : -x))
    end

    @constinferred Nothing iterate(1:5)
    @constinferred Nothing iterate(1:-1)
    @constinferred Tuple{Int,Int} iterate(1:-1)

    x = (2,3)
    @constinferred +(1, x...)

    @constinferred (1:3)[2]
    a = (0.1, 0.2, 0.3)
    b = (3, 2, 1)
    @constinferred a .+ b

    x = 3.
    @constinferred mysqrt(x; complex = false)
    @constinferred_broken mysqrt(x; complex = true)
end

# ensure constinferred only evaluates argument once
inferred_test_global = 0
function inferred_test_function()
    global inferred_test_global
    inferred_test_global += 1
    true
end
@constinferred inferred_test_function()
@test inferred_test_global == 1

f25835(;x=nothing) = _f25835(x)
_f25835(::Nothing) = ()
_f25835(x) = (x,)
# A keyword function that is never type stable
g25835(;x=1) = rand(Bool) ? 1.0 : 1
# A keyword function that is sometimes type stable
h25835(;x=1,y=1) = x isa Int ? x*y : (rand(Bool) ? 1.0 : 1)

@timedtestset begin
    @test @constinferred(f25835()) == ()
    @test @constinferred(f25835(x=nothing)) == ()
    @test @constinferred(f25835(x=$(rand() < 0 ? nothing : nothing))) == ()
    @test @constinferred(f25835(x=1)) == (1,)

    @test @constinferred(h25835()) == 1
    @test @constinferred(h25835(x=2,y=3)) == 6
    @test @constinferred(Union{Float64,Int64}, h25835(x=1.0,y=1.0)) == 1
end

@timedtestset "@test" begin
    @test true
    @test 1 == 1
    @test 1 != 2
    @test ==(1, 1)
    @test ==((1, 1)...)
    @test 1 ≈ 2 atol=1
    @test strip("\t  hi   \n") == "hi"
    @test strip("\t  this should fail   \n") != "hi"
    @test isequal(1, 1)
    @test isapprox(1, 1, atol=0.1)
    @test isapprox(1, 1; atol=0.1)
    @test isapprox(1, 1; [(:atol, 0)]...)
end
@timedtestset "@test keyword precedence" begin
    # post-semicolon keyword, suffix keyword, pre-semicolon keyword
    @test isapprox(1, 2, atol=0) atol=1
    @test isapprox(1, 3, atol=0; atol=2) atol=1
end
@timedtestset "@test should only evaluate the arguments once" begin
    g = Int[]
    f = (x) -> (push!(g, x); x)
    @test f(1) == 1
    @test g == [1]

    empty!(g)
    @test isequal(f(2), 2)
    @test g == [2]
end

@timedtestset "@test_broken with fail" begin
    @test_broken false
    @test_broken 1 == 2
    @test_broken 1 != 1
    @test_broken strip("\t  hi   \n") != "hi"
    @test_broken strip("\t  this should fail   \n") == "hi"
end
@timedtestset "@test_broken with errors" begin
    @test_broken error()
    @test_broken absolute_nonsense
end
@timedtestset "@test_skip" begin
    @test_skip error()
    @test_skip true
    @test_skip false
    @test_skip gobbeldygook
end
@timedtestset "@test_warn" begin
    @test 1234 === @test_nowarn(1234)
    @test 5678 === @test_warn("WARNING: foo", begin println(stderr, "WARNING: foo"); 5678; end)
    let a
        @test_throws UndefVarError(:a) a
        @test_nowarn a = 1
        @test a === 1
    end
end


text:
workspace()
using PDFParser

function add!(dict, value)
    if haskey(dict, value)
        dict[value] += 1
    else
        dict[value] = 1
    end
end

function test()
    root = "C:/Users/hshindo/Dropbox/research/pdfparsing/data/pmc"
    dicts = Dict{String,Dict{String,Int}}()
    for dir in readdir(root)
        println(dir)
        path = joinpath(root, dir)
        trees = map(x -> parse(Tree,x), open(readlines,path))
        for tree in trees
            front = tree[1]
            front.name == "front" || continue
            #body = tree[2]
            #body.value == "body" || continue
            topdown(front) do t
                haskey(target,t.name) || return
                if haskey(dicts, t.name)
                    d = dicts[t.name]
                else
                    d = Dict{String,Int}()
                    dicts[t.name] = d
                end
                for c in t.children
                    if haskey(d, c.name)
                        d[c.name] += 1
                    else
                        d[c.name] = 1
                    end
                end
            end
        end
        for (k,d) in dicts
            println(k)
            println(d)
        end
    end
end
test()


text:
"""
    mutable struct TemperedState
        states          :: Array{Any}
        Δ               :: Vector{<:Real}
        Δ_index         :: Vector{<:Integer}
        chain_index     :: Vector{<:Integer}
        step_counter    :: Integer
        total_steps     :: Integer
        Δ_history       :: Array{<:Real, 2}
        Δ_index_history :: Array{<:Integer, 2}
        Ρ               :: Vector{AdaptiveState}
    end

A `TemperedState` struct contains the `states` of each of the parallel chains
used throughout parallel tempering as pairs of `Transition`s and `VarInfo`s,
it also stores necessary information for tempering:
- `states` is an Array of pairs of `Transition`s and `VarInfo`s, one for each
  tempered chain
- `Δ` contains the ordered sequence of inverse temperatures
- `Δ_index` contains the current ordering to apply the temperatures to each chain, tracking swaps,
    i.e., contains the index `Δ_index[i] = j` of the temperature in `Δ`, `Δ[j]`, to apply to chain `i`
- `chain_index` contains the index `chain_index[i] = k` of the chain tempered by `Δ[i]`
    NOTE that to convert between this and `Δ_index` we simply use the `sortperm()` function
- `step_counter` maintains the number of steps taken since the last swap attempt
- `total_steps` maintains the count of the total number of steps taken
- `Δ_index_history` records the history of swaps that occur in sampling by recording the `Δ_index` at each step
- `Δ_history` records the values of the inverse temperatures, these will change if adaptation is being used
- `Ρ` contains all of the information required for adaptation of Δ

Example of swaps across 4 chains and the values of `chain_index` and `Δ_index`:

Chains:        chain_index:     Δ_index:
| | | |        1  2  3  4       1  2  3  4
| | | |
 V  | |        2  1  3  4       2  1  3  4
 Λ  | |
| | | |        2  1  3  4       2  1  3  4
| | | |
|  V  |        2  3  1  4       3  1  2  4
|  Λ  |
| | | |        2  3  1  4       3  1  2  4
| | | |
"""
mutable struct TemperedState
    states          :: Array{Any}
    Δ               :: Vector{<:Real}
    Δ_index         :: Vector{<:Integer}
    chain_index     :: Vector{<:Integer}
    step_counter    :: Integer
    total_steps     :: Integer
    Δ_history       :: Array{<:Real, 2}
    Δ_index_history :: Array{<:Integer, 2}
    Ρ               :: Vector{AdaptiveState}
end

"""
For each `β` in `Δ`, carry out a step with a tempered model at the corresponding `β` inverse temperature,
resulting in a list of transitions and states, the transition associated with `β₀ = 1` is then returned with the
rest of the information being stored in the state.
"""
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model,
    spl::TemperedSampler;
    kwargs...
)
    states = [
        AbstractMCMC.step(
            rng,
            make_tempered_model(model, spl.Δ[spl.Δ_init[i]]),
            spl.internal_sampler;
            kwargs...
        )
        for i in 1:length(spl.Δ)
    ]
    return (
        states[sortperm(spl.Δ_init)[1]][1],
        TemperedState(
            states,spl.Δ, spl.Δ_init, sortperm(spl.Δ_init), 1, 1, Array{Real, 2}(spl.Δ'), Array{Integer, 2}(spl.Δ_init'), spl.Ρ
        )
    )
end
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model,
    spl::TemperedSampler,
    ts::TemperedState;
    kwargs...
)
    if ts.step_counter == spl.N_swap
        ts = swap_step(rng, model, spl, ts)
        ts.step_counter = 0
    else
        ts.states = [
            AbstractMCMC.step(
                rng,
                make_tempered_model(model, ts.Δ[ts.Δ_index[i]]),
                spl.internal_sampler,
                ts.states[i][2];
                kwargs...
            )
            for i in 1:length(ts.Δ)
        ]
        ts.step_counter += 1
    end

    ts.Δ_history = vcat(ts.Δ_history, Array{Real, 2}(ts.Δ'))
    ts.Δ_index_history = vcat(ts.Δ_index_history, Array{Integer, 2}(ts.Δ_index'))
    ts.total_steps += 1
    return ts.states[ts.chain_index[1]][1], ts  # Use chain_index[1] to ensure the sample from the target is always returned for the step
end

"""
    swap_step(rng, model, spl, ts)

Uses the internals of the passed `TemperedSampler` - `spl` - and `TemperedState` -
`ts` - to perform a "swap step" between temperatures, in accordance with the relevant
swap strategy.
"""
function swap_step(
    rng::Random.AbstractRNG,
    model,
    spl::TemperedSampler,
    ts::TemperedState
)
    L = length(ts.Δ) - 1
    sampler = spl.internal_sampler

    if spl.swap_strategy == :standard

        k = rand(rng, Distributions.Categorical(L))  # Pick randomly from 1, 2, ..., k - 1
        ts = swap_attempt(model, sampler, ts, k, spl.adapt, ts.total_steps / L)

    else

        # Define a vector to populate with levels at which to propose swaps according to swap_strategy
        levels = Vector{Int}(undef, L)
        if spl.swap_strategy == :nonrev
            if ts.step_counter % (2 * spl.N_swap) == 0
                levels = 1:2:L
            else
                levels = 2:2:L
            end
        elseif spl.swap_strategy == :randperm
            randperm!(rng, levels)
        end

        for k in levels
            ts = swap_attempt(model, sampler, ts, k, spl.adapt, ts.total_steps)
        end

    end
    return ts
end


text:
##
# Sparse BroadcastStyle
##
BroadcastStyle(::StructuredMatrixStyle{<:Diagonal}, ::BandedBlockBandedStyle) =
    BandedBlockBandedStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::StructuredMatrixStyle{<:Diagonal}) =
    BandedBlockBandedStyle()

function blockbandwidths(P::PseudoBlockMatrix{<:Any,<:Diagonal})
    bs = blocksizes(P)
    cumulsizes(bs)[1] == cumulsizes(bs)[2] || throw(DimensionMismatch())
    (0,0)
end

bandeddata(P::PseudoBlockMatrix) = bandeddata(P.blocks)
bandwidths(P::PseudoBlockMatrix) = bandwidths(P.blocks)

BroadcastStyle(::Type{<:SubArray{<:Any,2,<:PseudoBlockMatrix{<:Any,<:Diagonal},
                                NTuple{2,BlockSlice1}}}) = BandedStyle()

const BlockDiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Diagonal{VT}}

BlockDiagonal(A) = mortar(Diagonal(A))

function sizes_from_blocks(A::Diagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    BlockSizes(size.(A.diag, 1), size.(A.diag,2))
end

# Block Tridiagonal
const BlockTridiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Tridiagonal{VT}}

BlockTridiagonal(A,B,C) = mortar(Tridiagonal(A,B,C))

function sizes_from_blocks(A::Tridiagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    BlockSizes(size.(A.d, 1), size.(A.d,2))
end

blockbandwidths(A::BlockArray) = bandwidths(A.blocks)
isblockbanded(A::BlockArray) = isbanded(A.blocks)

@inline function getblock(block_arr::BlockTridiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
    @boundscheck blockcheckbounds(block_arr, K, J)
    abs(J-K) ≥ 2 && return convert(VT, Zeros{T}(blocksize(block_arr,(K,J))))
    block_arr.blocks[K,J]
end

function checksquareblocks(A)
    m,n = cumulsizes(blocksizes(A))
    m == n || throw(DimensionMismatch("blocks are not square: block dimensions are $(blocksizes(A))"))
    m
end

for op in (:-, :+)
    @eval begin
        function $op(A::BlockTridiagonal, λ::UniformScaling)
            checksquareblocks(A)
            mortar(Tridiagonal(A.blocks.dl, broadcast($op, A.blocks.d, Ref(λ)), A.blocks.du))
        end
        function $op(λ::UniformScaling, A::BlockTridiagonal)
            checksquareblocks(A)
            mortar(Tridiagonal(A.blocks.dl, broadcast($op, Ref(λ), A.blocks.d), A.blocks.du))
        end
    end
end

function replace_in_print_matrix(A::BlockDiagonal, i::Integer, j::Integer, s::AbstractString)
    bi = global2blockindex(A.block_sizes, (i, j))
    I,J = bi.I
    i,j = bi.α
    J-I == 0 ? s : Base.replace_with_centered_mark(s)
end

function replace_in_print_matrix(A::BlockTridiagonal, i::Integer, j::Integer, s::AbstractString)
    bi = global2blockindex(A.block_sizes, (i, j))
    I,J = bi.I
    i,j = bi.α
    -1 ≤ J-I ≤ 1 ? s : Base.replace_with_centered_mark(s)
end


text:
using TransformVariables
const TV = TransformVariables

struct Const{T} <: Function
    value::T
end

export asConst
struct AsConst{T} <: TV.VectorTransform
    value::T
end

asConst(x) = AsConst(x)

TV.as(c::Const) = AsConst(c.value)

function Base.show(io::IO, c::Const)
    io = IOContext(io, :compact => true)
    print(io, "_ -> ")
    print(io, c.t)
end

(c::Const)(x) = c.value

TV.dimension(t::AsConst) = 0

function TV.transform_with(flag::TV.NoLogJac, t::AsConst, x, index)
    return (t.value, TV.NoLogJac(), index)
end

function TV.transform_with(flag::TV.LogJac, t::AsConst, x, index)
    return (t.value, 0.0, index)
end

TV.inverse_eltype(t::AsConst, x::AbstractVector) = Any

function TV.inverse_at!(x::AbstractVector, index, t::AsConst, y::AbstractVector)
    return index
end


text:
module ComposeTests

using InplaceArrays.TensorValues
using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField
using InplaceArrays.Fields: Comp
using FillArrays

import InplaceArrays.Fields: ∇

np = 4
p = Point(1,2)
x = fill(p,np)

fun(x,y) = 2*x
∇fun(x,y) = VectorValue(2*one(x[1]),2*one(x[1]))
∇(::typeof(fun)) = ∇fun

v = 3.0
d = 2
f = MockField{d}(v)
fx = fill(v,np)

k = Comp(fun)
test_kernel(k,(v,v),2*v)

g = compose(fun,f,f)
gx = 2*fx
∇gx = fill(∇fun(v,v),np)
test_field(g,x,gx,grad=∇gx)

l = 10
af = Fill(f,l)
ax = fill(x,l)
ag = compose(fun,af,af)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

end # module


text:
"""
    SupervisedLossContainer(loss_function)
    SupervisedLossContainer(loss_function, λ, λⱼ)

A container class for supervised loss functions.
"""
Base.@kwdef struct SupervisedLossContainer{L,T}
    loss_function::L
    λ::T = 0.0f0
    λⱼ::T = 0.0f0
end

function (lc::SupervisedLossContainer)(soln::DeepEquilibriumSolution)
    return lc.λ * mean(abs, soln.u₀ .- soln.z_star) + lc.λⱼ * soln.jacobian_loss
end

function (lc::SupervisedLossContainer)(soln::DeepEquilibriumSolution{T}) where {T<:Tuple}
    return lc.λ * mapreduce((x, y) -> mean(abs, x .- y), +, soln.u₀, soln.z_star) +
           lc.λⱼ * soln.jacobian_loss
end

function (lc::SupervisedLossContainer)(model::Union{DeepEquilibriumNetwork,SkipDeepEquilibriumNetwork,DEQChain}, x, y;
                                       kwargs...)
    ŷ, soln = model(x; kwargs...)
    return lc.loss_function(ŷ, y) + lc(soln)
end

function (lc::SupervisedLossContainer)(model::Union{MultiScaleDeepEquilibriumNetwork,
                                                    MultiScaleSkipDeepEquilibriumNetwork}, x, ys::Tuple; kwargs...)
    yŝ, soln = model(x; kwargs...)
    return mapreduce(lc.loss_function, +, ys, yŝ) + lc(soln)
end

function (lc::SupervisedLossContainer)(model::Union{MultiScaleDeepEquilibriumNetwork,
                                                    MultiScaleSkipDeepEquilibriumNetwork}, x, y; kwargs...)
    yŝ, soln = model(x; kwargs...)
    return sum(Base.Fix2(lc.loss_function, y), yŝ) + lc(soln)
end

# Default fallback
function (lc::SupervisedLossContainer)(model, x, y; kwargs...)
    return lc.loss_function(model(x; kwargs...), y)
end


text:
# Add `Primitive` type to BalanceLaws, AtmosModel

# Vars wrapper
function prognostic_to_primitive!(
    bl::BalanceLaw,
    prim::AbstractArray,
    prog::AbstractArray,
    aux::AbstractArray,
)
    FT = eltype(prim)
    prognostic_to_primitive!(
        bl,
        Vars{vars_state(bl, Primitive(), FT)}(prim),
        Vars{vars_state(bl, Prognostic(), FT)}(prog),
        Vars{vars_state(bl, Auxiliary(), FT)}(aux),
    )
end
function primitive_to_prognostic!(
    bl::BalanceLaw,
    prog::AbstractArray,
    prim::AbstractArray,
    aux::AbstractArray,
)
    FT = eltype(prog)
    primitive_to_prognostic!(
        bl,
        Vars{vars_state(bl, Prognostic(), FT)}(prog),
        Vars{vars_state(bl, Primitive(), FT)}(prim),
        Vars{vars_state(bl, Auxiliary(), FT)}(aux),
    )
end

# By default the primitive is the prognostic
vars_state(bl::BalanceLaw, ::Primitive, FT) = vars_state(bl, Prognostic(), FT)

function prognostic_to_primitive!(bl, prim::Vars, prog::Vars, aux)
    prim_arr = parent(prim)
    prog_arr = parent(prog)
    prim_arr .= prog_arr
end
function primitive_to_prognostic!(bl, prog::Vars, prim::Vars, aux)
    prim_arr = parent(prim)
    prog_arr = parent(prog)
    prog_arr .= prim_arr
end

"""
    construct_face_auxiliary_state!(
        bl::BalanceLaw,
        aux_face::AbstractArray,
        aux_cell::AbstractArray,
        Δz::FT
    )

Default constructor
 - `bl` balance law
 - `aux_face` face auxiliary variables to be constructed
 - `aux_cell` cell center auxiliary variable
 - `Δz` cell vertical size
"""
function construct_face_auxiliary_state!(
    bl,
    aux_face::AbstractArray,
    aux_cell::AbstractArray,
    Δz::FT,
) where {FT <: Real}
    # TODO: is this safe/correct?
    aux_face .= aux_cell
end


text:
export Options

"""
Options for the continuation process.
"""
const Options = Dict{String,Any}

"""
Return the default options for a particular toolbox or problem type.
"""
function default_options end


text:
mutable struct AutoSwitch{nAlg,sAlg,tolType,T}
  count::Int
  nonstiffalg::nAlg
  stiffalg::sAlg
  is_stiffalg::Bool
  maxstiffstep::Int
  maxnonstiffstep::Int
  nonstifftol::tolType
  stifftol::tolType
  dtfac::T
  stiffalgfirst::Bool
end
AutoSwitch(nonstiffalg, stiffalg;
           maxstiffstep=10, maxnonstiffstep=3,
           nonstifftol=9//10, stifftol=9//10,
           dtfac=2, stiffalgfirst=false) =
           AutoSwitch(0, nonstiffalg, stiffalg, stiffalgfirst,
                      maxstiffstep, maxnonstiffstep, promote(nonstifftol, stifftol)...,
                      dtfac, stiffalgfirst)

function is_stiff(integrator, alg, ntol, stol, is_stiffalg)
  eigen_est, dt = integrator.eigen_est, integrator.dt
  stiffness = abs(eigen_est*dt/alg_stability_size(alg)) # call `abs` here just to be sure
  tol = is_stiffalg ? stol : ntol
  os = oneunit(stiffness)
  stiffness > os * tol
end

function (AS::AutoSwitch)(integrator)
  integrator.iter == 0 && return Int(AS.stiffalgfirst) + 1
  eigen_est, dt = integrator.eigen_est, integrator.dt
  # Successive stiffness test positives are counted by a positive integer,
  # and successive stiffness test negatives are counted by a negative integer
  AS.count = is_stiff(integrator, AS.nonstiffalg, AS.nonstifftol, AS.stifftol, AS.is_stiffalg) ?
                      AS.count < 0 ? 1 : AS.count + 1 :
                      AS.count > 0 ? -1 : AS.count - 1
  if (!AS.is_stiffalg && AS.count > AS.maxstiffstep)
    integrator.dt = dt*AS.dtfac
    AS.is_stiffalg = true
  elseif (AS.is_stiffalg && AS.count < -AS.maxnonstiffstep)
    integrator.dt = dt/AS.dtfac
    AS.is_stiffalg = false
  end
  return Int(AS.is_stiffalg) + 1
end

function AutoAlgSwitch(nonstiffalg, stiffalg; kwargs...)
  AS = AutoSwitch(nonstiffalg, stiffalg; kwargs...)
  StochasticCompositeAlgorithm((nonstiffalg, stiffalg), AS)
end

AutoSOSRI2(alg; kwargs...) = AutoAlgSwitch(SOSRI2(), alg; kwargs...)
AutoSOSRA2(alg; kwargs...) = AutoAlgSwitch(SOSRA2(), alg; kwargs...)


text:
# Automatically generated by update_leap_seconds.jl - do not edit!

const leap_index_offset = 9 # For use with searchsortedlast()

const utc_offsets = [
    typemin(Date),  # For use with searchsortedlast()
    Date(1972, 7, 1),
    Date(1973, 1, 1),
    Date(1974, 1, 1),
    Date(1975, 1, 1),
    Date(1976, 1, 1),
    Date(1977, 1, 1),
    Date(1978, 1, 1),
    Date(1979, 1, 1),
    Date(1980, 1, 1),
    Date(1981, 7, 1),
    Date(1982, 7, 1),
    Date(1983, 7, 1),
    Date(1985, 7, 1),
    Date(1988, 1, 1),
    Date(1990, 1, 1),
    Date(1991, 1, 1),
    Date(1992, 7, 1),
    Date(1993, 7, 1),
    Date(1994, 7, 1),
    Date(1996, 1, 1),
    Date(1997, 7, 1),
    Date(1999, 1, 1),
    Date(2006, 1, 1),
    Date(2009, 1, 1),
    Date(2012, 7, 1),
    Date(2015, 7, 1),
    Date(2017, 1, 1)
]

const leap_secs_expiry_date = Date(2017, 12, 28)


text:
module OptimTests

using Photon, Test

function get_model()
    Sequential(
        Conv2D(16, (3,3), relu),
        Conv2D(32, (5,5), relu),
        MaxPool2D(),
        Dense(64, relu),
        Dense(10)
    )
end

function testOptim(opt)
    model = get_model()
    dtype = getContext().dtype

    X = [randn(dtype,50,50,3,4) for _=1:5]
    Y = [randn(dtype,10,4) for _=1:5]

    workout = Workout(model, MSELoss(), opt)
    fit!(workout, zip(X,Y), epochs=2)
    @test hasmetric(workout, :loss)
end

@testset "Optim" begin
    testOptim(Adam())
    testOptim(SGD())
    testOptim(Momentum())
end

end


text:
using Compat, ClassicalLaminateTheory

forcesAndMoments =
 @compat Dict(
    :Nx => 300,
    :Ny => 100,
    :Nxy => 200,
    :Mx => 5,
    :My => 3,
    :Mxy => 1)

materialProperties =@compat Dict(
  :E1 => [20.0*10^6],
  :E2 => [1.5*10^6],
  :G12 => [1.0*10^6],
  :ν12 => [.34],
  :alpha1 => [0.2e-6],
  :alpha2 => [20e-6],
  :beta1 => [0.2e-4],
  :beta2 => [20e-4],
  :thickness => [0.0075])

laminateProperties =@compat Dict(
  :nplies => 5,
  :materials => [1, 1, 1, 1, 1],
  :repeats => 1,
  :symmetric => false,
  :symmetricrepeats => 1,
  :orientation => [0, 45, 90, 45, -45],
  :deltaTemp => -270.0,
  :deltaMoisture => 1
)

qm = qmat(materialProperties)
materialProperties |> display

println("\nQ:")
qm |> display

qb0 = qbarmat(qm, 0.0)
qb90 = qbarmat(qm, 90.0)
qb45 = qbarmat(qm, 45.0)
qbm45 = qbarmat(qm, -45.0)
println("\nQbar 0 Degrees:")
qb0 |> display
println("\nQbar 90 Degrees:")
qb90 |> display
println("\nQbar 45 Degrees:")
qb45 |> display
println("\nQbar -45 Degrees:")
qbm45 |> display

println()
createLaminate!(laminateProperties, materialProperties, forcesAndMoments)
show(laminateProperties)

println("ABD matrix:")
round(laminateProperties[:ABD], 3) |> display

println()
println("Total forces and moments:")
round(laminateProperties[:fTotal], 7) |> display


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.7.0-beta3
#     language: julia
#     name: julia-1.7
# ---

# # `NCTiles.jl` : Create Files, Embed Meta Data
#
# [NCTiles.jl](https://gaelforget.github.io/NCTiles.jl/dev/) creates [NetCDF](https://en.wikipedia.org/wiki/NetCDF) files that follow the [CF Metadata Conventions](http://cfconventions.org). It can be used either (1) in stand-alone mode or (2) in combination with [MeshArrays.jl](https://juliaclimate.github.io/MeshArrays.jl/dev/). The examples below include:
#
# 1. Writing mapped model output, on a regular `lat-lon` grid, to a single `NetCDF` file
#   - 2D example
#   - 3D example
# 2. Writing tiled model output, on `C-grid` subdomains, to a collection of `NetCDF` files
#   - 2D surface example
#   - 3D temperature example
#   - 3D staggered vector example

# ### Packages & Helper Functions

# +
#import Pkg; Pkg.add("NCDatasets"); Pkg.add("NCTiles"); Pkg.test("NCTiles")

# +
using NCTiles
include("helper_functions.jl");

pth_notebooks=joinpath(tempdir(),"tmp_JuliaClimateNotebooks")
!isdir(pth_notebooks) ? mkdir(pth_notebooks) : nothing
# -

# ### File Paths & I/O Back-End
#
# _These will be used throughout the notebook_

# +
# File Paths
inputs = joinpath(pth_notebooks,"nctiles-testcases/")
get_testcases_if_needed(inputs)

outputs = joinpath(pth_notebooks,"nctiles-newfiles/")
if ~ispath(outputs); mkpath(outputs); end

# I/O Back-End
nc=NCTiles.NCDatasets
# -

# ## Interpolated Data Examples
#
# This example uses 2D and 3D model output that has been interpolated to a rectangular half-degree grid. It reads the data from binary files, adds meta data, and then writes it all to a single `NetCDF` file per model variable.

writedir = joinpath(outputs,"interp") #output files path
if ~ispath(writedir); mkpath(writedir); end

# ### 2D example

# Choose variable to process + specify file set name and precision

nam = "ETAN"
set = "state_2d_set1"
prc = Float32

# Get meta data for the chosen variable. Whithin `field` this defines:
#
# - a `NCvar` struct that sets up the subsequent `write` operation & incl. a `BinData` struct.
# - a `BinData` struct that contains the file names, precision, and array size.

    pth=input_file_paths(inputs)
    flddatadir = joinpath(pth["interp"],nam)
    fnames = joinpath.(Ref(flddatadir),filter(x -> occursin(".data",x), readdir(flddatadir)))
    diaginfo = read_available_diagnostics(nam, filename=pth["diaglist"])

(field,savename,readme)=prep_nctiles_interp(inputs,set,nam,prc);

# Create the NetCDF file and write data to it.

write(field,savename,README=readme);

# ### 3D example

(field,savename,readme)=prep_nctiles_interp(inputs,"WVELMASS","WVELMASS",Float32);
write(field,savename,README=readme);

# ## Tiled Data Examples
#
# This example reads in global variables defined over a collection of subdomain arrays ( _tiles_ ) using `MeshArrays.jl`, and writes them to a collection of `NetCDF` files ( _nctiles_ ) using `NCTiles.jl`.
#
# ### 2D & 3D examples

# +
#output folder name
writedir = joinpath(outputs,"tiled")
~ispath(writedir) ? mkpath(writedir) : nothing

#2D example
(flds,savename,readme)=prep_nctiles_native(inputs,"state_2d_set1","ETAN",Float32)
write(flds,savename,README=readme);

#3D example
(flds,savename,readme)=prep_nctiles_native(inputs,"state_3d_set1","THETA",Float32);
write(flds,savename,README=readme);
# -

# ### 3D vector example
#
# Here we process the three components of a vector field called `UVELMASS`, `VVELMASS` and `WVELMASS`. Note: on a `C-grid` these components are staggered in space.

# +
(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","UVELMASS",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","VVELMASS",Float32);
write(flds,savename,README=readme);

(flds,savename,readme)=prep_nctiles_native(inputs,"trsp_3d_set1","WVELMASS",Float32);
write(flds,savename,README=readme);
# -

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
"""
        topsis(decisionMat, weights, fns)

Apply TOPSIS (Technique for Order of Preference by Similarity to Ideal Solution) method
for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of function that are either minimize or maximize.

# Description
topsis() applies the TOPSIS method to rank n strategies subject to m criteria which are supposed to be either maximized or minimized.

# Output
- `::TopsisResult`: TopsisResult object that holds multiple outputs including scores and best index.

# Examples
```julia-repl
julia> df = DataFrame();
julia> df[:, :x] = Float64[9, 8, 7];
julia> df[:, :y] = Float64[7, 7, 8];
julia> df[:, :z] = Float64[6, 9, 6];
julia> df[:, :q] = Float64[7, 6, 6];

julia> w = Float64[4, 2, 6, 8];

julia> df
3×4 DataFrame
 Row │ x        y        z        q
     │ Float64  Float64  Float64  Float64
─────┼────────────────────────────────────
   1 │     9.0      7.0      6.0      7.0
   2 │     8.0      7.0      9.0      6.0
   3 │     7.0      8.0      6.0      6.0

julia> fns = makeminmax([maximum, maximum, maximum, maximum]);
julia> result = topsis(df, w, fns);

julia> result.bestIndex
2

julia> result.scores
3-element Array{Float64,1}:
 0.38768695492211824
 0.6503238218850163
 0.08347670030339041
```

# References
Hwang, C.L.; Yoon, K. (1981). Multiple Attribute Decision Making: Methods and Applications. New York: Springer-Verlag

Celikbilek Yakup, Cok Kriterli Karar Verme Yontemleri, Aciklamali ve Karsilastirmali
Saglik Bilimleri Uygulamalari ile. Editor: Muhlis Ozdemir, Nobel Kitabevi, Ankara, 2018

İşletmeciler, Mühendisler ve Yöneticiler için Operasyonel, Yönetsel ve Stratejik Problemlerin
Çözümünde Çok Kriterli Karar verme Yöntemleri, Editörler: Bahadır Fatih Yıldırım ve Emrah Önder,
Dora, 2. Basım, 2015, ISBN: 978-605-9929-44-8
"""
    function topsis(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::TopsisResult

    w = unitize(weights)
    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)

    weightednormalizedMat = w * normalizedMat

    # col_max = colmaxs(weightednormalizedMat)
    # col_min = colmins(weightednormalizedMat)
    col_max = apply_columns(fns, weightednormalizedMat)
    col_min = apply_columns(reverseminmax(fns), weightednormalizedMat)

    distances_plus  = zeros(Float64, nalternatives)
    distances_minus = zeros(Float64, nalternatives)

    scores = zeros(Float64, nalternatives)

    @inbounds for i in 1:nalternatives
        ithrow = weightednormalizedMat[i,:] |> Array{Float64,1}
		distances_plus[i]  = euclidean(col_max, ithrow)
		distances_minus[i] = euclidean(col_min, ithrow)
		scores[i] = distances_minus[i] / (distances_minus[i] + distances_plus[i])
    end

    best_index = sortperm(scores) |> last

    topsisresult = TopsisResult(
        decisionMat,
        w,
        normalizedMat,
        weightednormalizedMat,
        best_index,
        scores
    )

    return topsisresult
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
# discovering binary CUDA dependencies

using Pkg, Pkg.Artifacts
using Libdl

const libcublas = Ref{String}("cublas")
const libcusparse = Ref{String}("cusparse")
const libcusolver = Ref{String}("cusolver")
const libcufft = Ref{String}("cufft")
const libcurand = Ref{String}("curand")
const libcudnn = Ref{String}("cudnn")
const libcutensor = Ref{String}("cutensor")

## discovery

# NOTE: we don't use autogenerated JLLs, because we have multiple artifacts and need to
#       decide at run time (i.e. not via package dependencies) which one to use.
const cuda_artifacts = Dict(
    v"10.2" => ()->artifact"CUDA10.2",
    v"10.1" => ()->artifact"CUDA10.1",
    v"10.0" => ()->artifact"CUDA10.0",
    v"9.2"  => ()->artifact"CUDA9.2",
    v"9.0"  => ()->artifact"CUDA9.0",
)

# utilities to look up stuff in the artifact (at known locations, so not using CUDAapi)
function get_library(artifact, name)
    filename = if Sys.iswindows()
        "$name.dll"
    elseif Sys.isapple()
        "lib$name.dylib"
    else
        "lib$name.so"
    end
    joinpath(artifact, Sys.iswindows() ? "bin" : "lib", filename)
end

# try use CUDA from an artifact
function use_artifact_cuda()
    # select compatible artifacts
    if haskey(ENV, "JULIA_CUDA_VERSION")
        wanted_version = VersionNumber(ENV["JULIA_CUDA_VERSION"])
        filter!(((version,artifact),) -> version == wanted_version, cuda_artifacts)
    else
        driver_version = CUDAdrv.release()
        filter!(((version,artifact),) -> version <= driver_version, cuda_artifacts)
    end

    # download and install
    artifact = nothing
    release = nothing
    for version in sort(collect(keys(cuda_artifacts)); rev=true)
        try
            artifact = cuda_artifacts[version]()
            release = version
            break
        catch
        end
    end
    artifact == nothing && error("Could not find a compatible artifact.")

    # discover libraries
    for name in  ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CuArrays, Symbol("lib$name"))

        # on Windows, the library name is version dependent
        if Sys.iswindows()
            name *= release >= v"10.1" ? "64_$(release.major)" : "64_$(release.major)$(release.minor)"
        end

        handle[] = get_library(artifact, name)
        Libdl.dlopen(handle[])
    end

    @debug "Using CUDA $(release) from an artifact at $(artifact)"
    return release, [artifact]
end

# try to use CUDA from a local installation
function use_local_cuda()
    dirs = find_toolkit()

    tool = find_cuda_binary("nvdisasm")
    tool == nothing && error("Your CUDA installation does not provide the nvdisasm binary")
    version = parse_toolkit_version(tool)

    # discover libraries
    for name in  ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CuArrays, Symbol("lib$name"))

        path = find_cuda_library(name, dirs, [version])
        if path !== nothing
            handle[] = path
        end
    end

    @debug "Using local CUDA $(version) at $(join(dirs, ", "))"
    return version, dirs
end

const cudnn_artifacts = Dict(
    v"10.2" => ()->artifact"CUDNN+CUDA10.2",
    v"10.1" => ()->artifact"CUDNN+CUDA10.1",
    v"10.0" => ()->artifact"CUDNN+CUDA10.0",
    v"9.2"  => ()->artifact"CUDNN+CUDA9.2",
    v"9.0"  => ()->artifact"CUDNN+CUDA9.0",
)

function use_artifact_cudnn(cuda_release)
    artifact = try
        cudnn_artifacts[cuda_release]()
    catch ex
        @debug "Could not use CUDNN from artifacts" exception=(ex, catch_backtrace())
        return
    end

    libcudnn[] = get_library(artifact, Sys.iswindows() ? "cudnn64_7" : "cudnn")
    Libdl.dlopen(libcudnn[])
    @debug "Using CUDNN from an artifact at $(artifact)"
end

function use_local_cudnn(cuda_dirs)
    path = find_cuda_library("cudnn", cuda_dirs, [v"7"])
    if path !== nothing
        libcudnn[] = path
        @debug "Using local CUDNN at $(path)"
    end
end

const cutensor_artifacts = Dict(
    v"10.2" => ()->artifact"CUTENSOR+CUDA10.2",
    v"10.1" => ()->artifact"CUTENSOR+CUDA10.1",
)

function use_artifact_cutensor(cuda_release)
    artifact = try
        cutensor_artifacts[cuda_release]()
    catch ex
        @debug "Could not use CUTENSOR from artifacts" exception=(ex, catch_backtrace())
        return
    end

    libcutensor[] = get_library(artifact, "cutensor")
    Libdl.dlopen(libcutensor[])
    @debug "Using CUTENSOR from an artifact at $(artifact)"
end

function use_local_cutensor(cuda_dirs)
    path = find_cuda_library("cutensor", cuda_dirs, [v"1"])
    if path !== nothing
        libcutensor[] = path
        @debug "Using local CUTENSOR at $(path)"
    end
end

function __init_bindeps__(; silent=false, verbose=false)
    cuda = try
        parse(Bool, get(ENV, "JULIA_CUDA_USE_BINARYBUILDER", "true")) ||
            error("Use of CUDA artifacts not allowed by user")
        cuda_release, artifact = use_artifact_cuda()
        use_artifact_cudnn(cuda_release)
        use_artifact_cutensor(cuda_release)
        cuda_release
    catch ex
        @debug "Could not use CUDA from artifacts" exception=(ex, catch_backtrace())
        cuda_version, cuda_dirs = use_local_cuda()
        use_local_cudnn(cuda_dirs)
        use_local_cutensor(cuda_dirs)
        VersionNumber(cuda_version.major, cuda_version.minor)
    end

    # library dependencies
    CUBLAS.version()
    CUSPARSE.version()
    CUSOLVER.version()
    CUFFT.version()
    CURAND.version()
    # CUDNN and CUTENSOR are optional

    # library compatibility
    if has_cutensor()
        cutensor = CUTENSOR.version()
        if cutensor < v"1"
            silent || @warn("CuArrays.jl only supports CUTENSOR 1.0 or higher")
        end

        cutensor_cuda = CUTENSOR.cuda_version()
        if cutensor_cuda.major != cuda.major || cutensor_cuda.minor != cuda.minor
            silent || @warn("You are using CUTENSOR $cutensor for CUDA $cutensor_cuda with CUDA toolkit $cuda; these might be incompatible.")
        end
    end
    if has_cudnn()
        cudnn = CUDNN.version()
        if cudnn < v"7.6"
            silent || @warn("CuArrays.jl only supports CUDNN v7.6 or higher")
        end

        cudnn_cuda = CUDNN.cuda_version()
        if cudnn_cuda.major != cuda.major || cudnn_cuda.minor != cuda.minor
            silent || @warn("You are using CUDNN $cudnn for CUDA $cudnn_cuda with CUDA toolkit $cuda; these might be incompatible.")
        end
    end
end


text:
# identity operators: joDirac

export joDirac
"""
    joDirac(m;[DDT=...,][RDT=...,][name=...])

Dirac operator

# Signature

    joDirac(m::Integer;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joDirac")

# Arguments

# Arguments

- `m`: size
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    A=joDirac(3)

examples with DDT/RDT

    A=joDirac(3; DDT=Float32)
    A=joDirac(3; DDT=Float32,RDT=Float64)

"""
joDirac(m::Integer;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joDirac") =
    joMatrix{DDT,RDT}(name,m,m,
        v1->jo_convert(RDT,v1,false),
        v2->jo_convert(DDT,v2,false),
        v3->jo_convert(DDT,v3,false),
        v4->jo_convert(RDT,v4,false),
        v5->jo_convert(DDT,v5,false),
        v6->jo_convert(RDT,v6,false),
        v7->jo_convert(RDT,v7,false),
        v8->jo_convert(DDT,v8,false)
        )


text:
using LinearAlgebra
using ForneyLab
import ForneyLab: unsafeMean, unsafeCov, unsafePrecision
using NARMAX

function generate_data(ϕ; θ_scale=0.1, τ_true=1e3, degree=1, M1=1, M2=1, M3=1, fMin=0.8, fMax=1.0, fs=1.0, uStd=1., T=100, split_index=50, start_index=10, num_periods=num_periods, points_period=points_period)

    # Generate signal
    input, output = gensignalNARMAX(ϕ, θ_scale, τ_true; M1=M1, M2=M2, M3=M3, degree=degree, fMin=fMin, fMax=fMax, fs=fs, uStd=uStd, T=T, num_periods=num_periods, points_period=points_period)

    # Split data
    ix_trn, ix_val = split_data(T, split_index, start_index=start_index)

    return input, output, ix_trn, ix_val
end

function model_specification(ϕ; M1=M1, M2=M2, M3=M3, N=N)

    graph = FactorGraph()

    # Observed variables
    @RV u_kmin1; placeholder(u_kmin1, :u_kmin1, dims=(M1,))
    @RV y_kmin1; placeholder(y_kmin1, :y_kmin1, dims=(M2,))
    @RV e_kmin1; placeholder(e_kmin1, :e_kmin1, dims=(M3,))
    @RV u_k; placeholder(u_k, :u_k)

    # Time-invariant parameters
    @RV τ ~ Gamma(placeholder(:a_τ), placeholder(:b_τ))
    @RV θ ~ GaussianMeanPrecision(placeholder(:m_θ, dims=(N,)), placeholder(:w_θ, dims=(N,N)))

    # Likelihood
    @RV y_k ~ NAutoRegressiveMovingAverageX(θ, u_k, u_kmin1, y_kmin1, e_kmin1, τ, ϕ=ϕ)
    placeholder(y_k, :y_k)

    q = PosteriorFactorization(θ, τ, ids=[:θ :τ])
    algorithm = messagePassingAlgorithm([θ; τ], q, free_energy=true)
    return algorithmSourceCode(algorithm, free_energy=true)

end

function experiment_FEM(input_trn, output_trn, input_tst, output_tst, ϕ, priors; M1=1, M2=1, M3=1, N=3, num_iters=5, computeFE=false)

    # Maximum delay
    maxM = maximum([M1,M2,M3])

    # Zero-padding of signals
    input_trn = [zeros(maxM,); input_trn]
    input_tst = [zeros(maxM,); input_tst]
    output_trn = [zeros(maxM,); output_trn]
    output_tst = [zeros(maxM,); output_tst]

    # Signal splits
    T_trn = length(input_trn)
    T_tst = length(input_tst)

    "Inference execution"

    # Extract prior parameters
    θ_k = priors["θ"]
    τ_k = priors["τ"]

    # Initialize marginals
    marginals = Dict(:θ => ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=θ_k[1], w=θ_k[2]),
                     :τ => ProbabilityDistribution(Univariate, Gamma, a=τ_k[1], b=τ_k[2]))

    # Preallocate free energy array
    Fq = zeros(T_trn-maxM, num_iters)

    # Preallocate parameter arrays
    params_θ = (zeros(T_trn-maxM,N), zeros(T_trn-maxM,N,N))
    params_τ = (zeros(T_trn-maxM,1), zeros(T_trn-maxM,1))

    # Preallocate prediction arrays
    predictions = (zeros(T_trn,), zeros(T_trn,))
    errors = zeros(T_trn,)

    for (ii,k) in enumerate(maxM+1:T_trn)

        # Update history vectors
        u_kmin1 = input_trn[k-1:-1:k-M1]
        y_kmin1 = output_trn[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        ϕx = ϕ([input_trn[k]; u_kmin1; y_kmin1; e_kmin1])
        predictions[1][k] = θ_k[1]'*ϕx
        # predictions[2][k] = ϕx'*inv(θ_k[2])*ϕx + inv(τ_k[1]/τ_k[2])

        # Compute prediction error
        errors[k] = output_trn[k] - predictions[1][k]

        # Set data
        data = Dict(:y_k => output_trn[k],
                    :u_k => input_trn[k],
                    :u_kmin1 => u_kmin1,
                    :y_kmin1 => y_kmin1,
                    :e_kmin1 => e_kmin1,
                    :m_θ => θ_k[1],
                    :w_θ => θ_k[2],
                    :a_τ => τ_k[1],
                    :b_τ => τ_k[2])

        # Iterate updates
        for n = 1:num_iters

            # Compute free energy
            if computeFE
                Fq[ii,n] = freeEnergy(data,marginals)
            end

            # Update parameters
            stepθ!(data, marginals)
            stepτ!(data, marginals)

        end

        # Update params
        θ_k = (unsafeMean(marginals[:θ]), unsafePrecision(marginals[:θ]))
        τ_k = (marginals[:τ].params[:a], marginals[:τ].params[:b])

        # Store params
        params_θ[1][ii,:] = θ_k[1]
        params_θ[2][ii,:,:] = θ_k[2]
        params_τ[1][ii] = τ_k[1]
        params_τ[2][ii] = τ_k[2]

    end

    "1-step ahead prediction"

    # Prepare array
    predictions = (zeros(T_tst,), zeros(T_tst,))
    errors = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = output_tst[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        predictions[1][k] = θ_k[1]'*ϕx
        # predictions[2][k] = ϕx'*inv(θ_k[2])*ϕx + inv(τ_k[1] / τ_k[2])

        # Update error
        errors[k] = output_tst[k] - predictions[1][k]

    end

    "Simulation"

    # Prepare array
    simulations = (zeros(T_tst,), zeros(T_tst,))

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = simulations[1][k-1:-1:k-M2]
        e_kmin1 = zeros(M3,)

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        simulations[1][k] = θ_k[1]'*ϕx
        # simulations[2][k] = ϕx'*inv(θ_k[2])*ϕx + inv(τ_k[1] / τ_k[2])

    end

    "Evaluation"

    # Compute root mean square error
    RMS_prd = sqrt(mean((predictions[1][maxM+1:end] - output_tst[maxM+1:end]).^2))
    RMS_sim = sqrt(mean((simulations[1][maxM+1:end] - output_tst[maxM+1:end]).^2))

    if computeFE
        return RMS_sim, RMS_prd, Fq
    else
        return RMS_sim, RMS_prd
    end
end

function experiment_RLS(input_trn, output_trn, input_tst, output_tst, ϕ; M1=1, M2=1, M3=1, N=1, λ=1.00)

    # Maximum delay
    maxM = maximum([M1,M2,M3])

    # Zero-padding of signals
    input_trn = [zeros(maxM,); input_trn]
    input_tst = [zeros(maxM,); input_tst]
    output_trn = [zeros(maxM,); output_trn]
    output_tst = [zeros(maxM,); output_tst]

    # Signal splits
    T_trn = length(input_trn)
    T_tst = length(input_tst)

    # Parameters
    P = λ.*Matrix{Float64}(I,N,N)
    w_k = zeros(N,)

    # Preallocate prediction array
    predictions = zeros(T_trn,)
    errors = zeros(T_trn,)

    for (ii,k) in enumerate(maxM+1:T_trn)

        # Update history vectors
        u_kmin1 = input_trn[k-1:-1:k-M1]
        y_kmin1 = output_trn[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Update data vector
        ϕx = ϕ([input_trn[k]; u_kmin1; y_kmin1; e_kmin1])

        # Update weights
        α = output_trn[k] - w_k'*ϕx
        g = P*ϕx*inv(λ + ϕx'*P*ϕx)
        P = inv(λ)*P - g*ϕx'*inv(λ)*P
        w_k = w_k + α*g

        # Prediction
        predictions[k] = w_k'*ϕx
        errors[k] = output_trn[k] - predictions[k]

    end

    "1-step ahead prediction"

    # Prepare array
    predictions = zeros(T_tst,)
    errors = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = output_tst[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        predictions[k] = w_k'*ϕx

        # Update error
        errors[k] = output_tst[k] - predictions[k]

    end

    "Simulation"

    # Prepare array
    simulations = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = simulations[k-1:-1:k-M2]
        e_kmin1 = zeros(M3,)

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        simulations[k] = w_k'*ϕx

    end

    "Evaluation"

    # Compute root mean square error
    RMS_prd = sqrt(mean((predictions[maxM+1:end] - output_tst[maxM+1:end]).^2))
    RMS_sim = sqrt(mean((simulations[maxM+1:end] - output_tst[maxM+1:end]).^2))

    return RMS_sim, RMS_prd
end

function experiment_SYS(input_tst, output_tst, ϕ, θ; M1=1, M2=1, M3=1, N=1)

    # Maximum delay
    maxM = maximum([M1,M2,M3])

    # Zero-padding of signals
    input_tst = [zeros(maxM,); input_tst]
    output_tst = [zeros(maxM,); output_tst]

    # Length of validation signal
    T_tst = length(input_tst)

    "1-step ahead prediction"

    # Prepare array
    predictions = zeros(T_tst,)
    errors = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = output_tst[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        predictions[k] = θ'*ϕx

        # Update error
        errors[k] = output_tst[k] - predictions[k]

    end

    "Simulation"

    # Prepare array
    simulations = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = simulations[k-1:-1:k-M2]
        e_kmin1 = zeros(M3,)

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        simulations[k] = θ'*ϕx

    end

    "Evaluation"

    # Compute root mean square error
    RMS_prd = sqrt(mean((predictions[maxM+1:end] - output_tst[maxM+1:end]).^2))
    RMS_sim = sqrt(mean((simulations[maxM+1:end] - output_tst[maxM+1:end]).^2))

    return RMS_sim, RMS_prd
end


text:
export interpLocalToGlobal, interpGlobalToLocal
# Interpolate a vector x from the local mesh to the global mesh

"""
function jInv.ForwardShare.interpLocalToGlobal(mLocal,P,const)
function jInv.ForwardShare.interpGlobalToLocal(mGlobal,P,const)

convert a local (or global) model to a global (or local) mesh so that
    mGlobal = P * mLocal + const
    mLocal = P' * mGlobal + const

Examples:
    P = prepareMesh2Mesh(Msmall, Mbig, true)
    modelLoc =  interpGlobalToLocal(modelInv,P) # uncompresses the indices and does P'*modelInv

    P = prepareMesh2Mesh(Msmall, Mbig, true)
    modelInv =  interpGlobalToLocal(modelLoc,P) # uncompresses the indices and does P'*modelInv

Inputs:
    mLocal  - a model defined on the local (PDE) mesh
    mGlobal - a model defined on the global (inverse) mesh
    const   - a scalar constant shift (if not provided: const = 0.0)

Outputs:
    mGlobal - a model defined on the global (inverse) mesh
    mLocal  - a model defined on the local (PDE) mesh
"""
function interpLocalToGlobal(x::Vector{Float64}, P::AbstractFloat,y0::AbstractFloat=0.0)
	return P * x .+ y0
end

function interpGlobalToLocal(x::Vector{Float64}, P::AbstractFloat,y0::AbstractFloat=0.0)
	return P * x .+ y0
end

function interpGlobalToLocal(x,P)
    return P'*x
end

function interpGlobalToLocal(x,P,y0)
    return P'*x + y0
end

function interpLocalToGlobal(x,P)
    return P*x
end

function interpLocalToGlobal(x::Vector{Float64}, P::SparseMatrixCSC)

	if (eltype(P.nzval) == Int16) || (eltype(P.nzval) == Int8)
   		nzv = P.nzval
    	rv  = P.rowval
      	y   = zeros(P.m)
    	for col = 1 : P.n
	        xc = x[col]
	        @inbounds  for k = P.colptr[col] : (P.colptr[col+1]-1)
	            y[rv[k]] += xc / (1 << (-3 * nzv[k])) # = xc * 2 ^ (3 * nzv[k]); nzv[k] <= 0
        	end
		end
	else
		y = P * x
	end
	return y
end

# Interpolate a vector x from the global mesh to the local mesh
function interpGlobalToLocal(x::Vector{Float64}, P, y0::Vector{Float64})
	return y0 + interpGlobalToLocal(x, P)
end

# Interpolate a vector x from the global mesh to the local mesh
function interpGlobalToLocal(x::Vector{Float64}, P::SparseMatrixCSC)
	n = div(length(x), size(P,1))
	(length(x) == n * size(P,1)) || error("Incompatible input sizes")
	x = reshape(x,(size(P,1),n))
	if (eltype(P.nzval) == Int16) || (eltype(P.nzval) == Int8)
		nzv = P.nzval
		rv  = P.rowval
		y   = zeros(P.n,n)
		@inbounds begin
			for k = 1:n
				for i = 1 : P.n
					tmp = 0.0
					for j = P.colptr[i] : (P.colptr[i+1]-1)
						tmp += x[rv[j],k] / (1 << (-3 * nzv[j])) # = x[rv[j],k] * 2 ^ (3 * nzv[j]); nzv[j] <= 0
					end
					y[i,k] = tmp
				end
			end
		end
	else
		y = P' * x
	end
	return vec(y)
end


text:
using MiscBenchmarks
using Test

@testset "MiscBenchmarks.jl" begin
    # Write your own tests here.
end