text:
"""
Represents the processed spectral data necessary to efficiently filter-fit one or more unknown spectra.
A `FilterFitPacket{S<:Detector, T<:AbstractFloat}` contains the data necessary to filter the unknown and to
apply pre-filtered references.  The type `T <: AbstractFloat` allows you to specify the bit-depth used to
perform subsequent calculations using this `FilterFitPacket`.  `Float32` uses less memory and is a little
faster than `Float64`.  If there are duplicate `FilteredReference` for an elemental ROI, the preference
is for the first one.  This allows you to fill in unavailable "FilteredReference" elemental ROIs with more
general ones.
"""
struct FilterFitPacket{S<:Detector, T<:AbstractFloat}
    detector::S
    filter::TopHatFilter{T}
    references::Vector{FilteredReference{T}}

    function FilterFitPacket(det::S, filt::TopHatFilter{T}, refs::Vector{FilteredReference{T}})  where {S<:Detector, T<:AbstractFloat}
        # Only permit one FilteredReference for each type of label for a set of X-rays
        filtrefs = FilteredReference{T}[]
        for ref in refs
            exists = false
            for fr in filtrefs
                # If a reference already exists of this type for these X-rays reassign it
                if (typeof(fr.label) == typeof(ref.label))  && (fr.label.xrays==ref.label.xrays)
                    # filtrefs[i] = ref
                    exists = true
                    break
                end
            end
            if !exists
                push!(filtrefs, ref)
            end
        end
        new{S,T}(det, filt, filtrefs)
    end
end # struct

function Base.show(io::IO, ffp::FilterFitPacket)
    xrays = join(map(r -> "\t$(r.label),", ffp.references), "\n")
    print(io, "References[\n\t$(ffp.detector), \n$xrays\n]")
end

"""
    NeXLUncertainties.asa(::Type{DataFrame}, ffp::FilterFitPacket)

Summarize the `FilteredReference` structs within a `FilterFitPacket` as a `DataFrame`.
"""
function NeXLUncertainties.asa(::Type{DataFrame}, ffp::FilterFitPacket)
    function p2b(fref)
        # charonly over roi, data over ffroi
        croi = max(1,first(fref.roi)-first(fref.ffroi)+1):min(length(fref.data), last(fref.roi)-first(fref.ffroi)+1)
        @assert length(croi)==length(fref.roi)
        sum(fref.charonly) / (sum(fref.data[croi])-sum(fref.charonly))
    end
    function s2n(fref)
        croi = max(1,first(fref.roi)-first(fref.ffroi)+1):min(length(fref.data), last(fref.roi)-first(fref.ffroi)+1)
        @assert length(croi)==length(fref.roi)
        sum(fref.charonly) / sqrt(sum(fref.data[croi])-sum(fref.charonly))
    end
    DataFrame(
        :Spectrum => [ name(fr.label.spectrum) for fr in ffp.references ],
        Symbol("Beam Energy (keV)") => [ get(fr.label.spectrum, :BeamEnergy, missing)/1000.0 for fr in ffp.references ],
        Symbol("Probe Current (nA)") => [ get(fr.label.spectrum, :ProbeCurrent, missing) for fr in ffp.references ],
        Symbol("Live Time (s)") => [ get(fr.label.spectrum, :LiveTime, missing) for fr in ffp.references ],
        :Material => [ get(fr.label.spectrum, :Composition, nothing) for fr in ffp.references],
        :Lines => [ fr.label.xrays for fr in ffp.references],
        :ROI => [ fr.roi for fr in ffp.references],
        # :FullROI => [ fr.ffroi for fr in ffp.references],
        Symbol("P-to-B") => p2b.(ffp.references),
        Symbol("S-to-N") => s2n.(ffp.references)
    )
end

"""
    missingReferences(ffp::FilterFitPacket, elms::Vector{Element}, e0::Float64, ampl=1.0e-5)

Returns a `Vector{Tuple{Vector{CharXRay}, UnitRange{Int64}}}` containing the ROIs for which a
`FilteredReference` is missing from the `FilterFitPacket`.
"""
function missingReferences(ffp::FilterFitPacket, elms::Vector{Element}, e0::Float64, ampl=1.0e-5)
    # Find all collections of X-rays for which we'll need references
    allXrays = mapreduce(elm->NeXLSpectrum.labeledextents(elm, ffp.detector, ampl, e0), append!, elms)
    # Remove those for which there are references
    filter(allXrays) do xr
        isnothing(findfirst(ref->ref.label.xrays == xr[1], ffp.references))
    end
end

"""
Returns a tuple containing the unfiltered spectra associated with the references.
"""
spectra(ffp::FilterFitPacket) = (ref.label.spectrum for ref in ffp.references)

"""
A list of elements for which there are filtered references.
"""
NeXLCore.elms(ffp::FilterFitPacket) = union( [ element(ref.label) for ref in filter(r->r.label isa CharXRayLabel, ffp.references) ] )

struct ReferencePacket
    spectrum::Spectrum
    element::Element
    material::Material
end

Base.show(io::IO, rp::ReferencePacket) = print(io, "ReferencePacket[$(symbol(rp.element)), $(name(rp.material)), $(name(rp.spectrum))]")

"""
    reference( elm::Element, spec::Spectrum, mat::Material=spec[:Composition]; pc = nothing, lt = nothing, e0 = nothing, coating = nothing)::ReferencePacket
    reference(els::AbstractVector{Element}, spec::Spectrum, mat::Material = spec[:Composition]; pc = nothing, lt = nothing, e0 = nothing, coating = nothing)::Vector{ReferencePacket}

Construct a `ReferencePacket` from a `Spectrum` collected from the specified `Material` for the specified `Element`.
Often used with `references(...)` to build `FilterFitPacket`s.

Optional named arguments `pc`, `lt`, `e0`, `coating` allow you to specify the probe current, live time, beam energy and
sample coating.
"""
function reference(
    elm::Element,
    spec::Spectrum,
    mat::Material = spec[:Composition];
    pc = nothing,
    lt = nothing,
    e0 = nothing,
    coating = nothing,
)::ReferencePacket
    if !isnothing(lt)
        @assert lt > 0.0 "The live time must be larger than zero for $(spec[:Name])."
        spec[:LiveTime] = lt
    end
    if !isnothing(pc)
        @assert pc > 0.0 "The probe current must be larger than zero for $(spec[:Name])."
        spec[:ProbeCurrent] = pc
    end
    if !isnothing(e0)
        @assert e0 >= 1.0e3 "The beam energy ($(e0/1000.0)) keV is too low for $(spec[:Name])."
        spec[:BeamEnergy] = e0
    end
    if !isnothing(coating)
        @assert coating isa Film
        spec[:Coating] = coating
    end
    spec[:Composition] = mat
    @assert haskey(spec, :LiveTime) "The :LiveTime property must be defined for $(spec[:Name]).  (Use the `lt` keyword argument)"
    @assert haskey(spec, :ProbeCurrent) "The :ProbeCurrent property must be defined for $(spec[:Name]).  (Use the `pc` keyword argument)"
    @assert haskey(spec, :BeamEnergy) "The :BeamEnergy property must be defined for $(spec[:Name]).  (Use the `e0` keyword argument)"
    @assert haskey(mat, elm) "$(Symbol(elm)) is not present in $mat."
    return ReferencePacket(spec, elm, mat)
end

reference(elm::Element, filename::AbstractString, mat::Material; kwargs...) =
    reference(elm, loadspectrum(filename), mat; kwargs...)

reference(elm::Element, filename::AbstractString; kwargs...) =
    reference(elm, loadspectrum(filename); kwargs...)

function reference(
    els::AbstractVector{Element},
    spec::Spectrum,
    mat::Material = spec[:Composition];
    kwargs...)
    map(el->reference(el, spec, mat; kwargs...), els)
end
reference(elm::AbstractVector{Element}, filename::AbstractString, mat::Material; kwargs...) =
    reference(elm, loadspectrum(filename), mat; kwargs...)

"""
    references(refs::AbstractVector{ReferencePacket}, det::EDSDetector)::FilterFitPacket
    references(refs::AbstractVector{ReferencePacket}, fwhm::Float64)::FilterFitPacket

Constructs a FilterFitPacket from a vector of `ReferencePackets`.  Each `ReferencePacket` represents a
single ROI for an element.  It is possible more than one `ReferencePacket` might be defined for an
elemental ROI.  In this case, the `ReferencePacket` with the lower index will take preference over
later ones.  This allows you to fill in only the missing elemental ROIs using spectra collected from
alternative materials.  For example, a spectrum from F₂Fe is suitable for the Fe K-lines but not the
Fe L-lines. So we might specify F₂Fe first to specify the references for the Fe K-lines and then fill
in the L-lines with a spectrum from pure Fe.
"""
function references(
    refs::AbstractVector{ReferencePacket},
    det::EDSDetector;
    ftype::Type{<:AbstractFloat} = Float64
)::FilterFitPacket
    chcount = det.channelcount
    @assert all(length(r.spectrum) == chcount for r in refs) "The number of spectrum channels must match the detector for all spectra."
    @assert length(refs) > 0 "Please provide at least one ReferencePacket in references(...)"
    # Build the top-hat filter for det
    ff = buildfilter(ftype, det)
    # Apply the top-hat filter to all refs. Trying to thread this fails. :-(
    frefs = mapreduce(append!, refs) do ref
        frefs = filterreference(ff, ref.spectrum, ref.element, ref.material)
        length(frefs)==0 && @warn "Unable to create any filtered ROI references for $(ref.element) from $(name(ref.material))."
        frefs
    end
    return FilterFitPacket(det, ff, frefs)
end
references(refs::AbstractVector{ReferencePacket}, fwhm::Float64; ftype::Type{<:AbstractFloat}=Float64) =
    references(refs, matching(first(refs).spectrum, fwhm), ftype=ftype)

fit_spectrum(spec::Spectrum, ffp::FilterFitPacket{S, T}) where { S<:Detector, T<: AbstractFloat } =
    fit_spectrum(FilteredUnknownW{T}, spec, ffp.filter, ffp.references)

fit_spectrum(specs::AbstractVector{<:Spectrum}, ffp::FilterFitPacket{S, T}) where { S<:Detector, T<: AbstractFloat } =
    ThreadsX.map(spec -> fit_spectrum(FilteredUnknownW{T}, spec, ffp.filter, ffp.references), specs)

"""
    fit_spectrum(spec::Spectrum, ffp::FilterFitPacket)::FilterFitResult
    fit_spectrum(specs::AbstractVector{<:Spectrum}, ffp::FilterFitPacket)::Vector{FilterFitResult}

Fit a `Spectrum` or a vector of `Spectrum` using the specified `FilterFitPacket`.  The result is a
`FilterFitResult` structure which contains k-ratios, residuals, etc.

    fit_spectrum(hs::HyperSpectrum, ffp::FilterFitPacket; mode::Symbol=:Fast, zero = x -> max(0.0, x))::Array{KRatios}

  * `mode = :Fast` - Uses precomputed, filtered "vector" fit method.  No uncertainties are available.
  * `mode = :Intermediate` - Uses an optimized full fit without refits for negative k-ratios.
  * `mode = :Full` - Uses the full single spectrum fit algorithm including refitting when one or more k-ratio is less than zero.

Performs a filtered fit on a hyperspectrum returning an `Array{KRatios}`.

Selecting a mode:
  :Fast is good for generating k-ratio maps or exploratory analysis of a k-ratio map. :Full is best when a
  quantitative map of a high count hyperspectrum is desired.  Fit results for major elements are similar for
  all three but differ for minor and trace elements.  Particularly when a k-ratio is slightly negative. This
  negative k-ratio can effect the other k-ratios.  :Fast also works less well when many elements (>>10) (particularly
  interfering elements) are included in the fit. Unfortunately, :Intermediate and :Full slow down when many elements
  are fit - O(n(elements)²).

The following timing on a 512 x 512 x 2048 hyperspectrum fitting 21 elements with 33 distinct ROIs on a fast 6-core
laptop with 16 GiB memory give a relative feel for the speed of each algorithm.  Yes, :Fast is approximately 30 ×
faster than :Intermediate and used almost 80x less memory.  (64-bit references, 32-bit references use about 1/2
the memory and take about 4/5 the time.)

| mode=         | Threads  |Run time (s)   | Allocations  | Allocated (GiB) | GC time  |
|---------------|----------|---------------|--------------|-----------------|----------|
| :Fast         | 6        | 11.4          | 2.11 M       | 4.72            | 3.0%     |
| :Intermediate | 6        | 342.7         | 13.11 M      | 364.4           | 4.2%     |
| :Full         | 6        | 574.6         | 6.2 G        | 862.9           | 5.5%     |
| :Fast         | 1        | 25.5          | 2.1 M        | 4.72            | 1.8%     |
| :Intermediate | 1        | 1064.6        | 13.1 M       | 364.4           | 0.9%     |
| :Full         | 1        | 2186.2        | 6.2 G        | 862.1           | 2.8%     |
"""
function fit_spectrum(
    hs::HyperSpectrum,
    ffp::FilterFitPacket{S, T};
    mode::Symbol = :Fast,
    zero = x -> max(Base.zero(T), x),
    sigma = Base.zero(T)
)::Array{KRatios} where { S <: Detector, T <: AbstractFloat }
    @assert matches(hs[CartesianIndices(hs)[1]], ffp.detector) "The detector for the hyper-spectrum must match the detector for the filtered references."
    @assert (mode==:Fast) || (mode==:Intermediate) || (mode==:Full) "The mode argument must be :Fast, :Intermediate or :Full."
    if mode == :Fast
        vq = VectorQuant(ffp.references, ffp.filter)
        return fit_spectrum(hs, vq, zero)
    elseif mode == :Intermediate
        return fit_spectrum_int(hs, ffp, zero)
    elseif mode == :Full
        return fit_spectrum_full(hs, ffp, sigma)
    end
    # Should never get here...
    return  Array{KRatio}[]
end

function fit_spectrum_int(
    hs::HyperSpectrum,
    ffp::FilterFitPacket{S, T},
    zero::Function
)::Array{KRatios} where { S <: Detector, T <: AbstractFloat }
    unklabel = UnknownLabel(hs)
    function _tophatfilterhs(unklabel, data, thf, scale)
        @assert length(data) <= length(thf) "The reference spectra must have more channels than the hyperspectrum data."
        filtered = [filtereddatum(thf, data, i) for i in eachindex(data)]
        dp = [max(x, one(T)) for x in data] # To ensure covariance isn't zero or infinite precision
        covar = [filteredcovar(thf, dp, i, i) for i in eachindex(data)]
        FilteredUnknownW{T}(unklabel, scale, 1:length(data), data, filtered, covar)
    end
    fitcontiguousx(unk, ffs, chs) = #
        _buildscale(unk, ffs) * pinv(_buildmodel(ffs, chs), rtol = convert(T, 1.0e-6)) * extract(unk,chs)
    function _filterfitx(unk, ffs, fitrois)
        Iterators.flatten(map(fitrois) do fr
            fitcontiguousx(
                unk,
                filter(ff -> length(intersect(fr, ff.ffroi)) > 0, ffs),
                fr,
            )
        end)
    end
    @assert matches(hs[CartesianIndices(hs)[1]], ffp.detector) "The detector for the hyper-spectrum must match the detector for the filtered references."
    krs = zeros(Float32, length(ffp.references), size(hs)...)
    len = 1:min(depth(hs),length(ffp.filter))
    data = zeros(T, length(ffp.filter))
    fitrois = ascontiguous(map(fd -> fd.ffroi, ffp.references))
    cdata = counts(hs)
    ThreadsX.foreach(CartesianIndices(hs)) do ci
        data[len] .= T.(cdata[len,ci])
        unk = _tophatfilterhs(unklabel, data, ffp.filter, convert(T,1/dose(hs,ci)))
        krs[:, ci] .= zero.(_filterfitx(unk, ffp.references, fitrois))
    end
    return ThreadsX.map(filter(ii -> ffp.references[ii].label isa CharXRayLabel, eachindex(ffp.references))) do i
        lbl = ffp.references[i].label
        rprops = properties(spectrum(lbl))
        KRatios(xrays(lbl), properties(hs), rprops, rprops[:Composition], krs[i, :, :])
    end
end

function fit_spectrum_full(
    hs::HyperSpectrum,
    ffp::FilterFitPacket{S, T},
    sigma::T
)::Array{KRatios} where { S <: Detector, T <: AbstractFloat }
    unklabel = UnknownLabel(hs)
    function _tophatfilterhs(unklabel, data, thf, scale)
        @assert length(data) <= length(thf) "The reference spectra must have more channels than the hyperspectrum data."
        filtered = T[filtereddatum(thf, data, i) for i in eachindex(data)]
        dp = T[max(x, one(T)) for x in data] # To ensure covariance isn't zero or infinite precision
        covar = T[filteredcovar(thf, dp, i, i) for i in eachindex(data)]
        FilteredUnknownW{T}(unklabel, scale, 1:length(data), data, filtered, covar)
    end
    @assert matches(hs[CartesianIndices(hs)[1]], ffp.detector) "The detector for the hyper-spectrum must match the detector for the filtered references."
    krs = zeros(Float32, length(ffp.references), size(hs)...)
    len = 1:depth(hs)
    data = zeros(T, length(ffp.filter))
    cdata = counts(hs)
    ThreadsX.foreach(CartesianIndices(hs)) do ci
        data[len] .= view(cdata, len, ci)
        unk = _tophatfilterhs(unklabel, data, ffp.filter, convert(T,1/dose(hs,ci)))
        uvs = _filterfit(unk, ffp.references, true)
        krs[:, ci] = map(ffp.references) do ref
            id = ref.label
            isnan(uvs, id) || (value(uvs, id) < sigma*σ(uvs, id)) ? Base.zero(eltype(krs)) : convert(eltype(krs), value(uvs, id))
        end
    end
    return ThreadsX.map(filter(ii -> ffp.references[ii].label isa CharXRayLabel, eachindex(ffp.references))) do i
        lbl = ffp.references[i].label
        rprops = properties(spectrum(lbl))
        KRatios( xrays(lbl), properties(hs), rprops, rprops[:Composition], krs[i, :, :])
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
###############################################################################
#
# Generate density plot dataframe
#
###############################################################################
"""
    calculate_density(xs::Array, ys::Array)
    calculate_density(xs::Array, ys::Array, nthread::Int)

Plot density plot on axis, given
- `xs` Array of X
- `ys` Array of Y
- `nthread` Number of threads to run in parallel
"""
function calculate_density(xs::Array, ys::Array)
    mask = (.!isnan.(xs)) .* (.!isnan.(ys));
    xx   = xs[mask];
    yy   = ys[mask];
    # calculate the kernel density and sort the density of each data point
    ik = kde((xx,yy));
    cc = similar(xx);
    if length(xx) > 1000
        @showprogress for i in eachindex(xx)
            cc[i] = pdf(ik, xx[i], yy[i]);
        end
    else
        cc .= pdf.([ik], xx, yy);
    end
    df = DataFrame(X=xx, Y=yy, C=cc);
    sort!(df, [:C]);

    return df
end

function calculate_density(xs::Array, ys::Array, nthread::Int)
    # add workers
    if length(workers()) < nthread
        addprocs(nthread - length(workers()), exeflags="--project");
    end
    @everywhere Base.MainInclude.eval(using PlotPlants);

    # filter out NaN
    mask = (.!isnan.(xs)) .* (.!isnan.(ys));
    xx   = xs[mask];
    yy   = ys[mask];

    # calculate the kernel density and sort the density of each data point
    ik     = kde((xx,yy));
    params = [];
    for i in eachindex(xx)
        push!(params, [xx[i], yy[i]]);
    end

    # function to run in threads
    @inline pdf_thread(param) = return pdf(ik, param[1], param[2]);

    if length(xx) > 1000
        cc = @showprogress pmap(pdf_thread, params);
    else
        cc = pmap(pdf_thread, params);
    end

    # create and sort dataframe
    df = DataFrame(X=xx, Y=yy, C=cc);
    sort!(df, [:C]);

    return df
end

ext: jl


text:
using Random

# CVIT - Create itr and itst indeces for k-fold-cv
#
#    Description
#     [ITR,ITST]=CVITR(N,K) returns 1xK cell arrays ITR and ITST holding
#      cross-validation indeces for train and test sets respectively.
#      K-fold division is balanced with all sets having floor(N/K) or
#      ceil(N/K) elements.
#
#     [ITR,ITST]=CVITR(N,K,RS) with integer RS=true also makes random
#      permutation, using substream RS. This way different permutations
#      can be produced with different RS values, but same permutation is
#      obtained when called again with same RS. Function restores the
#      previous random stream before exiting.
#

# Copyright (c) 2010 Aki Vehtari

# This software is distributed under the GNU General Public
# License (version 2 or later); please refer to the file
# License.txt, included with the software, for details.

function cvit(n, k=10, rsubstream=false)

    a = k-rem(n,k)
    b = floor(Int, n/k);

    itst = Any[]
    itr = Any[]

    for cvi in 1:a
        push!(itst, collect(1:b) .+ (cvi-1) * b)
        push!(itr, setdiff(1:n,itst[cvi]))
    end
    for cvi in (a+1):k
        push!(itst, (a * b) + collect(1:(b + 1)) + (cvi - a - 1) * (b + 1))
        push!(itr, setdiff(1:n,itst[cvi]))
    end

    if rsubstream
        rng = MersenneTwister()
        rii = randperm(rng, n)
        for cvi in 1:k
            itst[cvi] = rii[itst[cvi]]
            itr[cvi] = rii[itr[cvi]]
        end
    end
    itr, itst
end


text:
function local_update!(model::Clock, T::Real, Js::AbstractArray)
    rng = model.rng
    nsites = numsites(model)
    nbonds = numbonds(model)
    mbeta = -1.0/T
    iQ2 = 2.0/model.Q

    @inbounds for site in 1:nsites
        center = model.spins[site]
        new_center = mod1(center+rand(rng, 1:(model.Q-1)), model.Q)
        de = 0.0
        for (n,b) in neighbors(model, site)
            de += model.cosines[mod1(model.spins[n]-center, model.Q)] * Js[bondtype(model,b)]
            de -= model.cosines[mod1(model.spins[n]-new_center, model.Q)] * Js[bondtype(model,b)]
        end
        if rand(rng) < exp(mbeta*de)
            model.spins[site] = new_center
        end
    end

    return nothing
end

function SW_update!(model::Clock, T::Real, Js::AbstractArray)
    rng = model.rng
    nsites = numsites(model)
    nbonds = numbonds(model)
    m2bJ = (-2/T).*Js
    m = rand(rng, 1:model.Q)-1
    rspins = zeros(Int, nsites)
    @inbounds for s in 1:nsites
        rspins[s] = mod1(model.spins[s]-m, model.Q)
    end
    uf = UnionFind(nsites)
    @inbounds for bond in bonds(model)
        s1,s2 = source(bond), target(bond)
        bt = bondtype(bond)
        if rand(rng) < -expm1(m2bJ[bt]*model.sines_sw[rspins[s1]]*model.sines_sw[rspins[s2]])
            unify!(uf, s1,s2)
        end
    end
    nc = clusterize!(uf)
    toflips = rand(rng, Bool, nc)
    clustersize = zeros(Int, nc)
    @inbounds for site in 1:nsites
        id = clusterid(uf, site)
        s = ifelse(toflips[id], model.Q-rspins[site]+1, rspins[site])
        clustersize[id] += 1
        model.spins[site] = mod1(s+m, model.Q)
    end

    return nothing
end

function Wolff_update!(model::Clock, T::Real, Js::AbstractArray)
    rng = model.rng
    nsites = numsites(model)
    b2J = (2.0/T).*Js

    clustersize = 0
    st = Stack(Deque{Int}())

    m = rand(rng, 1:model.Q)-1
    center = rand(rng, 1:nsites)
    push!(st, center)
    r = mod1(model.spins[center]-m, model.Q)
    model.spins[center] = mod1(m-r+1, model.Q)

    @inbounds while !isempty(st)
        clustersize += 1
        center = pop!(st)
        cr = mod1(model.spins[center]-m, model.Q)
        for (n,b) in neighbors(model, center)
            nr = mod1(model.spins[n]-m, model.Q)
            bt = bondtype(model,b)
            ## Note: center already flipped
            if rand(rng) < -expm1(b2J[bt]*model.sines_sw[cr]*model.sines_sw[nr])
                push!(st, n)
                model.spins[n] = mod1(m-nr+1,model.Q)
            end
        end
    end

    return nothing
end

@gen_convert_parameter(Clock, ("T", 1, 1.0), ("J", numbondtypes, 1.0))


text:
using SpecialSets

export rules

_strategy(::Val{S}) where {S} = string(S)
macro term(v::Val, ex)
    strategy = _strategy(v)
    :(throw(ArgumentError("Undefined @term strategy: " * $strategy)))
end
macro term(strategy::Symbol, expr)
    esc(:(@term $(Val(strategy)) $expr))
end

macro term(::Val{:AXIOMS}, ex)
    args = map(ex.args) do axiom
        @assert axiom.head == :tuple
        a, b = axiom.args
        a, b = Meta.quot(a), Meta.quot(b)
        :(convert.(Term, ($a, $b)))
    end
    Expr(ex.head, args...)
end
macro term(::Val{:RULES}, ex)
    args = map(ex.args) do rule
        if rule.head == :call
            @assert length(rule.args) == 3 && rule.args[1] == :(=>)
            pair = rule
            ps = []
        else
            @assert rule.head == :where
            pair = rule.args[1]
            ps = Expr(:vect, rule.args[2:end]...)
        end

        @assert pair.head == :call
        p, a, b = pair.args
        @assert p == :(=>)

        esc(:($PatternRule(@term($a), @term($b), $ps)))
    end
    :(Rules($Rule[$(args...)]))
end

rules(set::Symbol, args...; kwargs...) = rules(Val(set), args...; kwargs...)

rules() = [
    rules(:BASIC)
    rules(:ABSOLUTE_VALUE)
    rules(:BOOLEAN)
    rules(:CALCULUS)
    rules(:LOGARITHM)
    rules(:TRIGONOMETRY)
    rules(:TYPES)
]

function rules(::Val{:BASIC})
    @vars x y z
    @vars f

    [
        @term RULES [
            (f(f(x, y), z) => f(x, y, z)) where {σ -> isvalid(Associative(σ[f]))}
            (f(x, f(y, z)) => f(x, y, z)) where {σ -> isvalid(Associative(σ[f]))}

            x - y  => x + -y
            (x / y  => x * inv(y)) where {_image(y, Nonzero)}

            x + -x     => zero(x)
            (x * inv(x) => one(x)) where {_image(x, Nonzero)}

            -(x + y)   => -y + -x
            (inv(x * y) => inv(y) * inv(x)) where {_image(x, Nonzero)}

            -(-x)      => x
            (inv(inv(x)) => x) where {_image(x, Nonzero)}

            -1 * x     => -x
            -x * y     => -(x * y)
            (inv(-x)   => -inv(x)) where {_image(x, Nonzero)}

            x ^ 0      => one(x)
            x ^ 0.0    => one(x)
            x ^ 1      => x
            x ^ 1.0    => x
        ]
        Rules(
            OrderRule(x -> sprint(show, x)),
            EvalRule(+),
            EvalRule(-),
            EvalRule(*),
        )
    ]
end

function rules(::Val{:ABSOLUTE_VALUE})
    @vars x y

    [
        @term RULES [
            (abs(x)     =>  x)  where {_image(x, Nonnegative)}
            (abs(x)     => -x)  where {_image(x, Negative)   }
            abs(-x)     => abs(x)
            abs(x * y)  => abs(x) * abs(y)
            abs(inv(x)) => inv(abs(x))
        ]
        Rules(
            EvalRule(abs),
        )
    ]
end

function rules(::Val{:BOOLEAN}; and=&, or=|, neg=!)
    @vars x y z
    _bool(xs...) = σ -> all(x -> isvalid(Image(σ[x], TypeSet(Bool))), xs)

    [
        @term RULES [
            (and(x, and(y, z)) => and(x, y, z))  where {_bool(x, y, z)}
            (and(and(x, y), z) => and(x, y, z))  where {_bool(x, y, z)}

            (or(x, or(y, z)) => or(x, y, z))  where {_bool(x, y, z)}
            (or(or(x, y), z) => or(x, y, z))  where {_bool(x, y, z)}

            (or(x, false) => x)  where {_bool(x)}
            (and(x, true) => x)  where {_bool(x)}

            (or(x, true)   => true )  where {_bool(x)}
            (and(x, false) => false)  where {_bool(x)}

            (or(x, x)  => x)  where {_bool(x)}
            (and(x, x) => x)  where {_bool(x)}

            (or(x, and(x, y)) => x)  where {_bool(x, y)}
            (and(x, or(x, y)) => x)  where {_bool(x, y)}

            (or(x, neg(x))  => true )  where {_bool(x)}
            (and(x, neg(x)) => false)  where {_bool(x)}

            (neg(neg(x)) => x)  where {_bool(x)}
        ];
        Rules(
            EvalRule(and, &),
            EvalRule(or,  |),
            EvalRule(neg, !),
        );
    ]
end

function _diff(M, fn, arity)
    M === :Base || return
    args = [Variable(Symbol(:_, i)) for i ∈ 1:arity]
    x = Variable(:x)
    f = Expr(:call, getproperty(Base, fn), args...)

    partials = DiffRules.diffrule(M, fn, args...)

    arity == 1 && (partials = (partials,))
    partials = map(x -> get(eval(:(@term $x))), partials)

    ps = (:($*($p, $diff($a, $x))) for (p, a) ∈ zip(partials, args))
    rhs = arity == 1 ? first(ps) : Expr(:call, +, ps...)

    @term(diff(f, x)) => Term(rhs)
end
function rules(::Val{:CALCULUS})
    rules = []
    for (M, fn, arity) ∈ DiffRules.diffrules()
        try
            rule = _diff(M, fn, arity)
            rule === nothing && continue
            push!(rules, rule)
        catch
        end
    end

    @vars x

    Rules(
        rules...,
        @term(diff(x, x)) => @term(one(x)),
        DiffRule(),
    )
end

function rules(::Val{:LOGARITHM})
    n1  = NotEqual(1)
    n01 = NotEqual(0, 1)

    @vars x y r a b c

    @term RULES [
        (log(b, b) => 1)  where {_image(b, n01)}
        (log(b, 1) => 0)  where {_image(b,  n1)}

        (log(b, b ^ x) => x)  where {_image(b, n01)}
        b ^ log(b, x) => x

        log(b, x ^ r) => r * log(b, x)

        log(b, x * y) => log(b, x) + log(b, y)
        log(b, inv(x)) => -log(b, x)

        log(a, b) * log(b, c) => log(a, c)
    ]
end

function rules(::Val{:TRIGONOMETRY})
    @vars α β θ

    @term RULES [
        # Common angles
        sin(0) => 0
        cos(0) => 1
        tan(0) => 0

        sin(π * inv(6)) => 1 / 2
        cos(π * inv(6)) => √3 / 2
        tan(π * inv(6)) => √3 / 3

        sin(π * inv(4)) => √2 / 2
        cos(π * inv(4)) => √2 / 2
        tan(π * inv(4)) => 1

        sin(π * inv(3)) => √3 / 2
        cos(π * inv(3)) => 1 / 2
        tan(π * inv(3)) => √3

        sin(π * inv(2)) => 1
        cos(π * inv(2)) => 0
        # tan(π * inv(2)) => # TODO: infinite/undefined

        # Definitions of relations
        sin(θ) * sec(θ) => tan(θ)
        cos(θ) * csc(θ) => cot(θ)
        inv(cos(θ)) => sec(θ)
        inv(sec(θ)) => cos(θ)
        inv(sin(θ)) => csc(θ)
        inv(csc(θ)) => sin(θ)
        inv(tan(θ)) => cot(θ)
        inv(cot(θ)) => tan(θ)

        # Pythagorean identities
        sin(θ)^2 + cos(θ)^2 => one(θ)
        one(θ) + tan(θ)^2 => sec(θ)^2  # NOTE: will not match any one constants
        one(θ) + cot(θ)^2 => csc(θ)^2

        # Negative angles
        sin(-θ) => -sin(θ)
        cos(-θ) => cos(θ)
        tan(-θ) => -tan(θ)
        (csc(-θ) => -csc(θ))  where {_image(θ, Nonzero)}
        sec(-θ) => sec(θ)
        (cot(-θ) => -cot(θ))  where {_image(θ, Nonzero)}

        # Periodic formulae
        #=
        FIXME where clause requires predicates
        sin(θ + 2πn) where n isa Int => sin(θ)
        cos(θ + 2πn) where n isa Int => cos(θ)
        tan(θ + πn) where n isa Int => tan(θ)
        csc(θ + 2πn) where n isa Int => csc(θ)
        sec(θ + 2πn) where n isa Int => sec(θ)
        cot(θ + πn) where n isa Int => cot(θ)
        =#

        # Double-angle formulae
        2sin(θ)cos(θ) => sin(2θ)
        cos(θ)^2 + -sin(θ)^2 => cos(2θ)
        2cos(θ)^2 + -1 => cos(2θ)

        # Sum and difference formulae
        sin(α)cos(β) + cos(α)sin(β) => sin(α + β)
        sin(α)cos(β) + -(cos(α)sin(β)) => sin(α - β)
        cos(α)cos(β) + -(sin(α)sin(β)) => cos(α + β)
        cos(α)cos(β) + sin(α)sin(β) => cos(α - β)
        (tan(α) + tan(β)) * inv(1 + -tan(α)tan(β)) => tan(α + β)
        (tan(α) + -tan(β)) * inv(1 + tan(α)tan(β)) => tan(α - β)

        # Product to sum formulae
        cos(α + -β) + -cos(α + β) => 2sin(α)sin(β)
        cos(α + -β) + cos(α + β) => 2cos(α)cos(β)
        sin(α + β) + sin(α + -β) => 2sin(α)cos(β)
        sin(α + β) + -sin(α + -β) => 2cos(α)sin(β)

        # Sum to product formulae
        2sin((α + β) * inv(2))cos((α + -β) * inv(2)) => sin(α) + sin(β)
        2cos((α + β) * inv(2))sin((α + -β) * inv(2)) => sin(α) - sin(β)
        2cos((α + β) * inv(2))cos((α + -β) * inv(2)) => cos(α) + cos(β)
        -2sin((α + β) * inv(2))sin((α + -β) * inv(2)) => cos(α) - cos(β)

        # Cofunction formulae
        sin(π * inv(2) + -θ) => cos(θ)
        cos(π * inv(2) + -θ) => sin(θ)
        csc(π * inv(2) + -θ) => sec(θ)
        (sec(π * inv(2) + -θ) => csc(θ))  where {_image(θ, Nonzero)}
        (tan(π * inv(2) + -θ) => cot(θ))  where {_image(θ, Nonzero)}
        cot(π * inv(2) + -θ) => tan(θ)
    ]
end

function rules(::Val{:TYPES})
    @vars x

    rs = []
    types = [Number, Int, Float64]
    for T ∈ types
        T_ = TypeSet(T)

        push!(rs, @term RULES [
            (zero(x)        => zero(T)   ) where {_image(x, T_)}
            (x + zero(x)    => x         ) where {_image(x, T_)}
            (x + $(zero(T)) => x         ) where {_image(x, T_)}
            (-($(zero(T)))  => $(zero(T)))

            (one(x)         => one(T)    ) where {_image(x, T_)}
            (x * one(x)     => x         ) where {_image(x, T_)}
            (x * $(one(T))  => x         ) where {_image(x, T_)}
            (one(x) * x     => x         ) where {_image(x, T_)}
            ($(one(T)) * x  => x         ) where {_image(x, T_)}
            (inv($(one(T))) => $(one(T)) )

            (x * zero(x)    => zero(x)   ) where {_image(x, T_)}
            (x * $(zero(T)) => $(zero(T))) where {_image(x, T_)}
            (zero(x) * x    => zero(x)   ) where {_image(x, T_)}
            ($(zero(T)) * x => $(zero(T))) where {_image(x, T_)}
        ])
    end

    [
        rs...
        Rules(
            EvalRule(zero),
            EvalRule(one),
        )
    ]
end


text:
# build domain and true image
domain = [0.0 1.0 0.0 1.0]
n      = [32,32]
Minv   = getRegularMesh(domain,n)
xc     = getCellCenteredGrid(Minv)
xtrue = sin.(2*pi*xc[:,1]).*sin.(pi*xc[:,2])

# get noisy data
A     = sparse(1.0I,Minv.nc,Minv.nc)
ids   = sort(rand(1:Minv.nc,round(Int64,Minv.nc*.8)))
A     = A[ids,:]
btrue = A*xtrue
bdata = btrue + .1*randn(size(btrue))/norm(btrue,Inf)
Wd    = ones(length(btrue))

# generate misfit param
i1     = (1:round(Int64,size(A,1)/4));
i2     = (round(Int64,size(A,1)/4)+1:round(Int64,2*size(A,1)/4));
i3 	   = (round(Int64,2*size(A,1)/4)+1 : round(Int64,3*size(A,1)/4));
i4 	   = (round(Int64,3*size(A,1)/4)+1 : size(A,1));
pFor1  = LSparam(A[i1,:],[])
pFor2  = LSparam(A[i2,:],[])

sigmaBack = zeros(length(xtrue))
Iact         = sparse(1.0I,Minv.nc,Minv.nc);
gl1          = getGlobalToLocal(1.0,sigmaBack)
gl2          = getGlobalToLocal(Iact,sigmaBack)

bd1          = bdata[i1]
bd2          = bdata[i2]
bd3          = bdata[i3]
bd4          = bdata[i4]
Wd1          = Wd[i1]
Wd2          = Wd[i2]
Wd3          = Wd[i3]
Wd4          = Wd[i4]

pMisRefs    = Array{RemoteChannel}(undef, 6)
workerList  = workers()
nw          = nworkers()
pMisRefs[1] = initRemoteChannel(getMisfitParam,workerList[1%nw + 1],
                                pFor1,Wd1,bd1,SSDFun,fMod,gl1)
pMisRefs[2] = initRemoteChannel(getMisfitParam,workerList[2%nw + 1],
                                pFor2,Wd2,bd2,SSDFun,fMod,gl2)

pForRefs = Array{RemoteChannel}(undef, 4)
pForRefs[1] = initRemoteChannel(LSparam,workerList[3%nw+1],A[i3,:],[])
pForRefs[2] = initRemoteChannel(LSparam,workerList[4%nw+1],A[i4,:],[])
pForRefs[3] = initRemoteChannel(LSparam,workerList[5%nw+1],A[i3,:],[])
pForRefs[4] = initRemoteChannel(LSparam,workerList[6%nw+1],A[i4,:],[])

speye_t = (n)->SparseMatrixCSC(I, n, n);

Mesh2Mesh = Array{Future}(undef, 2)
for i=1:length(Mesh2Mesh)
	k = pForRefs[i].where
	if i==1
		Mesh2Mesh[i] = remotecall(speye_t,k,Minv.nc);
		wait(Mesh2Mesh[i]);
	elseif i==2
		Mesh2Mesh[i] = remotecall(identity,k,1.0);
		wait(Mesh2Mesh[i]);
	end
end

Wd = Array{Array{Float64}}(undef, 2);
dobs = Array{Array{Float64}}(undef, 2);
dobs[1] = bd3;
dobs[2] = bd4;
Wd[1] = Wd3;
Wd[2] = Wd4;
pMisRefs[3:4] = getMisfitParam(pForRefs[1:2],Wd,dobs,SSDFun,Iact,sigmaBack,Mesh2Mesh);
pMisRefs[5:6] = getMisfitParam(pForRefs[3:4],Wd,dobs,SSDFun,Iact,sigmaBack); # single mesh

M2M3		 = sparse(1.0I,Minv.nc,Minv.nc);
M2M4		 = sparse(1.0I,Minv.nc,Minv.nc);
pMis    = Array{MisfitParam}(undef, 6)
pMis[1] = getMisfitParam(pFor1,Wd1,bd1,SSDFun,fMod,gl1);
pMis[2] = getMisfitParam(pFor2,Wd2,bd2,SSDFun,fMod,gl2);
pMis[3] = fetch(pMisRefs[3]);
pMis[4] = fetch(pMisRefs[4]);
pMis[5] = fetch(pMisRefs[5]);
pMis[6] = fetch(pMisRefs[6]);

# setup pInv
alpha        = 1.
x0           = (sum(xtrue)./length(xtrue))*ones(Minv.nc)
boundsLow    = minimum(xtrue)*ones(Minv.nc)
boundsHigh   = maximum(xtrue)*ones(Minv.nc)
sigmaBack    = zeros(Minv.nc)

#  solve with automatic distribution
pInv         = getInverseParam(Minv,fMod,diffusionReg,alpha,x0,boundsLow,boundsHigh)
pInv.maxIter = 5
x1, = projGN(x0,pInv,pMis)

pInv.maxIter = 5
x2, = projGN(x0,pInv,pMisRefs)
@test norm(x1-x2)/norm(x1) .< 1e-12

clear!(pMisRefs);
for k=1:length(pMis)
	clear!(pMis[k],clearPFor=true, clearData=true,clearMesh2Mesh=true);
end

#Finalize pFor references
for pFo in pForRefs
	finalize(pFo)
end


text:
# BAM Writer
# ==========

"""
    BAMWriter(output::BGZFStream, header::SAMHeader)

Create a data writer of the BAM file format.

# Arguments
* `output`: data sink
* `header`: SAM header object
"""
type BAMWriter <: Bio.IO.AbstractWriter
    stream::BGZFStream
end

function BAMWriter(stream::BGZFStream, header::SAMHeader)
    refseqnames = String[]
    refseqlens = Int[]
    for rec in header["SQ"]
        push!(refseqnames, rec["SN"])
        len = rec["LN"]
        if isa(len, AbstractString)
            push!(refseqlens, parse(Int, len))
        else
            push!(refseqlens, len)
        end
    end
    write_header(stream, header, refseqnames, refseqlens)
    return BAMWriter(stream)
end

function Bio.IO.stream(writer::BAMWriter)
    return writer.stream
end

function write_header(stream, header, refseqnames, refseqlens)
    @assert length(refseqnames) == length(refseqlens)
    n = 0

    # magic bytes
    n += write(stream, "BAM\1")

    # SAM header
    buf = IOBuffer()
    l = write(SAMWriter(buf), header)
    n += write(stream, Int32(l))
    n += write(stream, takebuf_array(buf))

    # reference sequences
    n += write(stream, Int32(length(refseqnames)))
    for (seqname, seqlen) in zip(refseqnames, refseqlens)
        namelen = length(seqname)
        n += write(stream, Int32(namelen + 1))
        n += write(stream, seqname, '\0')
        n += write(stream, Int32(seqlen))
    end

    return n
end

function Base.write(writer::BAMWriter, record::BAMRecord)
    n = 0
    n += unsafe_write(writer.stream, pointer_from_objref(record), BAM_FIXED_FIELDS_BYTES)
    n += unsafe_write(writer.stream, pointer(record.data), data_size(record))
    return n
end


text:
#
# runtests.jl -
#
# Run all tests.
#
module TestingLazyAlgebra

using Test

@testset "Utilities                        " begin
    include("utils-tests.jl")
end
@testset "Rules                            " begin
    include("rules-tests.jl")
end
@testset "Generalized matrix multiplication" begin
    include("genmult-tests.jl")
end
@testset "Vectorized operations            " begin
    include("vect-tests.jl")
end
@testset "Mappings                         " begin
    include("map-tests.jl")
end
@testset "Finite differences               " begin
    include("diff-tests.jl")
end
@testset "Sparse operators                 " begin
    include("sparse-tests.jl")
end
@testset "Cropping and padding             " begin
    include("crop-tests.jl")
end
@testset "FFT methods                      " begin
    include("fft-tests.jl")
end
#@testset "Gram operators                   " begin
    include("gram-tests.jl")
#end
@testset "Conjugate gradient               " begin
    include("cg-tests.jl")
end

end # module


text:
# This file is maintained in a way to support CUDA-only test via
# `julia --project=test/cuda -e 'include("runtests.jl")'`
using ImageFiltering
using CUDA
using TestImages
using ImageBase
using ImageQualityIndexes
using Test
using Random

CUDA.allowscalar(false)

@testset "ImageFiltering" begin
    if CUDA.functional()
        include("models.jl")
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Detect division by zero

function isdefinite{T<:Number}(n::T)
    !isequal(n, NaN) && abs(n) != Inf
end

for n in {1, 1//1, 1.0, 1im, 0}
    d = n/0
    print("Divding ", n, " by 0 ")
    if isdefinite(d)
        println("results in ", d, ".")
    else
        println("yields an indefinite value (", d, ").")
    end
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# This file is a part of BinaryCommitteeMachineFBP.jl. License is MIT: http://github.com/carlobaldassi/BinaryCommitteeMachineFBP.jl/LICENCE.md

module BinaryCommitteeMachineFBP

export focusingBP, MagT64, MagP64,
       read_messages, write_messages,
       FocusingProtocol, StandardReinforcement, Scoping, PseudoReinforcement, FreeScoping

using StatsFuns
using GZip
using ExtractMacro
using SpecialFunctions
using Printf, LinearAlgebra, Random
import Base: copyto!

include("Magnetizations.jl")
using .Magnetizations

include("Util.jl")
using .Util

struct Messages{F<:Mag64}
    M::Int
    N::Int
    K::Int

    # 2 layers:
    #   0  --> external fields
    #   1  --> theta (perceptron) nodes
    #   2  --> theta (perceptron) nodes
    #
    # notation:
    #   m* --> total fields (includes external fields)
    #   u* --> messages directed down (node->variable)
    #   U* --> messages directed up (node->variable)
    #
    # variable names (quantity):
    #   w  --> weights (NxK)
    #   τ1 --> outputs from the first layer of perceptrons (KxM)
    #          also inputs to the second layer
    #   τ2 --> outputs from the second layer of perceptrons (M)
    #
    #                 DEPTH
    ux::MagVec2{F}    # 0
    mw::MagVec2{F}    # 0+1
    mτ1::MagVec2{F}   # 1+2
    uw::MagVec3{F}    # 1
    Uτ1::MagVec2{F}   # 1
    mτ2::MagVec{F}    # 2+
    uτ1::MagVec2{F}   # 2

    function Messages{F}(M::Int, N::Int, K::Int, ux::MagVec2{F}, mw::MagVec2{F}, mτ1::MagVec2{F}, uw::MagVec3{F},
                         Uτ1::MagVec2{F}, mτ2::MagVec{F}, uτ1::MagVec2{F}; check::Bool = true) where {F<:Mag64}
        if check
            checkdims(ux, K, N)
            checkdims(mw, K, N)
            checkdims(mτ1, M, K)
            checkdims(uw, M, K, N)
            checkdims(Uτ1, M, K)
            checkdims(mτ2, M)
            checkdims(uτ1, M, K)
        end
        new(M, N, K, ux, mw, mτ1, uw, Uτ1, mτ2, uτ1)
    end

end

function Messages(::Type{F}, M::Integer, N::Integer, K::Integer, x::Real) where {F<:Mag64}
    ux = [mflatp(F, N) for k = 1:K]
    mw = [mflatp(F, N) for k = 1:K]
    mτ1 = [mflatp(F, K) for a = 1:M]
    uw = [MagVec{F}[map(F, x*(2*rand(N) .- 1)) for k = 1:K] for a = 1:M]
    Uτ1 = [mflatp(F, K) for a = 1:M]
    mτ2 = mflatp(F, M)
    uτ1 = [mflatp(F, K) for a = 1:M]

    for k = 1:K, i = 1:N, a = 1:M
        mw[k][i] = mw[k][i] ⊗ uw[a][k][i]
    end
    for a = 1:M, k = 1:K
        mτ1[a][k] = mτ1[a][k] ⊗ Uτ1[a][k] ⊗ uτ1[a][k]
    end

    return Messages{F}(M, N, K, ux, mw, mτ1, uw, Uτ1, mτ2, uτ1, check=false)
end

Messages(::Type{F}, messages::Messages{F}) where {F<:Mag64} = messages
function Messages(::Type{F}, messages::Messages) where {F<:Mag64}
    @extract messages : M N K ux mw mτ1 uw Uτ1 mτ2 uτ1
    return Messages{F}(M, N, K,
                       chgeltype(ux, F), chgeltype(mw, F), chgeltype(mτ1, F), chgeltype(uw, F),
                       chgeltype(Uτ1, F), chgeltype(mτ2, F), chgeltype(uτ1, F),
                       check=false)
end

function read_messages(io::IO, ::Type{F}) where {F<:Mag64}
    l = split(readline(io))
    length(l) == 2 && l[1] == "fmt:" || error("invalid messages file")
    fmt = Val{Symbol(l[2])}
    l = split(readline(io))
    length(l) == 4 && l[1] == "N,K,M:" || error("invalid messgaes file")
    N, K, M = parse(Int, l[2]), parse(Int, l[3]), parse(Int, l[4])

    ux = [mflatp(F, N) for k = 1:K]
    mw = [mflatp(F, N) for k = 1:K]
    mτ1 = [mflatp(F, K) for a = 1:M]
    uw = [[mflatp(F, N) for k = 1:K] for a = 1:M]
    Uτ1 = [mflatp(F, K) for a = 1:M]
    mτ2 = mflatp(F, M)
    uτ1 = [mflatp(F, K) for a = 1:M]

    expected_lines = K + M + M*K + M + 1 + M + K
    for (i,l) in enumerate(eachline(io))
        i > expected_lines && (strip(l) == "END" || error("invalid messages file"); break)
        @readmagvec(l, fmt, ux, mw, mτ1, uw, Uτ1, mτ2, uτ1)
    end
    eof(io) || error("invalid messages file")
    return Messages{F}(M, N, K, ux, mw, mτ1, uw, Uτ1, mτ2, uτ1, check=false)
end

"""
    read_messages(filename, mag_type)

Reads messages from a file. `mag_type` is the internal storage format used in the resulting `Messages` object,
it can be either `MagT64` (uses tanhs, accurate but slower) or `MagP64` (plain format, faster but inaccurate).

The file format is the one produced by [`write_messages`](@ref).
"""
read_messages(filename::AbstractString, ::Type{F}) where {F<:Mag64} = gzopen(io->read_messages(io, F), filename, "r")

"""
    write_messages(filename, messages)

Writes messages to a file. The messages can be read back with [`read_messages`](@ref). Note that the output
is a plain text file compressed with gzip.
"""
function write_messages(filename::AbstractString, messages::Messages)
    gzopen(filename, "w") do f
        write_messages(f, messages)
    end
end

function write_messages(io::IO, messages::Messages{F}) where {F<:Mag64}
    @extract messages : N K M ux mw mτ1 uw Uτ1 mτ2 uτ1

    println(io, "fmt: ", magformat(F))
    println(io, "N,K,M: $N $K $M")
    @dumpmagvecs(io, ux, mw, mτ1, uw, Uτ1, mτ2, uτ1)
    println(io, "END")
end

Base.eltype(messages::Messages{F}) where {F<:Mag64} = F

function copyto!(dest::Messages{F}, src::Messages{F}) where {F<:Mag64}
    dest.N == src.N || throw(ArgumentError("incompatible arguments: dest.N=$(dest.N) src.N=$(src.N)"))
    dest.K == src.K || throw(ArgumentError("incompatible arguments: dest.K=$(dest.K) src.K=$(src.K)"))
    dest.M == src.M || throw(ArgumentError("incompatible arguments: dest.M=$(dest.M) src.M=$(src.M)"))
    for k = 1:dest.K
        copyto!(dest.ux[k], src.ux[k])
        copyto!(dest.mw[k], src.mw[k])
    end
    for a = 1:dest.M, k = 1:dest.K
        copyto!(dest.uw[a][k], src.uw[a][k])
    end
    for a = 1:dest.M
        copyto!(dest.mτ1[a], src.mτ1[a])
        copyto!(dest.Uτ1[a], src.Uτ1[a])
        copyto!(dest.uτ1[a], src.uτ1[a])
    end
    copyto!(dest.mτ2, src.mτ2)
    return dest
end

function set_outfields!(messages::Messages{F}, output::Vector, β::Float64) where {F<:Mag64}
    @extract messages : N K M mτ2
    @assert length(output) == M
    t = tanh(β / 2)
    for a = 1:M
        mτ2[a] = forcedmag(F, output[a] * t) # forced avoids clamping
    end
end

print_mags(messages::Messages) = print_mags(stdout, messages)
function print_mags(io::IO, messages::Messages)
    @extract messages : N K mw
    for k = 1:K, i = 1:N
        @printf(io, "%i %i %.15f\n", k, i, Float64(mw[k][i]))
    end
end

mutable struct Params{F<:Mag64}
    damping::Float64
    ϵ::Float64
    β::Float64
    max_iters::Int
    accuracy1::Symbol
    accuracy2::Symbol
    r::Float64
    pol::F
    dγ::Float64
    quiet::Bool
end

struct Patterns
    M::Int
    X::Vec2
    output::IVec
    function Patterns(X::AbstractVector, output::AbstractVector)
        M = length(X)
        length(output) == M || throw(ArgumentError("incompatible lengths of inputs and outputs: $M vs $(length(output))"))
        all(ξ->all(ξi->abs(ξi) == 1, ξ), X) || throw(ArgumentError("inputs must be ∈ {-1,1}"))
        all(o->abs(o) == 1, output) || throw(ArgumentError("outputs must be ∈ {-1,1}"))
        new(M, X, output)
    end
end
Patterns(Xo::Tuple{Vec2,Vec}) = Patterns(Xo...)

@doc """
    Patterns(inputs::AbstractVector, outputs::AbstractVector)

Construct a `Patterns` object. `inputs` and `outputs` must have the same length.
`outputs` entries must be ∈ {-1,1}. `intputs` entries must be Vectors in which
each element is ∈ {-1,1}, and all the vectors must have the same length.
""" Patterns(inputs::AbstractVector, outputs::AbstractVector)

"""
    Patterns(NM::Tuple{Integer,Integer}; teacher::Union{Bool,Int,Vector{Vector{Float64}}} = false)

Constructs random (unbiasad, unifom i.i.d.) binary patterns. The `NM` argument is a Tuple with
two items, the size of the inputs (`N`) and the number of patterns (`M`).
The keyword argument `teacher` controls how the outputs are generated:
* If `false`, they are random i.i.d. (the default)
* If `true`, a teacher unit with a single hidden layer (a perceptron) is generated randomly
  and it's used to compute the outputs
* If an `Int` is given, it's the number of hidden units of the teacher, which is generated
  randomly and used to compute the outputs
* If a `Vector{Vector{Float64}}` is given, it represents the weights of the teacher, with one
  `Vector{Float64}` for each hidden unit.
"""
function Patterns(NM::Tuple{Integer,Integer}; teacher::Union{Bool,Int,Vec2} = false)
    N, M = NM
    X = [rand(-1.0:2.0:1.0, N) for a = 1:M]
    if teacher == false
        output = ones(M)
    else
        if !isa(teacher, Vec2)
            K::Int = isa(teacher, Bool) ? 1 : teacher
            tw = [rand(-1.0:2.0:1.0, N) for k = 1:K]
        else
            K = length(teacher)
            K > 0 || throw(ArgumentError("empty teacher"))
            tw = teacher
            all(w -> length(w) == N, tw) || throw(ArgumentError("invalid teacher length, expected $N, given: $(sort!(unique(collect(w->length(w) for w in tw))))"))
            all(w -> all(x -> abs(x) == 1, w), tw) || throw(ArgumentError("invalid teacher, entries must be all ±1"))
        end
        r0 = Array{Float64}(undef, K)
        output = Float64[transf1!(r0, tw, ξ) for ξ in X]
        @assert all(o->abs(o) == 1, output)
    end
    Patterns(X, output)
end
Patterns(patterns::Patterns) = deepcopy(patterns)

"""
    Patterns(patternsfile::AbstractString)

Read patterns from a file. It only reads the inputs, one pattern per line, entries separated by
whitespace. All lines must have the same length. All outputs are assumed to be `1`. The file may
optionally be gzipped.
"""
function Patterns(patternsfile::AbstractString)
    X = Vec[]
    N = 0
    M = 0
    gzopen(patternsfile) do f
        M = 0
        for l in eachline(f)
            push!(X, map(x->parse(Float64, x), split(l)))
            M += 1
        end
    end
    o = Int[1.0 for a = 1:M]

    return Patterns(X, o)
end

function computeσ²(w::Vec)
    σ² = 0.0
    @inbounds @simd for wi in w
        # σ² += (1 - wi^2)
        σ² += (1 - wi*wi)
    end
    return σ²
end

# function computeσ²(w::Vec, ξ::Vec)
#     σ² = 0.0
#     @inbounds for (wi,ξi) in zip(w,ξ)
#         σ² += (1 - wi^2) * ξi^2
#     end
#     return σ²
# end

computeσ(σ²::Float64) = √(2σ²)
computeσ(w::Vec) = √(2computeσ²(w))
# computeσ(w::Vec, ξ::Vec) = √(2computeσ²(w, ξ))

# transf0(w::Vec, ξ::Vec) = transf0(w, computeσ(w, ξ), ξ)
transf0(w::Vec, ξ::Vec) = transf0(w, computeσ(w), ξ) # NOTE: assuming ξ = ±1 here!!!
transf0(w::Vec, σ::Float64, ξ::Vec) = erf((ξ ⋅ w) / σ)

function transf0!(r::Vec, w1::Vec2, ξ::Vec)
    for (i,w) in enumerate(w1)
        r[i] = transf0(w, ξ)
    end
    return r
end

gauss(x, dσ²) = ℯ^(-x^2 / dσ²) / √(π * dσ²)
gauss(x) = ℯ^(-x^2 / 2) / √(2π)

function subfield!(h::MagVec, m::MagVec, u::MagVec)
    @inbounds @simd for i = 1:length(m)
        h[i] = m[i] ⊘ u[i]
    end
end

function addfield!(m::MagVec, h::MagVec, u::MagVec)
    @inbounds @simd for i = 1:length(m)
        m[i] = h[i] ⊗ u[i]
    end
end

let hsT = Dict{Int,MagVec{MagT64}}(), hsP = Dict{Int,MagVec{MagP64}}(), vhs = Dict{Int,Vec}(),
    vHs = Dict{Int,Vec}(), leftCs = Dict{Int,Vec2}(), rightCs = Dict{Int,Vec2}()

    geth(::Type{MagT64}, N::Int) = Base.get(()->Array{MagT64}(undef, N), hsT, N)
    geth(::Type{MagP64}, N::Int) = Base.get(()->Array{MagP64}(undef, N), hsP, N, )

    global theta_node_update_approx!
    function theta_node_update_approx!(m::MagVec{F}, M::F, ξ::Vec, u::MagVec{F}, U::F, params::Params) where {F<:Mag64}
        @extract params : λ=damping

        N = length(m)
        h::MagVec{F} = geth(F, N)
        vh = Base.get(()->Array{Float64}(undef, N), vhs, N)

        subfield!(h, m, u)
        H = M ⊘ U

        @inbounds for i = 1:N
            vh[i] = h[i]
        end

        vH = Float64(H)
        # σ² = computeσ²(vh, ξ)
        σ² = computeσ²(vh) # assume ξ = ±1

        μ = vh ⋅ ξ

        dσ² = 2σ²
        newU = merf(F, μ / √dσ²)

        maxdiff = abs(U - newU)
        U = damp(newU, U, λ)
        newM = H ⊗ U
        M = newM

        g = gauss(μ, dσ²)

        p0 = 2vH * g / (1 + vH * U)

        pμ = p0 * (p0 + μ / σ²)

        pσ = p0 * (1 - μ / σ² - μ * p0) / dσ²

        @inbounds for i = 1:N
            ξi = ξ[i]
            hi = vh[i]
            newu = convert(F, clamp(ξi * (p0 + ξi * (hi * pμ + ξi * (1-hi^2) * pσ)), -1+3e-16, 1-3e-16)) # use mag-functions?
            d = conv_diff(newu, u[i])
            maxdiff = ifelse(maxdiff < d, d, maxdiff)
            u[i] = damp(newu, u[i], λ)
        end
        addfield!(m, h, u)

        return maxdiff, U, M
    end

    global theta_node_update_accurate!
    function theta_node_update_accurate!(m::MagVec{F}, M::F, ξ::Vec, u::MagVec{F}, U::F, params::Params) where {F<:Mag64}
        @extract params : λ=damping

        N = length(m)
        h::MagVec{F} = geth(F, N)
        vh = Base.get(()->Array{Float64}(undef, N), vhs, N)

        subfield!(h, m, u)
        H = M ⊘ U

        @inbounds for i = 1:N
            vh[i] = h[i]
        end

        # σ² = computeσ²(vh, ξ)
        σ² = computeσ²(vh) # assume ξ = ±1
        μ = vh ⋅ ξ

        dσ² = 2σ²
        newU = merf(F, μ / √dσ²)

        maxdiff = 0.0
        U = damp(newU, U, λ)
        M = H ⊗ U

        @inbounds for i = 1:N
            ξi = ξ[i]
            hi = vh[i]
            μ̄ = μ - ξi * hi
            σ̄² = σ² - (1-hi^2) # * ξi^2 # assume ξ = ±1
            sdσ̄² = √(2σ̄²)
            m₊ = (μ̄ + ξi) / sdσ̄²
            m₋ = (μ̄ - ξi) / sdσ̄²
            newu = erfmix(H, m₊, m₋)
            d = abs(newu - u[i])
            maxdiff = ifelse(maxdiff < d, d, maxdiff)
            u[i] = damp(newu, u[i], λ)
        end
        addfield!(m, h, u)
        return maxdiff, U, M
    end

    global theta_node_update_exact!
    function theta_node_update_exact!(m::MagVec{F}, M::F, ξ::Vec, u::MagVec{F}, U::F, params::Params) where {F<:Mag64}
        @extract params : λ=damping

        N = length(m)
        h::MagVec{F} = geth(F, N)
        vh = Base.get(()->Array{Float64}(undef, N), vhs, N)
        leftC = Base.get(()->[zeros(i+1) for i = 1:N], leftCs, N)
        rightC = Base.get(()->[zeros((N-i+1)+1) for i = 1:N], rightCs, N)

        subfield!(h, m, u)
        H = M ⊘ U

        @inbounds for i = 1:N
            vh[i] = h[i]
        end

        leftC[1][1] = (1-ξ[1]*vh[1])/2
        leftC[1][2] = (1+ξ[1]*vh[1])/2
        for i = 2:N
            lC0 = leftC[i-1]
            lC = leftC[i]
            hi = ξ[i] * vh[i]
            lC[1] = lC0[1] * (1-hi)/2
            for j = 2:i
                lC[j] = lC0[j-1] * (1+hi)/2 + lC0[j] * (1-hi)/2
            end
            lC[end] = lC0[end] * (1+hi)/2
        end

        rightC[end][1] = (1-ξ[end]*vh[end])/2
        rightC[end][2] = (1+ξ[end]*vh[end])/2
        for i = (N-1):-1:1
            rC0 = rightC[i+1]
            rC = rightC[i]
            hi = ξ[i] * vh[i]
            rC[1] = rC0[1] * (1-hi)/2
            for j = 2:(N-i+1)
                rC[j] = rC0[j-1] * (1+hi)/2 + rC0[j] * (1-hi)/2
            end
            rC[end] = rC0[end] * (1+hi)/2
        end

        @assert maximum(abs.(leftC[end] .- rightC[1])) ≤ 1e-10 (leftC[end], rightC[1])

        @assert isodd(N)
        z = (N+1) ÷ 2
        pm = sum(rightC[1][1:z])
        pp = sum(rightC[1][(z+1):end])

        newU = Mag64(F, pp, pm)

        @assert isfinite(newU)

        maxdiff = 0.0
        U = damp(newU, U, λ)
        newM = H ⊗ U
        M = newM

        @assert isfinite(newM) (H, U)

        u1 = ones(1)

        @inbounds for i = 1:N
            ξi = ξ[i]
            @assert ξi^2 == 1

            lC = i > 1 ? leftC[i-1] : u1
            rC = i < N ? rightC[i+1] : u1

            pm = 0.0
            pz = 0.0
            pp = 0.0
            for j = 1:N
                p = 0.0
                for k = max(1,j+i-N):min(j,i)
                    p += lC[k] * rC[j-k+1]
                end
                if j < z
                    pm += p
                elseif j == z
                    pz = p
                else
                    pp += p
                end
            end

            mp = convert(F, clamp(pp + ξi * pz - pm, -1.0, 1.0))
            mm = convert(F, clamp(pp - ξi * pz - pm, -1.0, 1.0))
            newu = exactmix(H, mp, mm)

            d = abs(newu - u[i])
            maxdiff = ifelse(maxdiff < d, d, maxdiff)
            u[i] = damp(newu, u[i], λ)
            @assert isfinite(u[i]) (u[i],)
        end
        addfield!(m, h, u)
        return maxdiff, U, M
    end

    global free_energy_theta
    function free_energy_theta(m::MagVec{F}, M::F, ξ::Vec, u::MagVec{F}, U::F) where {F<:Mag64}
        N = length(m)
        h::MagVec{F} = geth(F, N)
        vh = Base.get(()->Array{Float64}(undef, N), vhs, N)

        f = 0.0

        subfield!(h, m, u)
        H = M ⊘ U

        @inbounds for i = 1:N
            vh[i] = h[i]
        end

        # σ = computeσ(vh, ξ)
        σ = computeσ(vh) # assume ξ = ±1
        μ = vh ⋅ ξ

        b = merf(F, μ / σ)

        f -= log1pxy(H, b)
        @assert isfinite(f)

        for i = 1:N
            f += log1pxy(h[i], u[i])
        end
        return f
    end

    global free_energy_theta_exact
    function free_energy_theta_exact(m::MagVec{F}, M::F, ξ::Vec, u::MagVec{F}, U::F) where {F<:Mag64}

        N = length(m)
        h::MagVec{F} = geth(F, N)
        vh = Base.get(()->Array{Float64}(undef, N), vhs, N)
        leftC = Base.get(()->[zeros(i+1) for i = 1:N], leftCs, N)
        rightC = Base.get(()->[zeros((N-i+1)+1) for i = 1:N], rightCs, N)

        f = 0.0

        subfield!(h, m, u)
        H = M ⊘ U

        @inbounds for i = 1:N
            vh[i] = h[i]
        end
        vH = Float64(H)

        leftC[1][1] = (1-ξ[1]*vh[1])/2
        leftC[1][2] = (1+ξ[1]*vh[1])/2
        for i = 2:N
            lC0 = leftC[i-1]
            lC = leftC[i]
            hi = ξ[i] * vh[i]
            lC[1] = lC0[1] * (1-hi)/2
            for j = 2:i
                lC[j] = lC0[j-1] * (1+hi)/2 + lC0[j] * (1-hi)/2
            end
            lC[end] = lC0[end] * (1+hi)/2
        end

        rightC[end][1] = (1-ξ[end]*vh[end])/2
        rightC[end][2] = (1+ξ[end]*vh[end])/2
        for i = (N-1):-1:1
            rC0 = rightC[i+1]
            rC = rightC[i]
            hi = ξ[i] * vh[i]
            rC[1] = rC0[1] * (1-hi)/2
            for j = 2:(N-i+1)
                rC[j] = rC0[j-1] * (1+hi)/2 + rC0[j] * (1-hi)/2
            end
            rC[end] = rC0[end] * (1+hi)/2
        end

        @assert maximum(abs.(leftC[end] .- rightC[1])) ≤ 1e-10 (leftC[end], rightC[1])

        @assert isodd(N)
        z = (N+1) ÷ 2
        pm = sum(rightC[1][1:z])
        pp = sum(rightC[1][(z+1):end])

        b = Mag64(F, pp, pm)
        f -= log1pxy(H, b)
        @assert isfinite(f)

        for i = 1:N
            f += log1pxy(h[i], u[i])
        end
        return f
    end

end

function entro_node_update(m::F, u::F, params::Params{F}) where {F<:Mag64}
    @extract params : λ=damping r pol

    h = m ⊘ u
    if r == 0 || iszero(pol)
        newu = zero(F)
    elseif r == Inf
        newu = ifelse(iszero(h), zero(F), copysign(pol, h))
    else
        newu::F = ((h * pol) ↑ r) * pol
    end

    diff = abs(newu - u)
    newu = damp(newu, u, λ)
    newm = h ⊗ newu

    return diff, newu, newm
end

function iterate!(messages::Messages{F}, patterns::Patterns, params::Params) where {F<:Mag64}
    @extract messages : N M K ux mw mτ1 uw Uτ1 mτ2 uτ1
    @extract patterns : X output
    @extract params   : accuracy1 accuracy2
    maxdiff = 0.0
    tnu1! = accuracy1 == :exact ? theta_node_update_exact! :
            accuracy1 == :accurate ? theta_node_update_accurate! :
            accuracy1 == :none ? theta_node_update_approx! :
            error("accuracy must be one of :exact, :accurate, :none (was given $accuracy)")
    tnu2! = accuracy2 == :exact ? theta_node_update_exact! :
            accuracy2 == :accurate ? theta_node_update_accurate! :
            accuracy2 == :none ? theta_node_update_approx! :
            error("accuracy must be one of :exact, :accurate, :none (was given $accuracy)")
    for a = randperm(M + N*K)
        if a ≤ M
            ξ = X[a]
            for k = 1:K
                diff, Uτ1[a][k], mτ1[a][k] = tnu1!(mw[k], mτ1[a][k], ξ, uw[a][k], Uτ1[a][k], params)
                maxdiff = max(maxdiff, diff)
            end
            diff, _, mτ2[a] = tnu2!(mτ1[a], mτ2[a], ones(K), uτ1[a], zero(F), params)
            maxdiff = max(maxdiff, diff)
        else
            (params.r == 0 || iszero(params.pol)) && continue
            j = a - M
            k = (j-1) ÷ N + 1
            i = (j-1) % N + 1

            diff, ux[k][i], mw[k][i] = entro_node_update(mw[k][i], ux[k][i], params)
            maxdiff = max(diff, maxdiff)
        end
    end
    return maxdiff
end

function converge!(messages::Messages, patterns::Patterns, params::Params)
    @extract params : ϵ max_iters λ₀=damping quiet

    λ = λ₀
    ok = false
    strl = 0
    t = @elapsed for it = 1:max_iters
        diff = iterate!(messages, patterns, params)

        if !quiet
            str = "[it=$it Δ=$diff λ=$λ]"
            print("\r", " "^strl, "\r", str)
            strl = length(str)
            #println(str)
            flush(stdout)
            strl = length(str)
        end
        if diff < ϵ
            ok = true
            quiet || println("\nok")
            break
        end
    end
    if !quiet
        ok || println("\nfailed")
        println("elapsed time = $t seconds")
    end
    return ok
end

transf1(w::Vec) = sign(sum(w))

transf1!(r0::Vec, ws::Vec2, ξ::Vec) = transf1(transf0!(r0, ws, ξ))

function test!(r0::Vec, ws::Vec2, ξ::Vec, out::Int)
    o = transf1!(r0, ws, ξ)
    return o != out
end

let r0d = Dict{Int,Vec}()
    global test
    function test(ws::Vec2, ξs::Vec2, output::IVec)
        N = length(ws)
        r0 = get!(r0d, N) do; Vec(undef, N); end
        sum([test!(r0, ws, ξ, out) for (ξ,out) in zip(ξs, output)])
    end
end

function test(messages::Messages, patterns::Patterns)
    @extract messages : N K mw
    @extract patterns : X output
    ws = [Float64[mw[k][i] for i = 1:N] for k = 1:K]
    return test(ws, X, output)
end

function nonbayes_test(messages::Messages, patterns::Patterns)
    @extract messages : N K mw
    @extract patterns : X output
    ws = [Float64[sign0(mw[k][i]) for i = 1:N] for k = 1:K]
    return test(ws, X, output)
end

function free_energy(messages::Messages{F}, patterns::Patterns, params::Params{F}) where {F<:Mag64}
    @extract messages : M N K ux mw mτ1 uw Uτ1 mτ2 uτ1
    @extract patterns : X output
    @extract params   : r pol

    f = 0.0

    for a = 1:M
        ξ = X[a]
        for k = 1:K
            f += free_energy_theta(mw[k], mτ1[a][k], ξ, uw[a][k], Uτ1[a][k])
        end
        f += free_energy_theta_exact(mτ1[a], mτ2[a], ones(K), uτ1[a], zero(F))
    end

    for k = 1:K, i = 1:N
        f -= logZ(ux[k][i], F[uw[a][k][i] for a = 1:M])

        f -= logtwo / 2
        f += log1pxy(pol, -pol) / 2
        hkix = mw[k][i] ⊘ ux[k][i]
        f += log1pxy(hkix, ux[k][i])
        hpol = hkix * pol
        mx = hpol ↑ (r + 1)
        f += mcrossentropy(mx, hpol)
    end

    return f / (N * K)
end

function compute_S(messages::Messages{F}, params::Params{F}) where {F<:Mag64}
    @extract messages : N K ux mw
    @extract params   : r pol
    S = 0.0
    for k = 1:K, i = 1:N
        hkix = mw[k][i] ⊘ ux[k][i]
        hxki = (hkix * pol) ↑ r
        hh = hkix * hxki
        S += Float64(hh ⊗ pol)
    end
    return S / (N * K)
end

function compute_q̃(messages::Messages{F}, params::Params{F}) where {F<:Mag64}
    @extract messages : N K ux mw
    @extract params   : r pol
    q̃ = 0.0
    for k = 1:K, i = 1:N
        hkix = mw[k][i] ⊘ ux[k][i]
        mx = (hkix * pol) ↑ (r + 1)
        q̃ += mx^2
    end
    return q̃ / (N * K)
end

function compute_q(messages::Messages)
    @extract messages : N K mw
    q = 0.0
    for k = 1:K, i = 1:N
        q += Float64(mw[k][i])^2
    end
    return q / (N * K)
end

function mags_symmetry(messages::Messages)
    @extract messages : N K mw
    overlaps = Matrix(1.0I, K, K)
    qs = zeros(K)
    for k1 = 1:K
        z = 0.0
        for i = 1:N
            z += Float64(mw[k1][i])^2
        end
        qs[k1] = √z
    end
    for k1 = 1:K, k2 = k1+1:K
        s = 0.0
        for i = 1:N
            s += Float64(mw[k1][i]) * Float64(mw[k2][i])
        end
        s /= qs[k1] * qs[k2]
        overlaps[k1,k2] = s
        overlaps[k2,k1] = s
    end
    return overlaps, qs / N
end

"""
    FocusingProtocol

Abstract type representing a protocol for the focusing procedure, i.e. a way to produce
successive values for the quantities `γ`, `y` and `β`. Currently, however, only `β=Inf`
is supported. To be provided as an argument to [`focusingBP`](@ref).

Available protocols are: [`StandardReinforcement`](@ref), [`Scoping`](@ref), [`PseudoReinforcement`](@ref) and
[`FreeScoping`](@ref).
"""
abstract type FocusingProtocol end

struct StandardReinforcement <: FocusingProtocol
    r::AbstractRange{Float64}
    StandardReinforcement(r::AbstractRange{T}) where {T<:Real} = new(r)
end

@doc """
    StandardReinforcement(r::AbstractRange) <: FocusingProtocol

Standard reinforcement protocol, returns `γ=Inf` and `y=1/(1-x)`, where `x` is taken from the given range `r`.
""" StandardReinforcement(r::AbstractRange)

"""
    StandardReinforcement(dr::Float64) <: FocusingProtocol

Shorthand for [`StandardReinforcement`](@ref)`(0:dr:(1-dr))`.
"""
StandardReinforcement(dr::Float64) = StandardReinforcement(0.0:dr:(1-dr))

function Base.iterate(s::StandardReinforcement, i = 1)
    n = iterate(s.r, i)
    n ≡ nothing && return nothing
    return (Inf, 1/(1-n[1]), Inf), n[2]
end

"""
    Scoping(γr::AbstractRange, y) <: FocusingProtocol

Focusing protocol with fixed `y` and a varying `γ` taken from the given `γr` range.
"""
struct Scoping <: FocusingProtocol
    γr::AbstractRange{Float64}
    y::Float64
    Scoping(γr::AbstractRange, y) = new(γr, y)
end

function Base.iterate(s::Scoping, i = 1)
    n = iterate(s.γr, i)
    n ≡ nothing && return nothing
    return (n[1], s.y, Inf), n[2]
end

struct PseudoReinforcement <: FocusingProtocol
    r::Vector{Float64}
    x::Float64
    PseudoReinforcement(r::AbstractRange{T}...; x::Real=0.5) where {T<:Real} = new(vcat(map(collect, r)...), x)
end
@doc """
    PseudoReinforcement(r::AbstractRange...; x=0.5) <: FocusingProtocol

A focusing protocol in which both `γ` and `y` are progressively increased, according to
the formulas

```julia
γ = atanh(ρ^x)
y = 1+ρ^(1-2x)/(1-ρ)
```

where `ρ` is taken from the given range(s) `r`. With `x=0`, this is basically the same as
[`StandardReinforcement`](@ref).
""" PseudoReinforcement(r::AbstractRange...)

"""
    PseudoReinforcement(dr::Float64; x=0.5) <: FocusingProtocol

Shorthand for [`PseudoReinforcement`](@ref)`(0:dr:(1-dr); x=x)`.
""" PseudoReinforcement(dr::Float64; x::Real=0.5) = PseudoReinforcement(0.0:dr:(1-dr), x=x)

function Base.iterate(s::PseudoReinforcement, i = 1)
    n = iterate(s.r, i)
    n ≡ nothing && return nothing
    x = s.x
    ρ = n[1]
    # some special cases just to avoid possible 0^0
    if x == 0.5
        return (atanh(√ρ), (2-ρ)/(1-ρ), Inf), n[2]
    elseif x == 0
        return (Inf, 1/(1-ρ), Inf), n[2]
    else
        return (atanh(ρ^x), 1+ρ^(1-2x)/(1-ρ), Inf), n[2]
    end
end

"""
    FreeScoping(list::Vector{NTuple{2,Float64}}) <: FocusingProtocol

A focusing protocol which just returns the values of `(γ,y)` from the given `list`.

Example:

```julia
FreeScoping([(1/(1-x), (2-x)/(1-x)) for x = 0:0.01:0.99])
```
"""
struct FreeScoping <: FocusingProtocol
    list::Vector{NTuple{3,Float64}}
    FreeScoping(list::Vector{NTuple{3,Float64}}) = new(list)
end
FreeScoping(list::Vector{NTuple{2,Float64}}) = FreeScoping(NTuple{3,Float64}[(γ,y,Inf) for (γ,y) in list])

Base.iterate(s::FreeScoping, i = 1) = iterate(s.list, i)

"""
    focusingBP(N, K, patternspec; keywords...)

Run the Focusing Belief Propagation algorithm on a fully-connected committee machine with binary weights.
`N` is the input (first layer) size, `K` the number of the hidden units (second layer size), and
`patternspec` specifies how to build the patterns for the training set. Note that with the defult
settings `K` must be odd (see notes for the `accuracy1` and `accuracy2` arguments below).

Possible values of `patternspec` are:

* a `Float64` number: this is interpreted as the `α` parameter, and `M = α*N*K` random ±1 patterns are generated.
* a `Tuple` with `Vector{Vector{Float64}}` and a `Vector{Float64}`: these are the inputs and associated desired outputs.
* a string: the patterns are read from a file (one input pattern per line, entries separated by whitespace, outputs are
            assumed to be all 1); the file can be gzipped.
* a [`Patterns`](@ref) object (which could be the output of a previous run of the function).

*Note*: all inputs and outputs must be ∈ {-1,1}.

The keyword arguments are:

* `max_iters` (default = `1000`): maximum number of BP iterations per step. If convergence is not achieved in this many iterations,
                                  the algorithm proceeds to the next step.
* `max_steps` (default = `typemax(Int)`): maximum number of focusing steps.
* `seed` (default = `1`): random seed.
* `damping` (default = `0`): BP damping parameter (between `0` and `1`; `0` = no damping).
* `quiet` (default = `false`): whether to print on screen
* `accuracy1` (default = `:accurate`): accuracy of the messages computation at the hidden units level. Possible values are:
                                       `:accurate` (Gaussian approximation, good for large `N`, works in `O(N)` time),
                                       `:exact` (no approximation, uses convolutions, good for small `N`, works in `O(N³)` time, requires
                                       `N` to be odd),
                                       `:none` (a TAP-like approximation, fast but never very good, should probably be removed or done
                                       properly...).
* `accuracy2` (default = `:exact`): accuracy of the messages computation at the output node level.
                                    See `accuracy1` (and think of `K` instead of `N`).
* `randfact` (default = `0.01`): random factor used in the initialization of the messages. Must be between `0` and `1`. Large values are
                                 not a good idea.
* `fprotocol` (default = `StandardReinforcement(1e-2)`): focusing protocol specification. See [`FocusingProtocol`](@ref).
* `ϵ` (default = `1e-3`): convergence criterion: BP is assumed to have converged when the difference between messages in two successive
                          iterations is smaller than this value. Reduce it (e.g. to 1e-6) for more precise results, while increasing
                          `max_iters`.
* `messfmt` (default = `:tanh`): internal storage format for messages: it can be either `:tanh` or `:plain`; `:tanh` is much more
                                 precise but slower.
* `initmess` (default = `nothing`): how to initialize the messages. If `nothing`, they are initialized randomly; If a string is given,
                                    they are read from a file (see also [`read_messages`](@ref) and [`write_messages`](@ref)). If a
                                    `Messages` object is given (e.g. one returned from an earlier run of the function, or from
                                    [`read_messages`](@ref)) it is used (and overwritten).
* `outatzero` (default = `true`): if `true`, the algorithm exits as soon as a solution to the learning problem is found, without waiting
                                  for the focusing protocol to terminate.
* `writeoutfile` (default = `:auto`): whether to write results on an output file. Can be `:never`, `:always` or `:auto`. The latter means
                                      that a file is written when `outatzero` is set to `false` and only when BP converges. It can make sense
                                      setting this to `:always` even when `outfile` is `nothing` to force the computation of the local
                                      entropy and other thermodynamic quantities.
* `outfile` (default = `nothing`): the output file name. `nothing` means no output file is written. An empty string means using a default
                                   file name of the form: `"results_BPCR_N\$(N)_K\$(K)_M\$(M)_s\$(seed).txt"`.
* `outmessfiletmpl` (default = `nothing`): template file name for writing the messages at each focusing step. The file name should include a
                                           substring `"%gamma%"` which will be substituted with the value of the `γ` parameter at each step.
                                           If `nothing`, it is not used. If empty, the default will be used:
                                           `"messages_BPCR_N\$(N)_K\$(K)_M\$(M)_g%gamma%_s\$(seed).txt.gz"`.
                                           *Note*: this can produce a lot of fairly large files!.

The function returns three objects: the number of training errors, the messages and the patterns. The last two can be used as inputs to
successive runs of the algorithms, as the `initmessages` keyword argument and the `patternspec` argument, respectively.

Example of a run which solves a problem with `N * K = 1605` synapses with `K = 5` at `α = 0.3`:
```
julia> errs, messages, patterns = focusingBP(321, 5, 0.3, randfact=0.1, seed=135, max_iters=1, damping=0.5);
```
"""
function focusingBP(N::Integer, K::Integer,
                    initpatt::Union{AbstractString,Tuple{Vec2,Vec},Real,Patterns};

                    max_iters::Integer = 1000,
                    max_steps::Integer = typemax(Int),
                    seed::Integer = 1,
                    damping::Real = 0.0,
                    quiet::Bool = false,
                    accuracy1::Symbol = :accurate,
                    accuracy2::Symbol = :exact,
                    randfact::Real = 0.01,
                    fprotocol::FocusingProtocol = StandardReinforcement(1e-2),
                    ϵ::Real = 1e-3,
                    messfmt::Symbol = :tanh,
                    initmess::Union{Messages,Nothing,AbstractString} = nothing,
                    outatzero::Bool = true,
                    writeoutfile::Symbol = :auto, # note: ∈ [:auto, :always, :never]; auto => !outatzero && converged
                    outfile::Union{AbstractString,Nothing} = nothing, # note: "" => default, nothing => no output
                    outmessfiletmpl::Union{AbstractString,Nothing} = nothing) # note: same as outfile

    Random.seed!(seed)

    N > 0 || throw(ArgumentError("N must be positive; given: $N"))
    K > 0 || throw(ArgumentError("K must be positive; given: $K"))

    writeoutfile ∈ [:auto, :always, :never] || error("invalide writeoutfile, expected one of :auto, :always, :never; given: $writeoutfile")
    max_iters ≥ 0 || throw(ArgumentError("max_iters must be non-negative; given: $max_iters"))
    max_steps ≥ 0 || throw(ArgumentError("max_steps must be non-negative; given: $max_steps"))
    0 ≤ damping < 1 || throw(ArgumentError("damping must be ∈ [0,1); given: $damping"))
    0 ≤ randfact ≤ 1 || throw(ArgumentError("randfact must be ∈ [0,1]; given: $randfact"))
    accuracy1 ∈ [:exact, :accurate, :none] || error("accuracy1 must be one of :exact, :accurate, :none; given: $accuracy1")
    accuracy2 ∈ [:exact, :accurate, :none] || error("accuracy2 must be one of :exact, :accurate, :none; given: $accuracy2")

    accuracy1 == :exact && iseven(N) && throw(ArgumentError("when accuracy1==:exact N must be odd, given: $N"))
    accuracy2 == :exact && iseven(K) && throw(ArgumentError("when accuracy2==:exact K must be odd, given: $K"))

    messfmt ∈ [:tanh, :plain] || throw(ArgumentError("invalid messfmt, should be :tanh or :plain; given: $messfmt"))
    F = messfmt == :tanh ? MagT64 : MagP64

    if isa(initpatt, Real)
        initpatt ≥ 0 || throw(ArgumentError("invalide negative initpatt; given: $initpatt"))
        initpatt = (N, round(Int, K * N * initpatt))
    end

    patterns = Patterns(initpatt)

    M = patterns.M

    messages::Messages = initmess ≡ nothing ? Messages(F, M, N, K, randfact) :
                         isa(initmess, AbstractString) ? read_messages(initmess, F) :
                         initmess

    messages.N == N || throw(ArgumentError("wrong messages size, expected N=$N; given: $(messages.N)"))
    messages.K == K || throw(ArgumentError("wrong messages size, expected K=$K; given: $(messages.K)"))
    messages.M == M || throw(ArgumentError("wrong messages size, expected M=$M; given: $(messages.M)"))
    F = eltype(messages)

    params = Params{F}(damping, ϵ, NaN, max_iters, accuracy1, accuracy2, 0.0, 0.0, 0.0, quiet)

    outfile == "" && (outfile = "results_BPCR_N$(N)_K$(K)_M$(M)_s$(seed).txt")
    outmessfiletmpl == "" && (outmessfiletmpl = "messages_BPCR_N$(N)_K$(K)_M$(M)_g%gamma%_s$(seed).txt.gz")
    lockfile = "bpcomm.lock"
    if outfile ≢ nothing && writeoutfile ∈ [:always, :auto]
        !quiet && println("writing outfile $outfile")
        exclusive(lockfile) do
            !isfile(outfile) && open(outfile, "w") do f
                println(f, "#1=pol 2=y 3=β 4=S 5=q 6=q̃ 7=βF 8=𝓢ᵢₙₜ 9=Ẽ")
            end
        end
    end

    ok = true
    if initmess ≢ nothing
        errs = nonbayes_test(messages, patterns)
        !quiet && println("initial errors = $errs")

        outatzero && errs == 0 && return errs, messages, patterns
    end
    !quiet && K > 1 && (println("mags overlaps="); display(mags_symmetry(messages)[1]); println())

    for (it, (γ,y,β)) in enumerate(fprotocol)
        it > max_steps && break
        isfinite(β) && error("finite β not yet supported; given: $β")
        pol = mtanh(F, γ)
        params.pol = pol
        params.r = y - 1
        params.β = β
        set_outfields!(messages, patterns.output, params.β)
        ok = converge!(messages, patterns, params)
        !quiet && K > 1 && (println("mags overlaps="); display(mags_symmetry(messages)[1]); println())
        errs = nonbayes_test(messages, patterns)

        if writeoutfile == :always || (writeoutfile == :auto && !outatzero)
            S = compute_S(messages, params)
            q = compute_q(messages)
            q̃ = compute_q̃(messages, params)
            βF = free_energy(messages, patterns, params)
            Σint = -βF - γ * S

            !quiet && println("it=$it pol=$pol y=$y β=$β (ok=$ok) S=$S βF=$βF Σᵢ=$Σint q=$q q̃=$q̃ Ẽ=$errs\n")
            (ok || writeoutfile == :always) && outfile ≢ nothing && exclusive(lockfile) do
                open(outfile, "a") do f
                    println(f, "$pol $y $β $S $q $q̃ $βF $Σint $errs")
                end
            end
            if outmessfiletmpl ≢ nothing
                outmessfile = replace(outmessfiletmpl, "%gamma%" => γ)
                write_messages(outmessfile, messages)
            end
        else
            !quiet && println("it=$it pol=$pol y=$y β=$β (ok=$ok) Ẽ=$errs\n")
            errs == 0 && return 0, messages, patterns
        end
    end
    return errs, messages, patterns
end

end # module


text:
export
        get_actions!,
        tick!,
        simulate!

function get_actions!{A<:DriveAction, D<:DriverModel}(
    actions::Vector{A},
    scene::Scene,
    roadway::Roadway,
    models::Dict{Int, D}, # id → model
    )

    for (i,veh) in enumerate(scene)
        model = models[veh.def.id]
        observe!(model, scene, roadway, veh.def.id)
        actions[i] = rand(model)
    end

    actions
end

function tick!{A<:DriveAction, D<:DriverModel}(
    scene::Scene,
    roadway::Roadway,
    actions::Vector{A},
    models::Dict{Int, D}, # id → model
    )

    for (veh, action) in zip(scene, actions)
        model = models[veh.def.id]
        context = action_context(model)
        veh.state = propagate(veh, action, context, roadway)
    end

    scene
end
function tick!{A<:DriveAction, C<:ActionContext}(
    scene::Scene,
    roadway::Roadway,
    actions::Vector{A},
    contexts::Vector{C}
    )

    for (veh, action, context) in zip(scene, actions, contexts)
        veh.state = propagate(veh, action, context, roadway)
    end

    scene
end
function tick!{A<:DriveAction}(
    scene::Scene,
    roadway::Roadway,
    actions::Vector{A},
    context::ActionContext,
    )

    for (veh, action) in zip(scene, actions)
        veh.state = propagate(veh, action, context, roadway)
    end

    scene
end

"""
    Run a simulation and store the resulting scenes in the provided SceneRecord.
Only the ego vehicle is simulated; the other vehicles are as they were in the provided trajdata
Other vehicle states will be interpolated
"""
function simulate!(
    rec::SceneRecord,
    models::DriverModel,
    egoid::Int,
    trajdata::Trajdata,
    time_start::Float64,
    time_end::Float64;
    prime_history::Int=0, # no prime-ing
    scene::Scene = Scene(),
    )

    Δt = rec.timestep
    roadway = trajdata.roadway

    # prime with history
    empty!(rec)
    reset_hidden_state!(model)
    for h in 1:prime_history
        t = time_start - Δt * (prime_history - h + 1)
        update!(rec, get!(scene, trajdata, t))
        observe!(model, scene, roadway, egoid)
    end

    # add current frame
    update!(rec, get!(scene, trajdata, time_start))
    observe!(model, scene, roadway, egoid)

    # run simulation
    t = time_start
    context = action_context(model)
    prev_ego_state = get_by_id(scene, egoid).state
    while t < time_end

        # pull orig scene
        get!(scene, trajdata, t)

        # propagate ego vehicle and set
        ego_action = rand(model)
        ego_veh = get_by_id(scene, egoid)
        ego_veh.state = prev_ego_state
        prev_ego_state = ego_veh.state = propagate(ego_veh, ego_action, context, roadway)

        # update record
        update!(rec, scene)

        # observe
        observe!(model, scene, roadway, ego_veh.def.id)

        # update time
        t += Δt
    end

    rec
end


text:
using CUDA.CUBLAS
using CUDA.CUDNN

import CUDA: libcudnn, @checked, @argout
import CUDA.CUDNN:
    @check, @runtime_ccall, @workspace, handle,
    cudnnStatus_t, cudnnConvolutionDescriptor_t,
    ConvDesc, TensorDesc, FilterDesc,
    CuArray, CuVector, CUDNNFloat, cdsize,
    cudnnConvolutionForward,
    cudnnGetConvolutionForwardWorkspaceSize, cudnnConvolutionFwdAlgo_t,
    cudnnConvolutionBackwardData,
    cudnnGetConvolutionBackwardDataWorkspaceSize, cudnnConvolutionBwdDataAlgo_t,
    cudnnConvolutionBackwardFilter,
    cudnnGetConvolutionBackwardFilterWorkspaceSize, cudnnConvolutionBwdFilterAlgo_t
import NNlib: depthwiseconv!, ∇depthwiseconv_filter!, ∇depthwiseconv_data!

const CuOrVector = Union{CuVector, Vector}
const CuOrMatrix = Union{CuMatrix, Matrix}
const CuOrArray = Union{CuArray, Array}

# Implement conversion to dense, diagonal matrix.

diagonal(x::CuArray{T, 1}) where T<:Real = convert(CuArray, Diagonal(x))

# Implement GPU support for depthwise separable convolutions.

@checked function cudnnGetConvolutionGroupCount(convDesc, count)
    @runtime_ccall(
        (:cudnnGetConvolutionGroupCount, libcudnn()),
        cudnnStatus_t,
        (cudnnConvolutionDescriptor_t, Ptr{Cint}),
        convDesc,
        count
    )
end

@checked function cudnnSetConvolutionGroupCount(convDesc, count)
    @runtime_ccall(
        (:cudnnSetConvolutionGroupCount, libcudnn()),
        cudnnStatus_t,
        (cudnnConvolutionDescriptor_t, Cint),
        convDesc,
        count
    )
end

function ConvDesc(T, cdims::DepthwiseConvDims)
    cd = Ref{cudnnConvolutionDescriptor_t}()
    CUDNN.cudnnCreateConvolutionDescriptor(cd)
    N = NNlib.spatial_dims(cdims)
    CUDNN.cudnnSetConvolutionNdDescriptor(
        cd[],
        N,
        # Asymmetric padding is not supported.
        cdsize(NNlib.padding(cdims)[1:2:end], N),
        cdsize(NNlib.stride(cdims), N),
        cdsize(NNlib.dilation(cdims), N),
        NNlib.flipkernel(cdims),
        CUDNN.cudnnDataType(T)
    )
    # Set number of groups equal to number of channels to get a depthwise
    # separable convolution.
    cudnnSetConvolutionGroupCount(cd[], NNlib.channels_in(cdims))
    this = ConvDesc(cd[])
    CUDNN.finalizer(CUDNN.unsafe_free!, this)
    return this
end

function cudnnConvolutionForward(
    y::CuArray{T, N},
    x::CuArray{T, N},
    w::CuArray{T, N},
    cdims::DepthwiseConvDims;
    algo=0,
    alpha=1,
    beta=0
) where {T, N}
    @workspace size = @argout(
        cudnnGetConvolutionForwardWorkspaceSize(
            handle(),
            TensorDesc(x),
            FilterDesc(w),
            ConvDesc(T, cdims),
            TensorDesc(y),
            cudnnConvolutionFwdAlgo_t(algo),
            out(Ref{Csize_t}())
        )
    )[] workspace -> begin
        cudnnConvolutionForward(
            handle(),
            Ref(T(alpha)),
            TensorDesc(x),
            x,
            FilterDesc(w),
            w,
            ConvDesc(T,cdims),
            cudnnConvolutionFwdAlgo_t(algo),
            workspace,
            sizeof(workspace),
            Ref(T(beta)),
            TensorDesc(y),
            y
        )
    end
    return y
end

function cudnnConvolutionBackwardData(
    dx::CuArray{T, N},
    w::CuArray{T, N},
    dy::CuArray{T, N},
    cdims::DepthwiseConvDims;
    algo=0,
    alpha=1,
    beta=0
) where {T, N}
    @workspace size = @argout(
        cudnnGetConvolutionBackwardDataWorkspaceSize(
            handle(),
            FilterDesc(w),
            TensorDesc(dy),
            ConvDesc(T, cdims),
            TensorDesc(dx),
            cudnnConvolutionBwdDataAlgo_t(algo),
            out(Ref{Csize_t}())
        )
    )[] workspace -> begin
        cudnnConvolutionBackwardData(
            handle(),
            Ref(T(alpha)),
            FilterDesc(w),
            w,
            TensorDesc(dy),
            dy,
            ConvDesc(T, cdims),
            cudnnConvolutionBwdDataAlgo_t(algo),
            workspace,
            sizeof(workspace),
            Ref(T(beta)),
            TensorDesc(dx),
            dx
        )
    end
    return dx
end

function cudnnConvolutionBackwardFilter(
    dw::CuArray{T, N},
    x::CuArray{T, N},
    dy::CuArray{T, N},
    cdims::DepthwiseConvDims;
    algo=0,
    alpha=1,
    beta=0
) where {T, N}
    @workspace size = @argout(
        cudnnGetConvolutionBackwardFilterWorkspaceSize(
            handle(),
            TensorDesc(x),
            TensorDesc(dy),
            ConvDesc(T, cdims),
            FilterDesc(dw),
            cudnnConvolutionBwdFilterAlgo_t(algo),
            out(Ref{Csize_t}())
        )
    )[] workspace -> begin
        cudnnConvolutionBackwardFilter(
            handle(),
            Ref(T(alpha)),
            TensorDesc(x),
            x,
            TensorDesc(dy),
            dy,
            ConvDesc(T, cdims),
            cudnnConvolutionBwdFilterAlgo_t(algo),
            workspace,
            sizeof(workspace),
            Ref(T(beta)),
            FilterDesc(dw),
            dw
        )
    end
    return dw
end

function depthwiseconv!(
    y::CuArray{T},
    x::CuArray{T},
    w::CuArray{T},
    cdims::DepthwiseConvDims;
    alpha=1,
    algo=0
) where T<:CUDNNFloat
    cudnnConvolutionForward(
        y,
        x,
        w,
        cdims;
        alpha=alpha,
        algo=algo
    )
    return y
end

function ∇depthwiseconv_filter!(
    dw::CuArray{T},
    x::CuArray{T},
    dy::CuArray{T},
    cdims::DepthwiseConvDims;
    alpha=1,
    algo=0
) where T<:CUDNNFloat
    cudnnConvolutionBackwardFilter(
        dw,
        x,
        dy,
        cdims;
        alpha=alpha,
        algo=algo
    )
    return dw
end

function ∇depthwiseconv_data!(
    dx::CuArray{T},
    dy::CuArray{T},
    w::CuArray{T},
    cdims::DepthwiseConvDims;
    alpha=1,
    algo=0
) where T<:CUDNNFloat
    cudnnConvolutionBackwardData(
        dx,
        w,
        dy,
        cdims;
        alpha=alpha,
        algo=algo
    )
end


text:
function hole_correction_pair_correlation(x1::AbstractVector{T},s1::Specie{T}, x2::AbstractVector{T},s2::Specie{T}) where T
    overlapping = norm(x1 - x2) > outer_radius(s1) * s1.exclusion_distance + outer_radius(s2) * s2.exclusion_distance

    return  overlapping ? one(T) : zero(T)
end


text:
using LinearAlgebra
import StatsBase.sample

const Dic = AbstractDict
abstract type AbstractGBM{T, I} end

"""
    Response structure, including spin configuration (:states) and energies.
"""
mutable struct Response{T, I} <: AbstractGBM{T, I}
    W::Matrix{T}
    b::Matrix{T}
    states::Vector{Vector{I}}
    energies::Vector{T}
    indices::Vector
    sample::Dict
end

"""
    Response(n_unit::Int)

Construct response structure with number of unit `n_unit`.

 * T - type of parameters [interacts and biases;] (default: Float64)
 * I - type of spin states (default: Int64)
"""
function Response(n_unit::Int64)
    Response{Float64, Int64}(
        rand(n_unit, n_unit),
        rand(n_unit, 1),
        Vector{Vector}[],
        [],
        [],
        Dict()
    )
end

# initial setting

"""
    create_response(linear::AbstractDict, quad::AbstractDict)

Configure the response structure.
"""
function create_response(linear::Dic, quad::Dic)

    Jkeys = collect.(keys(quad))
    key_list = collect(Set(vcat(keys(linear)..., Jkeys...)))
    elements_match(key_list)
    sort!(key_list)
    linear, quad = relabel_dict(key_list, linear, quad)

    # create field & response
    n_unit = length(key_list)
    resp = Response(n_unit)
    field = zeros((n_unit, 1))
    field[collect(keys(linear))] .= values(linear)

    copy!(resp.W, trans(quad, n_unit))
    copy!(resp.b, field)
    resp.indices = key_list
    return resp
end

"""
    relabel_dict(key_lst::Vector, linear::AbstractDict, quad::AbstractDict)

Replace the keys of the dictionary to unit numbers.
"""
function relabel_dict(key_lst::Vector, linear::Dic, quad::Dic)

    num_nodes = collect(1:length(key_lst))
    tab = Dict(zip(sort(key_lst), num_nodes))
    bias, interact = Dict(), Dict()

    for key in keys(quad)
        key1, key2 = key
        interact[(tab[key1], tab[key2])] = quad[key]
    end

    for key in keys(linear)
        bias[tab[key]] = linear[key]
    end

    return bias, interact
end

# annealing algorithm

"""
    get_energy(resp::Response{T}, state::Vector)

Calculate the energy in the input `state`.
"""
function get_energy(resp::Response{T}, state::Vector) where T
    E = -dot(vec(resp.b), vec(state)) # linear term
    E -= state' * triu(resp.W, 1) * state * 0.5 # interaction term
    return E
end

"""
    accept_rate(dE::T, beta::T)

Return the accpetance rate from energy difference `dE`.
"""
function accept_rate(dE::T, beta::T) where T
    return min(1, exp(-beta * dE))
end

"""
    get_neighbor(old_state::Vector)

Return the proposal state calculated from `old_state`.
"""
function get_neighbor(old_state::Vector)
    new_state = copy(old_state)
    idx = sample(1:length(old_state), 1)
    new_state[idx] *= -1 # single spin flip
    return new_state
end

"""
    anneal(linear::AbstractDict, quad::AbstractDict; [, opts])

Solve forward ising problem expressed by `linear` and `quad` biases.

# Arguments
- `linear::Dict`: Linear biases meaning that magnetic field.
- `quad::Dict`: Quadratic biases meaning that spin-spin interaction.

Options that can be provided in the `opts` dictionary:

* `:beta_min` - Initial (minimal) value of inverse temperature. (default: 5.0)
* `:beta_max` - Final (maximal) value of inverse temperature. (default: 15.0)
* `:n_sweep` - Number of division between `beta_min` and `beta_max`. (default: 1000)
* `:n_read` - Number of run repetitions of Simulated Annealing. (default: 1)

# Examples
```jldoctest
julia> using InverseIsing

julia> h = Dict(:a => 1) # Setting field biases.

julia> J = Dict((:a, :b) => -1) # Setting interaction.

julia> response = anneal(h, J)

julia> response.sample
Dict{Symbol,Int64} with 2 entries:
  :a => 1
  :b => -1
```
"""
function anneal(linear::Dic, quad::Dic; opts...)

    ctx = copy(Dict{Symbol, Real}(opts))
    option_checker!(ctx)
    num_reads = ctx[:n_read]
    num_sweeps = ctx[:n_sweep] # number of division between beta_min & beta_max
    beta_min = ctx[:beta_min] # beta_min: initial beta
    beta_max = ctx[:beta_max] # beta_max: final beta

    resp = create_response(linear, quad)
    state = ones(size(resp.W, 1))
    energy = zero(eltype(resp.W))

    for step in 1:num_reads
        @simd for beta in range(beta_min, beta_max, length = num_sweeps)

            # get neighbor
            new_state = get_neighbor(state)

            # check energy level
            energy = get_energy(resp, state)
            dE = get_energy(resp, new_state) - energy

            # MCMC algorithm
            if rand() < accept_rate(dE, beta)
                state = new_state
            end
        end

        # append some info
        push!(resp.states, Int.(state))
        push!(resp.energies, energy)

    end
    optimal_idx = argmin(resp.energies)
    resp.sample = Dict(zip(resp.indices, resp.states[optimal_idx]))
    return resp

end

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
###############################################
#FEniCS tutorial demo program: Linear elastic problem.
#  -div(σ(u)) = f
#The model is used to simulate an elastic beam clamped at
#its left end and deformed under its own weight.
###############################################

module ft06

using FenicsPy

# Scaled variables
L = 1; W = 0.2
μ = 1
rho = 1
δ = W/L
γ = 0.4*δ^2
beta = 1.25
λ = beta

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, "P", 1)

# Define boundary condition
tol = 1e-14

clamped_boundary = string("on_boundary and x[0] < ", tol)

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

ϵ(u) =  0.5*(nabla_grad(u) + nabla_grad(u).T)

σ(u) = λ*nabla_div(u)*Identity(d) + 2*μ*ϵ(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*γ))
T = Constant((0, 0, 0))
a = inner(σ(u), ϵ(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = FeFunction(V)
solve(a == L, u, bc)

# Compute stress
s = σ(u) - (1.0/3)*tr(σ(u))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3.0/2*inner(s, s))
V = FunctionSpace(mesh, "P", 1)
von_Mises = project(von_Mises, V)

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
println("min/max u:",
      min(array(u_magnitude)...),
      max(array(u_magnitude)...)
)

# Save solution to file in VTK format
File("elasticity/displacement.pvd") << u
File("elasticity/von_mises.pvd") << von_Mises
File("elasticity/magnitude.pvd") << u_magnitude

# Plot solution
#plot(u, title="Displacement", mode="displacement")

# Plot stress
#plot(von_Mises, title="Stress intensity")

# Plot magnitude of displacement
#plot(u_magnitude, "Displacement magnitude")

end # module ft06


text:
"""
Obtain coefficients (A, B) of the function Se = B(x - A)^2/x for
Se(1.2) = B(1.2 - A)^2/1.2 and Se(1.0) = B(1.0 - A)^2/1.0 as:
Se(1.0) = (Se(1.2) * 1.2) /(1.2 - A)^2 * (1.0 - A)^2/1.0 that yields
(1.2 - A)^2 Se(1.0) = Se(1.2) * 1.2 * (1.0 - A)^2 or expanding:
(1.2 * Se(1.2) - Se(1.0)) A^2 + (2.4 Se(1.0) - 2 * 1.2 * Se(1.2)) A + (1.2 * Se(1.2) - 1.44 Se(1.0)) = 0
and uses the negative solution of the quadratic equation.
"""
function get_quadratic_saturation(Se::Tuple{Float64, Float64})
    return calculate_saturation_coefficients((1.0, 1.2), Se)
end

"""
Obtain coefficients (A, B) of the function Se = Bx^A for
Se(1.2) = B(1.2)^A and Se(1.0) = B(1.0)^A as:
B = Se(1.0) and hence
(1.2)^A = Se(1.2)/B -> A = log(Se(1.2)/B) / log(1.2)
"""
function get_exponential_saturation(Se::Tuple{Float64, Float64})
    if Se[1] == 0.0 || Se[2] == 0.0
        @warn "$(Se[1]) or $(Se[2]) equals to zero. Ignoring saturation."
        return (0.0, 0.0)
    end
    if Se[2] <= Se[1]
        throw(
            IS.ConflictingInputsError(
                "Se(1.2) <= Se(1.0). Saturation data is inconsistent.",
            ),
        )
    end
    B = Se[1]
    A = log(Se[2] / B) / log(1.2)

    return (A, B)
end


text:
using ITensorNetworkAD
using Test

@testset "ITensorNetworks.jl" begin
  for filename in ["basic_ops.jl", "subnetwork.jl"]
    println("Running $filename in ITensorNetworks.jl")
    include(filename)
  end
end


text:
export EGSParticle, Latch, BitRegions
export @BitRegions

include("particle.jl")
include("io.jl")


text:
module TestCategoricalAlgebra

using Test

@testset "ShapeDiagrams" begin
  include("ShapeDiagrams.jl")
end

@testset "Matrices" begin
  include("Matrices.jl")
end

@testset "FinSets" begin
  include("FinSets.jl")
end

@testset "FinRelations" begin
  include("FinRelations.jl")
end

@testset "Permutations" begin
  include("Permutations.jl")
end

@testset "CSets.jl" begin
  include("CSets.jl")
  include("Graphs.jl")
end

end


text:
# work in progress...

bessel1arg = [
(:airyai, :(airyaiprime(x)), (-Inf,Inf)),
(:airyaiprime, :(x.*airyai(x)), (-Inf,Inf)),
(:airybi, :(airybiprime(x)), (-Inf,Inf)),
(:airybiprime, :(x.*airybi(x)), (-Inf,Inf)),
(:airyprime, :(x.*airyai(x)), (-Inf,Inf)),
(:besselj0, :(-besselj1(x)), (-Inf,Inf)),
(:besselj1, :((besselj0(x)-besselj(2,x))/2), (-Inf,Inf)),
(:bessely0, :(-bessely1(x)), (0,Inf)),
(:bessely1, :((bessely0(x)-bessely(2,x))/2), (0,Inf)),
]

for (f,g,r) in bessel1arg
    @eval @primitive $f(x),dy,y  (dy.*($g))
    addtest1(f,r)
end

@primitive besselj(nu,x)  error("No gradient for besselj")  error("No gradient for besselj")
@primitive bessely(nu,x)  error("No gradient for besselj")  error("No gradient for besselj")

# The following defined for (nu,x) where nu can be any integer or
# positive real

# bessel2arg = [
# (besseli, :((besseli(nu-1,x)+besseli(nu+1,x))/2)), # x>=0 || isinteger(nu)
# (besselix, :((besselix(nu-1,x)+besselix(nu+1,x))/2)), # x>=0 || isinteger(nu)
# (besselj, :((besselj(nu-1,x)-besselj(nu+1,x))/2)), # x>=0 || isinteger(nu)
# (besseljx, :((besseljx(nu-1,x)-besseljx(nu+1,x))/2)), # x>=0 || isinteger(nu)
# #besselk, # x>=0
# #besselkx, # x>=0
# #bessely, # x>=0
# #besselyx, # x>=0
# ]

# for (f,g) in bessel2arg
#     @eval @primitive $f(nu,x),dy,y  nothing  unbroadcast(x, dy.*($g))
#     addtest2(f,(0,Inf),(-Inf,Inf))
# end

# bessel2arg2 = [
# airy,
# airyx,
# # besselh  # has (nu,k,z) and (nu,z) argtypes where k=int
# # besselhx # (nu,k,x), no vectorization
# # hankelh1 # (nu,x), vectorized
# # hankelh1x # (nu,x), vectorized
# # hankelh2
# # hankelh2x
# ]
# for f in bessel2arg2
#     @eval @zerograd $f(x1,x2)
# end


text:
using ConstraintDomains
using Documenter

makedocs(;
    modules=[ConstraintDomains],
    authors="Jean-François Baffier",
    repo="https://github.com/JuliaConstraints/ConstraintDomains.jl/blob/{commit}{path}#L{line}",
    sitename="ConstraintDomains.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://JuliaConstraints.github.io/ConstraintDomains.jl",
        assets = ["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(; repo="github.com/JuliaConstraints/ConstraintDomains.jl.git", devbranch = "main")