text:
@blasfun clblasSrot(N::Csize_t,
                         X::cl.CL_mem, offx::Csize_t, incx::Cint,
                         Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                         C::cl.CL_float, S::cl.CL_float,
                         n_queues::cl.CL_uint,
                         queues::Ptr{cl.CL_command_queue},
                         n_events_in_wait_list::cl.CL_uint,
                         event_wait_list::Ptr{cl.CL_event},
                         events::Ptr{cl.CL_event})

@blasfun2 clblasSrot(N::Csize_t,
                          X::cl.CL_mem, offx::Csize_t, incx::Cint,
                          Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                          C::cl.CL_float, S::cl.CL_float)

@blasfun clblasDrot(N::Csize_t,
                         X::cl.CL_mem, offx::Csize_t, incx::Cint,
                         Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                         C::cl.CL_double, S::cl.CL_double,
                         n_queues::cl.CL_uint,
                         queues::Ptr{cl.CL_command_queue},
                         n_events_in_wait_list::cl.CL_uint,
                         event_wait_list::Ptr{cl.CL_event},
                         events::Ptr{cl.CL_event})

@blasfun2 clblasDrot(N::Csize_t,
                          X::cl.CL_mem, offx::Csize_t, incx::Cint,
                          Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                          C::cl.CL_double, S::cl.CL_double)

@blasfun clblasCsrot(N::Csize_t,
                          X::cl.CL_mem, offx::Csize_t, incx::Cint,
                          Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                          C::cl.CL_float, S::cl.CL_float,
                          n_queues::cl.CL_uint,
                          queues::Ptr{cl.CL_command_queue},
                          n_events_in_wait_list::cl.CL_uint,
                          event_wait_list::Ptr{cl.CL_event},
                          events::Ptr{cl.CL_event})

@blasfun2 clblasCrot(N::Csize_t,
                          X::cl.CL_mem, offx::Csize_t, incx::Cint,
                          Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                          C::cl.CL_float, S::cl.CL_float)

@blasfun clblasZdrot(N::Csize_t,
                          X::cl.CL_mem, offx::Csize_t, incx::Cint,
                          Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                          C::cl.CL_double, S::cl.CL_double,
                          n_queues::cl.CL_uint,
                          queues::Ptr{cl.CL_command_queue},
                          n_events_in_wait_list::cl.CL_uint,
                          event_wait_list::Ptr{cl.CL_event},
                          events::Ptr{cl.CL_event})

@blasfun2 clblasZdrot(N::Csize_t,
                           X::cl.CL_mem, offx::Csize_t, incx::Cint,
                           Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                           C::cl.CL_double, S::cl.CL_double)


text:
using Test, Random, LinearAlgebra, StatsBase, StatsPlots, Distributions, Traceur

include("../src/types.jl")
include("../src/functions.jl")
include("../src/params.jl")
using .Types, .Func, .Params

function runmodel()

    # Set up types
    market = Types.MarketIndex(
        zeros(Params.bigt))
    stocks = Types.Equity(
        zeros(Params.bigm, Params.bigt),
        zeros(Params.bigm),
        zeros(Params.bigm),
        zeros(Params.bigm),
        zeros(Params.bigm, Params.bigt))
    investors = Types.RetailInvestor(
        zeros(Params.bign, Params.bigk + 1),
        zeros(Params.bign),
        zeros(Params.bign))
    funds = Func.Types.EquityFund(
        zeros(Params.bigk, Params.bigm),
        zeros(Params.bigk, Params.bign),
        zeros(Params.bigk, Params.bigt))

    Func.initialise(market, stocks, investors, funds)
    #Func.boundstest(market, stocks, investors, funds)
    # TODO: Think about how I should have tested the initialisation
    # TODO: Set test=true flag for testing bounds and other such things
    #  or is that bad practice?

    # TODO: Create logging set-up that records parameters, plots and outputs
    # TODO: Think about how I should have tested the running cycle
    # TODO: Write functions that collect activity of a) selling out and b) buying in
    # TODO: Tests for cases with empty reviewers, divestments
    # QUESTION: Need 'if reviewers not empty'/'divestments not empty' branches?
    # Looping over time after the end of the  perfwindow[end]+1 initialisation phase

    Func.modelrun(market, stocks, investors, funds)
    Func.boundstest(market, stocks, investors, funds)
    #Func.plot_stylisedfacts(market.value, stocks.value, stocks.volume)
end
runmodel()

ext: jl


text:
@testset "Testing CHSF Descriptor for dc Si" begin
@info("Testing CHSF Descriptor for dc Si.")
using DescriptorZoo, JuLIP, Test

at = bulk(:Si, cubic=true)
desc = chsf(at, 6.5, n=2, l=2)
chsf_ref = [10.3698237,1.4503467,-8.2118063,51.6882200,-53.0113716,69.8233316] #n=2,l=2 case
chsf_now = vcat(desc[1,:]...)
println(@test chsf_now  ≈  chsf_ref)

end


text:
using StatsBase: sample

mutable struct MCTSPlayer
  env :: GameEnv
  network
  num_readouts::Int
  two_player_mode::Bool
  τ_threshold::Int
  qs::Vector{Float32}
  searches_π::Vector{Vector{Float32}}
  result::Int
  result_string
  root
  resign_threshold
  position

  function MCTSPlayer(env::T, network; num_readouts = 800, two_player_mode = false,
               resign_threshold = -0.9) where T <: GameEnv
    τ_threshold = two_player_mode ? -1 : (env.N * env.N ÷ 12) ÷ 2 * 2
    new(env, network, num_readouts, two_player_mode, τ_threshold,
        Array{Float32, 1}(), Array{Array{Float32, 1}, 1}(), 0, "",
        nothing, resign_threshold, nothing)
  end
end

function play_move!(mcts_player::MCTSPlayer, c)
  #=
  Notable side effects:
    - finalizes the probability distribution according to
    this roots visit counts into the class' running tally, `searches_π`
    - Makes the node associated with this move the root, for future
      `inject_noise` calls.
  =#
  if !mcts_player.two_player_mode
    push!(mcts_player.searches_π, children_as_π(mcts_player.root,
    mcts_player.root.position.n <= mcts_player.τ_threshold))
  end
  push!(mcts_player.qs, Q(mcts_player.root))  # Save our resulting Q.
  try
    mcts_player.root = maybe_add_child!(mcts_player.root, to_flat(c, mcts_player.root.position))
  catch IllegalMove
    println("Illegal move")
    if !mcts_player.two_player_mode pop!(mcts_player.searches_π) end
    pop!(mcts_player.qs)
    return false
  end
  mcts_player.position = mcts_player.root.position  # for showboard
  mcts_player.root.parent.children = Dict()
  return true
end

function pick_move(mcts_player::MCTSPlayer)
  #= Picks a move to play, based on MCTS readout statistics.

  Highest N is most robust indicator. In the early stage of the game, pick
  a move weighted by visit count; later on, pick the absolute max. =#

  if mcts_player.root.position.n >= mcts_player.τ_threshold
    #max_val = maximum(mcts_player.root.child_N)
    #possible_moves = find(x -> x == max_val, mcts_player.root.child_N)
    #fcoord = sample(possible_moves)
    fcoord = findmax(mcts_player.root.child_N)[2]
  else
    cdf = cumsum(mcts_player.root.child_N)
    cdf /= cdf[end - 1]  # Prevents passing via softpick.
    selection = rand()
    fcoord = searchsortedfirst(cdf, selection)
    @assert mcts_player.root.child_N[fcoord] != 0
  end
  return from_flat(fcoord, mcts_player.root.position)
end

function tree_search!(mcts_player::MCTSPlayer, parallel_readouts = 8)
  leaves = Vector{MCTSNode}()
  failsafe = 0
  while length(leaves) < parallel_readouts && failsafe < 2parallel_readouts
    failsafe += 1
    leaf = select_leaf(mcts_player.root)
    # if game is over, override the value estimate with the true score
    if is_done(leaf)
      value = result(leaf.position)
      backup_value!(leaf, value, mcts_player.root)
      continue
    end
    add_virtual_loss!(leaf, mcts_player.root)
    push!(leaves, leaf)
  end
  if !isempty(leaves)
    move_probs, values = mcts_player.network([leaf.position for leaf in leaves])
    move_probs, values = move_probs.tracker.data, values.tracker.data
    move_probs = [move_probs[:, i] for i = 1:size(move_probs, 2)]
    for (leaf, move_prob, value) in zip(leaves, move_probs, values)
      revert_virtual_loss!(leaf, mcts_player.root)
      incorporate_results!(leaf, move_prob, value, mcts_player.root)
    end
  end
  return leaves
end

function set_result!(mcts_player::MCTSPlayer, winner, was_resign)
  mcts_player.result = winner
  if was_resign
    str = winner == BLACK ? "B+R" : "W+R"
  else
    str = result_string(mcts_player.root.position)
  end
  mcts_player.result_string = str
end

function initialize_game!(mcts_player::MCTSPlayer, pos = nothing)
  if pos == nothing
    pos = Position(mcts_player.env)
  end
  mcts_player.root = MCTSNode(pos)
  mcts_player.result = 0
  mcts_player.searches_π = Vector{Vector{Float32}}()
  mcts_player.qs = Vector{Float32}()
end

is_done(mcts_player::MCTSPlayer) = mcts_player.result != 0 || is_done(mcts_player.root)

# Returns true if the player resigned.  No further moves should be played

should_resign(mcts_player::MCTSPlayer) = Q_perspective(mcts_player.root) < mcts_player.resign_threshold

function extract_data(mcts_player::MCTSPlayer)
  @assert length(mcts_player.searches_π) == mcts_player.root.position.n
  @assert mcts_player.result != 0
  positions = Vector{Position}()
  πs = deepcopy.(mcts_player.searches_π)
  results = Vector{Int}()

  pwcs = replay_position(mcts_player.root.position, mcts_player.result)
  for pwc in pwcs
    push!(positions, pwc.position)
    push!(results, pwc.result)
  end
  return positions, πs, results
end

get_position(mcts_player::MCTSPlayer) = mcts_player.root != nothing ?
                                            mcts_player.root.position : nothing

function suggest_move(mcts_player::MCTSPlayer)
  current_readouts = N(mcts_player.root)
  while N(mcts_player.root) < current_readouts + mcts_player.num_readouts
    tree_search!(mcts_player)
  end

  pick_move(mcts_player)
end


text:
using IntcodeVM, OffsetArrays

const code = IntcodeVM.load_intcode("input.txt")

function get_puzzle()
    o = IntcodeVM.run(code)[2]
    pop!(o) # remove doubled newline at the end
    s = foldr(*, Char.(o))
    println(s)

    width = findfirst(==(Int('\n')), o)
    m = reshape(o, width, :)' # transpose because Julia is column-major
    v = view(m, :, 1:width-1) # take a view without newlines

    # zero-index array indices
    z = OffsetArray(v, 0:size(v)[1]-1, 0:size(v)[2]-1)
    return z
end

function get_example()
    v = Int.(first.(hcat((split.(readlines("example.txt"), ""))...)))'
    return OffsetArray(v, 0:size(v)[1]-1, 0:size(v)[2]-1)
end

function find_intersections(a, value)
    alignment_parameters = Array{Int,1}(undef, 0)

    for i=1:size(a)[1]-2
        for j=1:size(a)[2]-2
            if a[i,j] == a[i,j+1] == a[i,j-1] == a[i-1,j] == a[i+1,j] == value
                push!(alignment_parameters, i * j)
            end
        end
    end

    return alignment_parameters
end

maze = get_puzzle()
print("Day 17 Part 1: ")
println(sum(find_intersections(maze, Int('#'))))


text:
# This file contains a number of additions to BlockArrays.jl. These are completely
# independent of Stheno.jl, and will (hopefully) move over to BlockArrays.jl at some point.

"""
    blocksizes(X::AbstractBlockVecOrMat)
    blocksizes(X::AbstractArray, d::Int)

Get a vector containing the block sizes over the `d`th dimension of `X`. If `d` not
specified, returns a `Tuple` over blocksizes over each dimension.
"""
blocksizes(X::AbstractBlockArray, d::Int) = blocklengths(axes(X, d))
ChainRulesCore.@non_differentiable blocksizes(::Any, ::Any)

cumulsizes(X::AbstractBlockArray, d::Int) = vcat(1, cumsum(blocksizes(X, d)) .+ 1)
ChainRulesCore.@non_differentiable cumulsizes(::Any, ::Any)

nblocks(X::BlockArray, d::Int) = size(X.blocks, d)

"""
    are_conformal(A::BlockVecOrMat, B::BlockVecOrMat)

Test whether two block matrices (or vectors) are conformal. This criterion is stricter than
that for general matrices / vectors as we additionally require that each block be conformal
with block of the other matrix with which it will be multiplied. This ensures that the
result is itself straightforwardly representable as `BlockVecOrMat`.
"""
are_conformal(A::BlockVecOrMat, B::BlockVecOrMat) = blocksizes(A, 2) == blocksizes(B, 1)

ZygoteRules.@adjoint function BlockArrays.mortar(_blocks::AbstractArray)
    function mortar_pullback(Δ::NamedTuple{(:blocks, :axes)})
        return (Δ.blocks, )
    end
    function mortar_pullback(Δ::BlockArray)
        return mortar_pullback((blocks = Δ.blocks, axes=nothing))
    end
    return BlockArrays.mortar(_blocks), mortar_pullback
end

ZygoteRules.@adjoint function BlockArrays.Array(X::BlockArray)
    function Array_pullback(Δ::Array)
        ΔX = (blocks=BlockArray(Δ, axes(X)).blocks, axes=nothing)
        return (ΔX,)
    end
    return Array(X), Array_pullback
end


text:
#
# Qubit site type
#

# Define Qubit space in terms of
# S=1/2 space, but use different
# defaults for QN names

"""
    space(::SiteType"Qubit";
          conserve_qns = false,
          conserve_parity = conserve_qns,
          conserve_number = false,
          qnname_parity = "Parity",
          qnname_number = "Number")

Create the Hilbert space for a site of type "Qubit".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
  ::SiteType"Qubit";
  conserve_qns=false,
  conserve_parity=conserve_qns,
  conserve_number=false,
  qnname_parity="Parity",
  qnname_number="Number",
)
  if conserve_number && conserve_parity
    return [
      QN((qnname_number, 0), (qnname_parity, 0, 2)) => 1,
      QN((qnname_number, 1), (qnname_parity, 1, 2)) => 1,
    ]
  elseif conserve_number
    return [QN(qnname_number, 0) => 1, QN(qnname_number, 1) => 1]
  elseif conserve_parity
    return [QN(qnname_parity, 0, 2) => 1, QN(qnname_parity, 1, 2) => 1]
  end
  return 2
end

ITensors.val(::ValName"0", st::SiteType"Qubit") = 1
ITensors.val(::ValName"1", st::SiteType"Qubit") = 2

ITensors.state(::StateName"0", ::SiteType"Qubit") = [1.0, 0.0]
ITensors.state(::StateName"1", ::SiteType"Qubit") = [0.0, 1.0]

# Use S=1/2 definition of any operators
# called using Qubit SiteType
ITensors.op(o::OpName, ::SiteType"Qubit"; kwargs...) = op(o, SiteType("S=1/2"); kwargs...)

#
# 1-Qubit gates
#

ITensors.op(::OpName"σx", t::SiteType"S=1/2") = op("X", t)

ITensors.op(::OpName"σ1", t::SiteType"S=1/2") = op("X", t)

function ITensors.op(::OpName"√NOT", ::SiteType"S=1/2")
  return [
    (1 + im)/2 (1 - im)/2
    (1 - im)/2 (1 + im)/2
  ]
end

ITensors.op(::OpName"√X", t::SiteType"S=1/2") = op("√NOT", t)

ITensors.op(::OpName"σy", t::SiteType"S=1/2") = op("Y", t)

ITensors.op(::OpName"σ2", t::SiteType"S=1/2") = op("Y", t)

ITensors.op(::OpName"iσy", t::SiteType"S=1/2") = op("iY", t)

ITensors.op(::OpName"iσ2", t::SiteType"S=1/2") = op("iY", t)

ITensors.op(::OpName"σz", t::SiteType"S=1/2") = op("Z", t)

ITensors.op(::OpName"σ3", t::SiteType"S=1/2") = op("Z", t)

ITensors.op(::OpName"H", ::SiteType"S=1/2") = [
  1/sqrt(2) 1/sqrt(2)
  1/sqrt(2) -1/sqrt(2)
]

# Rϕ with ϕ = π/2
ITensors.op(::OpName"Phase", ::SiteType"S=1/2") = [
  1 0
  0 im
]

ITensors.op(::OpName"P", t::SiteType"S=1/2") = op("Phase", t)

ITensors.op(::OpName"S", t::SiteType"S=1/2") = op("Phase", t)

## Rϕ with ϕ = π/4
ITensors.op(::OpName"π/8", ::SiteType"S=1/2") = [
  1 0
  0 1 / sqrt(2)+im / sqrt(2)
]

ITensors.op(::OpName"T", t::SiteType"S=1/2") = op("π/8", t)

# Rotation around X-axis
function ITensors.op(::OpName"Rx", ::SiteType"S=1/2"; θ::Number)
  return [
    cos(θ / 2) -im*sin(θ / 2)
    -im*sin(θ / 2) cos(θ / 2)
  ]
end

# Rotation around Y-axis
function ITensors.op(::OpName"Ry", ::SiteType"S=1/2"; θ::Number)
  return [
    cos(θ / 2) -sin(θ / 2)
    sin(θ / 2) cos(θ / 2)
  ]
end

# Rotation around Z-axis
ITensors.op(::OpName"Rz", ::SiteType"S=1/2"; ϕ::Number) = [
  1 0
  0 exp(im * ϕ)
]

# Rotation around generic axis n̂
function ITensors.op(::OpName"Rn", ::SiteType"S=1/2"; θ::Real, ϕ::Real, λ::Real)
  return [
    cos(θ / 2) -exp(im * λ)*sin(θ / 2)
    exp(im * ϕ)*sin(θ / 2) exp(im * (ϕ + λ))*cos(θ / 2)
  ]
end

ITensors.op(::OpName"Rn̂", t::SiteType"S=1/2"; kwargs...) = ITensors.op("Rn", t; kwargs...)

#
# 2-Qubit gates
#

ITensors.op(::OpName"CNOT", ::SiteType"S=1/2") = [
  1 0 0 0
  0 1 0 0
  0 0 0 1
  0 0 1 0
]

ITensors.op(::OpName"CX", t::SiteType"S=1/2") = op("CNOT", t)

ITensors.op(::OpName"CY", ::SiteType"S=1/2") = [
  1 0 0 0
  0 1 0 0
  0 0 0 -im
  0 0 im 0
]

ITensors.op(::OpName"CZ", ::SiteType"S=1/2") = [
  1 0 0 0
  0 1 0 0
  0 0 1 0
  0 0 0 -1
]

ITensors.op(::OpName"SWAP", ::SiteType"S=1/2") = [
  1 0 0 0
  0 0 1 0
  0 1 0 0
  0 0 0 1
]

function ITensors.op(::OpName"√SWAP", ::SiteType"S=1/2")
  return [
    1 0 0 0
    0 (1 + im)/2 (1 - im)/2 0
    0 (1 - im)/2 (1 + im)/2 0
    0 0 0 1
  ]
end

# Ising (XX) coupling gate
function ITensors.op(::OpName"XX", ::SiteType"S=1/2"; ϕ::Number)
  return [
    cos(ϕ) 0 0 -im*sin(ϕ)
    0 cos(ϕ) -im*sin(ϕ) 0
    0 -im*sin(ϕ) cos(ϕ) 0
    -im*sin(ϕ) 0 0 cos(ϕ)
  ]
end

# Ising (YY) coupling gate
function ITensors.op(::OpName"YY", ::SiteType"S=1/2"; ϕ::Number)
  return [
    cos(ϕ) 0 0 im*sin(ϕ)
    0 cos(ϕ) -im*sin(ϕ) 0
    0 -im*sin(ϕ) cos(ϕ) 0
    im*sin(ϕ) 0 0 cos(ϕ)
  ]
end

# Ising (ZZ) coupling gate
function ITensors.op(::OpName"ZZ", ::SiteType"S=1/2"; ϕ::Number)
  return [
    exp(-im * ϕ) 0 0 0
    0 exp(im * ϕ) 0 0
    0 0 exp(im * ϕ) 0
    0 0 0 exp(-im * ϕ)
  ]
end

#
# 3-Qubit gates
#

function ITensors.op(::OpName"Toffoli", ::SiteType"S=1/2")
  return [
    1 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0
    0 0 1 0 0 0 0 0
    0 0 0 1 0 0 0 0
    0 0 0 0 1 0 0 0
    0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 1
    0 0 0 0 0 0 1 0
  ]
end

ITensors.op(::OpName"CCNOT", t::SiteType"S=1/2") = op("Toffoli", t)

ITensors.op(::OpName"CCX", t::SiteType"S=1/2") = op("Toffoli", t)

ITensors.op(::OpName"TOFF", t::SiteType"S=1/2") = op("Toffoli", t)

function ITensors.op(::OpName"Fredkin", ::SiteType"S=1/2")
  return [
    1 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0
    0 0 1 0 0 0 0 0
    0 0 0 1 0 0 0 0
    0 0 0 0 1 0 0 0
    0 0 0 0 0 0 1 0
    0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 1
  ]
end

ITensors.op(::OpName"CSWAP", t::SiteType"S=1/2") = op("Fredkin", t)

ITensors.op(::OpName"CS", t::SiteType"S=1/2") = op("Fredkin", t)

#
# 4-Qubit gates
#

function ITensors.op(::OpName"CCCNOT", ::SiteType"S=1/2")
  return [
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
  ]
end


text:
using Revise
using NNFEM
using PyPlot
using MAT
using ADCME

NT = 100
T = 0.5
Δt = T/NT

m, n =  15, 15
h = 1/m

ν = 0.35
E = 2.0
H = zeros(3,3)
H[1,1] = E*(1. -ν)/((1+ν)*(1. -2. *ν));
H[1,2] = H[1,1]*ν/(1-ν);
H[2,1] = H[1,2];
H[2,2] = H[1,1];
H[3,3] = H[1,1]*0.5*(1. -2. *ν)/(1. -ν);

# Create a very simple mesh
elements = []
prop = Dict("name"=> "PlaneStrain", "rho"=> 1.0, "E"=> 2.0, "nu"=> 0.35)
coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 3
        nodes = [
            (i-1)*h (j-1)*h
            i*h (j-1)*h
            i*h j*h
            (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)

for i = 2:m
    idx = n*(m+1)+i
    EBC[idx,:] .= -2
    EBC[i,:] .= -2
end

for j = 1:n+1
    idx = (j-1)*(m+1)+1
    EBC[idx,:] .= -2
    idx = (j-1)*(m+1)+m+1
    EBC[idx,:] .= -2
end

ndims = 2
domain = Domain(coords, elements, ndims, EBC, g, FBC, f)

EBC_idx = findall(EBC[:,1] .== -2)
EBC_X = domain.nodes[EBC_idx,:]
function gt(t)
    x, y = EBC_X[:,1], EBC_X[:,2]
    u1 = @. (x^2+y^2)*exp(-t) * 0.1
    u2 = @. (x^2-y^2)*exp(-t) * 0.1
    return [u1;u2], -[u1;u2], [u1;u2]
end

FBC_idx = findall(FBC[:,1] .== -2)
FBC_X = domain.nodes[FBC_idx,:]
# function ft(t)
#     x, y = FBC_X[:,1], FBC_X[:,2]
#     T1 = @. -1.48148148148148*x*exp(-t) - 1.48148148148148*y*exp(-t)
#     T2 = @. -3.45679012345679*x*exp(-t) - 6.41975308641975*y*exp(-t)
#     return [T1;T2] * h * 0.1
# end

function bt(x, y, t)
    f1 = @. 0.1*(x^2 + y^2)*exp(-t) - 0.790123456790123*exp(-t)
    f2 = @. 0.1*(x^2 - y^2)*exp(-t) + 0.493827160493827*exp(-t)
    [f1 f2]
end

c = domain.nodes[domain.dof_to_eq]
nn = div(length(c), 2)
x = c[1:nn]
y = c[nn+1:end]

Dstate =  @. [x^2+y^2;x^2-y^2]* 0.1
state = @. [x^2+y^2;x^2-y^2]* 0.1
velo =  @. -[x^2+y^2;x^2-y^2]* 0.1
acce =  @. [x^2+y^2;x^2-y^2]* 0.1
globdat = GlobalData(state, Dstate, velo, acce, domain.neqs, gt, nothing, bt)

assembleMassMatrix!(globdat, domain)
# updateTimeDependentEssentialBoundaryCondition!(domain,globdat)
# updateStates!(domain, globdat)

# for i = 1:NT
#     @info i
#     # global globdat, domain = GeneralizedAlphaSolverStep(globdat, domain, Δt)
#     global globdat, domain = ExplicitSolverStep(globdat, domain, Δt)
# end

# d_ = hcat(domain.history["state"]...)'
# i = 3
# j = 4
# plot(d_[:,(j-1)*(m+1)+i])
# x0 = (i-1)*h
# y0 = (j-1)*h
# ts = LinRange(0, 1, NT+1)
# plot((@. (x0^2+y0^2)*exp(-ts))*0.1,"--")

# error()
Hs = zeros(domain.neles*length(domain.elements[1].weights), 3, 3)
for i = 1:size(Hs,1)
    Hs[i,:,:] = H
end

# ts = GeneralizedAlphaSolverTime(Δt, NT)
ts = ExplicitSolverTime(Δt, NT)

ubd = zeros(NT,sum(domain.EBC[:].!=0))
abd = zeros(NT,sum(domain.EBC[:].!=0))
Fext = zeros(NT,domain.neqs)

idx = findall(domain.EBC[:,1].!=0)
x = domain.nodes[idx,1]
y = domain.nodes[idx,2]
for i = 1:NT
    t = ts[i]
    ubd[i,:] = [
        @. (x^2+y^2)*0.1*exp(-t)
        @. (x^2-y^2)*0.1*exp(-t)
    ]
    abd[i,:] = [
        @. (x^2+y^2)*0.1*exp(-t)
        @. (x^2-y^2)*0.1*exp(-t)
    ]
    globdat.time = t
    Fext[i,:] = getBodyForce(domain, globdat)
end

x = domain.nodes[:,1]
y = domain.nodes[:,2]
a0 = @. [x^2+y^2;x^2-y^2]* 0.1
v0 =  @. -[x^2+y^2;x^2-y^2]* 0.1
d0 =  @. [x^2+y^2;x^2-y^2]* 0.1

H = constant(Hs)

function linear_law(ε)
    @info H, ε
    batch_matmul(H, ε)
end

function nn_law(ε)
    ae(ε, [20,20,20,3])
end
# d, v, a= GeneralizedAlphaSolver(globdat, domain, d0, v0, a0, Δt, NT, H, Fext, ubd, abd)
# d, v, a= ExplicitSolver(globdat, domain, d0, v0, a0, Δt, NT, nn_law, Fext, ubd, abd)
d, v, a= ExplicitSolver(globdat, domain, d0, v0, a0, Δt, NT, nn_law, Fext, ubd, abd; strain_type="finite")

# # ...1... use `linear_law` for 8th argument in the last line to generate data
# sess = Session(); init(sess)
# d_, v_, a_ = run(sess, [d,v,a])
# matwrite("data.mat", Dict("d"=>d_, "v"=>v_, "a"=>a_))
# using Random; Random.seed!(233)
# k = 0
# for i = 1:5
#     i = rand(2:m)
#     j = rand(2:n)
#     x0 = (i-1)*h
#     y0 = (j-1)*h
#     global k += 1
#     plot(d_[:,(j-1)*(m+1)+i], "o", markersize=2,color = "C$k")
#     plot((@. (x0^2+y0^2)*exp(-ts))*0.1,"-", color="C$k")

#     global k += 1
#     plot(d_[:,(j-1)*(m+1)+i+(m+1)*(n+1)], "o", markersize=2,color = "C$i")
#     plot((@. (x0^2-y0^2)*exp(-ts))*0.1,"-", color="C$i")
# end
# error()

# ...2...comment the last block and replace `linear_law` with `nn_law` to train a neural network
dobs = matread("data.mat")["d"]
loss = sum((d - dobs)^2)
BFGS!(sess, loss)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct GraphDef <: ProtoType
    node::Base.Vector{NodeDef}
    versions::VersionDef
    version::Int32
    library::FunctionDefLibrary
    GraphDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GraphDef
const __fnum_GraphDef = Int[1,4,3,2]
meta(t::Type{GraphDef}) = meta(t, ProtoBuf.DEF_REQ, __fnum_GraphDef, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export GraphDef


text:
@testset "Utils" begin
  @testset "Select" begin
    @test true
  end;
end;

max_stars_repo_name: yollotltamayo/Covid

max_issues_repo_name: yollotltamayo/Covid

max_forks_repo_name: yollotltamayo/Covid


text:
function func(t,x)
    return exp(2*t)*x^2
end

function newtonFunc(yn1, tn1, yn)
    return yn1-yn-h*func(tn1, yn)
    # x_{n+1} - x_n -hfunc(t_{n+1}, x_{n+1}) =0
end

function fixedpointFunc(tn, yn)
    return yn+h*func(tn, yn)
    # x_{n+1} - x_n -hfunc(t_{n+1}, x_{n+1}) =0
end

function newton(a, b, index, previous)
    steps = 0
    mid = (a + b)/2
    while b-a>=2*tolerance
        #println("Step", " ", steps)
        mid = (a + b)/2
        if (abs(newtonFunc(mid, time[index],previous ))<tolerance)
            return mid
        elseif newtonFunc(mid, time[index], previous)*newtonFunc(b, time[index], previous) <0
            a = mid
        elseif newtonFunc(a, time[index], previous)*newtonFunc(mid, time[index], previous) <0
            b = mid
        end
        steps = steps+1
        if steps>100
            if index==N
                println("newton steps: ", steps)
            end
            return mid
        end
    if index==N
        println("newton steps: ", steps)
    end
    end
    if index==N
        println("newton steps: ", steps)
    end
    return mid
end

function fixedpoint(x0, index)

    ## x = hf(t,x) -x0

    x1 = fixedpointFunc(time[index], x0)
    x2 = fixedpointFunc(time[index], x1)
    k = 0
    shouldReturn=0
    while(abs(x2-x1)>=tolerance && abs(fixedpointFunc(time[index], x2)-x2)>=tolerance && k<=maxSteps-1)
        x2temp=x2;
        x2=fixedpointFunc(time[index],x1);
        x1=x2temp;
        k=k+1;
        if (abs(fixedpointFunc(time[index],x2)-x2)<tolerance)
            if index==N
                println("fp steps: ", k)
            end
            return x2
        end
    end
    if index==N
        println("fp steps: ", k)
    end
    return x2
end

function exactsolution(index)
    return 1/(10.5-exp(2*time[index]))
end

function main()
    newtonOut  = zeros(N)
    fixedpointOut = zeros(N)
    exactOut = zeros(N)
    newtonOut[1] = 0.1
    fixedpointOut[1] = 0.1
    exactOut[1] = exactsolution(1)
    ### Implicit backward euler
    ## Newton

    for i=2:N
        #println("Running", " ", i)
        newtonOut[i] = newton(newtonOut[i-1], exactsolution(i)+10*tolerance, i, newtonOut[i-1])
        fixedpointOut[i] = fixedpoint(fixedpointOut[i-1], i)
        exactOut[i] = exactsolution(i)
        #println(newtonOut[i], " ", fixedpointOut[i], " ", exactOut[i])

    end
    #println(2*newtonOut[N], " ", fixedpointOut[N], " ", exactOut[N])
    #println(exactOut[N]-2*newtonOut[N], " ", exactOut[N]-fixedpointOut[N])#, " ", exactOut[N])
end

x = [4,8,16,32,64, 128]
N=0
h=0
tolerance = 10e-6
maxSteps = 100
for i=1:6
    global N = x[i]
    global h = 1/N
    global time = linrange(0,1,N)

    main()
end

max_stars_repo_name: NeveIsa/hatex

max_issues_repo_name: NeveIsa/hatex

max_forks_repo_name: saketkc/hatex


text:
module SimradEK60TestData

export EK60_DATA, EK60_SAMPLE, ECS_SAMPLE

const EK60_DATA = joinpath(dirname(@__FILE__), "../data")
const EK60_SAMPLE = joinpath(EK60_DATA,"JR230-D20091215-T121917.raw")
const ECS_SAMPLE = joinpath(EK60_DATA, "JR230.ecs")

end # module


text:
"""
    calc_main_freq(ode_sol::ODESolution, sol_sampling::Int)

Compute the main Fourier frequency and its power for each node of an ODE solution.

Input sol_sampling determines how many equally-spaced points are used to interpolate
the ODE solution.
"""
function calc_main_freq(ode_sol::ODESolution, sol_sampling::Int)
  N = length(ode_sol.u[1])
  freq = zeros(N)
  power = zeros(N)
  t_eval = LinRange(ode_sol.t[1], ode_sol.t[end], sol_sampling)
  fs = round(Int, sol_sampling/ode_sol.t[end])
  for i in 1:N
    # Remove mean to avoid peak at low frequencies
    signal = ode_sol(t_eval)[i,:] .- mean(ode_sol(t_eval)[i,:])
    pgram = periodogram(signal; fs=fs)
    pks, ~ = findmaxima(pgram.power)
    pks, ~ = peakproms(pks, pgram.power, minprom=maximum(pgram.power)/3.0)
    if length(pks) > 1
      # From the top peaks get the one with the lowest frequency
      best_vals = sort(collect(zip(pgram.freq[pks], pgram.power[pks])); by=first)
      freq[i] = best_vals[1][1]
      power[i] = best_vals[1][2]
    else
      freq[i] = pgram.freq[argmax(pgram.power)]
      power[i] = pgram.power[argmax(pgram.power)]
    end
  end
  return [freq, power]
end

"""
    calc_amp(ode_sol::ODESolution, peak_num_thresh::Int)

Compute the peak-to-peak distance for all variables in an ODE solution and the
coefficient of variation for maxima and minima peaks.

Input peak_num_thresh is an integer used to filter non-oscillatory solutions.
"""
function calc_amp(ode_sol::ODESolution, peak_num_thresh::Int)
  N = length(ode_sol.u[1])
  amp_val = zeros(N)
  amp_cv_max = zeros(N)
  amp_cv_min = zeros(N)
  for i in 1:N
    rough_amp = maximum(ode_sol[i,:]) - minimum(ode_sol[i,:])
    pks_max, ~ = findmaxima(ode_sol[i,:])
    pks_max, ~ = peakproms(pks_max, ode_sol[i,:], minprom=rough_amp/3.0)
    vals_max = ode_sol[i, pks_max]
    pks_min, ~ = findminima(ode_sol[i,:])
    pks_min, ~ = peakproms(pks_min, ode_sol[i,:], minprom=rough_amp/3.0)
    vals_min = ode_sol[i, pks_min]
    if length(vals_max) > peak_num_thresh
      if length(vals_min) > peak_num_thresh
        if length(vals_max) < 3*peak_num_thresh
          if length(vals_min) < 3*peak_num_thresh
            amp_val[i] = maximum(vals_max) - minimum(vals_min)
            amp_cv_max[i] = std(vals_max)/mean(vals_max)
            amp_cv_min[i] = std(vals_min)/mean(vals_min)
          end
        end
      end
    end
  end
  return [amp_val, amp_cv_max, amp_cv_min]
end

"""
    ode_params_from_param_set(model::ReactionSystem, p_set::Vector)

Reorder and normalize an array of parameters to be suitable for remaking an
ODEProblem.

All β and γ parameters are normalized by β1 to achieve a dimensionless time
system. The returned parameter array has the right order to use on ODE solvers.
Input variable p_set is assumed to be ordered as:
[α1, α2, ..., β1, ..., γ1, ..., κ1, ..., η1, ..., ηN]
"""
function ode_params_from_param_set(model::ReactionSystem, p_set::Vector)
  N = length(species(model))
  P = length(parameters(model))
  E = Int((P - 2N)/3)
  new_map_vals = []
  for j in 1:N
    @nonamespace push!(new_map_vals, (model.α[j], p_set[j]))
    @nonamespace push!(new_map_vals, (model.β[j], p_set[N+j]./p_set[N+1]))
  end
  for j in 1:E
    @nonamespace push!(new_map_vals, (model.γ[j], p_set[2*N+j]./p_set[N+1]))
    @nonamespace push!(new_map_vals, (model.κ[j], p_set[2*N+E+j]))
    @nonamespace push!(new_map_vals, (model.η[j], p_set[2*N+2*E+j]))
  end
  return ModelingToolkit.varmap_to_vars(Dict(new_map_vals), parameters(model))
end

"""
    generate_LHC_sample(model::ReactionSystem, samples::Int, lim::Dict)

Create an array of parameter sets by LatinHypercubeSampling.

A total of `samples` parameter sets are generated. Sample values are scaled by
the quantities defined in `lim`.
"""
function generate_LHC_sample(model::ReactionSystem, samples::Int, lim::Dict)
  N = length(species(model))
  P = length(parameters(model))
  E = Int((P - 2N)/3)
  LHCplan = randomLHC(samples, P)
  scaling_plan = Tuple{Float64,Float64}[]
  scaling_plan = vcat(scaling_plan, [(lim["α"][1], lim["α"][2]) for i in 1:N])
  scaling_plan = vcat(scaling_plan, [(lim["β"][1], lim["β"][2]) for i in 1:N])
  scaling_plan = vcat(scaling_plan, [(lim["γ"][1], lim["γ"][2]) for i in 1:E])
  scaling_plan = vcat(scaling_plan, [(lim["κ"][1], lim["κ"][2]) for i in 1:E])
  scaling_plan = vcat(scaling_plan, [(lim["η"][1], lim["η"][2]) for i in 1:E])
  return scaleLHC(LHCplan, scaling_plan)
end

"""
    create_generic_ode_problem(model::ReactionSystem)

Create a generic ODE problem from a model.
"""
function create_generic_ode_problem(model::ReactionSystem)
  N = length(species(model))
  P = length(parameters(model))
  E = Int((P - 2N)/3)
  p_dict_vals = []
  @nonamespace p_dict_vals = vcat(p_dict_vals, [(model.α[i], 0.5) for i in 1:N])
  @nonamespace p_dict_vals = vcat(p_dict_vals, [(model.β[i], 1.0) for i in 1:N])
  @nonamespace p_dict_vals = vcat(p_dict_vals, [(model.γ[i], 1.0) for i in 1:E])
  @nonamespace p_dict_vals = vcat(p_dict_vals, [(model.κ[i], 0.5) for i in 1:E])
  @nonamespace p_dict_vals = vcat(p_dict_vals, [(model.η[i], 1.0) for i in 1:E])
  pmap = Dict(p_dict_vals)
  @nonamespace u0_dict_vals = [(model.a[i], 0.5) for i in 1:N]
  u₀map = Dict(u0_dict_vals)
  return ODEProblem(model, u₀map, (0.0, 1.0), pmap)
end

"""
    infer_timescale(model::ReactionSystem, p_set::Vector)

Calculate the longest timescale of a model as the minimum of the reaction rates.
"""
function infer_timescale(model::ReactionSystem, p_set::Vector)
  N = length(species(model))
  P = length(parameters(model))
  E = Int((P - 2N)/3)
  # Only β2 through γN affect the timescale of the solution. Normalize by β1.
  rates = [p_set[k] for k in (N+1):(2*N+E)]./p_set[N+1]
  timescale = 1.0/minimum(rates)
  return timescale
end

"""
    is_oscil(freq::Vector, power::Vector, amp::Vector, amp_cv_max::Vector,
             amp_cv_min::Vector, hparams::Dict)

Decide if the features of an ODE solution indicate oscillations.
"""
function is_oscil(freq::Vector, power::Vector, amp::Vector, amp_cv_max::Vector,
                  amp_cv_min::Vector, hparams::Dict)
    N = length(freq)
    oscillatory = false
    freqtol = hparams["freq_tolerance"]
    freq_diffs = [abs(freq[j] - freq[j+1])/freq[j] < freqtol for j in 1:N-1]
    if all(freq_diffs)
      if all([power[j] > hparams["power_threshold"] for j in 1:N])
        if any(amp_cv_max .< hparams["amp_cv_thresh"])
          if any(amp_cv_min .< hparams["amp_cv_thresh"])
            if all(amp .> 0.0)
              oscillatory = true
            end
          end
        end
      end
    end
    return oscillatory
end

"""
    find_oscillations(model::ReactionSystem, samples::Int, param_limits::Dict,
                      hparams::Dict)

Compute parameter sets that produce oscillations in a model.
"""
function find_oscillations(model::ReactionSystem, samples::Int,
                           param_limits::Dict, hparams::Dict)
  N = length(species(model))
  p_sample = generate_LHC_sample(model, samples, param_limits)
  prob = create_generic_ode_problem(model)

  """
      prob_func(prob::ODEProblem, i, ~)

  Equilibrate parameter set and define both simulation time and starting point
  for the model's simulation using EnsembleProblem.
  """
  function prob_func(prob::ODEProblem, i, ~)
    sample_p_set = p_sample[Int(i),:]
    ode_p_set = ode_params_from_param_set(model, sample_p_set)
    #  Equilibrate parameter set
    timescale = infer_timescale(model, sample_p_set)
    t_final_eq = hparams["equil_tscales"]*timescale
    equil_prob = remake(prob, p=ode_p_set, u0=ones(N)/2.0,
                        tspan=[0.0, t_final_eq])
    equil_sol = solve(equil_prob, lsoda(), abstol=hparams["abstol"],
                      reltol=hparams["reltol"])
    # Set parameters for simulation
    midpoint = round(Int, length(equil_sol.t)/2.0)
    final_half_eq = equil_sol[midpoint:end]
    sol_sampling = length(final_half_eq.t)
    freq, ~ = calc_main_freq(final_half_eq, sol_sampling)
    if minimum(freq) > 0.0
      new_timescale = 0.5/minimum(freq) # Because I only considered half the sim
    else
      new_timescale = timescale
    end
    t_final = hparams["sim_tscales"] * new_timescale
    remake(prob, p=ode_p_set, u0=equil_sol.u[end], tspan=[0.0, t_final])
  end

  """
      output_func(sol::ODESolution, ~)

  Calculate frequency and amplitude and determine if solution is oscillatory.
  """
  function output_func(sol::ODESolution, ~)
    sol_sampling = length(sol.t)
    freq, power = calc_main_freq(sol, sol_sampling)
    amp, amp_cv_max, amp_cv_min = calc_amp(sol, hparams["peak_num_thresh"])
    oscillatory = is_oscil(freq, power, amp, amp_cv_max, amp_cv_min, hparams)
    endpoint = sol.u[end]
    out = [oscillatory, endpoint, freq, amp]
    return (out, false)
  end

  # Setup EnsembleProblem and solve it
  ensemble_prob = EnsembleProblem(prob, prob_func=prob_func, safetycopy=false,
                                  output_func=output_func)
  sim = solve(ensemble_prob, lsoda(), EnsembleThreads(),
              trajectories=samples, abstol=1e-12)
  # Return only data from oscillatory solutions
  osci_idxs = findall(x->x[1], sim.u)
  return [sim.u[osci_idxs], p_sample[osci_idxs,:]]
end


text:
function inverse(Side::HorizonSide)
    return HorizonSide((mod(Int(Side)+2, 4)))
end

function movements!(r::Robot, side::HorizonSide, n::Int)
    z = n
    while z != 0
        move!(r, side)
        z -= 1
    end
end

function shut!(r::Robot)
    n=0 # число шагов от начального положения
    while isborder(r,Nord)
        n += 1
        movements!(r, Ost, n)
        if !isborder(r, Nord)
            break
        end
        movements!(r, West, n*2)
        if !isborder(r, Nord)
            break
        end
        movements!(r, Ost, n)
    end
end


text:
using CSFML.LibCSFML

t1 = sfMicroseconds(1000000)
t2 = sfMilliseconds(1000)
t3 = sfSeconds(1)

clock = sfClock_create()
println("Clock initialized.")
println("sleep()")
sfSleep(t1)

seconds_elapsed = sfTime_asSeconds(sfClock_getElapsedTime(clock))
println("elapsed time : $seconds_elapsed sec.")

function get_time(clock ::Ptr{Nothing})
   sfTime_asSeconds(sfClock_getElapsedTime(clock))
end

function restart(clock ::Ptr{Nothing})
    sfTime_asSeconds(sfClock_restart(clock))
end

t = get_time(clock)
println("reseting clock at $t seconds.")
println(restart(clock))

println("sleep 3 seconds.")
sleep(3)
t = get_time(clock)
println("$t seconds have elapsed.")

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
using CamelCup

game = CamelCup.Game(6)
p, st = CamelCup.setup!(game)
while !game.over
    p, st = CamelCup.taketurn!(game, p[1].action)
    show(p)
end
p
st
p, st = CamelCup.taketurn!(game, CamelCup.TakeBet(CamelCup.Green))
p, st = CamelCup.taketurn!(game, CamelCup.Oasis(5, true))
p, st = CamelCup.taketurn!(game, CamelCup.Roll())

p, st = CamelCup.taketurn!(game, CamelCup.FinalBet(CamelCup.White, true))

p, st = CamelCup.taketurn!(game, CamelCup.Roll())
p, st = CamelCup.taketurn!(game, CamelCup.Roll())
p, st = CamelCup.taketurn!(game, CamelCup.Roll())
CamelCup.taketurn!(game, CamelCup.Roll())


text:
using JLD2, Compat, Base.Test

macro read(fid, sym)
    if !isa(sym, Symbol)
        error("Second input to @read must be a symbol (i.e., a variable)")
    end
    esc(:($sym = read($fid, $(string(sym)))))
end
macro write(fid, sym)
    if !isa(sym, Symbol)
        error("Second input to @write must be a symbol (i.e., a variable)")
    end
    esc(:(write($fid, $(string(sym)), $sym)))
end

# Define variables of different types
x = 3.7
A = reshape(1:15, 3, 5)
A3 = reshape(1:30, 3, 5, 2)
A4 = reshape(1:120, 2, 3, 4, 5)
Aarray = Vector{Float64}[[1.2,1.3],[2.2,2.3,2.4]]
basic_types = Any[UInt8(42), UInt16(42), UInt32(42), UInt64(42), UInt128(42),
                  Int8(42), Int16(42), Int32(42), Int64(42), Int128(42),
                  Float16(42), Float32(42), Float64(42)]
str = "Hello"
str_unicode = "pandapanda🐼panda"
str_embedded_null = "there is a null\0 in the middle of this string"
strings = String["It", "was", "a", "dark", "and", "stormy", "night"]
empty_string = ""
empty_string_array = String[]
empty_array_of_strings = String[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
c = Complex64(3,7)
cint = 1+im  # issue 108
C = reinterpret(Complex128, B, (4,))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
try
    global MyStruct
    mutable struct MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
d = Dict([(syms[1],"aardvark"), (syms[2], "banana")])
oidd = ObjectIdDict([(syms[1],"aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
T = UInt8
Tarr = DataType[UInt8, UInt16, UInt32, UInt64, UInt128]
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]]  # issue #123
typevar = Array{Int}[[1]]
typevar_lb = (Vector{U} where U<:Integer)[[1]]
typevar_ub = (Vector{U} where U>:Int)[[1]]
typevar_lb_ub = (Vector{U} where Int<:U<:Real)[[1]]
undef = Vector{Any}(1)
undefs = Matrix{Any}(2, 2)
ms_undef = MyStruct(0)
# Unexported type:
cpus = Base.Sys.cpu_info()
# Immutable type:
rng = 1:5
# Type with a pointer field (#84)
struct ObjWithPointer
    a::Ptr{Void}
end
objwithpointer = ObjWithPointer(0)
# Custom BitsType (#99)
primitive type MyBT 64 end
bt = reinterpret(MyBT, @compat Int64(55))
btarray = fill(bt, 5, 7)
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = view([1:5;], 1:5)
# Array of empty tuples (to test tuple type params)
arr_empty_tuple = (@compat Tuple{})[]
struct EmptyImmutable end
emptyimmutable = EmptyImmutable()
arr_emptyimmutable = [emptyimmutable]
empty_arr_emptyimmutable = EmptyImmutable[]
mutable struct EmptyType end
emptytype = EmptyType()
arr_emptytype = [emptytype]
empty_arr_emptytype = EmptyImmutable[]
uninitialized_arr_emptytype = Vector{EmptyType}(1)
struct EmptyII
    x::EmptyImmutable
end
emptyii = EmptyII(EmptyImmutable())
struct EmptyIT
    x::EmptyType
end
emptyit = EmptyIT(EmptyType())
mutable struct EmptyTI
    x::EmptyImmutable
end
emptyti = EmptyTI(EmptyImmutable())
type EmptyTT
    x::EmptyType
end
emptytt = EmptyTT(EmptyType())
struct EmptyIIOtherField
    x::EmptyImmutable
    y::Float64
end
emptyiiotherfield = EmptyIIOtherField(EmptyImmutable(), 5.0)
struct EmptyIIType
    x::Type{Float64}
    y::EmptyImmutable
end
emptyiitype = EmptyIIType(Float64, EmptyImmutable())

# Unicode type field names (#118)
mutable struct MyUnicodeStruct☺{τ}
    α::τ
    ∂ₓα::τ
    MyUnicodeStruct☺{τ}(α::τ, ∂ₓα::τ) where τ = new(α, ∂ₓα)
end
unicodestruct☺ = MyUnicodeStruct☺{Float64}(1.0, -1.0)
# Arrays of matrices (#131)
array_of_matrices = Matrix{Int}[[1 2; 3 4], [5 6; 7 8]]
# Tuple of arrays and bitstype
tup = (1, 2, [1, 2], [1 2; 3 4], bt)
# Empty tuple
empty_tup = ()
# Non-pointer-free struct
struct MyImmutable{T}
    x::Int
    y::Vector{T}
    z::Bool
end
nonpointerfree_immutable_1 = MyImmutable(1, [1., 2., 3.], false)
nonpointerfree_immutable_2 = MyImmutable(2, Any[3., 4., 5.], true)
struct MyImmutable2
    x::Vector{Int}
    MyImmutable2() = new()
end
nonpointerfree_immutable_3 = MyImmutable2()
struct Vague
    name::String
end
vague = Vague("foo")
# Immutable with a union of BitsTypes
struct BitsUnion
    x::Union{Int64,Float64}
end
bitsunion = BitsUnion(5.0)
# Immutable with a union of Types
struct TypeUnionField
    x::Type{T} where T<:Union{Int64,Float64}
end
typeunionfield = TypeUnionField(Int64)
# Generic union type field
struct GenericUnionField
    x::Union{Vector{Int},Int}
end
genericunionfield = GenericUnionField(1)
# Array references
arr_contained = [1, 2, 3]
arr_ref = typeof(arr_contained)[]
push!(arr_ref, arr_contained, arr_contained)
# Object references
mutable struct ObjRefType
    x::ObjRefType
    y::ObjRefType
    ObjRefType() = new()
    ObjRefType(x, y) = new(x, y)
end
ref1 = ObjRefType()
obj_ref = ObjRefType(ObjRefType(ref1, ref1), ObjRefType(ref1, ref1))
# Immutable that requires padding between elements in array
struct PaddingTest
    x::Int64
    y::Int8
end
padding_test = PaddingTest[PaddingTest(i, i) for i = 1:8]
# Empty arrays of various types and sizes
empty_arr_1 = Int[]
empty_arr_2 = Matrix{Int}(56, 0)
empty_arr_3 = Any[]
empty_arr_4 = Matrix{Any}(0, 97)
# Moderately big dataset (which will be mmapped)
bigdata = [1:1000000;]
# BigFloats and BigInts
bigint = big(3)
bigfloat = big(3.2)
bigints = big(3).^(1:100)
bigfloats = big(3.2).^(1:100)
struct BigFloatIntObject
    bigfloat::BigFloat
    bigint::BigInt
end
bigfloatintobj = BigFloatIntObject(big(pi), big(typemax(UInt128))+1)
# None
none = Union{}
nonearr = Vector{Union{}}(5)
# nothing/Void
scalar_nothing = nothing
vector_nothing = Union{Int,Void}[1,nothing]

# some data big enough to ensure that compression is used:
Abig = kron(eye(10), rand(20,20))
Bbig = Any[i for i=1:3000]
Sbig = "A test string "^1000

# Bitstype type parameters
mutable struct BitsParams{x}; end
bitsparamfloat  = BitsParams{1.0}()
bitsparambool   = BitsParams{true}()
bitsparamsymbol = BitsParams{:x}()
bitsparamint    = BitsParams{1}()
bitsparamuint   = BitsParams{0x01}()
bitsparamint16  = BitsParams{@compat Int16(1)}()

# Tuple of tuples
tuple_of_tuples = (1, 2, (3, 4, [5, 6]), [7, 8])

# Zero-dimensional arrays
zerod = Array{Int}()
zerod[] = 1
zerod_any = Array{Any}()
zerod_any[] = 1.0+1.0im

# Type with None typed field
struct NoneTypedField{T}
    a::Int
    b::T

    NoneTypedField{T}(a) where T = new(a)
end
nonetypedfield = NoneTypedField{Union{}}(47)

# Cyclic object
mutable struct CyclicObject
    x::CyclicObject

    CyclicObject() = new()
    CyclicObject(x) = new(x)
end
cyclicobject = CyclicObject()
cyclicobject.x = cyclicobject

# SimpleVector
simplevec = Core.svec(1, 2, Int64, "foo")
iseq(x::SimpleVector, y::SimpleVector) = collect(x) == collect(y)

# Issue #243
# Type that overloads != so that it is not boolean
mutable struct NALikeType; end
Base.:!=(::NALikeType, ::NALikeType) = NALikeType()
Base.:!=(::NALikeType, ::Void) = NALikeType()
Base.:!=(::Void, ::NALikeType) = NALikeType()
natyperef = Any[NALikeType(), NALikeType()]

iseq(x,y) = isequal(x,y)
function iseq(x::Array{EmptyType}, y::Array{EmptyType})
    size(x) != size(y) && return false
    for i = 1:length(x)
        def = isassigned(x, i)
        def != isassigned(y, i) && return false
        if def
            x[i] != y[i] && return false
        end
    end
    return true
end
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(x::MyImmutable, y::MyImmutable) = (isequal(x.x, y.x) && isequal(x.y, y.y) && isequal(x.z, y.z))
iseq(x::Union{EmptyTI,EmptyTT}, y::Union{EmptyTI,EmptyTT}) = isequal(x.x, y.x)
iseq(x::NoneTypedField{Union{}}, y::NoneTypedField{Union{}}) = x.a === y.a
iseq(c1::Array{Base.Sys.CPUinfo}, c2::Array{Base.Sys.CPUinfo}) = length(c1) == length(c2) && all([iseq(c1[i], c2[i]) for i = 1:length(c1)])
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in fieldnames(Base.Sys.CPUinfo)
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
iseq(x::MyUnicodeStruct☺, y::MyUnicodeStruct☺) = (x.α == y.α && x.∂ₓα == y.∂ₓα)
iseq(x::Array{Union{}}, y::Array{Union{}}) = size(x) == size(y)
iseq(x::BigFloatIntObject, y::BigFloatIntObject) = (x.bigfloat == y.bigfloat && x.bigint == y.bigint)
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                Base.show_backtrace(STDOUT, catch_backtrace())
                error("error reading ", $(string(sym)))
            end
            if !iseq(tmp, $sym)
                written = $sym
                error("For ", $(string(sym)), ", read value $tmp does not agree with written value $written")
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
        end
    end
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @assert a == b
function checkexpr(a::Expr, b::Expr)
    @assert a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if isa(a.args[i], Expr) && a.args[i].head == :line
            i += 1
            continue
        end
        if isa(b.args[j], Expr) && b.args[j].head == :line
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(a.args) && j >= length(b.args)
end

fn = joinpath(tempdir(),"test.jld")

# Issue #106
module Mod106
primitive type Typ{T} 64 end
typ{T}(x::Int64, ::Type{T}) = Base.box(Typ{T}, Base.unbox(Int64,x))
abstract type UnexportedT end
end

fid = jldopen(fn, "w")
println(fn)
@time begin
@write fid x
@write fid A
@write fid A3
@write fid A4
@write fid Aarray
@write fid basic_types
@write fid str
@write fid str_unicode
@write fid str_embedded_null
@write fid strings
@write fid empty_string
@write fid empty_string_array
@write fid empty_array_of_strings
@write fid tf
@write fid TF
@write fid AB
@write fid t
@write fid c
@write fid cint
@write fid C
@write fid emptyA
@write fid emptyB
@write fid ms
@write fid msempty
@write fid sym
@write fid syms
@write fid d
@write fid oidd
@write fid ex
@write fid T
@write fid Tarr
@write fid char
@write fid unicode_char
@write fid α
@write fid β
@write fid vv
@write fid cpus
@write fid rng
@write fid typevar
@write fid typevar_lb
@write fid typevar_ub
@write fid typevar_lb_ub
@write fid undef
@write fid undefs
@write fid ms_undef
@test_throws JLD2.PointerException @write fid objwithpointer
@write fid bt
@write fid btarray
@write fid sa_asc
@write fid sa_utf8
@write fid subarray
@write fid arr_empty_tuple
@write fid emptyimmutable
@write fid arr_emptyimmutable
@write fid empty_arr_emptyimmutable
@write fid emptytype
@write fid arr_emptytype
@write fid empty_arr_emptytype
@write fid uninitialized_arr_emptytype
@write fid emptyii
@write fid emptyit
@write fid emptyti
@write fid emptytt
@write fid emptyiiotherfield
@write fid emptyiitype
@write fid unicodestruct☺
@write fid array_of_matrices
@write fid tup
@write fid empty_tup
@write fid nonpointerfree_immutable_1
@write fid nonpointerfree_immutable_2
@write fid nonpointerfree_immutable_3
@write fid vague
@write fid bitsunion
@write fid typeunionfield
@write fid genericunionfield
@write fid arr_ref
@write fid obj_ref
@write fid padding_test
@write fid empty_arr_1
@write fid empty_arr_2
@write fid empty_arr_3
@write fid empty_arr_4
@write fid bigdata
@write fid bigfloat
@write fid bigint
@write fid bigfloats
@write fid bigints
@write fid bigfloatintobj
@write fid none
@write fid nonearr
@write fid scalar_nothing
@write fid vector_nothing
@write fid Abig
@write fid Bbig
@write fid Sbig
@write fid bitsparamint16
@write fid bitsparamfloat
@write fid bitsparambool
@write fid bitsparamsymbol
@write fid bitsparamint
@write fid bitsparamuint
@write fid tuple_of_tuples
@write fid zerod
@write fid zerod_any
@write fid nonetypedfield
@write fid cyclicobject
@write fid simplevec
@write fid natyperef
end
close(fid)

fidr = jldopen(fn, "r")
@time begin
@check fidr x
@check fidr A
@check fidr A3
@check fidr A4
@check fidr Aarray
@check fidr basic_types
@check fidr str
@check fidr str_unicode
@check fidr str_embedded_null
@check fidr strings
@check fidr empty_string
@check fidr empty_string_array
@check fidr empty_array_of_strings
@check fidr tf
@check fidr TF
@check fidr AB
@check fidr t
@check fidr c
@check fidr cint
@check fidr C
@check fidr emptyA
@check fidr emptyB
@check fidr ms
@check fidr msempty
@check fidr sym
@check fidr syms
@check fidr d
@check fidr oidd
exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
checkexpr(ex, exr)
@check fidr T
@check fidr Tarr
@check fidr char
@check fidr unicode_char
@check fidr α
@check fidr β
@check fidr vv
@check fidr cpus
@check fidr rng
@check fidr typevar
@check fidr typevar_lb
@check fidr typevar_ub
@check fidr typevar_lb_ub

# Special cases for reading undefs
undef = read(fidr, "undef")
if !isa(undef, Array{Any, 1}) || length(undef) != 1 || isassigned(undef, 1)
    error("For undef, read value does not agree with written value")
end
undefs = read(fidr, "undefs")
if !isa(undefs, Array{Any, 2}) || length(undefs) != 4 || any(map(i->isassigned(undefs, i), 1:4))
    error("For undefs, read value does not agree with written value")
end
ms_undef = read(fidr, "ms_undef")
if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
    error("For ms_undef, read value does not agree with written value")
end

@check fidr bt
@check fidr btarray
@check fidr sa_asc
@check fidr sa_utf8
@check fidr subarray
@check fidr arr_empty_tuple
@check fidr emptyimmutable
@check fidr arr_emptyimmutable
@check fidr empty_arr_emptyimmutable
@check fidr emptytype
@check fidr arr_emptytype
@check fidr empty_arr_emptytype
@check fidr uninitialized_arr_emptytype
@check fidr emptyii
@check fidr emptyit
@check fidr emptyti
@check fidr emptytt
@check fidr emptyiiotherfield
@check fidr emptyiitype
@check fidr unicodestruct☺
@check fidr array_of_matrices
@check fidr tup
@check fidr empty_tup
@check fidr nonpointerfree_immutable_1
@check fidr nonpointerfree_immutable_2
@check fidr nonpointerfree_immutable_3
vaguer = read(fidr, "vague")
@test typeof(vaguer) == typeof(vague) && vaguer.name == vague.name
@check fidr bitsunion
@check fidr typeunionfield
@check fidr genericunionfield

arr = read(fidr, "arr_ref")
@test arr == arr_ref
@test arr[1] === arr[2]

obj = read(fidr, "obj_ref")
@test obj.x.x === obj.x.y == obj.y.x === obj.y.y
@test obj.x !== obj.y

@check fidr padding_test
@check fidr empty_arr_1
@check fidr empty_arr_2
@check fidr empty_arr_3
@check fidr empty_arr_4
@check fidr bigdata
@check fidr bigfloat
@check fidr bigint
@check fidr bigfloats
@check fidr bigints
@check fidr bigfloatintobj
@check fidr none
@check fidr nonearr
@check fidr scalar_nothing
@check fidr vector_nothing
@check fidr Abig
@check fidr Bbig
@check fidr Sbig
@check fidr bitsparamfloat
@check fidr bitsparambool
@check fidr bitsparamsymbol
@check fidr bitsparamint
@check fidr bitsparamuint
@check fidr tuple_of_tuples
@check fidr zerod
@check fidr zerod_any
@check fidr nonetypedfield

obj = read(fidr, "cyclicobject")
@test obj.x === obj

@check fidr simplevec
@check fidr natyperef
end
close(fidr)


text:
"""
VolumeSlices

    volumeslices(x, y, z, v)

Draws heatmap slices of the volume v

## Attributes
$(ATTRIBUTES)
"""
@recipe(VolumeSlices, x, y, z, volume) do scene
    Attributes(;
        default_theme(scene, Heatmap)...,
        bbox_visible = true,
        bbox_color = RGBAf(0.5, 0.5, 0.5, 0.5)
    )
end

function plot!(plot::VolumeSlices)
    @extract plot (x, y, z, volume)
    replace_automatic!(plot, :colorrange) do
        map(extrema, volume)
    end

    # heatmap will fail if we don't keep its attributes clean
    attr = copy(Attributes(plot))
    bbox_color = pop!(attr, :bbox_color)
    bbox_visible = pop!(attr, :bbox_visible)
    pop!(attr, :model) # stops `transform!()` from working

    bbox = map(x, y, z) do x, y, z
        mx, Mx = extrema(x)
        my, My = extrema(y)
        mz, Mz = extrema(z)
        Rect3D(mx, my, mz, Mx-mx, My-my, Mz-mz)
    end
    linesegments!(plot, bbox, color = bbox_color, visible = bbox_visible, inspectable = false)

    axes = :x, :y, :z
    for (ax, p, r, (X, Y)) ∈ zip(axes, (:yz, :xz, :xy), (x, y, z), ((y, z), (x, z), (x, y)))
        hmap = heatmap!(plot, attr, X, Y, zeros(length(X[]), length(Y[])))
        plot[Symbol(:update_, p)] = i -> begin
            transform!(hmap, (p, r[][i]))
            indices = ntuple(Val(3)) do j
                axes[j] == ax ? i : (:)
            end
            hmap[3][] = view(volume[], indices...)
        end
        # Trigger once to place heatmaps correctly
        plot[Symbol(:update_, p)][](1)
    end
    plot
end


text:
#=
    All iterable expectations.
=#

# Callable behavior for the object. Parameters because we cannot add methods to an abstract type.
"""
    function (e::IterableExpectation{NT, WT})(f::Function; kwargs...) where {NT, WT}

Implements callable behavior for `IterableExpectation` objects.
"""
function (e::IterableExpectation{NT, WT})(f::Function; kwargs...) where {NT, WT}
    applicable(f, rand(nodes(e))) || throw(ArgumentError("The function doesn't accept elements from the distribution's support."))
    return dot(f.(nodes(e)), weights(e))
end

# Getters for the object.
"""
    nodes(e::IterableExpectation)

Returns the nodes of an `IterableExpectation`.
"""
nodes(e::IterableExpectation) = e.nodes

"""
    weights(e::IterableExpectation)

Returns the weights of an `IterableExpectation`.
"""
weights(e::IterableExpectation) = e.weights

# Linear operator behavior.
import Base.*

# Right-multiplying an expectation by something.
"""
    *(e::IterableExpectation, h::AbstractArray) = dot(h, weights(e))

Implements the right-application of an `IterableExpectation` by a vector of values on its nodes.
"""
*(e::IterableExpectation, h::AbstractArray) = dot(h, weights(e))

# Left-multiplying an expectation by a scalar.
"""
    *(r::Real, e::IterableExpectation) =  IterableExpectation(nodes(e), r * weights(e))

Implements left-multiplication of an `IterableExpectation` by a real scalar.
"""
*(r::Real, e::IterableExpectation) =  IterableExpectation(nodes(e), r * weights(e)) # Necessary because, for example, multiplying UnitRange * 2 = StepRange

#=
    Discrete iterable expectations.
=#

# Constructors for the object.
"""
    expectation(dist::DiscreteUnivariateDistribution, alg::Type{FiniteDiscrete} = FiniteDiscrete; kwargs...) = _expectation(dist, alg; kwargs...)

Dispatcher for (finite) discrete univariate expectations.
"""
expectation(dist::DiscreteUnivariateDistribution, alg::Type{FiniteDiscrete} = FiniteDiscrete; kwargs...) = _expectation(dist, alg; kwargs...)

"""
    function _expectation(dist::DiscreteUnivariateDistribution, alg::Type{FiniteDiscrete}; kwargs...)

Auxiliary constructor for an `IterableExpectation` object.
"""
function _expectation(dist::DiscreteUnivariateDistribution, alg::Type{FiniteDiscrete}; kwargs...)
    hasfinitesupport(dist) || throw(ArgumentError("Countably infinite distributions are not supported."))
    ourSupport = support(dist)
    ourWeights = pdf.(Ref(dist), support(dist))
    sum(ourWeights) ≈ 1.0 || warn("The distribution supplied is not approximately equal to 1 in mass.")
    return IterableExpectation(ourSupport, ourWeights);
end

#=
    Continuous iterable expectations (no nodes supplied.)
=#

# General catchall behavior --> Gauss-Legendre quadrature.
"""
    expectation(dist::ContinuousUnivariateDistribution, alg::Type{<:QuadratureAlgorithm} = Gaussian; kwargs...) = _expectation(dist, alg; kwargs...)

Dispatcher for continuous univariate expectations.
"""
expectation(dist::ContinuousUnivariateDistribution, alg::Type{<:QuadratureAlgorithm} = Gaussian; kwargs...) = _expectation(dist, alg; kwargs...)

"""
    function _expectation(dist::ContinuousUnivariateDistribution, alg::Type{Gaussian}; n = 500, kwargs...)

Implements Gauss-Legendre quadrature for continuous univariate distributions for which no specialized method exists.
"""
function _expectation(dist::ContinuousUnivariateDistribution, alg::Type{Gaussian}; n = 500, kwargs...)
    a = minimum(dist)
    b = maximum(dist)
    (a > -Inf && b < Inf) || throw(ArgumentError("The distribution must be defined on a compact interval."))
    rawNodes, rawWeights = gausslegendre(n)
    # Transform nodes to proper interval.
    nodes = map(x -> (0.5(b-a))*x + (a+b)/2, rawNodes)
    # Add pdf to weights.
    compoundWeights = [rawWeights[i] * pdf(dist, nodes[i]) for i in 1:n]
    # Add scale factor to weights.
    weights = (b-a)/2 * compoundWeights
    return IterableExpectation(nodes, weights);
end

# QNWDist implementation.
"""
    function _expectation(dist::ContinuousUnivariateDistribution, alg::Type{QuantileRange}; n::Int = 50, q0::Real = 0.001, qN::Real = 0.999 kwargs...)

Implementation of the qnwdist() quadrature scheme written by Spencer Lyon (PhD. NYU), as part of the QuantEcon.jl library. Used with permission. For further details, see: https://github.com/QuantEcon/QuantEcon.jl/blob/be0a32ec17d1f5b04ed8f2e52604c70c69f416b2/src/quad.jl#L892.
"""
function _expectation(dist::ContinuousUnivariateDistribution, alg::Type{QuantileRange}; n::Int = 50, q0::Real = 0.001, qN::Real = 0.999, kwargs...)
    # check nondegeneracy.
    all(isfinite.(params(dist))) || throw(ArgumentError("Distribution must be nondegenerate."))
    # _quadnodes in the QuantEcon.
    quantiles = Compat.range(q0, stop=qN, length=n)
    nodes = quantile.(Ref(dist), quantiles)
    # qnwdist in the QuantEcon.
    weights = zeros(n)
    for i in 2:n-1
        weights[i] = cdf(dist, (nodes[i] + nodes[i+1])/2) - cdf(dist, (nodes[i] + nodes[i-1])/2)
    end
    weights[1] = cdf(dist, (nodes[1] + nodes[2])/2)
    weights[end] = 1 - cdf(dist, (nodes[end-1] + nodes[end])/2)
    return IterableExpectation(nodes, weights)
end

# Specific method for normal distributions.
# Number of points was calibrated by trial.
"""
    function _expectation(dist::Normal, alg::Type{Gaussian}; n = 30, kwargs...)

Implements Gauss-Hermite quadrature for normal distributions.
"""
function _expectation(dist::Normal, alg::Type{Gaussian}; n = 30, kwargs...)
    σ = std(dist)
    μ = mean(dist)
    (isfinite(σ) && isfinite(μ)) || throw(ArgumentError("Parameters σ, μ must be finite."))
    gh = gausshermite(n)
    nodes = gh[1].*(sqrt(2)*(σ)) .+ μ
    weights = gh[2]./sqrt(pi)
    return IterableExpectation(nodes, weights)
end

# Specific method for lognormal distributions.
"""
    function _expectation(dist::LogNormal, alg::Type{Gaussian}; n = 30, kwargs...)

Implements Gauss-Hermite quadrature for lognormal distributions.
"""
function _expectation(dist::LogNormal, alg::Type{Gaussian}; n = 30, kwargs...) # Same settings for the normal method.
    m = mean(dist)
    v = var(dist)
    (isfinite(m) && isfinite(v)) || throw(ArgumentError("Infinite μ or σ^2 are not supported."))
    # get normal nodes
    gh = gausshermite(n)
    μ = log(m^2/sqrt(v + m^2))
    σ = sqrt(log(v/m^2 + 1))
    nodes = gh[1].*sqrt(2).*(σ) .+ μ
    weights = gh[2]./sqrt(pi)
    # get new nodes
    map!(x -> exp(x), nodes, nodes)
    return IterableExpectation(nodes, weights) # Transform the output.
end

# Specific method for beta distributions.
"""
    function _expectation(dist::Beta, alg::Type{Gaussian}; n = 32, kwargs...)

Implements Gauss-Jacobi quadrature for beta distributions.
"""
function _expectation(dist::Beta, alg::Type{Gaussian}; n = 32, kwargs...)
    α, β = params(dist)
    (isfinite(α) && isfinite(β)) || throw(ArgumentError("The beta distribution supplied is malformed."))
    gj = FastGaussQuadrature.JacobiRec(n, α-1, β-1)
    G = gamma(α)*gamma(β)/gamma(α+β)
    nodes = (1 .- gj[1])/2
    weights = gj[2]/((2.0^(α+β-1.0))*G)
    return IterableExpectation(nodes, weights)
end

# Specific method for exponential distributions.
"""
    function _expectation(dist::Exponential, alg::Type{Gaussian}; n = 32, kwargs...)

Implements Gauss-Laguerre quadrature for Exponential distributions.
"""
function _expectation(dist::Exponential, alg::Type{Gaussian}; n = 32, kwargs...)
    θ = inv(dist.θ) # To correct for the Distributions parametrization.
    isfinite(θ) || throw(ArgumentError("The beta distribution supplied is malformed."))
    gl = gausslaguerre(n)
    nodes = gl[1]./θ
    weights = gl[2]
    return IterableExpectation(nodes, weights)
end

# Specific method for gamma distributions.
"""
    function _expectation(dist::Gamma, alg::Type{Gaussian}; n = 32, kwargs...)

Implements Gauss-Laguerre quadrature for Gamma distributions.
"""
function _expectation(dist::Gamma, alg::Type{Gaussian}; n = 32, kwargs...)
    α, θ = params(dist)
    (isfinite(θ) && isfinite(θ)) || throw(ArgumentError("The beta distribution supplied is malformed."))
    gl = gausslaguerre(n, α-1)
    nodes = gl[1]./θ
    weights = gl[2]./gamma(α)
    return IterableExpectation(nodes, weights)
end

#=
    Continuous iterable distributions (nodes supplied.)
=#

# Dispatcher.
"""
    expectation(dist::ContinuousUnivariateDistribution, nodes, alg::Type{<:ExplicitQuadratureAlgorithm} = Trapezoidal; kwargs...)  = _expectation(dist, nodes, alg; kwargs...)

Dispatcher for distributions with user-defined nodes.
"""
expectation(dist::ContinuousUnivariateDistribution, nodes, alg::Type{<:ExplicitQuadratureAlgorithm} = Trapezoidal; kwargs...)  = _expectation(dist, nodes, alg; kwargs...)

# Trapezoidal general behavior.
"""
    function _expectation(dist, nodes::AbstractArray, alg::Type{Trapezoidal}; kwargs...)

Implements trapezoidal integration for general distributions with user-defined nodes.
"""
function _expectation(dist, nodes::AbstractArray, alg::Type{Trapezoidal}; kwargs...)
    M = length(nodes)
    Δ = diff(nodes)
    prepend!(Δ, NaN) # To keep the indexing straight. Now, Δ[2] = Δ_2 = z_2 - z_1. And NaN will throw an error if we try to use it.
    f_vec = pdf.(Ref(dist), nodes)
    interiorWeights = [f_vec[i]/2 * (Δ[i] + Δ[i+1]) for i = 2:M-1]
    allWeights = [f_vec[1]/2 * Δ[2]; interiorWeights; f_vec[M]/2 * Δ[M]]
    return IterableExpectation(nodes, allWeights)
end

# Trapezoidal for regular.
"""
    @compat function _expectation(dist, nodes::AbstractRange, alg::Type{Trapezoidal}; kwargs...)

Overloads trapezoidal integration for cases where the user-defined grids are regular.
"""
@compat function _expectation(dist, nodes::AbstractRange, alg::Type{Trapezoidal}; kwargs...)
    (first(nodes) >= minimum(dist) && last(nodes) <= maximum(dist)) || throw(ArgumentError("The nodes exceed the distribution's support."))
    M = length(nodes)
    Δ = nodes[2] - nodes[1]
    f_vec = pdf.(Ref(dist), nodes)
    interiorWeights = [f_vec[i] * Δ for i = 2:M-1]
    allWeights = [f_vec[1]/2 * Δ; interiorWeights; f_vec[M]/2 * Δ]
    return IterableExpectation(nodes, allWeights)
end

#=
    Convenience functions.
=#
"""
    expectation(f::Function, dist::DiscreteUnivariateDistribution, alg::Type{FiniteDiscrete} = FiniteDiscrete; kwargs...) = expectation(dist, alg; kwargs...)(f)

Convenience function for (finite) discrete univariate distributions.
"""
expectation(f::Function, dist::DiscreteUnivariateDistribution, alg::Type{FiniteDiscrete} = FiniteDiscrete; kwargs...) = expectation(dist, alg; kwargs...)(f)

"""
    expectation(f::Function, dist::ContinuousUnivariateDistribution, alg::Type{<:QuadratureAlgorithm} = Gaussian; kwargs...) = expectation(dist, alg; kwargs...)(f)

Convenience function for continuous univariate distributions.
"""
expectation(f::Function, dist::ContinuousUnivariateDistribution, alg::Type{<:QuadratureAlgorithm} = Gaussian; kwargs...) = expectation(dist, alg; kwargs...)(f)

"""
    expectation(f::Function, dist::ContinuousUnivariateDistribution, nodes::AbstractArray, alg::Type{<:ExplicitQuadratureAlgorithm} = Trapezoidal; kwargs...) = expectation(dist, nodes, alg; kwargs...)(f)

Convenience function for continuous univariate distributions with user-supplied nodes.
"""
expectation(f::Function, dist::ContinuousUnivariateDistribution, nodes::AbstractArray, alg::Type{<:ExplicitQuadratureAlgorithm} = Trapezoidal; kwargs...) = expectation(dist, nodes, alg; kwargs...)(f)


text:
"""
    X = spiraldata(N)

Generates data points on a 2D spiral returned as a N×2 matrix X.
"""
function spiraldata(N, seed=1)

    rg = MersenneTwister(seed)

    t = rand(rg, Uniform(3, 15), N)

    x1 = 0.04 * t .* sin.(t) .+ randn(rg, N)*0.01

    x2 = 0.04 * t .* cos.(t) .+ randn(rg, N)*0.01

    [x1 x2]

end


text:
include("core.jl")

function ck_train()
    aam = load_aam()

    # idxs = sample(1:count, count, replace=false)
    imgs = load_images(CKDataset, DATA_DIR_CK,
                       resizeratio=RESIZERATIO)
    shapes = load_shapes(CKDataset, DATA_DIR_CK,
                         resizeratio=RESIZERATIO)

    dataset = to_dataset(aam.wparams, imgs, shapes)
    mask = aam.wparams.warp_map
    imgs = nothing

    # nviewall([to_image(dataset[:, i], mask) for i=1:36])

    rbms = Array(Any, 8)
    for (i, p) in enumerate([512, 1024, 2048])
        println("-------- param1 = $p -----------")
        rbm = RBM(Degenerate, Bernoulli, size(dataset, 1),
                  p, sigma=0.01)
        @time fit(rbm, dataset, n_gibbs=1, lr=0.01,
                  batch_size=1000, n_epochs=20,
                  weight_decay_kind=:l2, weight_decay_rate=0.5,
                  sparsity_cost=0.1, sparsity_target=0.01)
        nviewall([to_image(rbm.W'[:, i], mask) for i=1:36])
        # nviewall([to_image(rbm.W'[:, i], mask) for i=37:72])
        rbms[i] = rbm
    end
    h5open(joinpath(MODEL_DIR, "rbm_512.h5"), "w") do h5
        save_params(h5, rbms[1], "model")
    end
    h5open(joinpath(MODEL_DIR, "rbm_1024.h5"), "w") do h5
        save_params(h5, rbms[2], "model")
    end
    h5open(joinpath(MODEL_DIR, "rbm_2048.h5"), "w") do h5
        save_params(h5, rbms[3], "model")
    end
end

function put_train()
    aam = load(joinpath(MODEL_DIR, "aam_put.jld"))["aam"]

    # idxs = sample(1:count, count, replace=false)
    imgs = collect(Array{Float64,3},
                   load_images(PutFrontalDataset, DATA_DIR_PUT))
    shapes = collect(Array{Float64,2},
                     load_shapes(PutFrontalDataset, DATA_DIR_PUT))

    dataset = to_dataset(aam.wparams, imgs, shapes)
    mask = aam.wparams.warp_map
    # imgs = nothing

    nviewall([to_image(dataset[:, i], mask)
              for i=rand(1:size(dataset,2), 9)])

    rbms = Array(Any, 8)
    for (i, p) in enumerate([3072])
        println("-------- param1 = $p -----------")
        rbm = RBM(Degenerate, Bernoulli, size(dataset, 1),
                  p, sigma=0.01)
        @time fit(rbm, dataset, n_gibbs=1, lr=0.01,
                  batch_size=1000, n_epochs=20,
                  weight_decay_kind=:l2, weight_decay_rate=0.5,
                  sparsity_cost=0.1, sparsity_target=0.01)
        nviewall([to_image(rbm.W'[:, i], mask) for i=1:36])
        # nviewall([to_image(rbm.W'[:, i], mask) for i=37:72])
        rbms[i] = rbm
    end
    h5open(joinpath(MODEL_DIR, "rbm_512.h5"), "w") do h5
        save_params(h5, rbms[1], "model")
    end
    h5open(joinpath(MODEL_DIR, "rbm_1024.h5"), "w") do h5
        save_params(h5, rbms[2], "model")
    end
    h5open(joinpath(MODEL_DIR, "rbm_2048.h5"), "w") do h5
        save_params(h5, rbms[3], "model")
    end
end

function kaggle_fer_train()
    imgs = load_images(:kaggle_fer, datadir=joinpath(DATA_DIR, "fer2013"))
    labels = load_lables(:kaggle_fer, datadir=joinpath(DATA_DIR, "fer2013"))
    dataset = hcat([reshape(img, 48*48) for img in imgs]...)
    rbm = RBM(Degenerate, Bernoulli, 48*48, 500)
    @time fit(rbm, dataset, lr=0.1, n_epochs=10,
              weight_decay_kind=:l2, weight_decay_rate=0.5)
              # sparsity_cost=0.1, sparsity_target=0.01)
    nviewall([reshape(rbm.W'[:, i], 48, 48) for i in 1:36])
end


text:
using OpenBLAS_jll, OpenBLAS32_jll, MKL_jll, CompilerSupportLibraries_jll
using Pkg, Artifacts, Base.BinaryPlatforms, Libdl, Test

include("utils.jl")

# Compile `dgemm_test.c` and `sgesv_test.c` against the given BLAS/LAPACK
function run_test((test_name, test_expected_outputs), libblas_name, libdirs, interface, backing_libs)
    # We need to configure this C build a bit
    cflags = String[
        "-g",
    ]
    if interface == :ILP64
        push!(cflags, "-DILP64")
    end

   # If our GCC is 64-bit but Julia is 32-bit, pass -m32
   if needs_m32()
       push!(cflags, "-m32")
   end

    ldflags = String[
        # Teach it to find that libblas and its dependencies at build time
        ("-L$(pathesc(libdir))" for libdir in libdirs)...,
        "-l$(libblas_name)",
    ]

    if !Sys.iswindows()
        # Teach it to find that libblas and its dependencies at run time
        append!(ldflags, ("-Wl,-rpath,$(pathesc(libdir))" for libdir in libdirs))
    end

    mktempdir() do dir
        @info("Compiling `$(test_name)` against $(libblas_name) ($(backing_libs)) in $(dir)")
        srcdir = joinpath(@__DIR__, test_name)
        make_cmd = `$(make) -sC $(pathesc(srcdir)) prefix=$(pathesc(dir)) CFLAGS="$(join(cflags, " "))" LDFLAGS="$(join(ldflags, " "))"`
        p = run(ignorestatus(make_cmd))
        if !success(p)
            @error("compilation failed", srcdir, prefix=dir, cflags=join(cflags, " "), ldflags=join(ldflags, " "))
        end
        @test success(p)

        env = Dict(
            # We need to tell it how to find CSL at run-time
            LIBPATH_env => append_libpath(libdirs),
            "LBT_DEFAULT_LIBS" => backing_libs,
        )
        cmd = `$(dir)/$(test_name)`
        output = capture_output(addenv(cmd, env))

        # Test to make sure the test ran properly
        has_expected_output = all(occursin(expected, output) for expected in test_expected_outputs)
        if !has_expected_output
            # Uh-oh, we didn't get what we expected.  Time to debug!
            @error("Test failed, got output:")
            println(output)

            # If we're not on CI, launch `gdb`
            if isempty(get(ENV, "CI", ""))
                debugger = Sys.isbsd() ? "lldb" : "gdb"
                @warn("Launching $debugger")
                cmd = `$(debugger) $(cmd)`
                env["LBT_VERBOSE"] = "1"
                run(addenv(cmd, env))
            end
        end
        @test has_expected_output
    end
end

# our tests
dgemm = ("dgemm_test", ("||C||^2 is:  24.3384",))
sgesv = ("sgesv_test", ("||b||^2 is:   3.0000",))
sdot  = ("sdot_test",  ("C is:   1.9900"))

# Build version that links against vanilla OpenBLAS
openblas_interface = :LP64
if Sys.WORD_SIZE == 64 && Sys.ARCH != :aarch64
    openblas_interface = :ILP64
end
openblas_jll_libname = splitext(basename(OpenBLAS_jll.libopenblas_path)[4:end])[1]
@testset "Vanilla OpenBLAS_jll ($(openblas_interface))" begin
    run_test(dgemm, openblas_jll_libname, OpenBLAS_jll.LIBPATH_list, openblas_interface, "")
    run_test(sgesv, openblas_jll_libname, OpenBLAS_jll.LIBPATH_list, openblas_interface, "")
    run_test(sdot,  openblas_jll_libname, OpenBLAS_jll.LIBPATH_list, openblas_interface, "")
end

# Build version that links against vanilla OpenBLAS32
@testset "Vanilla OpenBLAS32_jll (LP64)" begin
    run_test(dgemm, "openblas", OpenBLAS32_jll.LIBPATH_list, :LP64, "")
    run_test(sgesv, "openblas", OpenBLAS32_jll.LIBPATH_list, :LP64, "")
    run_test(sdot,  "openblas", OpenBLAS32_jll.LIBPATH_list, :LP64, "")
end

# Next, build a version that links against `libblastrampoline`, and tell
# the trampoline to forwards calls to `OpenBLAS_jll`
lbt_dir = joinpath(get_blastrampoline_dir(), binlib)

@testset "LBT -> OpenBLAS_jll ($(openblas_interface))" begin
    libdirs = unique(vcat(OpenBLAS_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list..., lbt_dir))
    run_test(dgemm, "blastrampoline", libdirs, openblas_interface, OpenBLAS_jll.libopenblas_path)
    run_test(sgesv, "blastrampoline", libdirs, openblas_interface, OpenBLAS_jll.libopenblas_path)
    run_test(sdot,  "blastrampoline", libdirs, openblas_interface, OpenBLAS_jll.libopenblas_path)
end

# And again, but this time with OpenBLAS32_jll
@testset "LBT -> OpenBLAS32_jll (LP64)" begin
    libdirs = unique(vcat(OpenBLAS32_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list..., lbt_dir))
    run_test(dgemm, "blastrampoline", libdirs, :LP64, OpenBLAS32_jll.libopenblas_path)
    run_test(sgesv, "blastrampoline", libdirs, :LP64, OpenBLAS32_jll.libopenblas_path)
    run_test(sdot,  "blastrampoline", libdirs, :LP64, OpenBLAS32_jll.libopenblas_path)
end

# Test against MKL_jll using `libmkl_rt`, which is :LP64 by default
if MKL_jll.is_available()
    @testset "LBT -> MKL_jll (LP64)" begin
        libdirs = unique(vcat(MKL_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list..., lbt_dir))
        run_test(dgemm, "blastrampoline", libdirs, :LP64, MKL_jll.libmkl_rt_path)
        run_test(sgesv, "blastrampoline", libdirs, :LP64, MKL_jll.libmkl_rt_path)
        run_test(sdot,  "blastrampoline", libdirs, :LP64, MKL_jll.libmkl_rt_path)
    end

    # Test that we can set MKL's interface via an environment variable to select ILP64, and LBT detects it properly
    if Sys.WORD_SIZE == 64
        @testset "LBT -> MKL_jll (ILP64, via env)" begin
            withenv("MKL_INTERFACE_LAYER" => "ILP64") do
                libdirs = unique(vcat(MKL_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list..., lbt_dir))
                run_test(dgemm, "blastrampoline", libdirs, :ILP64, MKL_jll.libmkl_rt_path)
                run_test(sgesv, "blastrampoline", libdirs, :ILP64, MKL_jll.libmkl_rt_path)
                run_test(sdot,  "blastrampoline", libdirs, :ILP64, MKL_jll.libmkl_rt_path)
            end
        end
    end
end

# Do we have Accelerate available?
veclib_blas_path = "/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/libBLAS.dylib"
if isfile(veclib_blas_path)
    @testset "LBT -> vecLib/libBLAS" begin
        run_test(dgemm, "blastrampoline", [lbt_dir], :LP64, veclib_blas_path)
        run_test(sdot,  "blastrampoline", [lbt_dir], :LP64, veclib_blas_path)
    end

    veclib_lapack_path = "/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/libLAPACK.dylib"
    @testset "LBT -> vecLib/libLAPACK" begin
        run_test(dgemm, "blastrampoline", [lbt_dir], :LP64, string(veclib_blas_path, ";", veclib_lapack_path))
        run_test(sgesv, "blastrampoline", [lbt_dir], :LP64, string(veclib_blas_path, ";", veclib_lapack_path))
        run_test(sdot,  "blastrampoline", [lbt_dir], :LP64, string(veclib_blas_path, ";", veclib_lapack_path))
    end
end

# Do we have a `blas64.so` somewhere?  If so, test with that for fun
blas64 = dlopen("libblas64", throw_error=false)
if blas64 !== nothing
    @testset "LBT -> libblas64 (ILP64, BLAS)" begin
        run_test(dgemm, "blastrampoline", [lbt_dir], :ILP64, dlpath(blas64))
        run_test(sdot,  "blastrampoline", [lbt_dir], :ILP64, dlpath(blas64))
        # Can't run `sgesv` here as we don't have LAPACK symbols in `libblas64.so`
    end

    # Check if we have a `liblapack` and if we do, run again, this time including `sgesv`
    lapack = dlopen("liblapack64", throw_error=false)
    if lapack !== nothing
        @testset "LBT -> libblas64 + liblapack64 (ILP64, BLAS+LAPACK)" begin
            run_test(dgemm, "blastrampoline", [lbt_dir], :ILP64, "$(dlpath(blas64));$(dlpath(lapack))")
            run_test(sgesv, "blastrampoline", [lbt_dir], :ILP64, "$(dlpath(blas64));$(dlpath(lapack))")
            run_test(sdot,  "blastrampoline", [lbt_dir], :ILP64, "$(dlpath(blas64));$(dlpath(lapack))")
        end
    end
end

# Finally the super-crazy test: build a binary that links against BOTH sets of symbols!
if openblas_interface == :ILP64
    inconsolable = ("inconsolable_test", ("||C||^2 is:  24.3384", "||b||^2 is:   3.0000"))
    @testset "LBT -> OpenBLAS 32 + 64 (LP64 + ILP64)" begin
        libdirs = unique(vcat(OpenBLAS32_jll.LIBPATH_list..., OpenBLAS_jll.LIBPATH_list..., CompilerSupportLibraries_jll.LIBPATH_list..., lbt_dir))
        run_test(inconsolable, "blastrampoline", libdirs, :wild_sobbing, "$(OpenBLAS32_jll.libopenblas_path);$(OpenBLAS_jll.libopenblas_path)")
    end
end

# Run our "direct" tests within Julia
include("direct.jl")

max_stars_repo_name: Gnimuc/libblastrampoline

max_issues_repo_name: Gnimuc/libblastrampoline

max_forks_repo_name: Gnimuc/libblastrampoline


text:
import MelGeneralizedCepstrums: allpass_alpha, glog_gamma

abstract type Filter end
abstract type SynthesisFilter <: Filter end

### Generic interfaces ###

# Note that f has state and filt! should update state value in this function
# TODO(ryuichi): separate state and filt definition?
function filt!(f::SynthesisFilter, x::Real, coef::AbstractVector)
    error("should provide sample-by-sample filt!")
end

# A synthesis filter must provide conversion from spectral parameter
# (e.g. mel-cepstrum) to filter coefficients.
function to_filtcoef(f::SynthesisFilter, param::AbstractArray)
    error("cannot determine how to convert paramters")
end

### Interfaces for mel-generalized ceptrum synthesis filters ###

abstract type MelGeneralizedCepstrumSynthesisFilter <: SynthesisFilter end

function allpass_alpha(f::MelGeneralizedCepstrumSynthesisFilter)
    error("should provide get access to allpass constant")
end

function glog_gamma(f::MelGeneralizedCepstrumSynthesisFilter)
    error("should provide get access to glog gamma")
end

# LPC, PARCOR and LSP
abstract type LinearPredictionVariantsSynthesisFilter <: SynthesisFilter end

### Waveform synthesis implementation ###

# synthesis_one_frame! generates a one frame speech waveform given a excitation
# signal and successive filter coefficients of a synthesis filter.
function synthesis_one_frame!(f::SynthesisFilter,
                              y::AbstractVector, # result will be stored
                              excitation::AbstractVector,
                              bᵗ⁻¹::AbstractVector,
                              bᵗ::AbstractVector)
    @assert length(y) == length(excitation)
    slope = (bᵗ - bᵗ⁻¹) / length(excitation)

    interpolated_coef = copy(bᵗ⁻¹)

    for i in 1:length(excitation)
        scaled_excitation = excitation[i] * exp(interpolated_coef[1])
        y[i] = filt!(f, scaled_excitation, interpolated_coef)
        for j=1:length(slope)
            interpolated_coef[j] += slope[j]
        end
    end

    y
end

function synthesis_one_frame!(f::SynthesisFilter,
                              excitation::AbstractVector,
                              bᵗ⁻¹::AbstractVector,
                              bᵗ::AbstractVector)
    synthesis_one_frame!(f, similar(excitation), excitation, bᵗ⁻¹, bᵗ)
end

# synthesis! generates a speech waveform given a excitation signal and
# a sequence of spectral envelope parameter.
function synthesis!(f::SynthesisFilter, # filter states are will be updated
                    excitation::AbstractVector,
                    b::AbstractMatrix,
                    hopsize::Integer)
    synthesized = similar(excitation)
    fill!(synthesized, zero(eltype(synthesized)))

    bᵗ⁻¹ = b[:,1]
    bᵗ = similar(bᵗ⁻¹)
    buf = Vector{eltype(synthesized)}(hopsize)

    for i in 1:size(b, 2)
        if i > 1
            bᵗ⁻¹ = b[:,i-1]
        end
        bᵗ = b[:,i]

        s, e = (i-1)*hopsize+1, i*hopsize
        e > length(excitation) && break

        synthesis_one_frame!(f, buf, excitation[s:e], bᵗ⁻¹, bᵗ)
        copy!(synthesized, s, buf, 1, length(buf))
    end

    synthesized
end

### Dispatch on specific filter types ###

function synthesis!{T<:MelGeneralizedCepstrum}(f::MelGeneralizedCepstrumSynthesisFilter,
                    excitation::AbstractVector,
                    state::SpectralParamState{T},
                    hopsize::Integer)
    synthesis!(f, excitation, to_filtcoef(f, state), hopsize)
end

function synthesis!{T<:LinearPredictionCoefVariants}(f::LinearPredictionVariantsSynthesisFilter,
                    excitation::AbstractVector,
                    state::SpectralParamState{T},
                    hopsize::Integer)
    synthesis!(f, excitation, to_filtcoef(f, state), hopsize)
end


text:
# Expression Matrix
# =================

# A matrix wrapper.
struct ExpressionMatrix{T}
    data::AbstractMatrix{T}
    featuremap::Dict{String,Int}
    featurenames::Vector{String}

    function ExpressionMatrix(data::AbstractMatrix{T}, featurenames) where T
        m, n = size(data)
        if m != length(featurenames)
            throw(ArgumentError("wrong number of feature names"))
        end
        featuremap = Dict{String,Int}()
        for (i, feature) in enumerate(featurenames)
            if haskey(featuremap, feature)
                throw(ArgumentError("found duplicated feature names: '$(feature)'"))
            end
            featuremap[feature] = i
        end
        return new{T}(data, featuremap, featurenames)
    end
end

Base.size(M::ExpressionMatrix) = size(M.data)
Base.size(M::ExpressionMatrix, d::Integer) = size(M.data, d)
Base.convert(::Type{ExpressionMatrix{S}}, M::ExpressionMatrix{T}) where {S,T} =
    ExpressionMatrix(convert(Matrix{S}, M.data), M.featurenames)

function Base.getindex(M::ExpressionMatrix, featurenames::AbstractVector{<:AbstractString}, j::Integer)
    m = length(featurenames)
    v = Vector{eltype(M)}(undef, m)
    for i in 1:m
        ii = get(M.featuremap, featurenames[i], 0)
        v[i] = ii > 0 ? M.data[ii,j] : 0
    end
    return v
end


text:
# Eric S. Tellez <eric.tellez@infotec.mx>

#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

using ArgParse
import JSON
import Glob

include("textmodel.jl")
include("io.jl")

function voc(fun, input, config::TextConfig, key)
    # patlist = [[normtext(pat, config) for pat in patterns] for patterns in patlist]
    thesaurus = Dict{String,Int}()

    itertweets(input) do tweet
        tokens = tokenize(getkeypath(tweet, key), config)

        for tok in tokens
            thesaurus[tok] = get(thesaurus, tok, 0) + 1
        end
    end

    L = [(v, k) for (k,v) in thesaurus]
    sort!(L, by=(x) -> x[1])
    for x in L
        fun(x)
    end
end

function main()
    config = TextConfig()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--config"
        help = "a json-formatted dictionary with the TextModel configuration"
        arg_type = String
        default = "{}"
        "--stopwords"
        help = "the file containing the list of stopwords, one per line"
        arg_type = String
        default = ""
        "--key"
        help = "specifies the key pointing to the text"
        arg_type = String
        default = "text"
    end

    args = parse_args(ARGS, s)
    config = TextConfig()
    for (k, v) in JSON.parse(args["config"])
        setfield!(config, Symbol(k), v)
    end
    load(config, args["stopwords"])

    voc(STDIN, config, args["key"]) do tweet
        write(STDOUT, JSON.json(tweet), '\n')
    end
end

main()

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
# Code borrowed from NNlib

expand(::Type{Val{N}}, i::Integer) where N = ntuple(_ -> i, Val(N))
expand(::Type{Val{N}}, i::NTuple{N, Integer}) where N = i

function cdims(x::NTuple{N}, w::NTuple{N}, pad, stride) where N
    ntuple(Val(N)) do i
        if i < N-1
            1 + div(x[i] - w[i] + 2*pad[i], stride[i])
        elseif i == N-1
            w[N]
        else # i == N
            x[N]
        end
    end
end

head(x) = reverse(Base.tail(reverse(x)))
padtuple(x::Tuple,p::Integer) = map(_->p, head(head(x)))
padtuple(x::Tuple,p::Tuple) = p
padtuple(x::AbstractArray,p) = padtuple(size(x),p)

function psize(p, x)
    nd = ndims(x)-2
    if isa(p,Number)
        fill(Int(p),nd)
    elseif length(p)==nd
        collect(Int,p)
    else
        throw(DimensionMismatch("psize: $p $nd"))
    end
end

function dilation_dims(w, dilation = 1)
    N = ndims(w)
    dims_w = size(w)
    dil = psize(dilation, w)
    ntuple(N) do i
        if i < N - 1
            (dims_w[i] - 1) * dil[i] + 1
        else
            dims_w[i]
        end
    end
end

function pdims(dims::Dims{N}, window, padding, stride) where N
    ntuple(Val(N)) do i
        if i < N-1
            1 + (dims[i] + 2*padding[i] - window[i])÷stride[i]
        else
            dims[i]
        end
    end
end

function check_support(x, k, pad, stride, dilation = 1)
    supported = true
    dilation == 1 || dilation == (1, 1) || (supported = false)
    pad_, stride_ = expand(Val{length(k)}, pad), expand(Val{length(k)}, stride)
    ((size(x, 1) - k[1] + 2 * pad_[1]) % stride_[1] == 0 && (size(x, 2) - k[2] + 2 * pad_[2]) % stride_[2] == 0) || (supported = false)
    !supported && error("Operation Not Supported by NNPACK")
    return pad_, stride_
end


text:
module EiSCor

    using Base.LinAlg: Eigen, chkstride1
    using Compat.LinAlg: checksquare
    import Base: copy
    import Base.LinAlg: eig, eigfact!, eigvals

    const libeiscor = joinpath(Pkg.dir("EiSCor"), "deps", "libeiscor")

    type UnitaryHessenbergMatrix{T} <: AbstractMatrix{T}
        data::Matrix{T}
    end
    copy(UH::UnitaryHessenbergMatrix) = UnitaryHessenbergMatrix(copy(UH.data))

    function ohfqr!(compz::Char, H::Matrix{Float64})
        # Check
        chkstride1(H)
        n = checksquare(H)
        compz in ('N', 'I', 'V') || throw(ArgumentError("argument compz must be either be either 'N', 'I', or 'V'"))

        # Allocate
        work = Array(Float64, 4n)
        Z = compz == 'N' ? Array(Float64, 0, 0) : Array(Float64, n, n)
        its = Array(Int32, n - 1)
        info = Array(Int32, 1)

        ccall((:dohfqr_, libeiscor), Void,
            (Ptr{UInt8}, Ptr{Int32}, Ptr{Float64}, Ptr{Float64},
             Ptr{Int32}, Ptr{Float64}, Ptr{Int32}),
             &compz, &n, H, Z,
             its, work, info)

        if info[1] != 0
            info[1] == 2 && error("DOFFQR failed")
            info[1] == 1 && error("DOHRFF failed")
        end
        H, Z, its
    end

    function uhfqr!(compz::Char, H::Matrix{Complex{Float64}})
        # Check
        chkstride1(H)
        n = checksquare(H)
        compz in ('N', 'I', 'V') || throw(ArgumentError("argument compz must be either be either 'N', 'I', or 'V'"))

        # Allocate
        work = Array(Complex{Float64}, 5n)
        Z = compz == 'N' ? Array(Complex{Float64}, 0, 0) : Array(Complex{Float64}, n, n)
        its = Array(Int32, n - 1)
        info = Array(Int32, 1)

        ccall((:zuhfqr_, libeiscor), Void,
            (Ptr{UInt8}, Ptr{Int32}, Ptr{Complex{Float64}}, Ptr{Complex{Float64}},
             Ptr{Int32}, Ptr{Complex{Float64}}, Ptr{Int32}),
             &compz, &n, H, Z,
             its, work, info)

        if info[1] != 0
            info[1] == 2 && error("ZUFFQR failed")
            info[1] == 1 && error("ZUHRFF failed")
        end
        H, Z, its
    end

    function eigvals(UH::UnitaryHessenbergMatrix{Float64})
        H, Z, its = ohfqr!('N', copy(UH.data))
        n = size(H, 1)
        v = zeros(Complex{Float64}, n)
        i = 1
        while i < n
            if H[i + 1, i] == 0
                v[i] = H[i,i]
                i += 1
            else
                v[i]    = complex(H[i,i], H[i,i+1])
                v[i+1]  = complex(H[i,i], H[i+1,i])
                i += 2
            end
        end
        i == n && (v[n] = H[n,n])
        return v
    end
    eigvals(UH::UnitaryHessenbergMatrix{Complex{Float64}}) = diag(uhfqr!('N', copy(UH.data))[1])

    function eigfact!(UH::UnitaryHessenbergMatrix{Float64})
        H, Z, its = ohfqr!('I', UH.data)
        n = size(H, 1)
        v = zeros(Complex{Float64}, n)
        i = 1
        while i < n
            if H[i + 1, i] == 0
                v[i] = H[i,i]
                i += 1
            else
                v[i]    = complex(H[i,i], H[i,i+1])
                v[i+1]  = complex(H[i,i], H[i+1,i])
                i += 2
            end
        end
        i == n && (v[n] = H[n,n])
        return Eigen(v, Z)
    end

    function eigfact!(UH::UnitaryHessenbergMatrix{Complex{Float64}})
        H, Z, its = uhfqr!('I', UH.data)
        return Eigen(diag(H), Z)
    end

end # module


text:
module TestExamples2dP4est

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "P4estMesh" begin
  @trixi_testset "elixir_advection_basic_p4est.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_p4est.jl"),
      l2   = [9.14468177884088e-6],
      linf = [6.437440532947036e-5])
  end

  @trixi_testset "elixir_advection_p4est_non_conforming_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4est_non_conforming_flag.jl"),
      l2   = [2.58174252995893e-5],
      linf = [2.65115939055204e-4])
  end

  @trixi_testset "elixir_advection_p4est_non_conforming_flag_unstructured.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4est_non_conforming_flag_unstructured.jl"),
      l2   = [3.038384623519386e-3],
      linf = [6.324792487776842e-2])
  end

  @trixi_testset "elixir_advection_amr_solution_independent_p4est.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent_p4est.jl"),
      l2   = [7.945073295691892e-5],
      linf = [0.0007454287896710293])
  end

  @trixi_testset "elixir_advection_amr_p4est_unstructured_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_p4est_unstructured_flag.jl"),
      l2   = [1.29561551e-03],
      linf = [3.92061383e-02])
  end

  @trixi_testset "elixir_advection_restart_p4est.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_p4est.jl"),
      l2   = [6.398955192910044e-6],
      linf = [3.474337336717426e-5])
  end

  @trixi_testset "elixir_euler_nonperiodic_p4est.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_p4est.jl"),
      l2   = [2.3653424742684444e-6, 2.1388875095440695e-6, 2.1388875095548492e-6, 6.010896863397195e-6],
      linf = [1.4080465931654018e-5, 1.7579850587257084e-5, 1.7579850592586155e-5, 5.956893531156027e-5])
  end

  @trixi_testset "elixir_euler_nonperiodic_p4est_unstructured.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_p4est_unstructured.jl"),
      l2   = [0.005689496253354319, 0.004522481295923261, 0.004522481295922983, 0.009971628336802528],
      linf = [0.05125433503504517, 0.05343803272241532, 0.053438032722404216, 0.09032097668196482])
  end

  @trixi_testset "elixir_euler_free_stream_p4est.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_p4est.jl"),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 4.636291350834654e-13, 4.956035581926699e-13],
      atol = 5e-13, # required to make CI tests pass on macOS
    )
  end

  @trixi_testset "elixir_eulergravity_eoc_test_p4est.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test_p4est.jl"),
      l2   = [0.0002456507718735638, 0.0003374236900407321, 0.00033742369004074064, 0.0007198398465056179],
      linf = [0.0013653498212873316, 0.0018996803803537077, 0.0018996803803517093, 0.004384818684607161],
      tspan = (0.0, 0.1))
  end
end

end # module


text:
using Revise
using ProtoBuf
using V5RPC
import V5RPC:Field,EventArguments,destructure

struct MyStrategy <: Strategy
    robot_state
end

MyStrategy() = begin
    robot_state = map(1:5) do _
        Dict(:e_int => 0.0, :e_last => 0.0)
    end
    MyStrategy(robot_state)
end

V5RPC.get_team_info(::MyStrategy) = "Julia中文社区"

function V5RPC.get_instruction(s::MyStrategy, field::Field)
    map(1:5) do id
        instruction(s, V5RPC.destructure(field), id)
    end
end

function V5RPC.on_event(::MyStrategy, event_type::Int32, event_args::EventArguments)
    ty = lookup(V5RPC.EventType, event_type)
    println("Event: $ty $event_args")
end
#=
using DataStructures
using Makie
samples=Node(CircularBuffer{Float64}(100))
push!(samples[],0.0)
scene=lines(samples)
function addsample(x)
    global samples,scene
    push!(samples[],Float64(x))
    samples[]=samples[]
    center!(scene)
    update_limits!(scene)
    #display(scene)
end
=#
includet("impl.jl")

task1 = @async run(V5Server(MyStrategy(), 20000))
task2 = @async run(V5Server(EmptyStrategy(), 20001))

max_stars_repo_name: Runfa-Zhang/MeetUpMaterials

max_issues_repo_name: Runfa-Zhang/MeetUpMaterials

max_forks_repo_name: Runfa-Zhang/MeetUpMaterials


text:
"""
    Snapshot(snapbase, subbase, config::GadgetConfig)

Returns a Snapshot from a `snapbase` and a `subbase`, which can also be `nothing`.
The `config` includes information on the block properties of the files.
Both can be files or the base part of files ending with `.0`, `.1`, etc.
"""
struct Snapshot{T}
    snapbase::Union{T,Nothing}
    subbase::Union{T,Nothing}
    config::GadgetConfig
end

"""
    Snapshot(; snapbase=nothing, subbase=nothing, config=GadgetConfig())

Convenience method for creating a `Snapshot` where one of the base files does not exist.
"""
Snapshot(; snapbase=nothing, subbase=nothing, config=GadgetConfig()) = Snapshot(snapbase, subbase, config)

"""
    Snapshot(path_box, snapbase, subbase; config=GadgetConfig())

Convenience method for creating a `Snapshot` where `snapbase` and `subbase` are given as
relative paths from `path_box`.
"""
Snapshot(path_box, snapbase, subbase; config=GadgetConfig()) =
    Snapshot(joinpath(path_box, snapbase), joinpath(path_box, subbase), config)

"""
    Snapshot(path_box, snap::Integer; snapbase=true, subbase=true, config=GadgetConfig())

Convenience method for creating a `Snapshot` from the path to a simulated box and a snapshot number.
Setting `snapbase` or `subbase` to `false` sets the respective file base to `nothing`.

The formats used are the following, for `snapbase` and `subbase`, respectively:
- `snapdir_XXX/snap_XXX`
- `groups_XXX/sub_XXX`
"""
function Snapshot(path_box, snap::Integer; snapbase::Bool=true, subbase::Bool=true, config=GadgetConfig())
    _snapbase = snapbase ? get_snapbase(path_box, snap) : nothing
    _subbase = subbase ? get_subbase(path_box, snap) : nothing
    Snapshot(_snapbase, _subbase, config)
end

get_snapbase(path_box, snap::String) = joinpath(path_box, "snapdir_$snap", "snap_$snap")
get_snapbase(path_box, snap::Integer) = get_snapbase(path_box, lpad(snap, 3, '0'))
get_subbase(path_box, snap::String) = joinpath(path_box, "groups_$snap", "sub_$snap")
get_subbase(path_box, snap::Integer) = get_subbase(path_box, lpad(snap, 3, '0'))

function Base.show(io::IO, ::MIME"text/plain", obj::Snapshot)
    printstyled(io, Snapshot; bold=true)
    print(io)
    isnothing(obj.snapbase) || print(io, "\nSnap: $(obj.snapbase)")
    isnothing(obj.subbase) || print(io, "\nSub:  $(obj.subbase)")
end