text:
"""
    aislogimpweights(rbm; ...)
Computes the logarithmised importance weights for estimating the ratio of the
partition functions of the given `rbm` to the RBM with zero weights,
but same visible and hidden bias as the `rbm`.
This function implements the Annealed Importance Sampling algorithm (AIS)
like described in section 4.1.3 of [Salakhutdinov, 2008].

# Optional keyword arguments (for all types of Boltzmann Machines):
* `ntemperatures`: Number of temperatures for annealing from the starting model
  to the target model, defaults to 100
* `temperatures`: Vector of temperatures. By default `ntemperatures` ascending
  numbers, equally spaced from 0.0 to 1.0
* `nparticles`: Number of parallel chains and calculated weights, defaults to
   100
* `burnin`: Number of steps to sample for the Gibbs transition between models
"""
function aislogimpweights(rbm::AbstractXBernoulliRBM;
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   logimpweights = zeros(nparticles)
   mixrbm = deepcopy(rbm)

   # start with samples from model with zero weights
   hh = repeat(rbm.hidbias', nparticles)
   sigm_bernoulli!(hh)

   vv = Matrix{Float64}(undef, nparticles, length(rbm.visbias))

   for k = 2:length(temperatures)
      logimpweights .+= unnormalizedproblogratios(rbm, hh,
            temperatures[k], temperatures[k-1])

      # Gibbs transition
      copyannealed!(mixrbm, rbm, temperatures[k])
      for burn = 1:burnin
         samplevisible!(vv, mixrbm, hh)
         samplehidden!(hh, mixrbm, vv)
      end
   end

   logimpweights
end

"""
    aislogimpweights(rbm1, rbm2; ...)
Computes the logarithmised importance weights for estimating the log-ratio
log(Z2/Z1) for the partition functions Z1 and Z2 of `rbm1` and `rbm2`, respectively.
Implements the procedure described in section 4.1.2 of [Salakhutdinov, 2008].
This requires that `rbm1` and `rbm2` are of the same type and have the same
number of visible units.
"""
function aislogimpweights(rbm1::R, rbm2::R;
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5,
      initburnin::Int = 30) where {R <: AbstractRBM}

   if length(rbm2.visbias) != length(rbm1.visbias)
      error("The two RBMs must have the same numer of visible units.")
   end

   logimpweights = zeros(nparticles)

   vv = BMs.sampleparticles(rbm1, nparticles, initburnin)[1]
   nhidden1 = length(rbm1.hidbias)
   nhidden2 = length(rbm2.hidbias)
   hh = Matrix{Float64}(undef, nparticles, nhidden1 + nhidden2)

   prevmixedrbm = BMs.mixedrbm(rbm1, rbm2, temperatures[1])

   for k = 2:length(temperatures)
      curmixedrbm = BMs.mixedrbm(rbm1, rbm2, temperatures[k])

      # log(p*(x)) = -freeenergy(x)
      logimpweights .+= BMs.freeenergydiffs(prevmixedrbm, curmixedrbm, vv)

      # Gibbs transition for the (visible) nodes
      for burn = 1:burnin
         BMs.samplehidden!(hh, curmixedrbm, vv)
         BMs.samplevisible!(vv, curmixedrbm, hh)
      end

      prevmixedrbm = curmixedrbm
   end

   # account for different model size
   logimpweights .+= log(2) * (nhidden2 - nhidden1)

   logimpweights
end

"""
    aislogimpweights(dbm; ...)
Computes the logarithmised importance weights in the Annealed Importance Sampling
algorithm (AIS) for estimating the ratio of the partition functions of the given
DBM `dbm` to the base-rate DBM with all weights being zero and all biases equal
to the biases of the `dbm`.

Implements algorithm 4 in [Salakhutdinov+Hinton, 2012].
For DBMs with Bernoulli-distributed nodes only
(i. e. here DBMs of type `PartitionedBernoulliDBM`),
it is possible to calculate the importance weights by summing out either
the even layers (h1, h3, ...) or the odd layers (v, h2, h4, ...).
In the first case, the nodes' activations in the odd layers are used to
calculate the probability ratios, in the second case the even layer are used.
If `dbm` is of type `PartitionedBernoulliDBM`, the optional keyword argument
`sumout` can be used to choose by specifying the values `:odd` (default) or
`:even`.
In the case of `MultimodalDBM`s, it is not possible to choose and
the second case applies there.
"""
function aislogimpweights(dbm::PartitionedBernoulliDBM;
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5,
      sumout::Symbol = :even)

   if sumout == :odd
      # summing out even layers is done in the implementation for
      # MultimodalDBMs
      return invoke(aislogimpweights, Tuple{MultimodalDBM, }, dbm;
            ntemperatures = ntemperatures,
            temperatures = temperatures,
            nparticles = nparticles,
            burnin = burnin)
   elseif sumout != :even
      error("Invalid value for argument `sumout`.")
   end
   # in the following implementation, the even layers are summed out

   logimpweights = zeros(nparticles)
      # Todo: sample from null model, which has changed
   particles = initparticles(dbm, nparticles, biased = true)
   nlayers = length(particles)

   # for performance reasons: preallocate input and combine biases
   input1 = newparticleslike(particles)
   input2 = newparticleslike(particles)
   biases = combinedbiases(dbm)
   mixdbm = deepcopy(dbm)

   for k = 2:length(temperatures)
      # Calculate probability ratios for importance weights
      # according to activation of odd layers (v, h2, h4, ...)
      aisupdatelogimpweights!(logimpweights, input1, input2,
            temperatures[k], temperatures[k-1], dbm, biases, particles)

      # Gibbs transition
      copyannealed!(mixdbm, dbm, temperatures[k])
      gibbssample!(particles, mixdbm, burnin)
   end

   logimpweights
end

function aislogimpweights(mdbm::MultimodalDBM;
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   if length(mdbm) == 1
      # DBM with only one layer --> RBM
      return aislogimpweights(mdbm[1];
            ntemperatures = ntemperatures,
            temperatures = temperatures,
            nparticles = nparticles,
            burnin = burnin)
   end

   logimpweights = zeros(nparticles)
      # Todo: sample from null model, which has changed
   particles = initparticles(mdbm, nparticles, biased = true)
   nlayers = length(particles)

   # for performance reasons: preallocate input and combine biases
   hiddbminput1 = newparticleslike(particles[2:end])
   hiddbminput2 = newparticleslike(particles[2:end])

   mixdbm = deepcopy(mdbm)
   hiddbm = converttopartitionedbernoullidbm(mdbm[2:end])
   hidbiases = combinedbiases(hiddbm)

   for k = 2:length(temperatures)
      # Calculate probability ratios for importance weights
      # according to activation of odd hidden layers (h1, h3, ...)
      logimpweights += unnormalizedproblogratios(mdbm[1],
            particles[2], temperatures[k], temperatures[k-1])
      aisupdatelogimpweights!(logimpweights, hiddbminput1, hiddbminput2,
            temperatures[k], temperatures[k-1], hiddbm, hidbiases, particles[2:end])

      # Gibbs transition
      copyannealed!(mixdbm, mdbm, temperatures[k])
      gibbssample!(particles, mixdbm, burnin)
   end

   logimpweights
end

"""
    aisprecision(logr, aissd, sdrange)
Returns the differences of the estimated logratio `r` to the lower
and upper bound of the range defined by the multiple `sdrange`
of the standard deviation of the ratio's estimator `aissd`.
"""
function aisprecision(logr::Float64, aissd::Float64, sdrange::Float64 = 1.0)

   t = sdrange * aissd * exp(-logr)

   if 1 - t <= 0 # prevent domainerror
      diffbottom = -Inf
   else
      diffbottom = log(1 - t)
   end

   difftop = log(1 + t)

   diffbottom, difftop
end

"""
    aisprecision(logimpweights, sdrange)
"""
function aisprecision(logimpweights::Array{Float64,1}, sdrange::Float64 = 1.0)
   aisprecision(
         logmeanexp(logimpweights),
         aisstandarddeviation(logimpweights),
         sdrange)
end

"""
Computes the standard deviation of the AIS estimator (not logarithmised)
(eq 4.10 in [Salakhutdinov+Hinton, 2012]) given the logarithmised
importance weights.
"""
function aisstandarddeviation(logimpweights::Array{Float64,1})
   # Explanation: var(exp(x)) =
   #    = mean(exp(x)^3) - (mean(exp(x))^2)
   #    = exp(logmeanexp(2x)) - exp(2*logmeanexp(x))

   varimpweights =
         exp(logmeanexp(2 * logimpweights)) -
         exp(2 * logmeanexp(logimpweights))
   aissd = sqrt(varimpweights / length(logimpweights))
end

"""
Updates the logarithmized importance weights `logimpweights` in AIS
by adding the log ratio of unnormalized probabilities of the states
of the odd layers in the PartitionedBernoulliDBM `dbm`.
The activation states of the DBM's nodes are given by the `particles`.
For performance reasons, the biases are specified separately.
"""
function aisupdatelogimpweights!(logimpweights,
      input1::Particles, input2::Particles,
      temperature1::Float64, temperature2::Float64,
      dbm::PartitionedBernoulliDBM,
      biases::Particle,
      particles::Particles)

   nlayers = length(particles)
   weightsinput!(input1, input2, dbm, particles)

   # analytically sum out all even layers
   for i = 2:2:nlayers
      input2[i] .= input1[i]
      input2[i] .*= temperature2
      broadcast!(+, input2[i], input2[i], biases[i]')
      input1[i] .*= temperature1
      broadcast!(+, input1[i], input1[i], biases[i]')

      for n = 1:size(input1[i],2)
         for j in eachindex(logimpweights)
            logimpweights[j] +=
                  log((1 + exp.(input1[i][j,n])) / (1 + exp.(input2[i][j,n])))
         end
      end
   end
end

"""
    combinedbiases(dbm)
Returns a vector containing in the i'th element the bias vector for the i'th
layer of the `dbm`. For intermediate layers, visible and hidden biases are
combined to a single bias vector.
"""
function combinedbiases(dbm::MultimodalDBM)
   biases = Particle(undef, length(dbm) + 1)
   # Create copy to avoid accidental modification of dbm.
   # Use functions `visiblebias` and `hiddenbias` instead of
   # fields `visbias` and `hidbias` of RBMs to be able to
   # also respect PartitionedRBMs on an abstract level.
   biases[1] = copy(visiblebias(dbm[1]))
   for i = 2:length(dbm)
      biases[i] = visiblebias(dbm[i]) + hiddenbias(dbm[i-1])
   end
   biases[end] = copy(hiddenbias(dbm[end]))
   biases
end

"""
    copyannealed!(annealedrbm, rbm, temperature)
Copies all parameters that are to be annealed from the RBM `rbm` to the RBM
`annealedrbm` and anneals them with the given `temperature`.
"""
function copyannealed!(annealedrbm::BRBM, rbm::BRBM, temperature::Float64
      ) where {BRBM <: Union{BernoulliRBM, Binomial2BernoulliRBM}}

   annealedrbm.weights .= rbm.weights
   annealedrbm.weights .*= temperature
   nothing
end

function copyannealed!(annealedrbm::GRBM, gbrbm::GRBM, temperature::Float64
      ) where {GRBM <: Union{GaussianBernoulliRBM, GaussianBernoulliRBM2}}

   annealedrbm.weights .= gbrbm.weights
   annealedrbm.visbias .= gbrbm.visbias
   tempsq = sqrt(temperature)
   annealedrbm.weights .*= tempsq
   annealedrbm.visbias .*= tempsq
   nothing
end

function copyannealed!(annealedrbm::PartitionedRBM,
      prbm::PartitionedRBM, temperature::Float64)

   for i in eachindex(annealedrbm.rbms)
      copyannealed!(annealedrbm.rbms[i], prbm.rbms[i], temperature)
   end
   nothing
end

function copyannealed!(annealedrbms::Vector{T}, rbms::Vector{T},
      temperature::Float64) where {T<:AbstractRBM}

   for i in eachindex(rbms)
      copyannealed!(annealedrbms[i], rbms[i], temperature)
   end
   nothing
end

"""
    empiricalloglikelihood(x, xgen)
    empiricalloglikelihood(bm, x, nparticles)
    empiricalloglikelihood(bm, x, nparticles, burnin)
Computes the mean empirical loglikelihood for the data set `x`.
The probability of a sample is estimated to be the empirical probability of the
sample in a dataset generated by the model. This data set can be given as `xgen`
or it is generated by running a Gibbs sampler with `nparticles` for `burnin`
steps (default 5) in the Boltzmann Machine `bm`.
Throws an error if a sample in `x` is not contained in the generated data set.
"""
function empiricalloglikelihood(bm::AbstractBM, x::Matrix{Float64},
      nparticles::Int, burnin::Int = 5)

   empiricalloglikelihood(x, sampleparticles(bm, nparticles, burnin)[1])
end

function empiricalloglikelihood(x::Matrix{Float64}, xgen::Matrix{Float64})

   if size(x, 2) != size(xgen, 2)
      error("Number of variables does not match.")
   end

   genfreq = samplefrequencies(xgen)
   loglikelihood = 0.0
   norigsamples = size(x,1)
   for j = 1:norigsamples
      p = get(genfreq, vec(x[j,:]), 0)
      if p == 0
         error("Not enough samples")
      else
         loglikelihood += log(p)
      end
   end
   loglikelihood /= norigsamples
   loglikelihood
end

# Computes the energy of the combination of activations, given by the particle
# `u`, in the DBM.
function energy(dbm::BasicDBM, u::Particle)
   energy = [0.0]
   for i = 1:length(dbm)
      energy .-= u[i]'*dbm[i].weights*u[i+1] .+ dbm[i].visbias'*u[i] .+ dbm[i].hidbias'*u[i+1]
   end
   energy[1]
end

"""
    energy(rbm, v, h)
Computes the energy of the configuration of the visible nodes `v` and the
hidden nodes `h`, specified as vectors, in the `rbm`.
"""
function energy(rbm::BernoulliRBM, v::Vector{Float64}, h::Vector{Float64})
   - dot(v, rbm.weights*h) - dot(rbm.visbias, v) - dot(rbm.hidbias, h)
end

function energy(b2brbm::Binomial2BernoulliRBM, v::Vector{Float64}, h::Vector{Float64})
   # Add term to correct for 0/1/2-valued v.
   - sum(v .== 1.0) * log(2) - dot(v, b2brbm.weights*h) -
         dot(b2brbm.visbias, v) - dot(b2brbm.hidbias, h)
end

function energy(gbrbm::GaussianBernoulliRBM, v::Vector{Float64}, h::Vector{Float64})
   sum(((v - gbrbm.visbias) ./ gbrbm.sd).^2) / 2 - dot(gbrbm.hidbias, h) -
         dot(gbrbm.weights*h, v ./ gbrbm.sd)
end

function energy(gbrbm::GaussianBernoulliRBM2, v::Vector{Float64}, h::Vector{Float64})
   sum(((v - gbrbm.visbias) ./ gbrbm.sd).^2) / 2 - dot(gbrbm.hidbias, h) -
         dot(gbrbm.weights*h, v ./ (gbrbm.sd .^ 2))
end

function energy(prbm::PartitionedRBM, v::Vector{Float64}, h::Vector{Float64})
   ret = 0.0
   for i in eachindex(prbm.rbms)
      visrange = prbm.visranges[i]
      hidrange = prbm.hidranges[i]
      ret += energy(prbm.rbms[i], v[visrange], h[hidrange])
   end
   ret
end

"""
    energyzerohiddens(rbm, v)
Computes the energy for the visible activations `v` in the RBM `rbm`, if all
hidden nodes have zero activation, i. e. yields the same as
`energy´(rbm, v, zeros(rbm.hidbias))`.
"""
function energyzerohiddens(rbm::BernoulliRBM, v::Vector{Float64})
   - dot(rbm.visbias, v)
end

function energyzerohiddens(b2brbm::Binomial2BernoulliRBM, v::Vector{Float64})
   - sum(v .== 1.0) * log(2) - dot(b2brbm.visbias, v)
end

function energyzerohiddens(gbrbm::GaussianBernoulliRBM, v::Vector{Float64})
   sum(((v - gbrbm.visbias) ./ gbrbm.sd).^2) / 2
end

function energyzerohiddens(gbrbm::GaussianBernoulliRBM2, v::Vector{Float64})
   sum(((v - gbrbm.visbias) ./ gbrbm.sd).^2) / 2
end

function energyzerohiddens(prbm::PartitionedRBM, v::Vector{Float64})
   ret = 0.0
   for i in eachindex(prbm.rbms)
      visrange = prbm.visranges[i]
      ret += energyzerohiddens(prbm.rbms[i], v[visrange])
   end
   ret
end

"""
    exactloglikelihood(rbm, x)
Computes the mean log-likelihood for the given dataset `x` and the RBM `rbm`
exactly.
The log of the partition function is computed exactly by
`exactlogpartitionfunction(rbm)`.
Besides that, the function simply calls `loglikelihood(rbm, x)`.
"""
function exactloglikelihood(rbm::AbstractRBM, x::Matrix{Float64},
      logz = BMs.exactlogpartitionfunction(rbm))
   loglikelihood(rbm, x, logz)
end

"""
    exactloglikelihood(dbm, x)
    exactloglikelihood(dbm, x, logz)
Computes the mean log-likelihood for the given dataset `x` and the DBM `dbm`
exactly.
If the value of the log of the partition function of the `dbm` is not supplied
as argument `logz`, it will be computed by `exactlogpartitionfunction(dbm)`.
"""
function exactloglikelihood(mdbm::MultimodalDBM, x::Matrix{Float64},
      logz = exactlogpartitionfunction(mdbm))

   nsamples = size(x, 1)
   hiddbm::PartitionedBernoulliDBM =
         converttopartitionedbernoullidbm(mdbm[2:end])
   combinedbiases = BMs.combinedbiases(hiddbm)

   # combinations of hidden layers with odd index (i. e. h1, h3, ...)
   hodd = initcombinationoddlayersonly(hiddbm)

   logpx = 0.0
   for j = 1:nsamples
      v = vec(x[j,:])
      px = 0.0
      while true
         pun = unnormalizedproboddlayers(hiddbm, hodd, combinedbiases)
         pun *= exp(-energy(mdbm[1], v, hodd[1]))
         px += pun

         # next combination of hidden nodes' activations
         next!(hodd) || break
      end

      logpx += log(px)
   end

   logpx /= nsamples
   logpx -= logz
   logpx
end

"""
    exactlogpartitionfunction(rbm)
Calculates the log of the partition function of the BernoulliRBM `rbm` exactly.
The execution time grows exponentially with the minimum of
(number of visible nodes, number of hidden nodes).
"""
function exactlogpartitionfunction(rbm::BernoulliRBM)
   nvisible = length(rbm.visbias)
   nhidden = length(rbm.hidbias)
   z = 0.0
   if nvisible <= nhidden
      v = zeros(nvisible)
      while true
         z += exp.(-freeenergy(rbm, v))
         next!(v) || break
      end
   else
      h = zeros(nhidden)
      revrbm = reversedrbm(rbm)
      while true
         z += exp.(-freeenergy(revrbm, h))
         next!(h) || break
      end
   end
   log(z)
end

function exactlogpartitionfunction(rbm::Binomial2BernoulliRBM)
   nhidden = length(rbm.hidbias)
   h = zeros(nhidden)
   z = 0.0
   while true
      z += unnormalizedprobhidden(rbm, h)
      next!(h) || break
   end
   log(z)
end

"""
    exactlogpartitionfunction(gbrbm)
Calculates the log of the partition function of the GaussianBernoulliRBM `gbrbm`
exactly. The execution time grows exponentially with the number of hidden nodes.
"""
function exactlogpartitionfunction(gbrbm::GaussianBernoulliRBM)
   nvisible = length(gbrbm.visbias)
   nhidden = length(gbrbm.hidbias)

   h = zeros(nhidden)
   z = 0.0
   while true
      wh = gbrbm.weights * h
      z += exp(dot(gbrbm.hidbias, h) + sum(0.5*wh.^2 + gbrbm.visbias ./ gbrbm.sd .* wh))
      next!(h) || break
   end
   log(z) + nvisible/2 * log(2*pi) + sum(log.(gbrbm.sd))
end

function exactlogpartitionfunction(gbrbm::GaussianBernoulliRBM2)
   nvisible = length(gbrbm.visbias)
   nhidden = length(gbrbm.hidbias)

   h = zeros(nhidden)
   z = 0.0
   while true
      wh = gbrbm.weights * h
      z += exp(dot(gbrbm.hidbias, h) +
            sum((wh.^2 / 2 + gbrbm.visbias .* wh) ./ (gbrbm.sd .^2)))
      next!(h) || break
   end
   log(z) + nvisible/2 * log(2*pi) + sum(log.(gbrbm.sd))
end

"""
    exactlogpartitionfunction(bgrbm)
Calculates the log of the partition function of the BernoulliGaussianRBM `bgrbm`
exactly. The execution time grows exponentially with the number of visible nodes.
"""
function exactlogpartitionfunction(bgrbm::BernoulliGaussianRBM)
   exactlogpartitionfunction(reversedrbm(bgrbm))
end

"""
    exactlogpartitionfunction(dbm)
Calculates the log of the partition function of the DBM `dbm` exactly.
If the number of hidden layers is even, the execution time grows exponentially
with the total number of nodes in hidden layers with odd indexes
(i. e. h1, h3, ...).
If the number of hidden layers is odd, the execution time grows exponentially
with the minimum of
(number of nodes in layers with even index,
number of nodes in layers with odd index).
"""
function exactlogpartitionfunction(dbm::PartitionedBernoulliDBM)
   nunits = BMs.nunits(dbm)
   nhiddenlayers = length(dbm)

   # Apply algorithm on reversed DBM if it requires fewer iterations there
   if nhiddenlayers % 2 == 1
      noddlayersnodes = sum(nunits[1:2:end])
      nevenlayersnodes = sum(nunits[2:2:end])
      if noddlayersnodes < nevenlayersnodes
         dbm = reverseddbm(dbm)
         nunits = reverse(nunits)
      end
   end

   invoke(exactlogpartitionfunction, Tuple{MultimodalDBM,}, dbm)
end

"""
    exactlogpartitionfunction(mdbm)
Calculates the log of the partition function of the MultimodalDBM `mdbm`
exactly.
The execution time grows exponentially with the total number of nodes in hidden
layers with odd indexes (i. e. h1, h3, ...).
"""
function exactlogpartitionfunction(mdbm::MultimodalDBM)
   hodd = initcombinationoddlayersonly(mdbm[2:end])
   hiddbm::PartitionedBernoulliDBM =
         converttopartitionedbernoullidbm(mdbm[2:end])
   hiddbmbiases = combinedbiases(hiddbm)
   z = 0.0
   while true
      pun = unnormalizedproboddlayers(hiddbm, hodd, hiddbmbiases)
      pun *= unnormalizedprobhidden(mdbm[1], hodd[1])
      z += pun
      next!(hodd) || break
   end
   log(z)
end

"""
    MultivariateBernoulliDistribution(bm)
Calculates and stores the probabilities for all possible combinations of a
multivariate Bernoulli distribution defined by a Boltzmann machine model
with Bernoulli distributed visible nodes.
Can be used for sampling from this distribution, see `samples`.
"""
struct MultivariateBernoulliDistribution
   cumprobs::Vector{Float64}
   samples::Vector{Vector{Float64}}

   function MultivariateBernoulliDistribution(bm::AbstractBM)
      nvisible = nunits(bm)[1]

      # create samples, a vector of vectors, covering all
      # theoretically possible combinations of the visible nodes' states
      nviscombinations = 2^nvisible
      samples = Vector{Vector{Float64}}(undef, nviscombinations)
      v = zeros(nvisible)
      for i = 1:nviscombinations
         samples[i] = copy(v)
         next!(v)
      end

      # calculate unnormalized probabilities of all samples
      probs = unnormalizedprobs(bm, samples)

      # sort both the unnormalized probs and the samples
      sortingperm = sortperm(probs)
      permute!(probs, sortingperm)
      permute!(samples, sortingperm)

      # transform the unnormalized probabilities in a vector of cumulative
      # probabilities (corresponding to intervals) that can be used for sampling
      z = sum(probs)
      probs ./= z
      cumsum!(probs, probs)

      new(probs, samples)
   end
end

"""
    freeenergy(rbm, x)
Computes the average free energy of the samples in the dataset `x` for the
AbstractRBM `rbm`.
"""
function freeenergy(rbm::AbstractRBM, x::Matrix{Float64})
   nsamples = size(x, 1)
   freeenergy = 0.0
   for j = 1:nsamples
      v = vec(x[j, :])
      freeenergy += BMs.freeenergy(rbm, v)
   end
   freeenergy /= nsamples
   freeenergy
end

function freeenergy(bgrbm::BernoulliGaussianRBM, x::Matrix{Float64})
   nsamples = size(x, 1)
   nhidden = length(bgrbm.hidbias)

   freeenergy = 0.0
   for j = 1:nsamples
      v = vec(x[j,:])
      wtv = bgrbm.weights' * v
      freeenergy -= dot(bgrbm.hidbias, wtv) + 0.5 * sum(wtv.^2) + dot(bgrbm.visbias, v)
   end
   freeenergy /= nsamples
   freeenergy -= nhidden / 2 * log(2pi)
   freeenergy
end

"""
    freeenergy(rbm, v)
Computes the free energy of the sample `v` (a vector) for the `rbm`.
"""
function freeenergy(rbm::BernoulliRBM, v::Vector{Float64})
  - dot(rbm.visbias, v) - sum(log1p.(exp.(hiddeninput(rbm, v))))
end

function freeenergy(b2brbm::Binomial2BernoulliRBM, v::Vector{Float64})
   # To get probabilities for 0/1/2-valued v, multiply probabilities for v's
   # with 2^(number of 1s in v), because each v can be represented by
   # this number of combinations in the 00/01/10/11 space, all having equal
   # probability.
   - sum(v .== 1.0) * log(2) -
         dot(b2brbm.visbias, v) - sum(log1p.(exp.(hiddeninput(b2brbm, v))))
end

function freeenergy(gbrbm::GaussianBernoulliRBM, v::Vector{Float64})
   nhidden = length(gbrbm.hidbias)
   vscaled = v ./ gbrbm.sd
   freeenergy = 0.0
   for k = 1:nhidden
      freeenergy -= log1pexp(gbrbm.hidbias[k] + dot(gbrbm.weights[:,k], vscaled))
   end
   freeenergy += 0.5 * sum(((v - gbrbm.visbias) ./ gbrbm.sd).^2)
end

function freeenergy(gbrbm::GaussianBernoulliRBM2, v::Vector{Float64})
   nhidden = length(gbrbm.hidbias)
   vscaled = v ./ (gbrbm.sd .^ 2)
   freeenergy = 0.0
   for k = 1:nhidden
      freeenergy -= log1pexp(gbrbm.hidbias[k] + dot(gbrbm.weights[:,k], vscaled))
   end
   freeenergy += 0.5 * sum(((v - gbrbm.visbias) ./ gbrbm.sd).^2)
end

"""
    freeeenergydiffs(rbm1, rbm2, x)
Computes the differences of the free energy for the samples in the dataset `x`
regarding the RBM models `rbm1` and `rbm2`. Returns a vector of differences.
"""
function freeenergydiffs(rbm1::AbstractRBM, rbm2::AbstractRBM, x::Matrix{Float64})
   nsamples = size(x, 1)
   freeenergydiffs = Vector{Float64}(undef, nsamples)
   for j = 1:nsamples
      v = vec(x[j,:])
      freeenergydiffs[j] = freeenergy(rbm1, v) - freeenergy(rbm2, v)
   end
   freeenergydiffs
end

# TODO: optimize freenergydiffs for BernoulliRBM (Matrix operation)

function hiddenbias(rbm::AbstractRBM)
   rbm.hidbias
end

function hiddenbias(prbm::PartitionedRBM)
   vcat(map(rbm -> rbm.hidbias, prbm.rbms)...)
end

"
Returns particle for DBM, initialized with zeros.
"
function initcombination(dbm::BasicDBM)
   nunits = BMs.nunits(dbm)
   nlayers = length(nunits)
   u = Particle(undef, nlayers)
   for i = 1:nlayers
      u[i] = zeros(nunits[i])
   end
   u
end

"""
    initcombinationoddlayersonly(dbm)
Creates and zero-initializes a particle for layers with odd indexes
in the `dbm`.
"""
function initcombinationoddlayersonly(dbm::MultimodalDBM)
   nunits = BMs.nunits(dbm)
   nlayers = length(nunits)
   uodd = Particle(undef, round(Int, nlayers/2, RoundUp))
   for i = eachindex(uodd)
      uodd[i] = zeros(nunits[2i-1])
   end
   uodd
end

"""
    loglikelihood(rbm, x)
    loglikelihood(rbm, x, logz)
Computes the average log-likelihood of an RBM on a given dataset `x`.
Uses `logz` as value for the log of the partition function
or estimates the partition function with Annealed Importance Sampling.
"""
function loglikelihood(rbm::AbstractRBM, x::Array{Float64,2},
      logz::Float64 = logpartitionfunction(rbm))

   -freeenergy(rbm, x) - logz
end

"""
    loglikelihood(dbm, x; ...)
Estimates the mean log-likelihood of the DBM on the data set `x`
with Annealed Importance Sampling.
This requires a separate run of AIS for each sample.
"""
function loglikelihood(mdbm::MultimodalDBM, x::Matrix{Float64}, logz::Float64 = -Inf;
      parallelized = false,
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   nsamples = size(x, 1)

   if logz == -Inf
      logz = logpartitionfunction(mdbm;
            parallelized = parallelized, ntemperatures = ntemperatures,
            temperatures = temperatures, nparticles = nparticles, burnin = burnin)
   end

   if parallelized
      # divide data set x into batches and compute unnormalized probabilties
      batches = mostevenbatches(nsamples)
      batchranges = ranges(batches)
      logp = @sync @distributed (+) for i in 1:length(batches)
         (batches[i] / nsamples) * # (weighted mean)
            unnormalizedlogprob(mdbm, x[batchranges[i], :];
                  ntemperatures = ntemperatures, temperatures = temperatures,
                  nparticles = nparticles, burnin = burnin)
      end
   else
      logp = unnormalizedlogprob(mdbm, x;
            ntemperatures = ntemperatures,
            temperatures = temperatures, nparticles = nparticles, burnin = burnin)
   end
   logp -= logz
   logp
end

"""
    loglikelihooddiff(rbm1, rbm2, x)
    loglikelihooddiff(rbm1, rbm2, x, logzdiff)
    loglikelihooddiff(rbm1, rbm2, x, logimpweights)
Computes difference of the loglikelihood functions of the two RBMs on the data
matrix `x`, averaged over the samples.
For this purpose, the partition function ratio Z2/Z1 is estimated by AIS unless
the importance weights are specified the by parameter `logimpweights` or the
difference in the log partition functions is given by `logzdiff`.

The first model is better than the second if the returned value is positive.
"""
function loglikelihooddiff(rbm1::R, rbm2::R,
      x::Array{Float64,2},
      logzdiff::Float64,
      ) where {R<:AbstractRBM}

   lldiff = mean(freeenergydiffs(rbm2, rbm1, x))
   lldiff += logzdiff # estimator for log(Z2/Z1)
   lldiff
end

function loglikelihooddiff(rbm1::R, rbm2::R,
      x::Array{Float64,2},
      logimpweights::Array{Float64,1} = aislogimpweights(rbm1, rbm2)
      ) where {R<:AbstractRBM}

   loglikelihooddiff(rbm1, rbm2, x, BMs.logmeanexp(logimpweights))
end

""" Performs numerically stable computation of the mean on log-scale. """
function logmeanexp(v::Vector{Float64})
   logsumexp(v) - log(length(v))
end

""" Performs numerically stable summation on log-scale. """
function logsumexp(v::Vector{Float64})
   vmax = maximum(v)
   vmax + log(sum(exp.(v .- vmax)))
end

"""
    logpartitionfunction(bm; ...)
    logpartitionfunction(bm, logr)
Calculates or estimates the log of the partition function
of the Boltzmann Machine `bm`.

`r` is an estimator of the ratio of the `bm`'s partition function Z to the
partition function Z_0 of the reference BM with zero weights but same biases
as the given `bm`. In case of a GaussianBernoulliRBM, the reference model
also has the same standard deviation parameter.
The estimated partition function of the Boltzmann Machine is Z = r * Z_0
with `r` being the mean of the importance weights.
Therefore, the log of the estimated partition function is
log(Z) = log(r) + log(Z_0)

If the log of `r` is not given as argument `logr`, Annealed Importance Sampling
(AIS) is performed to get a value for it. In this case,
the optional arguments for AIS can be specified (see `aislogimpweights`),
and the optional boolean argument `parallelized` can
be used to turn on batch-parallelized computing of the importance weights.
"""
function logpartitionfunction(bm::AbstractBM, logr::Float64)
   logr + logpartitionfunctionzeroweights(bm)
end

function logpartitionfunction(bm::AbstractBM;
      parallelized::Bool = false,
      # optional arguments for AIS:
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   if parallelized
      logimpweights = BMs.batchparallelized(
            n -> BMs.aislogimpweights(bm;
                  ntemperatures = ntemperatures, temperatures = temperatures,
                  nparticles = n, burnin = burnin),
            nparticles, vcat)
   else
      logimpweights = BMs.aislogimpweights(bm;
            ntemperatures = ntemperatures, temperatures = temperatures,
            nparticles = nparticles, burnin = burnin)
   end

   logpartitionfunction(bm, logmeanexp(logimpweights))
end

"""
    logpartitionfunctionzeroweights(bm)
Returns the value of the log of the partition function of the Boltzmann Machine
that results when one sets the weights of `bm` to zero,
and leaves the other parameters (biases) unchanged.
"""
function logpartitionfunctionzeroweights(rbm::AbstractRBM)
   logpartitionfunctionzeroweights_visterm(rbm) +
         logpartitionfunctionzeroweights_hidterm(rbm)
end

function logpartitionfunctionzeroweights(prbm::PartitionedRBM)
   mapreduce(rbm -> logpartitionfunctionzeroweights(rbm), +, prbm.rbms)
end

function logpartitionfunctionzeroweights(dbm::PartitionedBernoulliDBM)
   logz0 = 0.0
   biases = combinedbiases(dbm)
   for i in eachindex(biases)
      logz0 += sum(log1p.(exp.(biases[i])))
   end
   logz0
end

function logpartitionfunctionzeroweights(mdbm::MultimodalDBM)
   if length(mdbm) == 1
      logpartitionfunctionzeroweights(mdbm[1])
   else
      hiddbm = mdbm[2:end]
      hiddbm[1] = deepcopy(hiddbm[1])
      setvisiblebias!(hiddbm[1], visiblebias(hiddbm[1]) + hiddenbias(mdbm[1]))
      logpartitionfunctionzeroweights_visterm(mdbm[1]) +
            invoke(logpartitionfunctionzeroweights,
                  Tuple{PartitionedBernoulliDBM,},
                  converttopartitionedbernoullidbm(hiddbm))
   end
end

function logpartitionfunctionzeroweights_visterm(rbm::BernoulliRBM)
   sum(log1p.(exp.(rbm.visbias)))
end

function logpartitionfunctionzeroweights_visterm(bgrbm::BernoulliGaussianRBM)
   sum(log1p.(exp.(bgrbm.visbias)))
end

function logpartitionfunctionzeroweights_visterm(b2brbm::Binomial2BernoulliRBM)
   2*sum(log1p.(exp.(b2brbm.visbias)))
end

function logpartitionfunctionzeroweights_visterm(
      gbrbm::Union{GaussianBernoulliRBM, GaussianBernoulliRBM2})

   nvisible = length(gbrbm.visbias)
   nvisible / 2 * log(2*pi) + sum(log.(gbrbm.sd))
end

function logpartitionfunctionzeroweights_visterm(prbm::PartitionedRBM)
   mapreduce(logpartitionfunctionzeroweights_visterm, +, prbm.rbms)
end

function logpartitionfunctionzeroweights_hidterm(rbm::AbstractXBernoulliRBM)
   sum(log1p.(exp.(rbm.hidbias)))
end

function logpartitionfunctionzeroweights_hidterm(bgrbm::BernoulliGaussianRBM)
   nhidden = length(bgrbm.hidbias)
   nhidden / 2 * log(2pi)
end

function logpartitionfunctionzeroweights_hidterm(prbm::PartitionedRBM)
   mapreduce(logpartitionfunctionzeroweights_hidterm, +, prbm.rbms)
end

"""
    logproblowerbound(dbm, x; ...)
    logproblowerbound(dbm, x, logimpweights; ...)
    logproblowerbound(dbm, x, logz; ...)
Estimates the mean of the variational lower bound for the log probability
of the DBM on a given dataset `x` like described in Equation 38
in [Salakhutdinov, 2015].
The logarithmized partition function can be specified directly as `logz`
or by giving the `logimpweights` from estimating the partition function
with the Annealed Importance Sampling algorithm (AIS).
(See `aislogimpweights`.)
If neither `logimpweights` or `logz` is given, the partition function will be
estimated by AIS with default parameters.

# Optional keyword argument:
* The approximate posterior distribution may be given as argument `mu`
  or is calculated by the mean-field method.
"""
function logproblowerbound(dbm::MultimodalDBM,
      x::Matrix{Float64},
      logimpweights::Array{Float64,1} = aislogimpweights(dbm);
      mu::Particles = meanfield(dbm, x))

   logproblowerbound(dbm, x,
         logpartitionfunction(dbm, logmeanexp(logimpweights));
         mu = mu)
end

function logproblowerbound(dbm::MultimodalDBM,
      x::Matrix{Float64}, logz::Float64;
      mu::Particles = meanfield(dbm, x))

   nsamples = size(mu[1], 1)
   nrbms  = length(dbm)

   lowerbound = 0.0
   for j = 1:nsamples # TODO parallelize

      for i = 1:nrbms
         v = vec(mu[i][j,:])   # visible units of i'th RBM
         h = vec(mu[i+1][j,:]) # hidden units of i'th RBM

         # add energy
         lowerbound -= energy(dbm[i], v, h)

         # add entropy of approximate posterior Q
         h = h[(h .> 0.0) .& (h .< 1.0)]
         lowerbound += - dot(h, log.(h)) - dot(1.0 .- h, log.(1.0 .- h))
      end
   end

   lowerbound /= nsamples
   lowerbound -= logz
   lowerbound
end

function mixedrbm(rbm1::BernoulliRBM, rbm2::BernoulliRBM, temperature::Float64)
   visbias = (1 - temperature) * rbm1.visbias + temperature * rbm2.visbias
   weights = hcat(
         (1 - temperature) * rbm1.weights,
         temperature * rbm2.weights)
   hidbias = vcat(
         (1 - temperature) * rbm1.hidbias,
         temperature * rbm2.hidbias)
   BernoulliRBM(weights, visbias, hidbias)
end

function mixedrbm(rbm1::GBRBM, rbm2::GBRBM, temperature::Float64
      ) where {GBRBM <: Union{GaussianBernoulliRBM, GaussianBernoulliRBM2}}

   sdf1 = (1 - temperature) ./ rbm1.sd.^2
   sdf2 = temperature ./ rbm2.sd.^2
   sdsq = 1 ./ (sdf1 + sdf2)
   sd = sqrt.(sdsq)
   visbias = (sdf1 .* rbm1.visbias + sdf2 .* rbm2.visbias) .* sdsq
   weights = hcat(
         (1-temperature) * rbm1.weights ./ (rbm1.sd ./ sd),
         temperature * rbm2.weights ./ (rbm2.sd ./ sd))
   hidbias = vcat(
         (1 - temperature) * rbm1.hidbias,
         temperature * rbm2.hidbias)
   GBRBM(weights, visbias, hidbias, sd)
end

"""
    next!(combination)
Sets the vector `combination`, containing a sequence of the values 0.0 and 1.0,
to the next combination of 0.0s and 1.0s.
Returns false if the new combination consists only of zeros; true otherwise.
"""
function next!(combination::T) where {T <: AbstractArray{Float64, 1}}
   i = 1
   while i <= length(combination)
      if combination[i] == 0.0
         combination[i] = 1.0
         break
      else
         combination[i] = 0.0
         i += 1
      end
   end

   i <= length(combination)
end

"""
    next!(particle)
Sets `particle` to the next combination of nodes' activations.
Returns false if the loop went through all combinations; true otherwise.
"""
function next!(particle::Particle)
   n = length(particle)
   i = 1
   while i <= n && !next!(particle[i])
      i += 1
   end
   i <= n
end

"""
    nunits(bm)
Returns an integer vector that contans in the i'th entry the number of nodes
in the i'th layer of the `bm`.
"""
function nunits(rbm::AbstractRBM)
   [length(rbm.visbias); length(rbm.hidbias)]
end

function nunits(prbm::PartitionedRBM)
   [
      sum(map(rbm -> length(rbm.visbias), prbm.rbms));
      sum(map(rbm -> length(rbm.hidbias), prbm.rbms))
   ]
end

function nunits(dbm::MultimodalDBM)
   nrbms = length(dbm)
   if nrbms == 0
      error("Nodes and layers not defined in empty DBM")
   end
   nlayers = length(dbm) + 1
   nu = Vector{Int}(undef, nlayers)
   for i = 1:nrbms
      nu[i] = nvisiblenodes(dbm[i])
   end
   nu[nlayers] = nhiddennodes(dbm[nlayers-1])
   nu
end

"""
    nmodelparameters(bm)
Returns the number of parameters in the Boltzmann Machine model `bm`.
"""
function nmodelparameters(bm::AbstractBM)
   nunits = BMs.nunits(bm)
   nweights = 0
   for i = 1:(length(nunits) - 1)
      nweights += nunits[i] * nunits[i+1]
   end
   nbiases = sum(nunits) # number of bias variables = number of nodes
   nweights + nbiases
end

function nmodelparameters(gbrbm::GaussianBernoulliRBM)
   invoke(nmodelparameters, (AbstractRBM,), gbrbm) + length(gbrbm.sd)
end

# TODO respect model parameters for MultimodalDBM with Gaussian nodes

"""
    reconstructionerror(rbm, x)
Computes the mean reconstruction error of the RBM on the dataset `x`.
"""
function reconstructionerror(rbm::AbstractRBM,
      x::Array{Float64,2},
      upfactor::Float64 = 1.0,
      downfactor::Float64 = 1.0)

   nsamples = size(x, 1)
   reconstructionerror = 0.0

   vmodel = Vector{Float64}(undef, nvisiblenodes(rbm))
   hmodel = Vector{Float64}(undef, nhiddennodes(rbm))

   for sample = 1:nsamples
      v = vec(x[sample,:])
      hiddenpotential!(hmodel, rbm, v, upfactor)
      visiblepotential!(vmodel, rbm, hmodel, downfactor)
      reconstructionerror += sum(abs.(v - vmodel))
   end

   reconstructionerror /= nsamples
end

"
Returns the GBRBM with weights such that hidden and visible of the given `bgrbm`
are switched and a visible standard deviation of 1.
"
function reversedrbm(bgrbm::BernoulliGaussianRBM)
   sd = ones(length(bgrbm.hidbias))
   GaussianBernoulliRBM(bgrbm.weights', bgrbm.hidbias, bgrbm.visbias, sd)
end

function reversedrbm(rbm::BernoulliRBM)
   BernoulliRBM(rbm.weights', rbm.hidbias, rbm.visbias)
end

function reversedrbm(prbm::PartitionedRBM{BernoulliRBM})
   PartitionedRBM{BernoulliRBM}(map(reversedrbm, prbm.rbms))
end

function reverseddbm(dbm::PartitionedBernoulliDBM)
   revdbm = reverse(dbm)
   map!(reversedrbm, revdbm, revdbm)
end

function unnormalizedproblogratios(rbm::BernoulliRBM,
      hh::Matrix{Float64},
      temperature1::Float64,
      temperature2::Float64)

   weightsinput = hh * rbm.weights'
   vec(sum(log.(
         (1.0 .+ exp.(temperature1 .* weightsinput .+ rbm.visbias')) ./
         (1.0 .+ exp.(temperature2 .* weightsinput .+ rbm.visbias'))), dims = 2))
end

function unnormalizedproblogratios(rbm::Binomial2BernoulliRBM,
      hh::Matrix{Float64},
      temperature1::Float64,
      temperature2::Float64)

   weightsinput = hh * rbm.weights'
   vec(sum(log.(
         (1.0 .+ exp.(temperature1 .* weightsinput .+ rbm.visbias')) ./
         (1.0 .+ exp.(temperature2 .* weightsinput .+ rbm.visbias'))), dims = 2) * 2)
end

function unnormalizedproblogratios(gbrbm::GaussianBernoulliRBM,
      hh::Matrix{Float64},
      temperature1::Float64,
      temperature2::Float64)

   wht = hh * gbrbm.weights'
   vec((temperature1 - temperature2) * sum(
         0.5 * wht.^2 + broadcast(*, wht, (gbrbm.visbias ./ gbrbm.sd)'), dims = 2))
end

function unnormalizedproblogratios(gbrbm::GaussianBernoulliRBM2,
   hh::Matrix{Float64},
   temperature1::Float64,
   temperature2::Float64)

   wht = hh * gbrbm.weights'
   vec((temperature1 - temperature2) * sum(
         (0.5 * wht.^2 + wht .* gbrbm.visbias') ./ (gbrbm.sd .^2)', dims = 2))
end

function unnormalizedproblogratios(prbm::PartitionedRBM,
      hh::Matrix{Float64},
      temperature1::Float64,
      temperature2::Float64)

   # TODO does not work with views
   mapreduce(
         i -> unnormalizedproblogratios(prbm.rbms[i],
               hh[:, prbm.hidranges[i]], temperature1, temperature2),
         (x,y) -> broadcast(+, x, y),
         eachindex(prbm.rbms))
end

"""
    samplefrequencies(x)
Returns a dictionary containing the rows of the data set `x` as keys and their
relative frequencies as values.
"""
function samplefrequencies(x::Array{T,2}) where T
   dict = Dict{Array{T}, Float64}()
   nsamples = size(x, 1)
   onefrac = 1/nsamples
   for i=1:nsamples
      sample = vec(x[i,:]);
      if haskey(dict, sample)
         dict[sample] += onefrac
      else
         dict[sample] = onefrac
      end
   end
   dict
end

function samples(mbdist::MultivariateBernoulliDistribution, nsamples::Int)
   nvisible = length(mbdist.samples[1])
   ret = Matrix{Float64}(undef, nsamples, nvisible)
   for i = 1:nsamples
      sampledindex = searchsortedfirst(mbdist.cumprobs, rand())
      ret[i, :] .= mbdist.samples[sampledindex]
   end
   ret
end

function setvisiblebias!(rbm::BernoulliRBM, v::Vector{Float64})
   rbm.visbias .= v
end

function setvisiblebias!(prbm::PartitionedRBM, v::Vector{Float64})
   for i in eachindex(prbm.rbms)
      visrange = prbm.visranges[i]
      prbm.rbms[i].visbias .= v[visrange]
   end
end

"""
    unnormalizedlogprob(mdbm, x; ...)
Estimates the mean unnormalized log probability of the samples (rows in `x`)
in the MultimodalDBM `mdbm` by running the Annealed Importance Sampling (AIS)
in a smaller modified DBM for each sample.

The named optional arguments for AIS can be specified here.
(See `aislogimpweights`)
"""
function unnormalizedlogprob(mdbm, x::Matrix{Float64};
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   nsamples = size(x, 1)
   logp = 0.0
   hiddbm = deepcopy(mdbm[2:end])
   h1 = Vector{Float64}(undef, nhiddennodes(mdbm[1]))
   visbias2 = visiblebias(mdbm[2]) # visible bias of second RBM

   for j = 1:nsamples
      v = vec(x[j,:])

      setvisiblebias!(hiddbm[1], visbias2) # reset to original bias

      logp -= energyzerohiddens(mdbm[1], v)

      # Integrate the energy of the sample v
      # into the visible bias of the first hidden layer.
      hiddeninput!(h1, mdbm[1], v)
      h1 .+= visbias2
      setvisiblebias!(hiddbm[1], h1)

      # Then compute the log partition function of the new DBM consisting
      # only of the hidden layer RBMs (with the first layer modified).
      logp += logpartitionfunction(hiddbm;
            ntemperatures = ntemperatures, temperatures = temperatures,
             nparticles = nparticles, burnin = burnin)
   end

   logp /= nsamples
   logp
end

"""
    unnormalizedprobs(bm, samples)
Calculates the unnormalized probabilities for all `samples` (vector of vectors),
in the Boltzmann Machine `bm`.

The visible nodes of the `bm` must be Bernoulli distributed.
"""
function unnormalizedprobs(rbm::Union{BernoulliRBM, BernoulliGaussianRBM},
      samples::Vector{Vector{Float64}})

   nsamples = length(samples)
   probs = Vector{Float64}(undef, nsamples)
   for i = 1:nsamples
      probs[i] = exp(-freeenergy(rbm, samples[i]))
   end
   probs
end

function unnormalizedprobs(dbm::PartitionedBernoulliDBM,
      samples::Vector{Vector{Float64}})

   nsamples = length(samples)
   biases = combinedbiases(dbm)
   oddparticle = initcombinationoddlayersonly(dbm)
   oddhiddenparticles = oddparticle[2:end]
   probs = Vector{Float64}(undef, nsamples)
   for i = 1:nsamples
      oddparticle[1] = samples[i]
      probs[i] = 0.0
      while true
         probs[i] += unnormalizedproboddlayers(dbm, oddparticle, biases)
         next!(oddhiddenparticles) || break
      end
   end
   probs
end

"
Computes the unnormalized probability of the nodes in layers with odd indexes,
i. e. p*(v, h2, h4, ...).
"
function unnormalizedproboddlayers(dbm::PartitionedBernoulliDBM, uodd::Particle,
      combinedbiases = BMs.combinedbiases(dbm))

   nlayers = length(dbm) + 1
   nintermediatelayerstobesummedout = div(nlayers - 1, 2)
   pun = 1.0
   for i = 1:nintermediatelayerstobesummedout
      pun *= prod(1.0 .+ exp.(
            hiddeninput(dbm[2i-1], uodd[i]) .+ visibleinput(dbm[2i], uodd[i+1])))
   end
   if nlayers % 2 == 0
      pun *= prod(1.0 .+ exp.(hiddeninput(dbm[end], uodd[end])))
   end
   for i = 1:length(uodd)
      pun *= exp(dot(combinedbiases[2i-1], uodd[i]))
   end

   pun
end

"""
    unnormalizedprobhidden(rbm, h)
    unnormalizedprobhidden(gbrbm, h)
Calculates the unnormalized probability of the `rbm`'s hidden nodes'
activations given by `h`.
"""
function unnormalizedprobhidden(rbm::BernoulliRBM, h::Vector{Float64})
   exp(dot(rbm.hidbias, h)) * prod(1.0 .+ exp.(visibleinput(rbm, h)))
end

function unnormalizedprobhidden(rbm::Binomial2BernoulliRBM, h::Vector{Float64})
   exp(dot(rbm.hidbias, h)) * prod(1.0 .+ exp.(visibleinput(rbm, h)))^2
end

const sqrt2pi = sqrt(2pi)

function unnormalizedprobhidden(gbrbm::GaussianBernoulliRBM, h::Vector{Float64})
   nvisible = length(gbrbm.visbias)
   wh = gbrbm.weights * h
   exp(dot(gbrbm.hidbias, h) + sum(0.5*wh.^2 + gbrbm.visbias ./ gbrbm.sd .* wh)) *
         prod(gbrbm.sd) * sqrt2pi^nvisible
end

function unnormalizedprobhidden(gbrbm::GaussianBernoulliRBM2, h::Vector{Float64})
   nvisible = length(gbrbm.visbias)
   wh = gbrbm.weights * h
   exp(dot(gbrbm.hidbias, h) +
         sum((0.5*wh.^2 + gbrbm.visbias .* wh) ./ gbrbm.sd .^ 2)) *
         prod(gbrbm.sd) * sqrt2pi^nvisible
end

function unnormalizedprobhidden(prbm::PartitionedRBM, h::Vector{Float64})
   prob = 1.0
   for i in eachindex(prbm.rbms)
      hidrange = prbm.hidranges[i]
      prob *= unnormalizedprobhidden(prbm.rbms[i], h[hidrange])
   end
   prob
end

function visiblebias(rbm::AbstractRBM)
   rbm.visbias
end

function visiblebias(prbm::PartitionedRBM)
   vcat(map(rbm -> rbm.visbias, prbm.rbms)...)
end

"""
    weightsinput!(input, input2, dbm, particles)
Computes the input that results only from the weights (without biases)
and the previous states in `particles` for all nodes in the DBM
`dbm` and stores it in `input`.
The state of the `particles` and the `dbm` is not altered.
`input2` must have the same size as `input` and `particle`.
For performance reasons, `input2` is used as preallocated space for storing
intermediate results.
"""
function weightsinput!(input::Particles, input2::Particles,
      dbm::PartitionedBernoulliDBM, particles::Particles)

   # first layer gets input only from layer above
   weightsvisibleinput!(input[1], dbm[1], particles[2])

   # intermediate layers get input from layers above and below
   for i = 2:(length(particles) - 1)
      weightsvisibleinput!(input2[i], dbm[i], particles[i+1])
      weightshiddeninput!(input[i], dbm[i-1], particles[i-1])
      input[i] .+= input2[i]
   end

   # last layer gets only input from layer below
   weightshiddeninput!(input[end], dbm[end], particles[end-1])
   input
end

function weightshiddeninput!(hh::M, rbm::BernoulliRBM,
      vv::M) where {M<:AbstractArray{Float64}}

   mul!(hh, vv, rbm.weights)
end

function weightshiddeninput!(hh::M, prbm::PartitionedRBM{BernoulliRBM},
      vv::M) where {M<:AbstractArray{Float64}}

   for i in eachindex(prbm.rbms)
      visrange = prbm.visranges[i]
      hidrange = prbm.hidranges[i]
      weightshiddeninput!(
            view(hh, :, hidrange), prbm.rbms[i], view(vv,:, visrange))
   end
   hh
end

function weightsvisibleinput!(vv::M, rbm::BernoulliRBM,
      hh::M) where {M<:AbstractArray{Float64}}

   mul!(vv, hh, transpose(rbm.weights))
end

function weightsvisibleinput!(vv::M, prbm::PartitionedRBM{BernoulliRBM},
      hh::M) where {M<:AbstractArray{Float64}}

   for i in eachindex(prbm.rbms)
      visrange = prbm.visranges[i]
      hidrange = prbm.hidranges[i]
      weightsvisibleinput!(
            view(vv, :, visrange), prbm.rbms[i], view(hh,:, hidrange))
   end
   vv
end

################################################################################
####################   Unsorted Functions (Rename?)         ####################
################################################################################

"""
    akaikeinformationcriterion(bm, loglikelihood)
Calculates the Akaike information criterion for a Boltzmann Machine, given its
`loglikelihood`.
"""
function akaikeinformationcriterion(bm::AbstractBM, loglikelihood::Float64)
   2*BMs.nmodelparameters(bm) - 2*loglikelihood
end

"""
    bayesianinformationcriterion(bm, nvariables, loglikelihood)
Calculates the Akaike information criterion for a Boltzmann machine, given its
`loglikelihood` and the number of samples `nsamples`.
"""
function bayesianinformationcriterion(bm::AbstractBM, nsamples::Int, loglikelihood::Float64)
   -2*loglikelihood + BMs.nmodelparameters(bm)*log(nsamples)
end

"""
    bernoulliloglikelihoodbaserate(nvariables)
Calculates the log-likelihood for a random sample in the "base-rate" BM with all
parameters being zero and thus all visible units being independent and
Bernoulli distributed.
"""
function bernoulliloglikelihoodbaserate(nvariables::Int)
   - nvariables * log(2)
end

"""
    bernoulliloglikelihoodbaserate(x)
Calculates the log-likelihood for the data set `x` in the "base-rate" BM with
all weights being zero and visible bias set to the empirical probability of the
samples' components in `x` being 1.
"""
function bernoulliloglikelihoodbaserate(x::Matrix{Float64})
   p = mean(x, dims = 1)
   nsamples, nvariables = size(x)
   loglikelihood = 0.0
   for i = 1:nvariables
      for j = 1:nsamples
         loglikelihood += log(x[j,i] == 1 ? p[i] : 1-p[i])
      end
   end
   loglikelihood /= nsamples
   loglikelihood
end

"""
    gaussianloglikelihoodbaserate(x)
Calculates the mean log-likelihood for the data set `x` with all variables and
components of the variables being independent and Gaussian distributed.
The standard deviation and the mean of the i'th variable is the mean and
standard deviation of values of the i'th component of the sample vectors.
"""
function gaussianloglikelihoodbaserate(x::Matrix{Float64})
   nsamples, nvariables = size(x)
   sigmasq = vec(var(x, dims = 1))
   mu = vec(mean(x, dims = 1))
   loglikelihood = 0.0
   for j = 1:nsamples
      loglikelihood -= sum((vec(x[j,:]) - mu).^2 ./ sigmasq)
   end
   loglikelihood /= nsamples
   loglikelihood -= log(2*pi) * nvariables + sum(log.(sigmasq))
   loglikelihood /= 2
   loglikelihood
end


text:
function hfun_bar(vname)
  val = Meta.parse(vname[1])
  return round(sqrt(val), digits=2)
end

function hfun_m1fill(vname)
  var = vname[1]
  return pagevar("index", var)
end

function lx_baz(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return uppercase(brace_content)
end

"""
    lx_readhtml(com, _)

Embed a Pluto notebook via:
https://github.com/rikhuijzer/PlutoStaticHTML.jl
"""
function lx_readhtml(com, _)
    file = string(Franklin.content(com.braces[1]))::String
    dir = joinpath("notebooks")
    html_path = joinpath(dir, "$file.html")
    jl_path = joinpath(dir, "$file.jl")

    return """
        ```julia:pluto
        # hideall

        using PlutoStaticHTML: notebook2html

        path = "$html_path"
        html = read(path, String)
        println("~~~\n\$html\n~~~\n")
        println("_To run this blog post locally, open [this notebook](/$jl_path) with Pluto.jl._")
        ```
        \\textoutput{pluto}
        """
end

ext: jl


text:
using Compose, Viznet
using Test

@testset "sq lattice" begin
    lt = SquareLattice(10, 20)
    @test size(lt, 1) == 10
    @test bravais_size(lt) == (10, 20)
    @test size(lt, 2) == 20
    @test size(lt) == (10, 20)
    @test all(lt[1,1] .≈ (unit(lt) * (0.5), unit(lt) * (0.5)))
    lt = SquareLattice(10, 10)
    @test all(lt[end,end] .≈ (1-unit(lt) * (0.5), 1-unit(lt) * (0.5)))
end

@testset "bonds" begin
    lt = SquareLattice(4, 5)
    b1s = bonds(lt; nth=1)
    @test length(b1s) == 31
    for b in b1s
        @test isconnected(lt, b...)
    end
    @test length(bonds(lt; nth=2)) == 24
    @test length(bonds(lt; nth=3)) == 22
    @test length(bonds(lt; nth=4)) == 34
end


text:
#############################
##  Numerical experiments associated with section 5.3
#############################
## This file contains the numerical experiments used for the highly oscillatory Bessel kernel example,
## which was originally featured by Nick Hale in https://doi.org/10.1093/imanum/dry042
##

using ApproxFun, MultivariateOrthogonalPolynomials, BandedMatrices, BlockBandedMatrices, SpecialFunctions, Plots
using SparseVolterraExamples

####
## First we define the Kernel and g(x) as in section 5.3.
## As mentioned in the paper, it is sensible to use pre-existing implementations to compute the Bessel functions
Kfun(x,y)=besselj(μ,ω*(x-y))
g(x) = besselj(μ+ν,ω*x)+1/(2*x^2)*((ν-1)*(ν-2)*besselj(ν-1,ω*x)+(ν+1)*(ν+2)*besselj(ν+1,ω*x))

#####
## Solver function with step-by-step explanation for given polynomial degree 'n'
## Note the factor 10^(-3) placed as discussed in the paper.
function solveSec54(n,g,Kfun,ν,μ,ω)
    gF = Fun(x->g(x),Jacobi(2,3, 0..1))
    V = triVolterraFullKernelOpP01(Kfun,n+2,true,339)
        V = reflectPabtoPba(n+2)*WLoweringP01P00(n+2)*V[1:n+2,1:n+2]
        V = Conversion(Jacobi(0,0,0..1),Jacobi(2,3,0..1))[1:n+2,1:n+2]*V[1:n+2,1:n+2]
        V = 10^(-3)*Derivative(Jacobi(0,1,0..1),2)[1:n+2,1:n+2]+ω^2*Conversion(Jacobi(0,1,0..1),Jacobi(2,3,0..1))[1:n+2,1:n+2]+ω*V
    coeff = [DirectEvalLHSP10at0(n+2);DirectEvalLHSP10atPrime0(n+2)*Derivative(Jacobi(0,1,0..1),1)[1:n+2,1:n+2];V[1:n,1:n+2]] \ [0.0;0.0;pad(gF.coefficients,n)]
    return Fun(Jacobi(0,1, 0..1), coeff)
end

####
## Compute the solution with same parameters as discussed by Hale and in our paper
## This correponds to computing the solution in Figure 7(a).
ν = 3
μ = 2
ω = 20
u = solveSec54(2000,g,Kfun,ν,μ,ω)
plot(u,color=:black, legend=:false, xlabel = "x", ylabel = "u(x)", legendfontsize=12, tickfontsize=10, thickness_scaling = 1 , grid=:none)

####
## A basic loop to compute errors between the above high order solution and increasingly accurate approximations
function errorvals(max,step,u)
    errorvec = []
    for n=10:step:max
        diff = maximum(abs.([maximum(solveSec54(n,g,Kfun,ν,μ,ω)-u),minimum(solveSec54(n,g,Kfun,ν,μ,ω)-u)]))
        errorvec = push!(errorvec,diff)
    end
    return errorvec
end
####
## Plot successive errors to produce something like Figure 7(b). The exact shape of the plot depends on how dense we make the steps in n.
## Note that for low polynomial orders, since the order is too low to resolve the oscillations, the error varies a lot.
## Only once the polynomial order can resolve the oscillations do we see sensible convergence.
errorvec = errorvals(450,40,u)
plot(Array(10:40:450),errorvec,color=:black, legend=:bottomleft , yscale=:log10 , xlabel = "n", label="sparse method" , ylabel = "error", legendfontsize=12, tickfontsize=10, thickness_scaling = 1.2 , grid=:none)

max_stars_repo_path: examples/Section 5-3 - Bessel kernels with highly oscillatory solutions.jl


text:
export getSystemMatrix, getSystemMatrixReshaped, calculateSystemMatrixSNR

"""
  getSystemMatrix(f, [neglectBGFrames]; kargs...) => Array{ComplexF32,4}

Load the system matrix in frequency domain

Supported keyword arguments:
* frequencies
* bgCorrection
* loadasreal
"""
function getSystemMatrix(f::MPIFile, frequencies=nothing;
                         bgCorrection=false, loadasreal=false,
                         tfCorrection=rxHasTransferFunction(f),
                         numPeriodAverages=1, numPeriodGrouping=1, kargs...)

  if frequencies == nothing
    frequencies = 1:rxNumFrequencies(f,numPeriodGrouping)*rxNumChannels(f)
  end

  if measIsFastFrameAxis(f) && measIsFourierTransformed(f)
    # This is the fast path if the SM lays on disk in the required format
    data = systemMatrix(f, frequencies, bgCorrection)
  else
    data = reshape(getMeasurementsFD(f, frequencies=frequencies, sortFrames=true, bgCorrection=bgCorrection,
              spectralLeakageCorrection=false, transposed=true, tfCorrection=false,
              numPeriodAverages=numPeriodAverages, numPeriodGrouping=numPeriodGrouping), Val(2))
  end

  S = map(ComplexF32, data)

  if tfCorrection && !measIsTFCorrected(f)
    tf = rxTransferFunction(f)
    if tf != nothing
      _corrTFSF(S,tf[rowsToSubsampledRows(f,frequencies)])
      @info "System matrix has been corrected with a Transfer Function. Name of TF: $(rxTransferFunctionFileName(f))"
    else
      error("TF not available")
    end
  end

  if loadasreal
    return converttoreal(S)
  else
    return S
  end
end

function _corrTFSF(S,tf)
  for l=1:size(S,2)
    S[:,l] ./= tf[l]
  end
end

function getSystemMatrixReshaped(f::MPIFile; kargs...)
  return reshape(getSystemMatrix(f;kargs...),:,rxNumFrequencies(f),
                            rxNumChannels(f),acqNumPeriodsPerFrame(f))
end

function calculateSystemMatrixSNR(f::MPIFile)
  data = systemMatrixWithBG(f)
  SNR = calculateSystemMatrixSNR(f, data)
  return SNR
end

function calculateSystemMatrixSNR(f::MPIFile, S::Array; numPeriodAverages=1, numPeriodGrouping=1)
  J = div(acqNumPeriodsPerFrame(f),numPeriodAverages*numPeriodGrouping)
  R = rxNumChannels(f)
  K = rxNumFrequencies(f,numPeriodGrouping)
  N = acqNumFGFrames(f)

  SNR = zeros(K, R, J)

  calculateSystemMatrixSNRInner(S, SNR, J, R, K, N)
  return SNR
end

function calculateSystemMatrixSNRInner(S, SNR, J, R, K, N)
  for j=1:J
    for r=1:R
      for k=1:K
        SBG = S[(N+1):end,k,r,j]
        SFG = S[1:N,k,r,j]
        diffBG = diff(SBG)
        meanBG = mean(SBG)
        signal = maximum(abs.(SFG.-meanBG))
        #noise = mean(abs.(SFG.-meanBG))
        noise = mean(abs.(diffBG))
        SNR[k,r,j] = signal / noise
      end
    end
  end
  SNR[:,:,:] .= mean(SNR,dims=3)
  return
end
#=
function calculateSNRCustomSF(f::BrukerFile,fgFrames::Array,bgFramesFull::Array,bgFramesHalf::Array)
  SNR = zeros(rxNumFrequencies(f),rxNumChannels(f),1)
  for j=1:1
    for r=1:rxNumChannels(f)
      for k=1:rxNumFrequencies(f)
        meanBG = mean(bgFramesFull[:,k,r,j])
        signal = maximum(abs.(fgFrames[:,k,r,j].-meanBG))
	meanBGHalf = mean(bgFramesHalf[:,k,r,j])
        noise = sqrt(var(bgFramesFull[:,k,r,j].-meanBGHalf))#[:,k,r,j]))
        SNR[k,r,j] = signal / noise
      end
    end
  end
  SNR[:,:,:] .= mean(SNR,dims=3)
  return SNR
end
=#
function calculateSNRCustomSF(f::BrukerFile,fgFrames::Array,bgFramesFull::Array,bgFramesHalf::Array)
  SNR = zeros(rxNumFrequencies(f),rxNumChannels(f),1)
  for j=1:1
    for r=1:rxNumChannels(f)
      for k=1:rxNumFrequencies(f)
        diffBG = diff(bgFramesFull[:,k,r,j])
        meanBG = mean(bgFramesFull[:,k,r,j])
        signal = maximum(abs.(fgFrames[:,k,r,j].-meanBG))
        noise = mean(abs.(diffBG))
        SNR[k,r,j] = signal / noise
      end
    end
  end
  SNR[:,:,:] .= mean(SNR,dims=3)
  return SNR
end

function converttoreal(S::AbstractArray{Complex{T},2}) where {T}
  N = size(S,1)
  M = size(S,2)
  S = reshape(reinterpret(T,vec(S)),(2*N,M))
  for l=1:M
    tmp = S[:,l]
    S[1:N,l] = tmp[1:2:end]
    S[N+1:end,l] = tmp[2:2:end]
  end
  return reshape(S,(N,2*M))
end


text:
using HDF5,JLD
const ROOT="ROOT"

# Construct mapping from strings to positive integer values
typealias StrSet Set{UTF8String}
form_hash = StrSet()
postag_hash = StrSet()
deprel_hash = StrSet()

for l in eachline(STDIN)
    l == "\n" && continue
    (form, postag, head, deprel, wvec) = split(chomp(l), '\t')
    push!(form_hash, form)
    push!(postag_hash, postag)
    deprel != ROOT && push!(deprel_hash, deprel)
end

save(ARGS[1],
     "form", sort(collect(form_hash)),
     "postag", sort(collect(postag_hash)),
     "deprel", sort(collect(deprel_hash)))


text:
@testset "core" begin
    @testset "name" begin
        @system SSystemName(Controller)
        s = instance(SSystemName)
        @test Cropbox.namefor(s) == :SSystemName
        @test Cropbox.namefor(SSystemName) == :SSystemName
        S = typeof(s)
        @test startswith(string(S), "var\"##_SSystemName#")
        @test S <: SSystemName
        @test Cropbox.namefor(S) == :SSystemName
    end

    @testset "names" begin
        @test names(System) == [:System]
        @test names(Controller) == [:Controller]
        @test names(Context) == [:Context]
        @test names(Clock) == [:Clock]
    end

    @testset "iteration" begin
        @system SSystemCollect(Controller) begin
            a => 1 ~ preserve
            b(a) => 2a ~ track
        end
        s = instance(SSystemCollect)
        @test length(s) == 4
        @test collect(s) == [s.context, s.config, s.a, s.b]
        @test s[:a] == s.a
        @test s["a"] == s.a
        @test s."a" == s.a
        @test s["context.clock"] == s.context.clock
        @test s."context.clock" == s.context.clock
        @test s[nothing] == s
        @test s[""] == s
    end

    @testset "setvar!" begin
        @system SSystemSetVar(Controller) begin
            a => 1 ~ preserve(ref)
            b => 2 ~ preserve(ref)
        end
        s = instance(SSystemSetVar)
        @test s.a' == 1 && s.b' == 2
        a, b = s.a[], s.b[]
        Cropbox.setvar!(s, :b, a)
        Cropbox.setvar!(s, :a, b)
        @test s.a[] === b && s.b[] === a
        @test s.a' == 2 && s.b' == 1
    end

    @testset "value" begin
        @system SSystemValue(Controller) begin
            a => 1 ~ preserve
            b(a) => 2a ~ preserve
            c(a, b) => a + b ~ track
            d(c) ~ accumulate
        end
        s = instance(SSystemValue)
        @test s.a' == 1 && s.b' == 2 && s.c' == 3
        @test Cropbox.value(s, :b; a=2) == 4
        @test Cropbox.value(s, :c; a=0, b=1) == 1
        @test_throws ErrorException Cropbox.value(s, :c; a=0)
        @test_throws AssertionError Cropbox.value(s, :d; c=0)
    end
end


text:
# This file contains high level Measurement routines based on the low level
# send/receive routines

export measureTransferFunction

# TODO: Currently only one channel is recorded -> measure input channel and
#       relate the output to that

function measureTransferFunction(rp::RedPitaya, freqs)
  tf = zeros(Complex128, length(freqs))

  for (i,freq) in enumerate(freqs)
    dec = optimalDecimation(rp,freq)

    numPeriods = 5
    freqR = roundFreq(rp,dec,freq)
    numSampPerPeriod = numSamplesPerPeriod(rp,dec,freqR)
    numSamp = numSampPerPeriod*numPeriods

    println("Frequency = $freqR Hz")
    println("Number Sampling Points per Period: $numSampPerPeriod")

    # start sending
    send(rp,"GEN:RST")
    sendAnalogSignal(rp,1,"SINE",freqR,0.4)#, numPeriods*2)

    # receive data
    #u1 = receiveAnalogSignal(rp, 1, 0, numSamp, dec=dec, delay=0.2) # NO TRIGGER VERSION
    trigger = "NOW" #"AWG_PE" #"CH1_PE"
    uMeas, uRef = receiveAnalogSignalWithTrigger(rp, 0, 0, numSamp, dec=dec, delay=0.1, typ="OLD",
                                trigger=trigger, triggerLevel=-0.1,
                                triggerDelay=numSampPerPeriod)

    #tf[i] = rfft(uMeas)[numPeriods+1] / length(uMeas)
    tf[i] = rfft(uMeas)[numPeriods+1] / rfft(uRef)[numPeriods+1]
  end
  return tf
end


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : true
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 64 - 1
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 64 - 1
const nz          = haskey(ENV, "NZ"         ) ? parse(Int , ENV["NZ"]         ) : 64 - 1
###
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3)
else
    @init_parallel_stencil(Threads, Float64, 3)
end
using ImplicitGlobalGrid, Plots, Printf, Statistics, LinearAlgebra, MAT
import MPI
# Global reductions
mean_g(A)    = (mean_l = mean(A); MPI.Allreduce(mean_l, MPI.SUM, MPI.COMM_WORLD)/MPI.Comm_size(MPI.COMM_WORLD))
norm_g(A)    = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))
maximum_g(A) = (max_l  = maximum(A); MPI.Allreduce(max_l,  MPI.MAX, MPI.COMM_WORLD))
minimum_g(A) = (min_l  = minimum(A); MPI.Allreduce(min_l,  MPI.MIN, MPI.COMM_WORLD))
# CPU functions
@views av_xza(A) = (A[1:end-1,:,1:end-1] .+ A[1:end-1,:,2:end] .+ A[2:end,:,1:end-1] .+ A[2:end,:,2:end]).*0.25
@views av_zi(A)  = (A[2:end-1,2:end-1,2:end-2] .+ A[2:end-1,2:end-1,3:end-1]).*0.5
@views av_za(A)  = (A[:,:,1:end-1] .+ A[:,:,2:end]).*0.5
@views inn(A)    =  A[2:end-1,2:end-1,2:end-1]

@parallel function smooth!(A2::Data.Array, A::Data.Array, fact::Data.Number)
    @inn(A2)   = @inn(A) + 1.0/6.1/fact*(@d2_xi(A) + @d2_yi(A) + @d2_zi(A))
    return
end

@parallel function compute_maxloc!(Musτ2::Data.Array, Musτ::Data.Array)
    @inn(Musτ2) = @maxloc(Musτ)
    return
end

macro Mu_eff() esc(:(1.0/(1.0/@all(Musτ) + 1.0/(G*dt)))) end
@parallel function compute_iter_params!(dτ_Rho::Data.Array, Gdτ::Data.Array, Musτ::Data.Array, Vpdτ::Data.Number, G::Data.Number, dt::Data.Number, Re::Data.Number, r::Data.Number, max_lxyz::Data.Number)
    @all(dτ_Rho) = Vpdτ*max_lxyz/Re/@Mu_eff()
    @all(Gdτ)    = Vpdτ^2/@all(dτ_Rho)/(r+2)
    return
end

@parallel function assign_τ!(τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array,  τxz::Data.Array,  τyz::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τzz_o::Data.Array, τxy_o::Data.Array, τxz_o::Data.Array, τyz_o::Data.Array)
    @all(τxx_o) = @all(τxx)
    @all(τyy_o) = @all(τyy)
    @all(τzz_o) = @all(τzz)
    @all(τxy_o) = @all(τxy)
    @all(τxz_o) = @all(τxz)
    @all(τyz_o) = @all(τyz)
    return
end

@parallel function compute_P!(∇V::Data.Array, Pt::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Gdτ::Data.Array, r::Data.Number, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(∇V)  = @d_xa(Vx)/dx + @d_ya(Vy)/dy + @d_za(Vz)/dz
    @all(Pt)  = @all(Pt) - r*@all(Gdτ)*@all(∇V)
    return
end

macro inn_yz_Gr() esc(:( @inn_yz(Gdτ)/(G*dt) )) end
macro inn_xz_Gr() esc(:( @inn_xz(Gdτ)/(G*dt) )) end
macro inn_xy_Gr() esc(:( @inn_xy(Gdτ)/(G*dt) )) end
macro av_xyi_Gr() esc(:( @av_xyi(Gdτ)/(G*dt) )) end
macro av_xzi_Gr() esc(:( @av_xzi(Gdτ)/(G*dt) )) end
macro av_yzi_Gr() esc(:( @av_yzi(Gdτ)/(G*dt) )) end
@parallel function compute_τ!(τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, τxx_o::Data.Array, τyy_o::Data.Array, τzz_o::Data.Array, τxy_o::Data.Array, τxz_o::Data.Array, τyz_o::Data.Array, Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, Mus::Data.Array, Gdτ::Data.Array, G::Data.Number, dt::Data.Number, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(τxx) = (@all(τxx) + @all(τxx_o)*@inn_yz_Gr() + 2.0*@inn_yz(Gdτ)*(@d_xi(Vx)/dx))/(1.0 + @inn_yz(Gdτ)/@inn_yz(Mus) + @inn_yz_Gr())
    @all(τyy) = (@all(τyy) + @all(τyy_o)*@inn_xz_Gr() + 2.0*@inn_xz(Gdτ)*(@d_yi(Vy)/dy))/(1.0 + @inn_xz(Gdτ)/@inn_xz(Mus) + @inn_xz_Gr())
    @all(τzz) = (@all(τzz) + @all(τzz_o)*@inn_xy_Gr() + 2.0*@inn_xy(Gdτ)*(@d_zi(Vz)/dz))/(1.0 + @inn_xy(Gdτ)/@inn_xy(Mus) + @inn_xy_Gr())
    @all(τxy) = (@all(τxy) + @all(τxy_o)*@av_xyi_Gr() + 2.0*@av_xyi(Gdτ)*(0.5*(@d_yi(Vx)/dy + @d_xi(Vy)/dx))) /(1.0 + @av_xyi(Gdτ)/@av_xyi(Mus) + @av_xyi_Gr())
    @all(τxz) = (@all(τxz) + @all(τxz_o)*@av_xzi_Gr() + 2.0*@av_xzi(Gdτ)*(0.5*(@d_zi(Vx)/dz + @d_xi(Vz)/dx))) /(1.0 + @av_xzi(Gdτ)/@av_xzi(Mus) + @av_xzi_Gr())
    @all(τyz) = (@all(τyz) + @all(τyz_o)*@av_yzi_Gr() + 2.0*@av_yzi(Gdτ)*(0.5*(@d_zi(Vy)/dz + @d_yi(Vz)/dy))) /(1.0 + @av_yzi(Gdτ)/@av_yzi(Mus) + @av_yzi_Gr())
    return
end

@parallel function compute_dV!(dVx::Data.Array, dVy::Data.Array, dVz::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, dτ_Rho::Data.Array, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(dVx) = (@d_xa(τxx)/dx + @d_ya(τxy)/dy + @d_za(τxz)/dz - @d_xi(Pt)/dx)*@av_xi(dτ_Rho)
    @all(dVy) = (@d_ya(τyy)/dy + @d_xa(τxy)/dx + @d_za(τyz)/dz - @d_yi(Pt)/dy)*@av_yi(dτ_Rho)
    @all(dVz) = (@d_za(τzz)/dz + @d_xa(τxz)/dx + @d_ya(τyz)/dy - @d_zi(Pt)/dz)*@av_zi(dτ_Rho)
    return
end

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, Vz::Data.Array, dVx::Data.Array, dVy::Data.Array, dVz::Data.Array)
    @inn(Vx) = @inn(Vx) + @all(dVx)
    @inn(Vy) = @inn(Vy) + @all(dVy)
    @inn(Vz) = @inn(Vz) + @all(dVz)
    return
end

@parallel function compute_Res!(Rx::Data.Array, Ry::Data.Array, Rz::Data.Array, Pt::Data.Array, τxx::Data.Array, τyy::Data.Array, τzz::Data.Array, τxy::Data.Array, τxz::Data.Array, τyz::Data.Array, dx::Data.Number, dy::Data.Number, dz::Data.Number)
    @all(Rx)    = @d_xa(τxx)/dx + @d_ya(τxy)/dy + @d_za(τxz)/dz - @d_xi(Pt)/dx
    @all(Ry)    = @d_ya(τyy)/dy + @d_xa(τxy)/dx + @d_za(τyz)/dz - @d_yi(Pt)/dy
    @all(Rz)    = @d_za(τzz)/dz + @d_xa(τxz)/dx + @d_ya(τyz)/dy - @d_zi(Pt)/dz
    return
end

@parallel_indices (iy,iz) function bc_x!(A::Data.Array)
    A[  1, iy,  iz] = A[    2,   iy,   iz]
    A[end, iy,  iz] = A[end-1,   iy,   iz]
    return
end

@parallel_indices (ix,iz) function bc_y!(A::Data.Array)
    A[ ix,  1,  iz] = A[   ix,    2,   iz]
    A[ ix,end,  iz] = A[   ix,end-1,   iz]
    return
end

@parallel_indices (ix,iy) function bc_z!(A::Data.Array)
    A[ ix,  iy,  1] = A[   ix,   iy,    2]
    A[ ix,  iy,end] = A[   ix,   iy,end-1]
    return
end

@views function Stokes3D_()
    # Physics
    lx, ly, lz = 10.0, 10.0, 10.0  # domain extends
    ξ          = 1.0               # Maxwell relaxation time
    μs0        = 1.0               # matrix viscosity
    μsi        = 1e-3              # inclusion viscosity
    G          = 1.0               # elastic shear modulus
    εbg        = 1.0               # background strain-rate
    dt         = μs0/(G*ξ)
    # Numerics
    nt         = 5                 # number of time steps
    iterMax    = 2e5               # maximum number of pseudo-transient iterations
    nout       = 200               # error checking frequency
    Re         = 6.0π              # Reynolds number
    r          = 1.0               # Bulk to shear elastic modulus ratio
    CFL        = 0.8/sqrt(3)       # CFL number
    ε          = 1e-8              # nonlinear absolute tolerence
    # nx, ny, nz = 127, 127, 127     # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    b_width    = (8, 4, 4)         # boundary width for comm/comp overlap
    # Derived numerics
    me, dims   = init_global_grid(nx, ny, nz) # MPI initialisation
    @static if USE_GPU select_device() end    # select one GPU per MPI local rank (if >1 GPU per node)
    dx, dy, dz = lx/nx_g(), ly/ny_g(), lz/nz_g()      # cell sizes
    max_lxyz   = max(lx,ly,lz)
    Vpdτ       = min(dx,dy,dz)*CFL
    xc, yc, zc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny), LinRange(dz/2, lz-dz/2, nz)
    # Array allocations
    Pt         = @zeros(nx  ,ny  ,nz  )
    RhoG       = @zeros(nx  ,ny  ,nz  )
    dτ_Rho     = @zeros(nx  ,ny  ,nz  )
    Gdτ        = @zeros(nx  ,ny  ,nz  )
    ∇V         = @zeros(nx  ,ny  ,nz  )
    Vy         = @zeros(nx  ,ny+1,nz  )
    τxx        = @zeros(nx  ,ny-2,nz-2)
    τyy        = @zeros(nx-2,ny  ,nz-2)
    τzz        = @zeros(nx-2,ny-2,nz  )
    τxy        = @zeros(nx-1,ny-1,nz-2)
    τxz        = @zeros(nx-1,ny-2,nz-1)
    τyz        = @zeros(nx-2,ny-1,nz-1)
    τxx_o      = @zeros(nx  ,ny-2,nz-2)
    τyy_o      = @zeros(nx-2,ny  ,nz-2)
    τzz_o      = @zeros(nx-2,ny-2,nz  )
    τxy_o      = @zeros(nx-1,ny-1,nz-2)
    τxz_o      = @zeros(nx-1,ny-2,nz-1)
    τyz_o      = @zeros(nx-2,ny-1,nz-1)
    dVx        = @zeros(nx-1,ny-2,nz-2)
    dVy        = @zeros(nx-2,ny-1,nz-2)
    dVz        = @zeros(nx-2,ny-2,nz-1)
    Rx         = @zeros(nx-1,ny-2,nz-2)
    Ry         = @zeros(nx-2,ny-1,nz-2)
    Rz         = @zeros(nx-2,ny-2,nz-1)
    Mus2       = @zeros(nx  ,ny  ,nz  )
    Musτ       = @zeros(nx  ,ny  ,nz  )
    # Initial conditions
    Rad2       =  zeros(nx  ,ny  ,nz  )
    Vx         =  zeros(nx+1,ny  ,nz  )
    Vz         =  zeros(nx  ,ny  ,nz+1)
    Rad2      .= [(x_g(ix,dx,Rad2) +0.5*dx -0.5*lx)^2 + (y_g(iy,dy,Rad2) +0.5*dy -0.5*ly)^2 + (z_g(iz,dz,Rad2) +0.5*dz -0.5*lz)^2 for ix=1:size(Rad2,1), iy=1:size(Rad2,2), iz=1:size(Rad2,3)]
    Vx         = Data.Array( -εbg.*[(x_g(ix,dx,Vx) +0.5*dx -0.5*lx) for ix=1:size(Vx,1), iy=1:size(Vx,2), iz=1:size(Vx,3)] )
    Vz         = Data.Array(  εbg.*[(z_g(iz,dz,Vz) +0.5*dz -0.5*lz) for ix=1:size(Vz,1), iy=1:size(Vz,2), iz=1:size(Vz,3)] )
    Mus        = μs0*ones(nx,ny,nz)
    Mus[Rad2.<1.0] .= μsi
    Mus        = Data.Array( Mus )
    Mus2      .= Mus
    for ism=1:10#15
        @hide_communication b_width begin # communication/computation overlap
            @parallel smooth!(Mus2, Mus, 1.0)
            Mus, Mus2 = Mus2, Mus
            update_halo!(Mus)
        end
    end
    Musτ      .= Mus
    @hide_communication b_width begin # communication/computation overlap
        @parallel compute_maxloc!(Musτ, Mus)
        @parallel (1:size(Musτ,2), 1:size(Musτ,3)) bc_x!(Musτ)
        @parallel (1:size(Musτ,1), 1:size(Musτ,3)) bc_y!(Musτ)
        @parallel (1:size(Musτ,1), 1:size(Musτ,2)) bc_z!(Musτ)
        update_halo!(Musτ)
    end
    len_Rx_g   = ((nx-2-1)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    len_Ry_g   = ((nx-2-2)*dims[1]+2)*((ny-2-1)*dims[2]+2)*((nz-2-2)*dims[3]+2)
    len_Rz_g   = ((nx-2-2)*dims[1]+2)*((ny-2-2)*dims[2]+2)*((nz-2-1)*dims[3]+2)
    len_∇V_g   = ((nx-2  )*dims[1]+2)*((ny-2  )*dims[2]+2)*((nz-2  )*dims[3]+2)
    # Preparation of visualisation
    if do_viz || do_save_viz
        if (me==0) ENV["GKSwstype"]="nul"; if do_viz !ispath("../../figures") && mkdir("../../figures") end; end
        nx_v, ny_v, nz_v = (nx-2)*dims[1], (ny-2)*dims[2], (nz-2)*dims[3]
        if (nx_v*ny_v*nz_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        Pt_v   = zeros(nx_v, ny_v, nz_v) # global array for visu
        Vz_v   = zeros(nx_v, ny_v, nz_v)
        Rz_v   = zeros(nx_v, ny_v, nz_v)
        Mus_v  = zeros(nx_v, ny_v, nz_v)
        τxz_v  = zeros(nx_v, ny_v, nz_v)
        Pt_inn = zeros(nx-2, ny-2, nz-2) # no halo local array for visu
        Vz_inn = zeros(nx-2, ny-2, nz-2)
        Rz_inn = zeros(nx-2, ny-2, nz-2)
        Mus_inn= zeros(nx-2, ny-2, nz-2)
        τxz_inn= zeros(nx-2, ny-2, nz-2)
        y_sl2, y_sl = Int(ceil((ny_g()-2)/2)), Int(ceil(ny_g()/2))
        Xi_g, Zi_g  = dx+dx/2:dx:(lx-dx-dx/2), dz+dz/2:dz:(lz-dz-dz/2) # inner points only
    end
    # Time loop
    @parallel compute_iter_params!(dτ_Rho, Gdτ, Musτ, Vpdτ, G, dt, Re, r, max_lxyz)
    t=0.0; ittot=0; evo_t=[]; evo_τzz=[]
    for it = 1:nt
        err=2*ε; iter=0; err_evo1=[]; err_evo2=[]
        @parallel assign_τ!(τxx, τyy, τzz, τxy, τxz, τyz, τxx_o, τyy_o, τzz_o, τxy_o, τxz_o, τyz_o)
        # Pseudo-transient iteration
        while err > ε && iter <= iterMax
            if (it==1 && iter==11)  tic()  end
            @parallel compute_P!(∇V, Pt, Vx, Vy, Vz, Gdτ, r, dx, dy, dz)
            @parallel compute_τ!(τxx, τyy, τzz, τxy, τxz, τyz, τxx_o, τyy_o, τzz_o, τxy_o, τxz_o, τyz_o, Vx, Vy, Vz, Mus, Gdτ, G, dt, dx, dy, dz)
            @parallel compute_dV!(dVx, dVy, dVz, Pt, τxx, τyy, τzz, τxy, τxz, τyz, dτ_Rho, dx, dy, dz)
            @hide_communication b_width begin # communication/computation overlap
                @parallel compute_V!(Vx, Vy, Vz, dVx, dVy, dVz)
                @parallel (1:size(Vy,2), 1:size(Vy,3)) bc_x!(Vy)
                @parallel (1:size(Vz,2), 1:size(Vz,3)) bc_x!(Vz)
                @parallel (1:size(Vx,1), 1:size(Vx,3)) bc_y!(Vx)
                @parallel (1:size(Vz,1), 1:size(Vz,3)) bc_y!(Vz)
                @parallel (1:size(Vx,1), 1:size(Vx,2)) bc_z!(Vx)
                @parallel (1:size(Vy,1), 1:size(Vy,2)) bc_z!(Vy)
                update_halo!(Vx, Vy, Vz)
            end
            iter += 1
            if iter % nout == 0
                @parallel compute_Res!(Rx, Ry, Rz, Pt, τxx, τyy, τzz, τxy, τxz, τyz, dx, dy, dz)
                Vmin, Vmax = minimum_g(Vx), maximum_g(Vx)
                Pmin, Pmax = minimum_g(Pt), maximum_g(Pt)
                norm_Rx    = norm_g(Rx)/(Pmax-Pmin)*lx/sqrt(len_Rx_g)
                norm_Ry    = norm_g(Ry)/(Pmax-Pmin)*lx/sqrt(len_Ry_g)
                norm_Rz    = norm_g(Rz)/(Pmax-Pmin)*lx/sqrt(len_Rz_g)
                norm_∇V    = norm_g(∇V)/(Vmax-Vmin)*lx/sqrt(len_∇V_g)
                # norm_Rx = norm_g(Rx)/len_Rx_g; norm_Ry = norm_g(Ry)/len_Ry_g; norm_Rz = norm_g(Rz)/len_Rz_g; norm_∇V = norm_g(∇V)/len_∇V_g
                err = maximum([norm_Rx, norm_Ry, norm_Rz, norm_∇V])
                if isnan(err) error("NaN") end
                push!(err_evo1,maximum([norm_Rx, norm_Ry, norm_Rz, norm_∇V])); push!(err_evo2,iter)
                if (me==0) @printf("Step = %d, iter = %d, err = %1.3e [norm_Rx=%1.3e, norm_Ry=%1.3e, norm_Rz=%1.3e, norm_∇V=%1.3e] \n", it, iter, err, norm_Rx, norm_Ry, norm_Rz, norm_∇V) end
            end
        end
        ittot += iter; t += dt
        push!(evo_t, t); push!(evo_τzz, maximum(τzz))
    end
    # Performance
    wtime    = toc()
    A_eff    = (10*2 +1)/1e9*nx*ny*nz*sizeof(Data.Number) # Effective main memory access per iteration [GB] (Lower bound of required memory access: Te has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    wtime_it = wtime/(ittot-10)                           # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                             # Effective memory throughput [GB/s]
    if (me==0) @printf("Total iters = %d (%d steps), time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", ittot, nt, wtime, round(T_eff, sigdigits=2)) end
    # Visualisation
    if do_viz || do_save_viz
        Pt_inn .= inn(Pt);   gather!(Pt_inn, Pt_v)
        Vz_inn .= av_zi(Vz); gather!(Vz_inn, Vz_v)
        Rz_inn .= av_za(Rz); gather!(Rz_inn, Rz_v)
        Mus_inn.= inn(Mus);  gather!(Mus_inn, Mus_v)
        τxz_inn.= av_xza(τxz); gather!(τxz_inn, τxz_v)
        if me==0 && do_viz
            p1 = heatmap(Xi_g, Zi_g, Pt_v[:,y_sl,:]', aspect_ratio=1, xlims=(Xi_g[1],Xi_g[end]), zlims=(Zi_g[1],Zi_g[end]), c=:viridis, title="Pressure")
            p2 = heatmap(Xi_g, Zi_g, Vz_v[:,y_sl,:]', aspect_ratio=1, xlims=(Xi_g[1],Xi_g[end]), zlims=(Zi_g[1],Zi_g[end]), c=:viridis, title="Vz")
            p4 = heatmap(Xi_g, Zi_g, log10.(abs.(Rz_v[:,y_sl2,:]')), aspect_ratio=1,  xlims=(Xi_g[1],Xi_g[end]), zlims=(Zi_g[1],Zi_g[end]), c=:viridis, title="log10(Rz)")
            #p5 = plot(err_evo2,err_evo1, legend=false, xlabel="# iterations", ylabel="log10(error)", linewidth=2, markershape=:circle, markersize=3, labels="max(error)", yaxis=:log10)
            p3 = plot(evo_t, evo_τzz, legend=false, xlabel="time", ylabel="max(τzz)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
               #plot!(evo_t, 2.0.*εbg.*μs0.*(1.0.-exp.(.-evo_t.*G./μs0)), linewidth=2.0) # analytical solution
            plot(p1, p2, p4, p3)
            savefig("../../figures/Stokes_3D_ve_$(nx_g()).png")
        end
    end
    if me==0 && do_save
        !ispath("../../output") && mkdir("../../output")
        open("../../output/out_Stokes3D_ve.txt","a") do io
            println(io, "$(nx_g()) $(ny_g()) $(nz_g()) $(ittot) $(nt)")
        end
    end
    if me==0 && do_save_viz
        !ispath("../../out_visu") && mkdir("../../out_visu")
        matwrite("../../out_visu/Stokes_3D_ve.mat", Dict("Pt_3D"=> Pt_v, "Mus_3D"=> Mus_v, "Txz_3D"=> τxz_v, "Vz_3D"=> Vz_v, "dx_3D"=> dx, "dy_3D"=> dy, "dz_3D"=> dz); compress = true)
    end
    finalize_global_grid()
    return xc, yc, zc, Pt
end

if use_return
    xc, yc, zc, P = Stokes3D_();
else
    Stokes3D = begin Stokes3D_(); return; end
end


text:
module BenchParallelFindMinMax

using BenchmarkTools
using FLoops

function findminmax(exc, xs)
    @floop exc for (i, x) in pairs(xs)
        @reduce() do (imin; i), (xmin; x)
            if isless(x, xmin)
                xmin = x
                imin = i
            end
        end
        @reduce() do (imax; i), (xmax; x)
            if isless(xmax, x)
                xmax = x
                imax = i
            end
        end
    end
    return ((xmin, imin), (xmax, imax))
end

const SUITE = BenchmarkGroup()

xs = rand(100_000)
SUITE["threaded"] = @benchmarkable findminmax(ThreadedEx(), $xs)
SUITE["sequential"] = @benchmarkable findminmax(SequentialEx(), $xs)

end  # module
BenchParallelFindMinMax.SUITE


text:
function test_pairwise(F::Function, a, b, s::Integer)
    s1 = minhash(F, a, s)
    s2 = minhash(F, b, s)
    @test intersectionlength(s1, s2) == length(intersect(s1.hashes, s2.hashes))
end

function test_multi(F::Function, x, s)
    sketches = [minhash(F, i, s) for i in x]
    m = intersectionlength(sketches)
    allgood = true
    for i in 1:length(sketches) - 1
        for j in i+1:length(sketches)
            tp = intersectionlength(sketches[i], sketches[j])
            allgood &= tp == m[j, i]
        end
    end
    @test allgood
end

@testset "Pairwise" begin
    test_pairwise(hash, 1:100, 60:200, 40)
    test_pairwise(hash, 1:100, 'A':'z', 40)
    test_pairwise(hash, 1:100, 'A':'z', 1)
    test_pairwise(hash, 'A':'g', 'W':'s', 10)
end

@testset "Multi" begin
    xs = [collect(Set(rand(1:250, 100))) for i in 1:5]
    test_multi(hash, xs, 50)
end

@testset "Misc" begin
    s = minhash('A':'z', 38)
    @test length(s) == 38
    s = minhash(5:55, 250)
    @test length(s) == 51
    @test !isempty(s)
    @test isempty(minhash(hash, 1:0, 100))
end


text:
@llvmfunc foo(x::Int32) -> Int32 begin
    @block test begin
        @ret Int32 999
    end
end

ast2 = @LLVM.DSL begin
    @mod begin
        @func foo(bar::Int32) -> Int32 begin
            @block test begin
                @ret @const Int32 999
            end
        end
    end
end

ast1 = Ast.Mod("Factorial", nothing, nothing, [
    Ast.Func(:fac, Int32, (:n::Int32, :y::Int32),
        Ast.BasicBlock(:entry,
            :n <= Ast.ICmp{:EQ}(:n, 1, :test),
            Ast.CondBr(test, :result, :recur)
        ),
        Ast.BasicBlock(:recur,
            :n_1 <= Ast.ISub(n, Int32,),
            Ast.Br(:result)
        ),
        Ast.BasicBlock(:result,
            :fac <= Ast.Phi(Int32, 1, :n_fac_n_1, :fac)
            Ast.Ret(:fac)
        ),
    )
])

@llvmfunc fac(n::Int32, y::Int32) -> Int32 begin

    @r(0) = @type {Int32, Ptr{@r(1)}, Ptr{@r(0)}}
    @r(1) = @type opaque

    @v(0) = @global Int32 1
    @v(1) = @external protected addrspace(3) global Int32, section "foo"
    @v(3) = @unnamed_addr global @r(1)
    @v(4) = @external global Int32[4393948493]

    @v(".test") = @thread_local global Int32 0

    @v(three) = @alias private Int32 addrspace(3) Ptr{@v(1)}
    @v(two)   = @alias Int32 addrspace(3) Ptr{@v(three)}

    @define bar() -> Int32 begin
        @v(1) = @call zeroext Int32 foo(@inreg(1::Int32), @signext(4::Int32)) @m(0)
        @ret Int32 @v(1)
    end

    @block entry begin
        test = @icmp(eq, n, 1)
        @condbr(test, recur, result)
    end

    @block recur begin
        n_1 = @sub(n, 1)
        fac_n_1 = @call(fac, n_1)
        n_fac_n_1 = @mul(n, fac_n_1)
        @br(result)
    end

    @block result begin
        fac = @phi Int32 entry=>1 recur=>n_fac_n_1
        @ret fac
    end

end


text:
using Base.Test
using QuickTest

import Base: ==
import QuickTest: generate_test_value

@testset "unconditioned property test" begin
  @testprop a::Int == a
  @testprop a::Int == a && b::Int == b
  @testprop a::Int == a && b::Float64 == b
  @testprop a == a a::Int
  @testprop a::Int == a a::Int
  @testprop a::Int == a a::Int b::Float64
  @testprop 10 a::Int == a

  @test_throws ArgumentError eval(macroexpand(:( @testprop a::Int == a::Int128 )))
  @test_throws ArgumentError eval(macroexpand(:( @testprop a::Int == a a::Int128 )))
end

@testset "conditioned property test" begin
  @testprop -(a::Int) > 0  a < 0
  @testprop -a > 0  a::Int < 0
  @testprop -a > 0  a < 0  a::Int
  @testprop -a > 0  a::Int < 0  a::Int
  @testprop -a > 0  a::Int < 0  b::Float64
  @testprop 10 -a > 0  a::Int < 0
  @testprop 10 100 -a > 0  a::Int < 0

  @testprop isa(a::Int,Float64) a != a

  @test_throws ArgumentError eval(macroexpand(:( @testprop -a::Int128 > 0  a::Int < 0 )))
  @test_throws ArgumentError eval(macroexpand(:( @testprop -a > 0  a::Int < 0  a::Int128 )))
end

@testset "expressions" begin
  @testprop a+1 == b  a::Int b::Expression{a}(a+1)
  @testprop a+b == c  a::Int b::Int c::Exp{a,b}(a+b)
  @testprop 10 a+1 == b  a::Int b::Expression{a}(a+1)
  @testprop 10 a+b == c  a::Int b::Int c::Exp{a,b}(a+b)
end

type TestType
  a::Int
end
function generate_test_value(::Type{TestType}, size)
  TestType(generate_test_value(Int,size))
end

type ZZmod
  modulus::Int
end
type zz
  parent::ZZmod
  value::Int
end

==(a::ZZmod, b::ZZmod) = a.modulus == b.modulus
==(a::zz, b::zz) = a.parent == b.parent && a.value == b.value

generate_test_value(::Type{ZZmod}, gsize) = ZZmod(rand(1:gsize))
generate_test_value(a::ZZmod, gsize) = zz(a, rand(0:(a.modulus-1)))

@testset "custom types in test properties" begin
  @testprop isa(a::TestType,TestType)
  @testprop 10 (isa(a::TestType,TestType))

  @testprop a == a  p::ZZmod a::Element{p}
  @testprop isa(a,zz) == true  p::ZZmod a::Elem{p}
  @testprop isa(a,Vector{zz}) == true  p::ZZmod a::ElementVector{p}
  @testprop isa(a,Vector{zz}) == true  p::ZZmod a::ElemV{p}

  let
    p = ZZmod(5)
    @testprop a == a  a::Element{p}
    @testprop isa(a,zz) == true  a::Elem{p}
    @testprop isa(a,Vector{zz}) == true  a::ElementVector{p}
    @testprop isa(a,Vector{zz}) == true  a::ElemV{p}
  end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Qt6Declarative"
version = v"6.3.0"

# Set this to true first when updating the version. It will build only for the host (linux musl).
# After that JLL is in the registyry, set this to false to build for the other platforms, using
# this same package as host build dependency.
const host_build = false

# Collection of sources required to build qt6
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtdeclarative-everywhere-src-$version.tar.xz",
                  "d294b029dc2b2d4f65da516fdc3b8088d32643eb7ff77db135a8b9ce904caa37"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/11.0-11.1/MacOSX11.1.sdk.tar.xz",
                  "9b86eab03176c56bb526de30daa50fa819937c54b280364784ce431885341bf6"),
    ArchiveSource("https://sourceforge.net/projects/mingw-w64/files/mingw-w64/mingw-w64-release/mingw-w64-v10.0.0.tar.bz2",
                  "ba6b430aed72c63a3768531f6a3ffc2b0fde2c57a3b251450dcf489a894f0894")
]

script = raw"""
cd $WORKSPACE/srcdir

mkdir build
cd build/
qtsrcdir=`ls -d ../qtdeclarative-*`

case "$target" in

    x86_64-linux-musl*)
        cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_BUILD_TYPE=Release $qtsrcdir
    ;;

    *mingw*)
        cd $WORKSPACE/srcdir/mingw*/mingw-w64-headers
        ./configure --prefix=/opt/$target/$target/sys-root --enable-sdk=all --host=$target
        make install

        cd ../mingw-w64-crt/
        if [ ${target} == "i686-w64-mingw32" ]; then
            _crt_configure_args="--disable-lib64 --enable-lib32"
        elif [ ${target} == "x86_64-w64-mingw32" ]; then
            _crt_configure_args="--disable-lib32 --enable-lib64"
        fi
        ./configure --prefix=/opt/$target/$target/sys-root --enable-sdk=all --host=$target --enable-wildcard ${_crt_configure_args}
        make -j${nproc}
        make install

        cd ../mingw-w64-libraries/winpthreads
        ./configure --prefix=/opt/$target/$target/sys-root --host=$target --enable-static --enable-shared
        make -j${nproc}
        make install

        cd $WORKSPACE/srcdir/build
        cmake -DQT_HOST_PATH=$host_prefix \
            -DPython_ROOT_DIR=/usr \
            -DCMAKE_INSTALL_PREFIX=${prefix} \
            -DCMAKE_PREFIX_PATH=$host_prefix \
            -DCMAKE_FIND_ROOT_PATH=$prefix \
            -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
            -DCMAKE_BUILD_TYPE=Release \
        $qtsrcdir
    ;;

    *apple-darwin*)
        apple_sdk_root=$WORKSPACE/srcdir/MacOSX11.1.sdk
        sed -i "s!/opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root!$apple_sdk_root!" $CMAKE_TARGET_TOOLCHAIN
        cmake -DQT_HOST_PATH=$host_prefix \
        -DPython_ROOT_DIR=/usr \
        -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DCMAKE_PREFIX_PATH=$host_prefix \
        -DCMAKE_FIND_ROOT_PATH=$prefix \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_SYSROOT=$apple_sdk_root -DCMAKE_FRAMEWORK_PATH=$apple_sdk_root/System/Library/Frameworks -DCMAKE_OSX_DEPLOYMENT_TARGET=10.15 \
        -DCMAKE_BUILD_TYPE=Release \
        $qtsrcdir
    ;;

    *)
        cmake -DQT_HOST_PATH=$host_prefix \
        -DPython_ROOT_DIR=/usr \
        -DCMAKE_INSTALL_PREFIX=${prefix} \
        -DCMAKE_PREFIX_PATH=$host_prefix \
        -DCMAKE_FIND_ROOT_PATH=$prefix \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
        -DCMAKE_BUILD_TYPE=Release \
        $qtsrcdir
    ;;

esac

cmake --build . --parallel ${nproc}
cmake --install .
install_license $WORKSPACE/srcdir/qt*-src-*/LICENSE.LGPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
if host_build
    platforms = [Platform("x86_64", "linux",cxxstring_abi=:cxx11,libc="musl")]
else
    platforms = expand_cxxstring_abis(filter(!Sys.isapple, supported_platforms()))
    filter!(p -> arch(p) != "armv6l", platforms) # No OpenGL on armv6
    platforms_macos = [ Platform("x86_64", "macos"), Platform("aarch64", "macos") ]
end

# The products that we will ensure are always built
products = [
    LibraryProduct(["Qt6Qml", "libQt6Qml", "QtQml"], :libqt6qml),
    LibraryProduct(["Qt6QmlModels", "libQt6QmlModels", "QtQmlModels"], :libqt6qmlmodels),
    LibraryProduct(["Qt6QmlCore", "libQt6QmlCore", "QtQmlCore"], :libqt6qmlcore),
    LibraryProduct(["Qt6QmlWorkerScript", "libQt6QmlWorkerScript", "QtQmlWorkerScript"], :libqt6qmlworkerscript),
    LibraryProduct(["Qt6QmlLocalStorage", "libQt6QmlLocalStorage", "QtQmlLocalStorage"], :libqt6qmllocalstorage),
    LibraryProduct(["Qt6QmlXmlListModel", "libQt6QmlXmlListModel", "QtQmlXmlListModel"], :libqt6qmlxmllistmodel),
    LibraryProduct(["Qt6Quick", "libQt6Quick", "QtQuick"], :libqt6quick),
    LibraryProduct(["Qt6QuickLayouts", "libQt6QuickLayouts", "QtQuickLayouts"], :libqt6quicklayouts),
    LibraryProduct(["Qt6QuickTest", "libQt6QuickTest", "QtQuickTest"], :libqt6quicktest),
    LibraryProduct(["Qt6QuickParticles", "libQt6QuickParticles", "QtQuickParticles"], :libqt6quickparticles),
    LibraryProduct(["Qt6QuickShapes", "libQt6QuickShapes", "QtQuickShapes"], :libqt6quickshapes),
    LibraryProduct(["Qt6QuickWidgets", "libQt6QuickWidgets", "QtQuickWidgets"], :libqt6quickwidgets),
    LibraryProduct(["Qt6QuickTemplates2", "libQt6QuickTemplates2", "QtQuickTemplates2"], :libqt6quicktemplates2),
    LibraryProduct(["Qt6QuickControls2Impl", "libQt6QuickControls2Impl", "QtQuickControls2Impl"], :libqt6quickcontrols2impl),
    LibraryProduct(["Qt6QuickControls2", "libQt6QuickControls2", "QtQuickControls2"], :libqt6quickcontrols2),
    LibraryProduct(["Qt6QuickDialogs2Utils", "libQt6QuickDialogs2Utils", "QtQuickDialogs2Utils"], :libqt6quickdialogs2utils),
    LibraryProduct(["Qt6QuickDialogs2QuickImpl", "libQt6QuickDialogs2QuickImpl", "QtQuickDialogs2QuickImpl"], :libqt6quickdialogs2quickimpl),
    LibraryProduct(["Qt6QuickDialogs2", "libQt6QuickDialogs2", "QtQuickDialogs2"], :libqt6quickdialogs2),
    LibraryProduct(["Qt6LabsSettings", "libQt6LabsSettings", "QtLabsSettings"], :libqt6labssettings),
    LibraryProduct(["Qt6LabsQmlModels", "libQt6LabsQmlModels", "QtLabsQmlModels"], :libqt6labsqmlmodels),
    LibraryProduct(["Qt6LabsFolderListModel", "libQt6LabsFolderListModel", "QtLabsFolderListModel"], :libqt6labsfolderlistmodel),
    LibraryProduct(["Qt6LabsAnimation", "libQt6LabsAnimation", "QtLabsAnimation"], :libqt6labsanimation),
    LibraryProduct(["Qt6LabsWavefrontMesh", "libQt6LabsWavefrontMesh", "QtLabsWavefrontMesh"], :libqt6labswavefrontmesh),
    LibraryProduct(["Qt6LabsSharedImage", "libQt6LabsSharedImage", "QtLabsSharedImage"], :libqt6labssharedimage),
]

products_macos = [
    FrameworkProduct("QtQml", :libqt6qml),
    FrameworkProduct("QtQmlModels", :libqt6qmlmodels),
    FrameworkProduct("QtQmlCore", :libqt6qmlcore),
    FrameworkProduct("QtQmlWorkerScript", :libqt6qmlworkerscript),
    FrameworkProduct("QtQmlLocalStorage", :libqt6qmllocalstorage),
    FrameworkProduct("QtQmlXmlListModel", :libqt6qmlxmllistmodel),
    FrameworkProduct("QtQuick", :libqt6quick),
    FrameworkProduct("QtQuickLayouts", :libqt6quicklayouts),
    FrameworkProduct("QtQuickTest", :libqt6quicktest),
    FrameworkProduct("QtQuickParticles", :libqt6quickparticles),
    FrameworkProduct("QtQuickShapes", :libqt6quickshapes),
    FrameworkProduct("QtQuickWidgets", :libqt6quickwidgets),
    FrameworkProduct("QtQuickTemplates2", :libqt6quicktemplates2),
    FrameworkProduct("QtQuickControls2Impl", :libqt6quickcontrols2impl),
    FrameworkProduct("QtQuickControls2", :libqt6quickcontrols2),
    FrameworkProduct("QtQuickDialogs2Utils", :libqt6quickdialogs2utils),
    FrameworkProduct("QtQuickDialogs2QuickImpl", :libqt6quickdialogs2quickimpl),
    FrameworkProduct("QtQuickDialogs2", :libqt6quickdialogs2),
    FrameworkProduct("QtLabsSettings", :libqt6labssettings),
    FrameworkProduct("QtLabsQmlModels", :libqt6labsqmlmodels),
    FrameworkProduct("QtLabsFolderListModel", :libqt6labsfolderlistmodel),
    FrameworkProduct("QtLabsAnimation", :libqt6labsanimation),
    FrameworkProduct("QtLabsWavefrontMesh", :libqt6labswavefrontmesh),
    FrameworkProduct("QtLabsSharedImage", :libqt6labssharedimage),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    HostBuildDependency("Qt6Base_jll"),
    HostBuildDependency("Qt6ShaderTools_jll"),
    Dependency("Qt6Base_jll"),
    Dependency("Qt6ShaderTools_jll"),
]

if !host_build
    push!(dependencies, HostBuildDependency("Qt6Declarative_jll"))
end

include("../../fancy_toys.jl")

@static if !host_build
    if any(should_build_platform.(triplet.(platforms_macos)))
        build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
    end
end
if any(should_build_platform.(triplet.(platforms)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
end