text:
# GOKU-net model
#
# Based on
# https://arxiv.org/abs/2003.10775

struct GOKU <: LatentDE end

apply_feature_extractor(encoder::Encoder{GOKU}, x) = encoder.feature_extractor.(x)

@doc raw"""
    apply_pattern_extractor(encoder::Encoder{GOKU}, fe_out)
Passes extracted features through the pattern_extractor layer contained in the `encoder`, returning a tuple corresponding to the initial latent states and parameters, respectively.
# Arguments
`fe_out`: Output of feature extractor layer
"""
function apply_pattern_extractor(encoder::Encoder{GOKU}, fe_out)
    pe_z₀, pe_θ_forward, pe_θ_backward = encoder.pattern_extractor

    # reverse sequence
    fe_out_rev = reverse(fe_out)

    # pass it through the recurrent layers
    pe_z₀_out = map(pe_z₀, fe_out_rev)[end]
    pe_θ_out_f = map(pe_θ_forward, fe_out)[end]
    pe_θ_out_b = map(pe_θ_backward, fe_out_rev)[end]
    pe_θ_out = vcat(pe_θ_out_f, pe_θ_out_b)

    # reset hidden states
    Flux.reset!(pe_z₀)
    Flux.reset!(pe_θ_forward)
    Flux.reset!(pe_θ_backward)

    return pe_z₀_out, pe_θ_out
end

function apply_latent_in(encoder::Encoder{GOKU}, pe_out)
    pe_z₀_out, pe_θ_out = pe_out
    li_μ_z₀, li_logσ²_z₀, li_μ_θ, li_logσ²_θ = encoder.latent_in

    z₀_μ = li_μ_z₀(pe_z₀_out)
    z₀_logσ² = li_logσ²_z₀(pe_z₀_out)

    θ_μ = li_μ_θ(pe_θ_out)
    θ_logσ² = li_logσ²_θ(pe_θ_out)

    return (z₀_μ, θ_μ), (z₀_logσ², θ_logσ²)
end

@doc raw"""
    apply_latent_out(decoder::Decoder{GOKU}, l̃)
Applies the `decoder`'s `latent_out` layer to `l̃`, obtaining the inferred initial conditions and parameters for use in the differential equation layer.
# Arguments
`l̃`: Tuple containing sampled abstract representations of the initial conditions and parameters, respectively.
"""
function apply_latent_out(decoder::Decoder{GOKU}, l̃)
    z̃₀, θ̃ = l̃
    lo_z₀, lo_θ = decoder.latent_out

    ẑ₀ = lo_z₀(z̃₀)
    θ̂ = lo_θ(θ̃)

    return ẑ₀, θ̂
end

@doc raw"""
    diffeq_layer(decoder::Decoder{GOKU}, l̂, t)
Uses decoder.diffeq's DE solver to extrapolate the latent states (from the initial states and parameters l̂) to time t.
"""
function diffeq_layer(decoder::Decoder{GOKU}, l̂, t)
    ẑ₀, θ̂ = l̂
    prob = decoder.diffeq.prob
    solver = decoder.diffeq.solver
    sensealg = decoder.diffeq.sensealg
    kwargs = decoder.diffeq.kwargs

    # Function definition for ensemble problem
    prob_func(prob,i,repeat) = remake(prob, u0=ẑ₀[:,i], p = θ̂[:,i])

    # Check if solve was successful, if not, return NaNs to avoid problems with dimensions matches
    output_func(sol, i) = sol.retcode == :Success ? (Array(sol), false) : (fill(NaN32,(size(ẑ₀, 1), length(t))), false)

    ## Adapt problem to given time span and create ensemble problem definition
    prob = remake(prob; tspan = (t[1],t[end]))
    ens_prob = EnsembleProblem(prob, prob_func = prob_func, output_func = output_func)

    ## Solve
    ẑ = solve(ens_prob, solver, EnsembleThreads(); sensealg = sensealg, trajectories = size(θ̂, 2), saveat = t, kwargs...)

    # Transform the resulting output (mainly used for Kuramoto-like systems)
    ẑ = transform_after_diffeq(ẑ, decoder.diffeq)
    ẑ = Flux.unstack(ẑ, 2)
    return ẑ
end

# Identity by default
transform_after_diffeq(x, diffeq) = x

apply_reconstructor(decoder::Decoder{GOKU}, ẑ) = decoder.reconstructor.(ẑ)

@doc raw"""
    sample(μ::T, logσ²::T, model::LatentDiffEqModel{GOKU}) where T <: Tuple{Array, Array}
Samples abstract representations of the parameters and initial state from the normal distribution with mean μ and variance exp(logσ²).
"""
function sample(μ::T, logσ²::T, model::LatentDiffEqModel{GOKU}) where T <: Tuple{Array, Array}
    z₀_μ, θ_μ = μ
    z₀_logσ², θ_logσ² = logσ²

    ẑ₀ = z₀_μ + randn(Float32, size(z₀_logσ²)) .* exp.(z₀_logσ²/2f0)
    θ̂ =  θ_μ + randn(Float32, size( θ_logσ²)) .* exp.(θ_logσ²/2f0)

    return ẑ₀, θ̂
end

function sample(μ::T, logσ²::T, model::LatentDiffEqModel{GOKU}) where T <: Tuple{Flux.CUDA.CuArray, Flux.CUDA.CuArray}
    z₀_μ, θ_μ = μ
    z₀_logσ², θ_logσ² = logσ²

    ẑ₀ = z₀_μ + gpu(randn(Float32, size(z₀_logσ²))) .* exp.(z₀_logσ²/2f0)
    θ̂ =  θ_μ + gpu(randn(Float32, size( θ_logσ²))) .* exp.(θ_logσ²/2f0)

    return ẑ₀, θ̂
end

function default_layers(model_type::GOKU, input_dim, diffeq; device=cpu,
                            hidden_dim_resnet = 200, rnn_input_dim = 32,
                            rnn_output_dim = 16, latent_dim = 16,
                            latent_to_diffeq_dim = 200, θ_activation = softplus,
                            output_activation = σ, init = Flux.kaiming_uniform(gain = 1/sqrt(3)))

    z_dim = length(diffeq.prob.u0)
    θ_dim = length(diffeq.prob.p)

    ######################
    ### Encoder layers ###
    ######################
    # Resnet
    l1 = Dense(input_dim, hidden_dim_resnet, relu, init = init)
    l2 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l3 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l4 = Dense(hidden_dim_resnet, rnn_input_dim, relu, init = init)
    feature_extractor = Chain(l1,
                                SkipConnection(l2, +),
                                SkipConnection(l3, +),
                                l4) |> device

    # RNN
    pe_z₀ = Chain(RNN(rnn_input_dim, rnn_output_dim, relu, init = init),
                       RNN(rnn_output_dim, rnn_output_dim, relu, init = init)) |> device

    # Bidirectional LSTM
    pe_θ_forward = Chain(LSTM(rnn_input_dim, rnn_output_dim, init = init),
                       LSTM(rnn_output_dim, rnn_output_dim, init = init)) |> device

    pe_θ_backward = Chain(LSTM(rnn_input_dim, rnn_output_dim, init = init),
                        LSTM(rnn_output_dim, rnn_output_dim, init = init)) |> device

    pattern_extractor = (pe_z₀, pe_θ_forward, pe_θ_backward)

    # final fully connected layers before sampling
    li_μ_z₀ = Dense(rnn_output_dim, latent_dim, init = init) |> device
    li_logσ²_z₀ = Dense(rnn_output_dim, latent_dim, init = init) |> device

    li_μ_θ = Dense(rnn_output_dim*2, latent_dim, init = init) |> device
    li_logσ²_θ = Dense(rnn_output_dim*2, latent_dim, init = init) |> device

    latent_in = (li_μ_z₀, li_logσ²_z₀, li_μ_θ, li_logσ²_θ)

    encoder_layers = (feature_extractor, pattern_extractor, latent_in)

    ######################
    ### Decoder layers ###
    ######################

    # after sampling in the latent space but before the differential equation layer
    lo_z₀ = Chain(Dense(latent_dim, latent_to_diffeq_dim, relu, init = init),
                        Dense(latent_to_diffeq_dim, z_dim, init = init)) |> device

    lo_θ = Chain(Dense(latent_dim, latent_to_diffeq_dim, relu, init = init),
                        Dense(latent_to_diffeq_dim, θ_dim, θ_activation, init = init)) |> device

    latent_out = (lo_z₀, lo_θ)

    # going back to the input space
    # Resnet
    l1 = Dense(z_dim, hidden_dim_resnet, relu, init = init)
    l2 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l3 = Dense(hidden_dim_resnet, hidden_dim_resnet, relu, init = init)
    l4 = Dense(hidden_dim_resnet, input_dim, output_activation, init = init)
    reconstructor = Chain(l1,
                            SkipConnection(l2, +),
                            SkipConnection(l3, +),
                            l4)  |> device

    decoder_layers = (latent_out, diffeq, reconstructor)

    return encoder_layers, decoder_layers
end


text:
using DocStringExtensions

"""
Represents a rule scheduler for the equality saturation process

"""
abstract type AbstractScheduler end

"""
Should return true if the e-graph can be said to be saturated
"""
cansaturate(s::AbstractScheduler)::Bool = error("not implemented")

"""
Should return true if the rule `r` should be skipped
"""
shouldskip(s::AbstractScheduler, r::Rule) = error("not implemented")

"""
This function is called **before** pattern matching on the e-graph
"""
readstep!(s::AbstractScheduler) = nothing

"""
This function is called **after** pattern matching on the e-graph
"""
writestep!(s::AbstractScheduler, r::Rule) = nothing

ext: jl


text:
"""
ISORANK
-------
    solve an overlap matching problem with IsoRank
"""

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                 alpha::Float64,rtype::Int64,tol::Float64,
                 maxit::Int64,verbose::Bool,P::SparseMatrixCSC{Float64,Int64})
  csum = sum_kbn(w)
  v = w./csum
  assert(all(v.>=0))
  n = size(P,1)

  allstats = !isempty(li) && !isempty(lj)
  if allstats
    rhistsize = 5
    M_setup = MatrixNetworks.bipartite_matching_setup(w,li,lj,maximum(li),maximum(lj))
    tripi = M_setup.tripi
    matm = M_setup.m
    matn = M_setup.n
    rp = M_setup.rp
    ci = M_setup.ci
    mperm = tripi[tripi.>0] # a permutation for the matching problem
  else
    rhistsize = 1
  end
  r = alpha
  x = zeros(Float64,n) + v
  delta = 2
  iter = 0
  reshist = zeros(Float64,maxit,rhistsize)
  xbest = x
  fbest = 0
  fbestiter = 0
  if verbose && allstats #print the header
    @printf("%5s   %4s   %8s   %7s %7s %7s %7s\n","best","iter","pr-delta","obj","weight","card","overlap")
  elseif verbose
    @printf("%4s   %8s\n","iter","delta")
  end
  while iter < maxit && delta > tol
    y = r * (P' * x)
    omega = sum_kbn(x) - sum_kbn(y)
    y = y + omega * v
    delta = vecnorm(x-y,1)
    reshist[iter+1] = delta
    iter = iter + 1
    x = y * (1/sum_kbn(y))
    if allstats
      if rtype == 1
        xf = x
      elseif rtype == 2
        xf = a*v + b/2*(S*x) #add v to preserve scale
      end
      ai = zeros(Float64,length(tripi))
      ai[tripi.>0] = xf[mperm]
      M_output = MatrixNetworks.bipartite_matching_primal_dual(rp,ci,ai,matm,matn)
      ma = M_output.cardinality
      #get match:
      mi_int = MatrixNetworks.edge_indicator(M_output,li,lj)
      val = dot(w,mi_int)
      overlap = dot(mi_int,(S*mi_int)/2)
      f = a*val + b*overlap
      if f > fbest
        xbest = x
        fbest = f
        fbestiter = iter
        itermark = "*"
      else
        itermark = " "
      end
      if verbose && allstats
        @printf("%5s   %4i   %8.1e   %5.2f %7.2f %7i %7i\n", itermark, iter, delta, f, val, ma, overlap)
        reshist[iter,2:end] = [a*val + b*overlap, val, ma, overlap]
      elseif verbose
        @printf("%4i    %8.1e\n", iter, delta)
      end
    end
  end
  flag = delta>tol
  reshist = reshist[1:iter,:]
  if allstats
    x=xbest
  end
  return (x,flag,reshist)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                 alpha::Float64,rtype::Int64,tol::Float64,
                 maxit::Int64,verbose::Bool)
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                 alpha::Float64,rtype::Int64,tol::Float64,
                 maxit::Int64)
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                 alpha::Float64,rtype::Int64,tol::Float64)
  maxit = 100
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                 alpha::Float64,rtype::Int64)
  tol = 1e-12
  maxit = 100
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64},
                 alpha::Float64)
  rtype = 2
  tol = 1e-12
  maxit = 100
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64,li::Vector{Int64},lj::Vector{Int64})
  alpha = b/(a+b)
  rtype = 2
  tol = 1e-12
  maxit = 100
  verbose = true
  P = normout(S)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64,b::Int64)
  li = []
  lj = []
  alpha = b/(a+b)
  rtype = 2
  tol = 1e-12
  maxit = 100
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64},
                 a::Int64)
  b = 1
  li = []
  lj = []
  alpha = b/(a+b)
  rtype = 2
  tol = 1e-12
  maxit = 100
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

function isorank{T}(S::SparseMatrixCSC{T,Int64},w::Vector{Float64})
  a = 0.5
  b = 1
  li = []
  lj = []
  alpha = b/(a+b)
  rtype = 2
  tol = 1e-12
  maxit = 100
  verbose = true
  ss = sum(S,2)
  P = sparse(S./ss)
  return isorank(S,w,a,b,li,lj,alpha,rtype,tol,maxit,verbose,P)
end

max_stars_repo_name: nassarhuda/NetworkAlign

max_issues_repo_name: nassarhuda/NetworkAlign

max_forks_repo_name: nassarhuda/NetworkAlign


text:
# Bogumił Kamiński, 2022

# Codes for chapter 8

# Code for section 8.1

import Downloads
if isfile("puzzles.csv.bz2")
    @info "file already present"
else
    @info "fetching file"
    Downloads.download("https://database.lichess.org/" *
                       "lichess_db_puzzle.csv.bz2",
                       "puzzles.csv.bz2")
end

using CodecBzip2
compressed = read("puzzles.csv.bz2")
plain = transcode(Bzip2Decompressor, compressed)

open("puzzles.csv", "w") do io
    println(io, "PuzzleId,FEN,Moves,Rating,RatingDeviation," *
                "Popularity,NbPlays,Themes,GameUrl")
    write(io, plain)
end

readlines("puzzles.csv")

# Code for section 8.2

using CSV
using DataFrames
puzzles = CSV.read("puzzles.csv", DataFrame);

puzzles2 = CSV.read(plain, DataFrame;
                    header=["PuzzleId", "FEN", "Moves",
                            "Rating","RatingDeviation",
                            "Popularity", "NbPlays",
                            "Themes","GameUrl"]);
puzzles == puzzles2

compressed = nothing
plain = nothing

# Code for listing 8.1

puzzles

# Code for listing 8.2

describe(puzzles)

# Code for getting basic information about a data frame

ncol(puzzles)

nrow(puzzles)

names(puzzles)

CSV.write("puzzles2.csv", puzzles)

read("puzzles2.csv")

read("puzzles2.csv") == read("puzzles.csv")

# Code for section 8.3

puzzles.Rating

using BenchmarkTools
@btime $puzzles.Rating;

puzzles.Rating == copy(puzzles.Rating)

puzzles.Rating === copy(puzzles.Rating)

puzzles.Rating === puzzles.Rating

copy(puzzles.Rating) === copy(puzzles.Rating)

puzzles."Rating"

col = "Rating"

# data_frame_name[selected_rows, selected_columns]

puzzles[:, "Rating"]
puzzles[:, :Rating]
puzzles[:, 4]
puzzles[:, col]

columnindex(puzzles, "Rating")

columnindex(puzzles, "Some fancy column name")

hasproperty(puzzles, "Rating")
hasproperty(puzzles, "Some fancy column name")

@btime $puzzles[:, :Rating];

puzzles[!, "Rating"]
puzzles[!, :Rating]
puzzles[!, 4]
puzzles[!, col]

using Plots
plot(histogram(puzzles.Rating, label="Rating"),
     histogram(puzzles.RatingDeviation, label="RatingDeviation"),
     histogram(puzzles.Popularity, label="Popularity"),
     histogram(puzzles.NbPlays, label="NbPlays"))

plot([histogram(puzzles[!, col]; label=col) for
      col in ["Rating", "RatingDeviation",
              "Popularity", "NbPlays"]]...)

# Code for section 8.4

# Codes for Arrow examples
using Arrow
Arrow.write("puzzles.arrow", puzzles)

arrow_table = Arrow.Table("puzzles.arrow")
puzzles_arrow = DataFrame(arrow_table);
puzzles_arrow == puzzles

puzzles_arrow.PuzzleId
puzzles_arrow.PuzzleId[1] = "newID"

puzzles_arrow = copy(puzzles_arrow);
puzzles_arrow.PuzzleId

# Codes for SQLite examples
using SQLite
db = SQLite.DB("puzzles.db")

SQLite.load!(puzzles, db, "puzzles")

SQLite.tables(db)
SQLite.columns(db, "puzzles")

query = DBInterface.execute(db, "SELECT * FROM puzzles")
puzzles_db = DataFrame(query);
puzzles_db == puzzles

puzzles_db.PuzzleId

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
import PredictMD
import Test

Test.@test_throws ErrorException PredictMD.simple_moving_average([], -1)


text:
module BackwardAutodiff
import Base: exp, log, sin, cos, tan, +, -, *, /, sqrt, convert, promote_rule, zero

export D

# This is the structure that stores the autodiff information.  The name is
# `Infinitesimal` spelled backward.  The elements are `x`, the result of a
# computation; `dfdy` a storage location where the derivative of some overall
# computation (`f(x)`) with respect to the output of this element is stored;
# `parent` storing the "parent" locations that have been used to produce the
# result `x`; and `bp!` a function that propogates the derivative `dfdy` to the
# parents via `dfdx = dydx * dfdy` where `dydx` is the derivative of this piece
# of the computation, `dfdy` is the derivative of the overall function with
# respect to the output, and `dfdx` is the same for the next level up.
mutable struct Lamisetinifni{T <: Number} <: Number
    x::T
    dfdy::T
    parent::Union{Lamisetinifni{T}, Array{Lamisetinifni{T},1}, Int, Nothing}
    bp!
end

function convert(::Type{Lamisetinifni{T}}, x::T) where T <: Number
    Lamisetinifni(x, zero(T), nothing, (dfdy, parents) -> nothing)
end

function convert(::Type{Lamisetinifni{T}}, x::S) where {T, S <: Number}
    Lamisetinifni(T(x), zero(T), nothing, (dfdy, parents) -> nothing)
end

function convert(::Type{Lamisetinifni{T}}, x::Lamisetinifni{T}) where T
    x
end

function zero(x::Lamisetinifni{T}) where T
    Lamisetinifni(zero(T), zero(T), nothing, (dfdy, parents) -> nothing)
end

function promote_rule(::Type{Lamisetinifni{T}}, ::Type{S}) where {T, S <: Number}
    Lamisetinifni{promote_type(T,S)}
end

function promote_rule(::Type{T}, ::Type{Lamisetinifni{S}}) where {T <: Number, S}
    Lamisetinifni{promote_type(T,S)}
end

function promote_rule(::Type{S}, ::Type{Lamisetinifni{T}}) where {S <: AbstractIrrational, T}
    Lamisetinifni{promote_type(S, T)}
end

function push_parents!(queue::Array{Lamisetinifni{T}, 1}, ::Nothing) where T
    # Do nothing
end
function push_parents!(queue::Array{Lamisetinifni{T}, 1}, i::Int) where T
    # Do nothing
end
function push_parents!(queue::Array{Lamisetinifni{T}, 1}, ls::Array{Lamisetinifni{T}, 1}) where T
    append!(queue, ls)
end
function push_parents!(queue::Array{Lamisetinifni{T}, 1}, l::Lamisetinifni{T}) where T
    push!(queue, l)
end

function backprop!(l::Lamisetinifni{T}) where T
    # Apparently we need this construction because otherwise l gets copied when
    # put into the array.
    backprop!([l])
end
function backprop!(queue::Array{Lamisetinifni{T},1}) where T
    while length(queue) > 0
        l = popfirst!(queue)
        l.bp!(l.dfdy, l.parent)
        push_parents!(queue, l.parent)
    end
end

function collect_outputs(l::Lamisetinifni{T}) where T
    queue = Lamisetinifni{T}[l]

    outputs = Lamisetinifni{T}[]

    while length(queue) > 0
        l = popfirst!(queue)
        if typeof(l.parent) <: Int
            push!(outputs, l)
        elseif typeof(l.parent) == Lamisetinifni{T}
            push!(queue, l.parent)
        elseif typeof(l.parent) == Array{Lamisetinifni{T}, 1}
            append!(queue, l.parent)
        else # Nothing
            # Do nothing
        end
    end

    outputs
end

"""    D([i], f)

Returns a function that computes the derivative of `f` (if it is
single-argument) or the gradient of `f` (if it is multi-argument or takes a
structured argument).  If `i` is given returns the `i`th component of the
gradient (though this does not reduce the cost with backprop).

Currently works only for `f` with scalar outputs.  Also, note that the autodiff
will fail unless the output type is identical to the input type (the code
automatically converts constants and non-differentiated expressions to the
appropriate type---only the input and output type of the function needs to
match).

So, for example, `D(cos)(3)` will fail (because the output type is `Float64` not
`Int`), while `D(cos)(3.0)` will work fine and return `-sin(3.0)`.

"""
function D(f)
    function dfdx(x::T) where T <: Number
        # Pass the function a backward infinitesimal whose backprop function
        # stores the backprop derivative in dfdx_store

        x = Lamisetinifni(x, zero(x), 1, (dfdy, parents) -> nothing)

        result = f(x)

        result.dfdy = one(result.x)
        backprop!(result)

        y = collect_outputs(result)[1]

        return y.dfdy
    end

    function dfdx(x::Array{T}) where T <: Number
        fargs = [Lamisetinifni(xelt, zero(xelt), i, (dfdy, parents) -> nothing) for (i, xelt) in enumerate(x)]
        result = f(fargs)
        result.dfdy = one(result.x)
        backprop!(result)
        y = collect_outputs(result)

        grad = zeros(typeof(result.x), length(x))
        for yelt in y
            grad[yelt.parent] = yelt.dfdy
        end

        return grad
    end

    function dfdx(x...)
        fargs = [Lamisetinifni(xelt, zero(xelt), i, (dfdy, parents) -> nothing) for (i, xelt) in enumerate(x)]
        result = f(fargs...)
        result.dfdy = one(result.x)
        backprop!(result)
        y = collect_outputs(result)

        grad = zeros(typeof(result.x), length(x))
        for yelt in y
            grad[yelt.parent] = yelt.dfdy
        end

        return grad
    end

    return dfdx
end

function D(i::Integer, f)
    df = D(f)
    function df_wrapper(x...)
        g = df(x...)
        return g[i]
    end
    return df_wrapper
end

function bpp!(dfdy, xy)
    x, y = xy
    x.dfdy += dfdy
    y.dfdy += dfdy
end

function +(x::Lamisetinifni{T}, y::Lamisetinifni{T}) where T
    Lamisetinifni(x.x + y.x, zero(T), [x, y], bpp!)
end

function bpm!(dfdy, xy)
    x, y = xy
    x.dfdy += dfdy
    y.dfdy -= dfdy
end
function -(x::Lamisetinifni{T}, y::Lamisetinifni{T}) where T
    Lamisetinifni(x.x - y.x, zero(T), [x, y], bpm!)
end

function bpum!(dfdy, x)
    x.dfdy -= dfdy
end
function -(x::Lamisetinifni{T}) where T
    Lamisetinifni(-x.x, zero(T), x, bpum!)
end

function bpt!(dfdy, xy)
    x,y = xy
    x.dfdy += dfdy*y.x
    y.dfdy += x.x*dfdy
end
function *(x::Lamisetinifni{T}, y::Lamisetinifni{T}) where T
    Lamisetinifni(x.x*y.x, zero(T), [x,y], bpt!)
end

function /(x::Lamisetinifni{T}, y::Lamisetinifni{T}) where T
    yinv = one(T)/y.x

    function bp!(dfdy, xy)
        a,b = xy
        a.dfdy += dfdy*yinv
        b.dfdy -= a.x*dfdy*(yinv*yinv)
    end

    Lamisetinifni(x.x*yinv, zero(T), [x,y], bp!)
end

function exp(x::Lamisetinifni{T}) where T
    expx = exp(x.x)

    function bp!(dfdy, p)
        p.dfdy += dfdy*expx
    end

    Lamisetinifni(expx, zero(expx), x, bp!)
end

function sin(x::Lamisetinifni{T}) where T
    function bp!(dfdy, p)
        p.dfdy += cos(x.x)*dfdy
    end

    sx = sin(x.x)
    Lamisetinifni(sx, zero(sx), x, bp!)
end

function cos(x::Lamisetinifni{T}) where T
    function bp!(dfdy, p)
        p.dfdy -= sin(x.x)*dfdy
    end

    cx = cos(x.x)
    Lamisetinifni(cx, zero(cx), x, bp!)
end

function tan(x::Lamisetinifni{T}) where T
    c = cos(x.x)
    function bp!(dfdy, p)
        p.dfdy += dfdy/(c*c)
    end

    tx = tan(x.x)
    Lamisetinifni(tx, zero(tx), x, bp!)
end

function sqrt(x::Lamisetinifni{T}) where T
    sqrtx = sqrt(x.x)

    function bp!(dfdy, p)
        p.dfdy += dfdy/(2*sqrtx)
    end

    Lamisetinifni(sqrtx, zero(sqrtx), x, bp!)
end

function log(x::Lamisetinifni{T}) where T
    function bp!(dfdy, p)
        p.dfdy += dfdy/x.x
    end

    Lamisetinifni(log(x.x), zero(T), x, bp!)
end

end


text:
export CuDeviceArray

struct CuDeviceArray{T,N}
    ptr::Ptr{T}
    dims::NTuple{N,Cint}
    strides::NTuple{N,Cint}
end

function CuDeviceArray(x::AbstractCuArray)
    CuDeviceArray(pointer(x), map(Cint,size(x)), map(Cint,strides(x)))
end

Base.length(x::CuDeviceArray) = Int(prod(x.dims))

@generated function Base.copyto!(dest::CuDeviceArray{T,N}, src::CuDeviceArray{T,N}) where {T,N}
    Ct = cstring(T)
    k = Kernel("""
    __global__ void copy(Array<$Ct,$N> dest, Array<$Ct,$N> src) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= src.length()) return;

        int sub[$N];
        dest.ind2sub(sub, idx);
        dest(sub) = src(sub);
    }""")
    quote
        gdims, bdims = cudims(length(src))
        $k(gdims, bdims, dest, src)
        dest
    end
end

@generated function addto!(dest::CuDeviceArray{T,N}, src::CuDeviceArray{T,N}) where {T,N}
    Ct = cstring(T)
    k = Kernel("""
    __global__ void add(Array<$Ct,$N> dest, Array<$Ct,$N> src) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= src.length()) return;

        int sub[$N];
        dest.ind2sub(sub, idx);
        dest(sub) += src(sub);
    }
    """)
    quote
        gdims, bdims = cudims(length(src))
        $k(gdims, bdims, dest, src)
        dest
    end
end


text:
ENV["GKS_ENCODING"] = "utf-8"

atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch e
        @warn "Error initializing Revise" exception=(e, catch_backtrace())
    end

    try
        @eval using OhMyREPL
    catch e
        @warn "error while importing OhMyREPL" e
    end

    try
        @eval using ClearStacktrace
    catch e
        @warn e.msg
    end
end

max_stars_repo_name: briochemc/dotfiles

max_issues_repo_name: briochemc/dotfiles

max_forks_repo_name: briochemc/dotfiles


text:
using SLEEFInline
using Test
using Printf

using Base.Math: significand_bits

isnzero(x::T) where {T <: AbstractFloat} = signbit(x)
ispzero(x::T) where {T <: AbstractFloat} = !signbit(x)

function cmpdenorm(x::Tx, y::Ty) where {Tx <: AbstractFloat, Ty <: AbstractFloat}
    sizeof(Tx) < sizeof(Ty) ? y = Tx(y) : x = Ty(x) # cast larger type to smaller type
    (isnan(x) && isnan(y)) && return true
    (isnan(x) || isnan(y)) && return false
    (isinf(x) != isinf(y)) && return false
    (x == Tx(Inf)  && y == Ty(Inf))  && return true
    (x == Tx(-Inf) && y == Ty(-Inf)) && return true
    if y == 0
        (ispzero(x) && ispzero(y)) && return true
        (isnzero(x) && isnzero(y)) && return true
        return false
    end
    (!isnan(x) && !isnan(y) && !isinf(x) && !isinf(y)) && return sign(x) == sign(y)
    return false
end

# the following compares the ulp between x and y.
# First it promotes them to the larger of the two types x,y
const infh(::Type{Float64}) = 1e300
const infh(::Type{Float32}) = 1e37
function countulp(T, x::AbstractFloat, y::AbstractFloat)
    X, Y = promote(x, y)
    x, y = T(X), T(Y) # Cast to smaller type
    (isnan(x) && isnan(y)) && return 0
    (isnan(x) || isnan(y)) && return 10000
    if isinf(x)
        (sign(x) == sign(y) && abs(y) > infh(T)) && return 0 # relaxed infinity handling
        return 10001
    end
    (x ==  Inf && y ==  Inf) && return 0
    (x == -Inf && y == -Inf) && return 0
    if y == 0
        x == 0 && return 0
        return 10002
    end
    if isfinite(x) && isfinite(y)
        return T(abs(X - Y) / ulp(y))
    end
    return 10003
end

const DENORMAL_MIN(::Type{Float64}) = 2.0^-1074
const DENORMAL_MIN(::Type{Float32}) = 2f0^-149

function ulp(x::T) where {T<:AbstractFloat}
    x = abs(x)
    x == T(0.0) && return DENORMAL_MIN(T)
    val, e = frexp(x)
    return max(ldexp(T(1.0), e - significand_bits(T) - 1), DENORMAL_MIN(T))
end

countulp(x::T, y::T) where {T <: AbstractFloat} = countulp(T, x, y)

# get rid off annoying warnings from overwritten function
macro nowarn(expr)
    quote
        _stderr = stderr
        tmp = tempname()
        stream = open(tmp, "w")
        redirect_stderr(stream)
        result = $(esc(expr))
        redirect_stderr(_stderr)
        close(stream)
        result
    end
end

# overide domain checking that base adheres to
using Base.MPFR: ROUNDING_MODE
for f in (:sin, :cos, :tan, :asin, :acos, :atan, :asinh, :acosh, :atanh, :log, :log10, :log2, :log1p)
    @eval begin
        import Base.$f
        @nowarn function ($f)(x::BigFloat)
            z = BigFloat()
            ccall($(string(:mpfr_, f), :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
            return z
        end
    end
end

strip_module_name(f::Function) = last(split(string(f), '.')) # strip module name from function f

# test the accuracy of a function where fun_table is a Dict mapping the function you want
# to test to a reference function
# xx is an array of values (which may be tuples for multiple arugment functions)
# tol is the acceptable tolerance to test against
function test_acc(T, fun_table, xx, tol; debug = false, tol_debug = 5)
    @testset "accuracy $(strip_module_name(xfun))" for (xfun, fun) in fun_table
        rmax = 0.0
        rmean = 0.0
        xmax = map(zero, first(xx))
        for x in xx
            q = xfun(x...)
            c = fun(map(BigFloat, x)...)
            u = countulp(T, q, c)
            rmax = max(rmax, u)
            xmax = rmax == u ? x : xmax
            rmean += u
            if debug && u > tol_debug
                @printf("%s = %.20g\n%s  = %.20g\nx = %.20g\nulp = %g\n", strip_module_name(xfun), q, strip_module_name(fun), T(c), x, ulp(T(c)))
            end
        end
        rmean = rmean / length(xx)

        t = @test trunc(rmax, digits=1) <= tol

        fmtxloc = isa(xmax, Tuple) ? string('(', join((@sprintf("%.5f", x) for x in xmax), ", "), ')') : @sprintf("%.5f", xmax)
        println(rpad(strip_module_name(xfun), 18, " "), ": max ", @sprintf("%f", rmax),
            rpad(" at x = " * fmtxloc, 40, " "),
            ": mean ", @sprintf("%f", rmean))
    end
end

function runtests()
    @testset "SLEEF" begin
        include("exceptional.jl")
        include("accuracy.jl")
    end
end

runtests()


text:
module BMTest

using DelimitedFiles
using LinearAlgebra
using Random
using Test
using Statistics

import BoltzmannMachines
const BMs = BoltzmannMachines

macro check(expr)
   msg = string(expr) # TODO: print vales
   quote
      if !($(esc(expr)))
         @warn "Check failed:\n\t" * $msg
         return false
      else
         true # no return
      end
   end
end

function createsamples(nsamples::Int, nvariables::Int, samerate::Float64 = 0.7)

   x = round.(rand(nsamples, nvariables))
   samerows = ((x[:,1] .== x[:,2]) .& (x[:,2] .== x[:,3]))
   sameratestart = count(samerows) / nsamples

   sameratediff = samerate - sameratestart

   if sameratediff > 0
      addedsamerange = findall(.!samerows)
      shuffle!(addedsamerange)
      nsameadded = round(Int, sameratediff * nsamples)
      addedsamerange = addedsamerange[1:nsameadded]
      x[addedsamerange, 3] = x[addedsamerange, 2] = x[addedsamerange, 1]
   else
      subtractedsamerange = findall(samerows)
      shuffle!(subtractedsamerange)
      nsamesubtracted = round(Int, - sameratediff * nsamples)
      subtractedsamerange = subtractedsamerange[1:nsamesubtracted]
      for row in subtractedsamerange
         while ((x[row,1] == x[row,2]) & (x[row,2] == x[row,3]))
            xrow = view(x, row, 1:3)
            rand!(xrow)
            map!(round, xrow, xrow)
         end
      end
   end
   x
end

function createsamples_categorical(nsamples::Int, nvariables::Int,
      samerates::Vector{Float64} = [0.2; 0.2; 0.4])

   x = zeros(nsamples, nvariables)
   nsames = Int.(floor.(nsamples .* samerates))
   sampleranges = BMs.ranges(nsames)
   for i in eachindex(samerates)
      samplerange = sampleranges[i]
      x[samplerange, 2] .= x[samplerange, 3] .= x[samplerange, 1] .= i - 1
   end

   notsamerange = (sampleranges[end][end] + 1):nsamples
   ncategories = length(samerates) - 1
   for i in 1:ncategories
      x[notsamerange, 1:3] .+= round.(rand(length(notsamerange), 3))
      x[:, 4:end] .+= round.(rand(nsamples, nvariables - 3))
   end
   x = x[shuffle!(collect(1:size(x,1))), :]
   x
end

function randrbm(nvisible, nhidden, factorw = 1.0, factora = 1.0, factorb = 1.0)
   w = factorw*randn(nvisible, nhidden)
   a = factora*rand(nvisible)
   b = factorb*(0.5 .- rand(nhidden))
   BMs.BernoulliRBM(w, a, b)
end

function randgbrbm(nvisible, nhidden, factorw = 1.0, factora = 1.0, factorb = 1.0, factorsd = 1.0)
   w = factorw*randn(nvisible, nhidden)
   a = factora*rand(nvisible)
   b = factorb*ones(nhidden)
   sd = factorsd*ones(nvisible)
   BMs.GaussianBernoulliRBM(w, a, b, sd)
end

function randbgrbm(nvisible, nhidden, factorw = 1.0, factora = 1.0, factorb = 1.0)
   w = factorw*randn(nvisible, nhidden)
   a = factora*rand(nvisible)
   b = factorb*ones(nhidden)
   BMs.BernoulliGaussianRBM(w, a, b)
end

function randdbm(nunits)
   nrbms = length(nunits) - 1
   dbm = BMs.BasicDBM(undef, nrbms)
   for i = 1:nrbms
      dbm[i] = randrbm(nunits[i], nunits[i+1])
   end
   dbm
end

"""
Tests the functions for computing the activation potentials
"""
function test_potentials()
   nvisible = 3
   nhidden = 2
   nsamples = 5
   vv = rand(nsamples, nvisible)

   # Test Test activation potential of hidden nodes for BernoulliRBM
   hh = rand(nsamples, nhidden)
   rbm = BMTest.randrbm(nvisible, nhidden)
   BMs.hiddenpotential!(hh, rbm, vv)
   @test isapprox(sum(abs.(hh - BMs.hiddenpotential(rbm, vv))), 0)

   # Test that activation potential of hidden nodes of a
   # Binomial2BernoulliRBM is the same as that of a BernoulliRBM
   b2brbm = BMs.Binomial2BernoulliRBM(rbm.weights, rbm.visbias, rbm.hidbias)
   hh2 = BMs.hiddenpotential(b2brbm, vv)
   @test sum(abs.(hh - hh2)) == 0
   BMs.hiddenpotential!(hh, b2brbm, vv)
   @test isapprox(sum(abs.(hh - hh2)), 0)

   # Test activation potential of visible nodes for BernoulliRBM
   BMs.visiblepotential!(vv, rbm, hh)
   @test isapprox(sum(abs.(vv - BMs.visiblepotential(rbm, hh))), 0)

   # Test activation potential of the visible nodes of a Binomial2BernoulliRBM
   BMs.visiblepotential!(vv, b2brbm, hh)
   @test isapprox(sum(abs.(vv- BMs.visiblepotential(b2brbm, hh))), 0)

   # Test activation potential of hidden nodes for GBRBM
   hh = rand(nsamples, nhidden)
   gbrbm = BMTest.randgbrbm(nvisible, nhidden)
   BMs.hiddenpotential!(hh, gbrbm, vv)
   @test isapprox(sum(abs.(hh - BMs.hiddenpotential(gbrbm, vv))), 0)

   # Test activation potential of visible nodes for GBRBM
   BMs.visiblepotential!(vv, gbrbm, hh)
   @test isapprox(sum(abs.(vv - BMs.visiblepotential(gbrbm, hh))), 0)

   # Test activation potential of hidden nodes for PartitionedRBM.
   # For a PartitionedRBM consisting of BernoulliRBMs, the
   # activation potential has to be the same as the activation potential
   # of the BernoulliRBM resulting from joining the RBMs.
   nvisible2 = 7
   nhidden2 = 4
   vv = rand(nsamples, nvisible + nvisible2)
   hhpartitioned = rand(nsamples, nhidden + nhidden2)
   rbm2 = BMTest.randrbm(nvisible2, nhidden2)
   prbm = BMs.PartitionedRBM{BMs.BernoulliRBM}([rbm; rbm2])
   joinedrbm = BMs.joinrbms([rbm; rbm2])
   BMs.hiddenpotential!(hhpartitioned, joinedrbm, vv)
   @test all(isapprox.(hhpartitioned, BMs.hiddenpotential(joinedrbm, vv)))

   # Test activation potential of visible nodes for PartitionedRBM
   BMs.visiblepotential!(vv, prbm, hhpartitioned)
   @test all(isapprox.(vv, BMs.visiblepotential(joinedrbm, hhpartitioned)))
end

function rbmexactloglikelihoodvsbaserate(x::Matrix{Float64}, nhidden::Int)
   a = BMs.logit.(vec(mean(x, dims = 1)))
   nvisible = length(a)
   rbm = BMs.BernoulliRBM(zeros(nvisible, nhidden), a, ones(nhidden))
   baserate = BMs.bernoulliloglikelihoodbaserate(x)
   exactloglik = BMs.loglikelihood(rbm, x, BMs.exactlogpartitionfunction(rbm))
   baserate - exactloglik
end

function bgrbmexactloglikelihoodvsbaserate(x::Matrix{Float64}, nhidden::Int)
   a = BMs.logit.(vec(mean(x, dims = 1)))
   nvisible = length(a)
   bgrbm = BMs.BernoulliGaussianRBM(zeros(nvisible, nhidden), a, ones(nhidden))
   baserate = BMs.bernoulliloglikelihoodbaserate(x)
   exactloglik = BMs.loglikelihood(bgrbm, x, BMs.exactlogpartitionfunction(bgrbm))
   baserate - exactloglik
end

function gbrbmexactloglikelihoodvsbaserate(x::Matrix{Float64}, nhidden::Int)
   a = vec(mean(x, dims = 1))
   nvisible = length(a)
   sd = vec(std(x, dims = 1))
   gbrbm = BMs.GaussianBernoulliRBM(zeros(nvisible, nhidden), a, ones(nhidden), sd)
   baserate = BMs.gaussianloglikelihoodbaserate(x)
   exactloglik = BMs.loglikelihood(gbrbm, x, BMs.exactlogpartitionfunction(gbrbm))
   baserate - exactloglik
end

function test_likelihoodbaserates()
   @test isapprox(BMs.bernoulliloglikelihoodbaserate(10),
         BMs.bernoulliloglikelihoodbaserate(rand([0.0 1.0], 10000, 10)),
         atol = 0.002)

   # Test exact loglikelihood of a baserate RBM
   # vs the calculated loglikelihoodbaserate
   x = rand([0.0 0.0 0.0 1.0 1.0], 100, 10);
   @test isapprox(BMTest.bgrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   @test isapprox(BMTest.rbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   x = rand(100, 10) + randn(100, 10);
   @test isapprox(BMTest.gbrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)

   @test isapprox(BMs.bernoulliloglikelihoodbaserate(10),
         BMs.bernoulliloglikelihoodbaserate(rand([0.0 1.0], 10000, 10)),
         atol = 0.002)

   # Test exact loglikelihood of a baserate RBM
   # vs the calculated loglikelihoodbaserate
   x = rand([0.0 0.0 0.0 1.0 1.0], 100, 10);
   @test isapprox(BMTest.bgrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   @test isapprox(BMTest.rbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)
   x = rand(100, 10) + randn(100, 10);
   @test isapprox(BMTest.gbrbmexactloglikelihoodvsbaserate(x, 10), 0, atol = 1e-10)

   nvisible = 10
   nhidden = 10
   bgrbm = BMs.BernoulliGaussianRBM(zeros(nvisible, nhidden), rand(nvisible), ones(nhidden));
   @test isapprox(BMs.logpartitionfunction(bgrbm, 0.0),
         BMs.exactlogpartitionfunction(bgrbm))
end

function test_stackrbms_preparetrainlayers()
   x = Matrix{Float64}(undef, 22, 10);
   epochs = 17
   learningrate = 0.007
   nhiddens = [7;6;5]
   learningrate1 = 0.006
   learningrate2 = 0.017
   epochs1 = 27
   epochs2 = 17
   batchsize = 5
   batchsize1 = 2
   batchsize2 = 4
   optimizer = BMs.loglikelihoodoptimizer(learningrate = learningrate)
   optimizer1 = BMs.loglikelihoodoptimizer(learningrate = learningrate1)
   optimizer2 = BMs.loglikelihoodoptimizer(learningrate = learningrate2)

   # no layer-specific instructions
   trainlayers = Vector{BMs.TrainLayer}()
   trainlayers = BMs.stackrbms_preparetrainlayers(x;
         trainlayers = trainlayers, epochs = epochs,
         learningrate = learningrate, nhiddens = nhiddens,
         batchsize = batchsize, optimizer = optimizer)
   @test length(trainlayers) == length(nhiddens)
   @test all(map(t -> t.learningrate == learningrate, trainlayers))
   @test all(map(t -> t.epochs == epochs, trainlayers))
   @test map(t-> t.nhidden, trainlayers) == nhiddens
   @test all(map(t-> t.batchsize == batchsize, trainlayers))
   @test all(map(t-> t.optimizer.learningrate == learningrate, trainlayers))

   # layer-specific instructions without partitioning
   trainlayers = [
         BMs.TrainLayer(learningrate = learningrate1, epochs = epochs1,
               batchsize = batchsize1, nhidden = nhiddens[1],
               optimizer = optimizer1);
         BMs.TrainLayer(learningrate = learningrate2, epochs = epochs2,
               batchsize = batchsize2, nhidden = nhiddens[2],
               optimizer = optimizer2);
         BMs.TrainLayer(nhidden = nhiddens[3])
   ]
   trainlayers = BMs.stackrbms_preparetrainlayers(x;
         trainlayers = trainlayers, epochs = epochs,
         learningrate = learningrate, nhiddens = Vector{Int}(),
         batchsize = batchsize, optimizer = optimizer)
   @test length(trainlayers) == length(nhiddens)
   @test trainlayers[1].learningrate == learningrate1
   @test trainlayers[2].learningrate == learningrate2
   @test trainlayers[3].learningrate == learningrate
   @test trainlayers[1].epochs == epochs1
   @test trainlayers[2].epochs == epochs2
   @test trainlayers[3].epochs == epochs
   @test map(t-> t.nhidden, trainlayers) == nhiddens
   @test trainlayers[1].batchsize == batchsize1
   @test trainlayers[2].batchsize == batchsize2
   @test trainlayers[3].batchsize == batchsize
   @test trainlayers[1].optimizer.learningrate == learningrate1
   @test trainlayers[2].optimizer.learningrate == learningrate2
   @test trainlayers[3].optimizer.learningrate == learningrate

   # wrong/misleading specification of number of visible nodes
   badtrainlayers = deepcopy(trainlayers)
   badtrainlayers[1].nvisible = 100
   @test_throws ErrorException BMs.stackrbms_preparetrainlayers(x;
         trainlayers = badtrainlayers)

   # partitioning in first layer
   trainlayers = [
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(learningrate = learningrate1, epochs = epochs1,
                  rbmtype = BMs.GaussianBernoulliRBM, nhidden = 3,
                  nvisible = 3, batchsize = batchsize2);
            BMs.TrainLayer(nhidden = 4, nvisible = 7)
         ]);
         BMs.TrainLayer(nhidden = 6);
         BMs.TrainLayer(nhidden = 5)
   ]
   trainlayers = BMs.stackrbms_preparetrainlayers(x;
         trainlayers = trainlayers,
         epochs = epochs, learningrate = learningrate,
         batchsize = batchsize1)
   @test length(trainlayers) == length(nhiddens)
   @test trainlayers[1].parts[1].learningrate == learningrate1
   @test trainlayers[1].parts[2].learningrate == learningrate
   @test trainlayers[3].learningrate == learningrate
   @test trainlayers[1].parts[1].epochs == epochs1
   @test trainlayers[1].parts[2].epochs == epochs
   @test trainlayers[2].epochs == epochs
   @test trainlayers[3].epochs == epochs
   @test [trainlayers[1].parts[1].nhidden + trainlayers[1].parts[2].nhidden;
         trainlayers[2].nhidden; trainlayers[3].nhidden] == nhiddens
   @test trainlayers[1].parts[1].batchsize == batchsize2
   @test trainlayers[1].parts[2].batchsize == batchsize1

   # incorrect partitioning must not be allowed
   badtrainlayers = deepcopy(trainlayers)
   badtrainlayers[1].parts[1].nvisible = 20
   @test_throws ErrorException BMs.stackrbms_preparetrainlayers(x;
         trainlayers = badtrainlayers)

   # too many visible nodes in layer 1
   badtrainlayers = deepcopy(trainlayers)
   badtrainlayers[1].parts[1].nvisible = 100
   @test_throws ErrorException BMs.stackrbms_preparetrainlayers(x;
         trainlayers = badtrainlayers)

   nothing
end

function test_summing_out()
   # Test exact computation of log partition function of DBM:
   # Compare inefficient simple implementation with more efficient one that
   # sums out layers analytically.
   nhiddensvecs = Array[[5;4;6;3], [5;6;7], [4;4;2;4;3], [4;3;2;2;3;2]];
   for i = 1:length(nhiddensvecs)
      dbm = BMTest.randdbm(nhiddensvecs[i]);
      @test isapprox(BMs.exactlogpartitionfunction(dbm),
            BMTest.exactlogpartitionfunctionwithoutsummingout(dbm))
   end
end

function test_oneornone_encoding()
   x = [1 2 3; 0 1 0; 3 0 1]
   ncategories = 4
   @test all(x .== BMs.oneornone_decode(BMs.oneornone_encode(x, ncategories), ncategories))
   nscategories = [4;3;4]
   @test all(x .== BMs.oneornone_decode(BMs.oneornone_encode(x, nscategories), nscategories))
end

function test_softmaxvssigm()
   nsamples = 5
   nvariables = 6
   x1 = 6. .* rand(nsamples, nvariables) .- 3.
   x2 = copy(x1)
   # softmax transfomation to each of the columns
   for i in 1:nvariables
      BMs.softmax0!(view(x1, :, i:i))
   end
   @test all(isapprox.(x1, BMs.sigm!(x2)))
end

"""
    testaisvsexact(bm, percentalloweddiff)
Tests whether the exact log partition function is approximated by the value
estimated by AIS for the given Boltzmann Machine.
The test is successful, if the difference in percent of the log of the exact value
between the exact log partition function and AIS-estimated one is less than
`percentalloweddiff`.
"""
function checkaisvsexact(bm::BMs.AbstractBM, percentalloweddiff::Float64;
      ntemperatures = 100, nparticles = 100)

   exact = BMs.exactlogpartitionfunction(bm)
   estimated = BMs.logpartitionfunction(bm;
         ntemperatures = ntemperatures, nparticles = nparticles)

   @check abs((exact - estimated)/exact) < percentalloweddiff / 100
end

function exactlogpartitionfunctionwithoutsummingout(dbm::BMs.BasicDBM)
   nlayers = length(dbm) + 1
   u = BMs.initcombination(dbm)
   z = 0.0
   while true
      z += exp(-BMs.energy(dbm, u))
      BMs.next!(u) || break
   end
   log(z)
end

function testsummingoutforexactloglikelihood(nunits::Vector{Int})
   x = BMTest.createsamples(1000, nunits[1]);
   dbm = BMs.stackrbms(x, nhiddens = nunits[2:end], batchsize = 5,
         epochs = 50, predbm = true, learningrate = 0.001);
   dbm = BMs.traindbm!(dbm, x,
         learningrates = [0.02*ones(10); 0.01*ones(10); 0.001*ones(10)],
         epochs = 30);
   logz = BMs.exactlogpartitionfunction(dbm)
   @test isapprox(BMs.exactloglikelihood(dbm, x, logz),
         BMTest.exactloglikelihoodwithoutsummingout(dbm, x, logz))
end

function exactloglikelihoodwithoutsummingout(dbm::BMs.BasicDBM, x::Array{Float64,2},
      logz = BMs.exactlogpartitionfunction(dbm))

   nsamples = size(x,1)
   nlayers = length(dbm) + 1

   u = BMs.initcombination(dbm)
   logp = 0.0
   for j = 1:nsamples
      u[1] = vec(x[j,:])

      p = 0.0
      while true
         p += exp(-BMs.energy(dbm, u))

         # next combination of hidden nodes' activations
         BMs.next!(u[2:end]) || break
      end

      logp += log(p)
   end

   logp /= nsamples
   logp -= logz
   logp
end

function test_dbmjoining()
   dbm1 = BMTest.randdbm([5;4;3])
   dbm2 = BMTest.randdbm([4;5;2])
   dbm3 = BMTest.randdbm([6;7;8])

   exactlogpartitionfunction1 = BMs.exactlogpartitionfunction(dbm1)
   exactlogpartitionfunction2 = BMs.exactlogpartitionfunction(dbm2)
   exactlogpartitionfunction3 = BMs.exactlogpartitionfunction(dbm3)

   jointdbm1 = BMs.joindbms(BMs.BasicDBM[dbm1, dbm2])
   # Test use of visibleindexes
   indexes = randperm(15)

   jointdbm2 = BMs.joindbms(BMs.BasicDBM[dbm1, dbm2, dbm3],
            [indexes[1:5], indexes[6:9], indexes[10:15]])

   @test isapprox(
               exactlogpartitionfunction1 + exactlogpartitionfunction2,
               BMs.exactlogpartitionfunction(jointdbm1))

   @test isapprox(
            exactlogpartitionfunction1 + exactlogpartitionfunction2 +
                  exactlogpartitionfunction3,
            BMs.exactlogpartitionfunction(jointdbm2))
end

"""
Tests whether the empirical loglikelihood from samples that are generated from
the exact distribution of Boltzmann machines is equal to the exact loglikelihood.
"""
function check_exactsampling()
   nsamples = 200
   x = BMs.barsandstripes(nsamples, 9)[:, 1:6]

   dbm1 = BMs.fitdbm(x; epochspretraining = 20, epochs = 15, nhiddens = [8; 5])
   BMTest.check_exactsampling(dbm1, x)

   dbm1 = BMs.fitdbm(x; epochspretraining = 20, epochs = 15, nhiddens = [6; 5; 4])
   BMTest.check_exactsampling(dbm1, x)

   rbm = BMs.fitrbm(x; epochs = 100, nhidden = 20)
   BMTest.check_exactsampling(rbm, x)
end

function check_exactsampling(bm::BMs.AbstractBM, x::Matrix{Float64};
      percentalloweddiff = 2.0, nsamples::Int = 35000)

   mbdist = BMs.MultivariateBernoulliDistribution(bm)
   exactsamples = BMs.samples(mbdist, nsamples)

   emploglikexactsamples = BMs.empiricalloglikelihood(x, exactsamples)
   exactloglik = BMs.exactloglikelihood(bm, x)
   @check abs((emploglikexactsamples - exactloglik) / exactloglik) <
         percentalloweddiff / 100
end

"""
Fits a larger RBM and a DBM and tests whether the AIS estimation
yields approximately the same results, if performed parallel or not.
"""
function check_likelihoodconsistency()
   nsamples = 1000
   nvariables = 64
   x = BMs.barsandstripes(nsamples, nvariables)
   dbm = BMs.fitdbm(x; epochspretraining = 20, epochs = 15,
         nhiddens = [50; 25; 20])

   logp1 = BMs.loglikelihood(dbm, x[1:10,:];
         parallelized = false, nparticles = 300, ntemperatures = 200)
   logp2 = BMs.loglikelihood(dbm, x[1:10,:];
         parallelized = true, nparticles = 300, ntemperatures = 200)
   @check isapprox(logp1, logp2, rtol = 0.055)

   rbm1 = BMs.fitrbm(x; epochs = 20)
   logz1 = BMs.logpartitionfunction(rbm1; parallelized = false)
   logz1_2 = BMs.logpartitionfunction(rbm1; parallelized = true)
   @check isapprox(logz1, logz1_2, rtol = 0.02)

   logp1 = BMs.loglikelihood(rbm1, x, logz1)
   rbm2 = BMs.fitrbm(x; epochs = 20, nhidden = 30)
   logz2 = BMs.logpartitionfunction(rbm2)
   logp2 = BMs.loglikelihood(rbm2, x, logz2)

   samplemeanrbm =
         BMs.BernoulliRBM(fill(0.0, nvariables, 1), vec(mean(x, dims = 1)), fill(0.0, 1))
   samplemeanrbm2 =
         BMs.BernoulliRBM(fill(0.0, size(rbm1.weights)), vec(mean(x, dims = 1)),
               fill(0.0, length(rbm1.hidbias)))

   # Annealing between two RBMs with zero weights
   @check isapprox(
         BMs.logmeanexp(BMs.aislogimpweights(samplemeanrbm, samplemeanrbm2)),
         (BMs.logpartitionfunctionzeroweights(samplemeanrbm2) -
               BMs.logpartitionfunctionzeroweights(samplemeanrbm)))

   # Annealing from RBM to RBM with sample mean and zero weights
   @check abs(logz1 - BMs.logpartitionfunctionzeroweights(samplemeanrbm) -
         BMs.logmeanexp(BMs.aislogimpweights(samplemeanrbm, rbm1))) <
               logz1 * 0.013

   # Annealing from one RBM to another RBM: Compare partition functions ...
   logimpweights = BMs.aislogimpweights(rbm1, rbm2;
         ntemperatures = 1000, nparticles = 200, burnin = 10)
   @check abs(BMs.logmeanexp(logimpweights) - (logz2 - logz1)) <
         max(logz2, logz1) * 0.03
   # ... and loglikelihood
   lldiff2 = BMs.loglikelihooddiff(rbm1, rbm2, x, logimpweights)
   lldiff1 = logp1 - logp2
   # very weak test, TODO: check and improve
   @check sign(lldiff1) == sign(lldiff2)
end

"""
Tests the likelihood estimation approaches for GaussianBernoulliRBMs
and alternative GaussianBernoulliRBMs.
"""
function check_likelihoodconsistency_gaussian(gbrbmtype::Type{GBRBM}
      ) where GBRBM <: Union{BMs.GaussianBernoulliRBM, BMs.GaussianBernoulliRBM2}

   x = BMs.curvebundles(nperbundle = 300, nvariables = 10,
         nbundles = 4)

   gbrbm1 = BMs.fitrbm(x; nhidden = 30, rbmtype = gbrbmtype,
         cdsteps = 10,
         pcd = false,
         learningrate = 0.001, epochs = 50,
         sdlearningrate = 0.00001)

   datainitrbm = GBRBM(zeros(size(x,2), 1),
         vec(mean(x, dims = 1)), [0.0], vec(std(x, dims = 1)))

   # Compare estimated difference of log partition functions to difference
   # between estimation and exact log partition function of GBRBM with zero weights
   logz1 = BMs.logpartitionfunction(gbrbm1)
   @check abs(logz1 - BMs.logpartitionfunctionzeroweights(datainitrbm) -
         BMs.logmeanexp(BMs.aislogimpweights(datainitrbm, gbrbm1,
               nparticles = 200, ntemperatures = 200))) < 0.03 * logz1

   gbrbm2 = BMs.fitrbm(x; nhidden = 15, rbmtype = gbrbmtype,
         cdsteps = 10,
         learningrate = 0.001, epochs = 100,
         sdlearningrates = [fill(0.0, 30); fill(0.00001, 300)])

   logz2 = BMs.logpartitionfunction(gbrbm2)

   logdiffimpweights = BMs.aislogimpweights(gbrbm1, gbrbm2)
   logzdiff = BMs.logmeanexp(logdiffimpweights)

   # Annealing from one RBM to another RBM: Compare partition functions ...
   @check abs(logzdiff - (logz2 - logz1)) < max(logz2, logz1) * 0.03
   # ... and loglikelihood
   logp1 = BMs.loglikelihood(gbrbm1, x, logz1)
   logp2 = BMs.loglikelihood(gbrbm2, x, logz2)
   lldiff2 = BMs.loglikelihooddiff(gbrbm1, gbrbm2, x, logdiffimpweights)
   lldiff1 = logp1 - logp2
   # very weak test, TODO: check and improve
   @check sign(lldiff1) == sign(lldiff2)
end

"
Tests whether the log-likelihood of Binomial2BernoulliRBMs is approximately
equal to the empirical loglikelihood of data generated by the models, and whether
the partition function estimated by AIS is near the exact value.
"
function check_b2brbm()
   x = BMTest.createsamples(100, 4) + BMTest.createsamples(100, 4)
   b2brbm = BMs.fitrbm(x, rbmtype = BMs.Binomial2BernoulliRBM, epochs = 150,
         nhidden = 4, learningrate = 0.001)
   @check checklikelihoodempirically(b2brbm, x; percentalloweddiff = 1.0)
end

function check_rbm()
   x = BMTest.createsamples(100, 4)
   rbm = BMs.fitrbm(x, epochs = 30,
         nhidden = 4, learningrate = 0.001)
   @check checklikelihoodempirically(rbm, x)
end

function check_sampling()
   realsamerate = 0.8
   x = BMTest.createsamples(100, 4, realsamerate)
   rbm = BMs.fitrbm(x, epochs = 400, nhidden = 5, learningrate = 0.002)

   samerate(x) = sum((x[:,1] .== x[:,2]) .& (x[:,2] .== x[:,3])) ./ size(x, 1)
   ngensamples = 5000
   @check abs(samerate(BMs.samples(rbm, ngensamples)) - realsamerate) < 0.1
   @test samerate(BMs.samples(rbm, 5, conditions = [1 => 1.0, 2 => 1.0, 3 => 1.0])) == 1.0
   @check samerate(BMs.samples(rbm, 5, samplelast = false)) == 0.0
   @test samerate(BMs.samples(rbm, 5, conditions = [1 => 0.0, 2 => 0.0, 3 => 0.0],
         samplelast = true)) == 1.0
   true
end

function check_softmax0rbm()
   samerates = [0.2; 0.3; 0.4]
   categories = length(samerates)
   x = BMs.oneornone_encode(BMTest.createsamples_categorical(100, 5, samerates), categories)
   datadict = BMs.DataDict("x" => x)
   monitor = BMs.Monitor()
   rbm = BMs.fitrbm(x, epochs = 150, rbmtype = BMs.Softmax0BernoulliRBM,
         categories = 3, nhidden = 4, learningrate = 0.01,
         monitoring = (rbm, epoch) -> begin
            if epoch % 50 == 0
               BMs.monitorreconstructionerror!(monitor, rbm, epoch, datadict)
               BMs.monitorexactloglikelihood!(monitor, rbm, epoch, datadict)
               BMs.monitorloglikelihood!(monitor, rbm, epoch, datadict)
            end
         end)

   xtest = x[1:5, :]

   # Softmax0BernoulliRBM with two categories must be equivalent to BernoulliRBM
   rbm2_softmax = BMs.Softmax0BernoulliRBM(rbm.weights, rbm.visbias, rbm.hidbias, 2)
   rbm2_bernoulli = BMs.BernoulliRBM(rbm.weights, rbm.visbias, rbm.hidbias)
   @test isapprox(BMs.exactloglikelihood(rbm2_softmax, xtest),
         BMs.exactloglikelihood(rbm2_bernoulli, xtest))

   zeroweightsrbm = BMs.Softmax0BernoulliRBM(
         zeros(size(rbm.weights)), rbm.visbias, rbm.hidbias, 3)
   @test isapprox(BMs.logpartitionfunctionzeroweights(zeroweightsrbm),
         BMs.exactlogpartitionfunction(zeroweightsrbm))

   @check BMTest.check_exactsampling(rbm, xtest; nsamples = 200000)

   # TODO evaluate quality of sample generation
   # sampled = BMs.oneornone_decode(BMs.samples(rbm, 1500), categories)
   # sampled[:, 1] .== sampled[:, 2] .== sampled[:, 3]

   # Note: Needs a little bit more burnin than the default. Problem?
   @check isapprox(BMs.exactloglikelihood(rbm, x),
         BMs.empiricalloglikelihood(rbm, x, 500000, 100), rtol = 0.03)

   # Test AIS for estimation of partition function
   logz_exact = BMs.exactlogpartitionfunction(rbm)
   logz_ais = BMs.logpartitionfunction(rbm)
   @check isapprox(logz_exact, logz_ais, rtol = 0.02)
end

function check_softmax0dbm()
   samerates = [0.2; 0.3; 0.4]
   categories = length(samerates)
   x = BMs.oneornone_encode(BMTest.createsamples_categorical(100, 5, samerates), categories)
   dbm = BMs.fitdbm(x, epochs = 150,
         pretraining = [
               BMs.TrainLayer(rbmtype = BMs.Softmax0BernoulliRBM, categories = 3,
                     nhidden = 5);
               BMs.TrainLayer(nhidden = 4)])

   xtest = x[1:5, :]

   @check BMTest.check_exactsampling(dbm, xtest; nsamples = 200000)

   # Test exact vs empirical loglikelihood
   @check isapprox(BMs.exactloglikelihood(dbm, x),
         BMs.empiricalloglikelihood(dbm, x, 50000, 150), rtol = 0.03)

   # Test AIS for estimation of partition function
   logz_exact = BMs.exactlogpartitionfunction(dbm)
   logz_ais = BMs.logpartitionfunction(dbm)
   @check isapprox(logz_exact, logz_ais, rtol = 0.02)

   # Test AIS for likelihood estimation vs exact calculation
   @check isapprox(BMs.loglikelihood(dbm, xtest; nparticles = 200, ntemperatures = 200),
         BMs.exactloglikelihood(dbm, xtest), rtol = 0.05)
end

function check_gbrbm(gbrbmtype::Type{GBRBM}
      ) where GBRBM <:Union{BMs.GaussianBernoulliRBM, BMs.GaussianBernoulliRBM2}

   x = BMs.curvebundles(nperbundle = 300, nvariables = 10,
         nbundles = 4)
   gbrbm = BMs.fitrbm(x; nhidden = 15, rbmtype = gbrbmtype,
         cdsteps = 10,
         pcd = false,
         learningrate = 0.001, epochs = 50,
         sdlearningrate = 0.00001)

   logz = BMs.logpartitionfunction(gbrbm; parallelized = true,
         ntemperatures = 500, nparticles = 500)
   estloglik = BMs.loglikelihood(gbrbm, x, logz)
   exactloglik = BMs.exactloglikelihood(gbrbm, x)
   @check abs((exactloglik - estloglik) / exactloglik) < 2.5 / 100
end

function test_intensities()
   x = randn(5, 3) .* [1 2 3]
   transformed, transformation = BMs.intensities_encode(x)
   @test isapprox(BMs.intensities_decode(transformed, transformation), x)
end

function test_rbm_monitoring(gbrbmtype::Type{GBRBM}
      ) where GBRBM <:Union{BMs.GaussianBernoulliRBM, BMs.GaussianBernoulliRBM2}

   x = BMs.curvebundles(nperbundle = 20, nvariables = 5,
      nbundles = 4)

   # Test some more options and monitoring
   monitor = BMs.Monitor()
   datadict = BMs.DataDict("Small subset" => x[1:3, :])
   gbrbm = BMs.fitrbm(x; nhidden = 5, rbmtype = gbrbmtype,
         learningrate = 0.000001, epochs = 20,
         batchsize = 10,
         pcd = false,
         sdlearningrate = 0.000000001,
         monitoring = (rbm, epoch) -> begin
            if epoch == 10 || epoch == 20
               BoltzmannMachines.monitorexactloglikelihood!(monitor, rbm, epoch, datadict)
               BoltzmannMachines.monitorreconstructionerror!(monitor, rbm, epoch, datadict)
            end
         end)
   nothing
end

function test_mdbm_monitoring()
   x1 = BMTest.createsamples(80, 4)
   x2 = BMs.curvebundles(nperbundle = 20, nvariables = 5, nbundles = 4)
   data = hcat(x1, x2)
   traindata, testdata = BMs.splitdata(data, 0.2)
   monitor = BMs.Monitor(); monitor1 = BMs.Monitor()
   datadict = BMs.DataDict("Training data" => traindata, "Test data" => testdata)
   dbm = BMs.fitdbm(traindata,
      epochs = 2,
      monitoringdatapretraining = datadict,
      batchsizepretraining = 10,
      pretraining = [
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(nvisible = 4, nhidden = 5,
                  monitoring = (rbm, epoch, datadict) ->
                        BMs.monitorreconstructionerror!(monitor1, rbm, epoch, datadict))
            BMs.TrainLayer(nhidden = 6, rbmtype = BMs.GaussianBernoulliRBM2)
         ]),
         BMs.TrainLayer(nhidden = 5)
      ],
      monitoring = (dbm, epoch) -> BMs.monitorlogproblowerbound!(monitor, dbm, epoch, datadict))
   nothing
end

function test_monitored_fitrbm()
   nvisible1 = 4
   x = BMTest.createsamples(5, nvisible1)

   # without monitoring:
   # Test whether monitored_fitrbm and fitrbm do the same thing
   randomseed = abs(rand(Int))
   Random.seed!(randomseed)
   monitor, rbm1 = BMs.monitored_fitrbm(x;
         pcd = false, cdsteps = 3, nhidden = 5, epochs = 2)
   @test isempty(monitor)
   Random.seed!(randomseed)
   rbm2 = BMs.fitrbm(x; pcd = false, cdsteps = 3, nhidden = 5, epochs = 2)
   @test all(isapprox.(rbm1.weights, rbm2.weights))

   # one monitoring function
   monitor, rbm1 = BMs.monitored_fitrbm(x;
         monitoring = BMs.monitorreconstructionerror!, batchsize = 4, epochs = 3)
   @test length(monitor) == 3

   # two monitoring functions
   # (also test uneven batch size)
   monitor, rbm1 = BMs.monitored_fitrbm(x; epochs = 3, batchsize = 3,
         monitoring = [BMs.monitorfreeenergy!;
               BMs.monitorexactloglikelihood!])
   @test length(monitor) == 6

   nothing
end

function test_monitored_stackrbms()
   nvisible1 = 4
   nvisible2 = 3
   x = BMTest.createsamples(20, nvisible1)
   xpart = hcat(x, rand(20, nvisible2))
   nhiddens = [4; 3]

   # Test whether monitored_stackrbm and stackrbm do the same thing
   randomseed = abs(rand(Int))
   Random.seed!(randomseed)
   monitors1, dbm1 = BMs.monitored_stackrbms(x; predbm = true, nhiddens = nhiddens)
   @test length(monitors1) == 2
   @test all(isempty.(monitors1))
   Random.seed!(randomseed)
   dbm2 = BMs.stackrbms(x; predbm = true, nhiddens = nhiddens)
   @test all([all(isapprox.(dbm1[i].weights, dbm2[i].weights)) for i in eachindex(dbm1)])

   # partitioned
   monitor1 = BMs.Monitor()
   trainlayers = [
      BMs.TrainPartitionedLayer([
         BMs.TrainLayer(nvisible = nvisible1, nhidden = 2)
         BMs.TrainLayer(nhidden = 2, rbmtype = BMs.GaussianBernoulliRBM2)
      ]),
      BMs.TrainLayer(nhidden = 4)
   ]
   epochs = 5
   monitors1, dbm1 = BMs.monitored_stackrbms(xpart;
         trainlayers = trainlayers,
         monitoring = BMs.monitorexactloglikelihood!,
         predbm = true, epochs = epochs)
   @test length(monitors1) == 2
   @test length(monitors1[1]) == 2 # layers with two partitions
   @test length(monitors1[1][1]) == epochs
   @test length(monitors1[1][2]) == epochs
   @test length(monitors1[2]) == epochs

   # Vector of monitoring functions
   monitors, dbn = BMs.monitored_stackrbms(x;
         nhiddens = nhiddens,
         monitoring = [BMs.monitorexactloglikelihood!, BMs.monitorreconstructionerror!],
         predbm = false, epochs = epochs)
   @test length(monitors) == length(nhiddens) == length(dbn)
   @test all(map(length, monitors) .== 2*epochs)

   nothing
end

function dbms_are_approx_equal(dbm1, dbm2)
   all([all(isapprox.(dbm1[i].weights, dbm2[i].weights)) for i in eachindex(dbm1)])
end

function test_monitored_fitdbm()
   nvisible1 = 4
   nvisible2 = 3
   x = BMTest.createsamples(20, nvisible1)
   xpart = hcat(x, rand(20, nvisible2))

   kwargs = Dict{Symbol, Any}(:epochs=> 5, :epochspretraining => 10,
         :nhiddens => [4; 3], :batchsizepretraining => 5,
         :learningrate => 0.001, :nparticles => 17)
   # Test whether monitored_fitdbm and fitdbm do the same thing
   randomseed = abs(rand(Int))
   Random.seed!(randomseed)
   monitors1, dbm1 = BMs.monitored_fitdbm(x; kwargs...)
   @test length(monitors1) == 3
   @test all(isempty.(monitors1))
   Random.seed!(randomseed)
   dbm2 = BMs.fitdbm(x; kwargs...)
   @test dbms_are_approx_equal(dbm1, dbm2)

   # ensure real types are not unified
   BMs.monitored_fitdbm(x; epochs = 1, learningrate = 0.05)

   # partitioned
   monitoringdata = BMs.DataDict("Only part" => BMs.splitdata(xpart, 0.5)[1])
   epochspretraining = 4
   epochs = 5
   trainlayers = [
      BMs.TrainPartitionedLayer([
         BMs.TrainLayer(nvisible = nvisible1, nhidden = 2)
         BMs.TrainLayer(nhidden = 2, rbmtype = BMs.GaussianBernoulliRBM2)
      ]),
      BMs.TrainLayer(nhidden = 4)
   ]
   monitors, dbm = BMs.monitored_fitdbm(xpart;
         pretraining = trainlayers,
         monitoringpretraining = BMs.monitorfreeenergy!,
         monitoring = [BMs.monitorexactloglikelihood!,
               BMs.monitorlogproblowerbound!],
         epochs = epochs,
         epochspretraining = epochspretraining)
   @test length(monitors) == 3
   @test length(monitors[1]) == 2 # layers with two partitions
   @test length(monitors[1][1]) == epochspretraining
   @test length(monitors[1][2]) == epochspretraining
   @test length(monitors[2]) == epochspretraining
   lowerbounditems = filter(m -> m.evaluation == BMs.monitorlogproblowerbound, monitors[3])
   exactloglikitems = filter(m -> m.evaluation == BMs.monitorexactloglikelihood, monitors[3])
   @test length(lowerbounditems) == epochs
   @test length(exactloglikitems) == epochs

   # passing an unknown argument must throw an exception
   @test_throws ErrorException BMs.monitored_fitdbm(x; bla = 10, epochs = 5)

   nothing
end

function test_batchsize_dbm_training()
   x = BMTest.createsamples(20, 4)

   randomseed = abs(rand(Int))
   # Test correct defaults for batchsize
   Random.seed!(randomseed)
   dbm1 = BMs.fitdbm(x; epochs = 2)
   Random.seed!(randomseed)
   _, dbm2 = BMs.monitored_fitdbm(x; epochs = 2, batchsize = 1, batchsizefinetuning = size(x,1))
   @test BMTest.dbms_are_approx_equal(dbm1, dbm2)

   # Test specifying different batchsizefinetuning
   Random.seed!(randomseed)
   dbm2 = BMs.fitdbm(x; epochs = 2, batchsizefinetuning = 2)
   Random.seed!(randomseed)
   _, dbm3 = BMs.monitored_fitdbm(x; epochs = 2, batchsize = 1, batchsizefinetuning = 2)
   @test BMTest.dbms_are_approx_equal(dbm2, dbm3)
   @test !BMTest.dbms_are_approx_equal(dbm1, dbm3)

   # Test that batchsize argument determines both fine-tuning and pretraining
   Random.seed!(randomseed)
   dbm4 = BMs.fitdbm(x; epochs = 2, batchsize = 2)
   Random.seed!(randomseed)
   _, dbm5 = BMs.monitored_fitdbm(x; epochs = 2, batchsizepretraining = 2,  batchsizefinetuning = 2)
   @test BMTest.dbms_are_approx_equal(dbm4, dbm5)
   @test !BMTest.dbms_are_approx_equal(dbm3, dbm5)

   nothing
end

function test_epochs_dbm_training()
   x = BMTest.createsamples(20, 4)

   randomseed = abs(rand(Int))
   # Test correct defaults for epochs
   Random.seed!(randomseed)
   dbm1 = BMs.fitdbm(x; epochs = 2)
   Random.seed!(randomseed)
   _, dbm2 = BMs.monitored_fitdbm(x; epochspretraining = 2, epochsfinetuning = 2)
   @test BMTest.dbms_are_approx_equal(dbm1, dbm2)

   # Test specifying different epochs for finetuning
   Random.seed!(randomseed)
   dbm2 = BMs.fitdbm(x; epochs = 1, epochsfinetuning = 2)
   Random.seed!(randomseed)
   _, dbm3 = BMs.monitored_fitdbm(x; epochspretraining = 1, epochsfinetuning = 2)
   @test BMTest.dbms_are_approx_equal(dbm2, dbm3)
   @test !BMTest.dbms_are_approx_equal(dbm1, dbm3)

   nothing
end

function test_learningrates_dbm_training()
   x = BMTest.createsamples(20, 4)

   randomseed = abs(rand(Int))

   Random.seed!(randomseed)
   dbm1 = BMs.fitdbm(x; learningratefinetuning = 0.1, learningratepretraining = 0.2, epochs = 1)
   Random.seed!(randomseed)
   # arguments learningratefinetuning and learningratepretraining work
   _, dbm2 = BMs.monitored_fitdbm(x; learningratefinetuning = 0.1, learningratepretraining = 0.2,
         epochs = 1)
   @test BMTest.dbms_are_approx_equal(dbm1, dbm2)

   Random.seed!(randomseed)
   # learningratefinetuning used as starting value for fine tuning learning rates
   _, dbm3 = BMs.monitored_fitdbm(x; learningratesfinetuning = [0.1], learningratepretraining = 0.2,
         epochs = 1)
   @test BMTest.dbms_are_approx_equal(dbm1, dbm3)

   nothing
end

function test_top2latentdims()
   nsamples = 10
   x = BMTest.createsamples(nsamples, 4)

   dbm = BMs.fitdbm(x, nhiddens = [3, 2], epochs = 1)
   dimred = BMs.top2latentdims(dbm, x)
   @test size(dimred) == (nsamples, 2)

   dbm = BMs.fitdbm(x, nhiddens = [4, 4], epochs = 1)
   dimred = BMs.top2latentdims(dbm, x)
   @test size(dimred) == (nsamples, 2)

   dbm = BMs.fitdbm(x, nhiddens = [2, 1], epochs = 1)
   dimred = BMs.top2latentdims(dbm, x)
   @test size(dimred) == (nsamples, 2)

   nothing
end

function check_mdbm_rbm_b2brbm()
   nsamples = 100
   nvariables = 4

   x1 = BMTest.createsamples(nsamples, nvariables) +
         BMTest.createsamples(nsamples, nvariables);
   x2 = BMTest.createsamples(nsamples, nvariables);
   x = hcat(x1, x2);
   trainlayers1 = [
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(rbmtype = BMs.Binomial2BernoulliRBM,
                  nvisible = nvariables, nhidden = nvariables - 1,
                  learningrate = 0.0015);
            BMs.TrainLayer(rbmtype = BMs.BernoulliRBM,
                  nvisible = nvariables, nhidden = nvariables)
         ]);
         BMs.TrainLayer(nhidden = nvariables);
         BMs.TrainLayer(nhidden = 3)]

   dbm1 = BMs.fitdbm(x, epochs = 15,
         epochspretraining = 15,
         learningratepretraining = 0.001,
         pretraining = trainlayers1)

   @check BMTest.checklikelihoodempirically(dbm1, x; percentalloweddiff = 5.5,
         ntemperatures = 300, nparticles = 300)

   # partitioned second layer
   trainlayers2 = [
         trainlayers1[1],
         BMs.TrainPartitionedLayer([
            BMs.TrainLayer(nhidden = 3);
            BMs.TrainLayer(nhidden = 2)
         ]),
         BMs.TrainLayer(nhidden = 4),
         BMs.TrainLayer(nhidden = 3)
   ]

   dbm2 = BMs.fitdbm(x, pretraining = trainlayers2)
   @check BMTest.checklikelihoodempirically(dbm2, x; percentalloweddiff = 6.3,
         ntemperatures = 300, nparticles = 300)
end

function checklikelihoodempirically(rbm::BMs.AbstractRBM, x::Matrix{Float64};
      percentalloweddiff = 0.5, ntemperatures::Int = 100, nparticles::Int = 100)

   logz = BMs.logpartitionfunction(rbm; parallelized = true,
      nparticles = nparticles, ntemperatures = ntemperatures)
   estloglik = BMs.loglikelihood(rbm, x, logz)
   exactloglik = BMs.exactloglikelihood(rbm, x)
   @check abs((exactloglik - estloglik) / exactloglik) < percentalloweddiff / 100

   emploglik = BMs.empiricalloglikelihood(rbm, x, 1000000)
   @check abs((exactloglik - emploglik) / exactloglik) < percentalloweddiff / 100
end

function checklikelihoodempirically(dbm::BMs.MultimodalDBM, x::Matrix{Float64};
      percentalloweddiff = 0.5, ntemperatures::Int = 100, nparticles::Int = 100)

   emploglik = BMs.empiricalloglikelihood(dbm, x, 1000000)
   estloglik = BMs.loglikelihood(dbm, x; parallelized = true,
         ntemperatures = ntemperatures, nparticles = nparticles)
   exactloglik = BMs.exactloglikelihood(dbm, x)
   @check abs((exactloglik - emploglik) / exactloglik) < percentalloweddiff / 100
   @check abs((exactloglik - estloglik) / exactloglik) < percentalloweddiff / 100
end

"""
Test DBMs with Gaussian visible nodes.
"""
function check_mdbm_gaussianvisibles()

   x = convert(Matrix{Float64}, readdlm(
         joinpath(dirname(pathof(BoltzmannMachines)), "..", "test/data/iris.csv"), ',',
         header = true)[1][:,1:4]);

   datadict = BMs.DataDict("x" => x)

   monitor1 = BMs.Monitor()
   trainlayers = [
         BMs.TrainLayer(rbmtype = BMs.GaussianBernoulliRBM,
               nhidden = 4,
               sdlearningrate = 0.000001,
               monitoring = (rbm, epoch) -> begin
                  BMs.monitorexactloglikelihood!(monitor1, rbm, epoch, datadict)
               end);
         BMs.TrainLayer(nhidden = 4);
         BMs.TrainLayer(nhidden = 4)]

   learningrates = [0.02*ones(10); 0.01*ones(10); 0.001*ones(10)]

   seed = round(Int, rand()*typemax(Int), RoundDown)
   Random.seed!(seed)
   dbm1 = BMs.stackrbms(x, epochs = 20, predbm = true, learningrate = 0.001,
         trainlayers = trainlayers)
   # BoltzmannMachinesPlots.plotevaluation(monitor1, BMs.monitorexactloglikelihood)
   dbm1 = BMs.traindbm!(dbm1, x,
         learningrates = learningrates,
         epochs = 30);

   Random.seed!(seed)
   dbm2 = BMs.fitdbm(x, epochs = 30,
         epochspretraining = 20,
         learningratepretraining = 0.001,
         pretraining = trainlayers,
         learningratesfinetuning = learningrates)

   # first and second dbm must be equal
   @test length(dbm1) == length(dbm2)
   @test isapprox(dbm1[1].sd, dbm2[1].sd)
   for i in 1:length(dbm1)
      @test isapprox(dbm1[i].weights, dbm2[i].weights)
      @test isapprox(dbm1[i].visbias, dbm2[i].visbias)
      @test isapprox(dbm1[i].hidbias, dbm2[i].hidbias)
   end

   # Test AIS
   @check BMTest.checkaisvsexact(dbm1, 0.6)

   # Test exact likelihood vs estimated likelihood
   estloglik = BMs.loglikelihood(dbm1, x; nparticles = 400, ntemperatures = 200)
   exactloglik = BMs.exactloglikelihood(dbm1, x)
   @check abs((exactloglik - estloglik) / exactloglik) < 3.5 / 100
end

function check_softmaxsampling()
   probs = [0.5 0.3 0.7 0.2; 0.2 0.25 0.1 0.15]
   sampled = zeros(size(probs))
   s0brbm = BMs.Softmax0BernoulliRBM(zeros(4,4), zeros(4), zeros(2), 3)
   nsampled = 1000
   for i = 1:nsampled
      newsamples = BMs.samplevisiblepotential!(copy(probs), s0brbm)
      @test all(map(s -> s in [1.0 0.0], newsamples))
      sampled .+= newsamples
   end
   @check all(1.12 .> (sampled ./ nsampled ./ probs) .> 0.88)
end

function test_mdbm_architectures()
   nsamples = 5
   data = hcat(float.(rand(nsamples, 3) .< 0.5), randn(nsamples, 3))

   BMs.fitdbm(data; epochs = 2, pretraining = [
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.GaussianBernoulliRBM)]);
               BMs.TrainLayer(nhidden = 4);
               BMs.TrainLayer(nhidden = 3)])

   BMs.fitdbm(data; epochs = 2, pretraining = [
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.GaussianBernoulliRBM2)]);
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nhidden = 2)
               BMs.TrainLayer(nhidden = 2)]);
         BMs.TrainLayer(nhidden = 4)])

   data2 = hcat(data,
         float.(rand(nsamples, 3) .< 0.5) .+ float.(rand(nsamples, 3) .< 0.5))
   BMs.fitdbm(data2; epochs = 2, pretraining = [
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.GaussianBernoulliRBM2);
               BMs.TrainLayer(nvisible = 3, nhidden = 3,
                     rbmtype = BMs.Binomial2BernoulliRBM)]);
         BMs.TrainPartitionedLayer([
               BMs.TrainLayer(nvisible = 3, nhidden = 3)
               BMs.TrainLayer(nvisible = 6, nhidden = 3)]);
         BMs.TrainLayer(nhidden = 3)])
   nothing
end

function test_unevenbatches()
   x = BMs.barsandstripes(20, 4)
   rbm = BMs.fitrbm(x, batchsize = 3, epochs = 2) # 20 % 3 !== 0
   @test_throws ErrorException BMs.fitrbm(x, batchsize = 21, epochs = 2)
   nothing
end

function test_beam()
   x = randn(100, 10)
   rbm = BMs.fitrbm(x, rbmtype = BMs.GaussianBernoulliRBM2,
         batchsize = 20,
         cdsteps = 50,
         nhidden = 20, epochs = 30,
         sdinitfactor = 0.1,
         optimizers = [
               fill(BMs.loglikelihoodoptimizer(learningrate = 0.0001,
                     sdlearningrate = 0.000001), 10);
               fill(BMs.beamoptimizer(learningrate = 0.0001,
                     sdlearningrate = 0.000001,
                     adversarialweight = 0.5), 20)
         ])
   nothing
end

end # of module BMTest


text:
function computemeasure(;
    MeasureType::Val = Val(PARAM[:measuretype]),
    regressionmethod=REGRESSION_METHOD[PARAM[:regressionmethod]],
    testdate::Date = PARAM[:testdate],
    includecomp::Bool = PARAM[:compincludecomp]
  )

  #acquire the panel. be careful about short-circuiting refresh panel and changing progrma parameters
  panel = formpanel()

  #debugmsg("dbgfinal", panel)
  #throw("stop")

  (includecomp == (:gvkey ∈ propertynames(panel))) || error("includecomp = $includecomp
  yet :gvkey ∈ panel.propertynames = $(:gvkey ∈ propertynames(panel))!!! This implies that
  crsp or crsp and comp must be refreshed.")

  #this assigns the primary weights and sets parameters for estimating the measure
  #see Measures.jl for details
  (panel, ms::MeasureSpec) = measure!(panel, MeasureType)

  #used in calculations for the controls
  formmarketweights!(panel)

  #get the strategy returns for post estimation purposes
  returns = formreturns(panel, ms)
  funds = prepfunds(returns, ms)

  #form the controls seperate from the focal characteristics
  formcontrols!(panel, ms)

  #estimate the measure
  regressionmethod(panel, ms; returns, funds)

  @info "Measure $(PARAM[:refreshmeasure] ? "" : "NOT") computed.
    Panel size: $(size(panel)), uniquepermno: $(unique(panel.permno) |> length)"

  testregressioninputs(panel)

  return nothing
end

function testregressioninputs(panel::DataFrame)
  if PARAM[:testmeasureinputs]

    #time series test output
    testpermno::Int = PARAM[:testpermnomult]
    testview::SubDataFrame = view(panel, (panel.permno .% testpermno) .== 0, :)
    if size(testview,1)>20_000 #prevents scenarios where the file is larger than needed
      testview=view(testview, 1:20_000,:)
    end
    testview |> CSV.write("$(PARAM[:testpath])\\$(
      PARAM[:crspdaily])_$(PARAM[:crspfrequencysuffix])_measuredat.csv")

    #cross-sectional test output
    testdate::Date = PARAM[:testdate]
    testview = view(panel, panel.date .== testdate,:)
    testview |> CSV.write("$(PARAM[:testpath])\\$(PARAM[:crspdaily])_$(
        PARAM[:crspfrequencysuffix])_measuredat_$testdate.csv")
  end

  #special case where we replace the actuald ata with simulated data
  if PARAM[:measuretype] == :testsimulated
    @warn "Estimated simulated data only!!!"
    panel |> CSV.write("$(PARAM[:testpath])\\completesimulated_measuredat_$testdate.csv")
  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# ---
# title: 1408. String Matching in an Array
# id: problem1408
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/string-matching-in-an-array/description/>
# hidden: true
# ---
#
# Given an array of string `words`. Return all strings in `words` which is
# substring of another word in **any** order.
#
# String `words[i]` is substring of `words[j]`, if can be obtained removing some
# characters to left and/or right side of `words[j]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["mass","as","hero","superhero"]
#     Output: ["as","hero"]
#     Explanation: "as" is substring of "mass" and "hero" is substring of "superhero".
#     ["hero","as"] is also a valid answer.
#
#
# **Example 2:**
#
#
#
#     Input: words = ["leetcode","et","code"]
#     Output: ["et","code"]
#     Explanation: "et", "code" are substring of "leetcode".
#
#
# **Example 3:**
#
#
#
#     Input: words = ["blue","green","bu"]
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 100`
#   * `1 <= words[i].length <= 30`
#   * `words[i]` contains only lowercase English letters.
#   * It's **guaranteed**  that `words[i]` will be unique.
#
#
## @lc code=start
using LeetCode

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


text:
using HDF5
using Test
using Pkg

println("HDF5 version ", HDF5.h5_get_libversion())

include("plain.jl")
include("readremote.jl")
include("extend_test.jl")
include("gc.jl")
include("external.jl")
include("swmr.jl")
if get(Pkg.installed(), "MPI", nothing) !== nothing
  # basic MPI tests, for actual parallel tests we need to run in MPI mode
  include("mpio.jl")
end


text:
"""
```
residuename(residue::String)
```

Function to return the long residue name from other residue codes. The function is case-insensitive.

### Examples

```julia-repl
julia> residuename("A")
"Alanine"

julia> residuename("Glu")
"Glutamic Acid"

```

"""
function residuename(residue::String)
    l = length(residue)
    code = uppercase(residue)
    if l == 1
        ires = findfirst(r -> r.one_letter_code == code, natural_aminoacids)
    elseif l == 3
        ires = findfirst(r -> r.three_letter_code == code, natural_aminoacids)
    else
        ires = findfirst(r -> uppercase(r.name) == code, natural_aminoacids)
    end
    code = (ires == nothing ? nothing : natural_aminoacids[ires].name)
    return code
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
# f(x) = x^2 - 2x + 2

using LaTeXStrings, Measures
gr()

xf(x) = x^2 - 2x + 2

Plots.plot(xf,
    -3, 3,
    marker=:circle,
    c=:darkorange3,
    linewidth= 2,
    xlabel=L"n",
    ylabel=L"F_n",
    annotation=(2, 10, L"F_n = x^2 - 2x + 2"))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module AcousticFWI

# Julia modules
using Seismic

# Inversion
export FWI
include("Inversion/FWI.jl")

# Modeling
export Helmholtz,
Source
include("Modeling/Helmholtz.jl")
include("Modeling/Source.jl")

# Operators
export Attenuation,
Laplacian,
MassMatrix,
Restriction,
SensitivityMultiShot
include("Operators/Attenuation.jl")
include("Operators/Laplacian.jl")
include("Operators/MassMatrix.jl")
include("Operators/Restriction.jl")
include("Operators/SensitivityMultiShot.jl")

end