text:
using MathOptInterface
const MOI = MathOptInterface

# Test, whether the callback of the separator is called.
@testset "DummySepa (no separation)" begin
    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add the separator
    sepa = DummySepa.Sepa()
    SCIP.include_sepa(mscip, sepa)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # the separator is called
    @test sepa.called >= 1

    # free the problem
    finalize(mscip)
end

# Test, whether adding cuts in `exec_lp` via `add_cut_sepa` works [1/2].
@testset "AddSingleCut (cut off one optimal solution)" begin
    atol, rtol = 1e-6, 1e-6

    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add the separator
    varrefs = [SCIP.VarRef(x.value)]
    coefs = [1.0]
    sepa = AddSingleCut.Sepa(mscip, varrefs, coefs, 0.0, 0.0)
    SCIP.include_sepa(mscip, sepa)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # SCIP found the single remaining optimal solution
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol

    # free the problem
    finalize(mscip)
end

# Test, whether adding cuts in `exec_lp` via `add_cut_sepa` works [2/2].
@testset "AddSingleCut (cut off another optimal solution)" begin
    atol, rtol = 1e-6, 1e-6

    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add the separator
    varrefs = [SCIP.VarRef(y.value)]
    coefs = [1.0]
    sepa = AddSingleCut.Sepa(mscip, varrefs, coefs, 0.0, 0.0)
    SCIP.include_sepa(mscip, sepa)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # SCIP found the single remaining optimal solution
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol

    # free the problem
    finalize(mscip)
end

# Test, whether we can cut the optimal solution.
@testset "AddSingleCut (too strong cut)" begin
    atol, rtol = 1e-6, 1e-6

    # create an empty problem
    optimizer = SCIP.Optimizer()
    mscip = optimizer.mscip
    sepa_set_scip_parameters((par,val) -> SCIP.set_parameter(mscip, par, val))

    # add variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # add constraint: x + y ≤ 1.5
    MOI.add_constraint(optimizer,
                       MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
                       MOI.LessThan(1.5))

    # maximize x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add the separator
    varrefs = [SCIP.VarRef(x.value), SCIP.VarRef(y.value)]
    coefs = [1.0, 1.0]
    sepa = AddSingleCut.Sepa(mscip, varrefs, coefs, 0.0, 0.0)
    SCIP.include_sepa(mscip, sepa)

    # solve the problem
    SCIP.@SC SCIP.SCIPsolve(mscip.scip[])

    # SCIP found the non-optimal solution, that remains after the cut.
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol

    # free the problem
    finalize(mscip)
end


text:
module TestCluster
using Test

include("../src/cluster.jl")

@testset "Partition 1112                 " begin
    nums = [1, 1, 1, 2]
    @test all(nums .> 0)
    p = Cluster.PartitionIndex(nums)
    @test length(p) == 2
    @test Set(p[1]) == Set([1, 2, 3])
    @test Set(p[2]) == Set([4])
end

@testset "Test Partition5                " begin
    nums = [2, 3, 1, 1, 3]
    @test all(nums .> 0)
    p = Cluster.PartitionIndex(nums)
    @test length(p) == 3
    for i in 1:3
        @test Set(findall(nums .== i)) == Set(p[i])
    end
end

@testset "Test Cluster                   " begin
    neigh = Dict(1 => [2, 3],
                 2 => [1, 4],
                 3 => [1, 4],
                 4 => [2, 3])
    x = [ 0.0  0.1
         -0.1  1.1]
    @test size(x) == (2,2)
    p = Cluster.cluster(x, neigh, eps=0.2)
    @test length(p) == 2
    @test isperm(p.V)
    @test Set(p[1]) == Set([4])
    @test Set(p[2]) == Set([1,2,3])
end

end

ext: jl

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
## ---------  Table that shows how each parameter affects each deviation

"""
	$(SIGNATURES)

Table that shows how each parameter affects each deviation.
Rows are deviations. Columns are parameters.
"""
mutable struct ChangeTable{T1 <: AbstractFloat}
    # Name of each parameter
    paramNameV
    # Name of each deviation
    devNameV
    # Vector of intial deviations
    dev0V :: Vector{T1}
    # Matrix of deviations by [deviation, parameter]
    devM :: Matrix{T1}
    # Scalar deviations for all parameters
    scalarDevV :: Vector{T1}
    scalarDev0 :: T1
end

"""
	$(SIGNATURES)

Initialize a `ChangeTable` from a [`DevVector`](@ref).

# Arguments
- `d`: the `DevVector`
- `paramNameV`: the names of all parameters to be displayed
- `scalarDev`: the scalar deviation for the initial point
"""
function ChangeTable(d :: DevVector, paramNameV; scalarDev = nothing)
    nParams = length(paramNameV);
    devNameV = collect(keys(scalar_dev_dict(d)));
    nDev = length(devNameV);
    dev0V = zeros(nParams);
    devM = zeros(nDev, nParams);
    scalarDevV = zeros(nParams);
    scalarDev0 = 0.0;
    ct = ChangeTable(paramNameV, devNameV, dev0V, devM, scalarDevV, scalarDev0);
    ct.scalarDev0, ct.dev0V = changes_one_param(ct, d; scalarDev = nothing);
    return ct
end

dev_names(ct :: ChangeTable) = ct.devNameV;
param_names(ct :: ChangeTable) = ct.paramNameV;
n_params(ct :: ChangeTable) = length(param_names(ct));
n_devs(ct :: ChangeTable) = length(dev_names(ct));

function validate_change_table(ct :: ChangeTable)
    isValid = true;
    if !isequal(size(ct.devM),  (n_devs(ct), n_params(ct)))
        isValid = false;
        @warn "Invalid size of devation matrix: $(size(ct.devM))"
    end
    return isValid
end

# Compute changes for one parameter
function changes_one_param(ct :: ChangeTable, d :: DevVector;
    scalarDev = nothing)

    devDict = scalar_dev_dict(d);
    devV = zeros(n_devs(ct))
    for (j, key) in enumerate(dev_names(ct))
        devV[j] = devDict[key];
    end

    if isnothing(scalarDev)
        scalarDev = scalar_deviation(d);
        @assert isapprox(scalarDev, sum(devV), atol = 1e-3)
    end

    return scalarDev, devV
end

"""
	$(SIGNATURES)

Set values for one parameter. After the `ChangeTable` has been initialized, the model is solved with a different values for paramete `iParam`. This results in the deviations given by `d`.
"""
function set_param_values!(ct :: ChangeTable, iParam :: Integer, d :: DevVector;
    scalarDev = nothing)

    ct.scalarDevV[iParam], devV = changes_one_param(ct, d; scalarDev = scalarDev);
    ct.devM[:, iParam] = devV;
    return nothing
end

"""
	$(SIGNATURES)

Make a formatted table that shows how each parameter changes each deviation.
Returns header row, header column, table with values (as strings).
"""
function make_table(ct :: ChangeTable)
    @assert validate_change_table(ct)
    # The `ChangeTable` stores levels. We want differences relative to initial point.
    devChangeM = ct.devM .- ct.dev0V;
    dScalarDevV = ct.scalarDevV .- ct.scalarDev0;
    # Col 1: initial point
    dev0V = [ct.scalarDev0; ct.dev0V];
    # Row 1: scalar devs
    valueM = [dev0V  [dScalarDevV';  devChangeM]];
    @assert size(valueM) == (n_devs(ct) + 1,  n_params(ct) + 1)
    valueM = string.(round.(valueM, digits = 3));

    hdColumn = ["Scalar"; string.(dev_names(ct))];

    hdRow = [" ";  "Baseline";  string.(param_names(ct))];
    return  hdRow,  hdColumn,  valueM
end

"""
	$(SIGNATURES)

Show a `ChangeTable` as a `PrettyTable`.
"""
function show_table(ct :: ChangeTable, io = stdout)
    hdRow, hdColumn, valueM = make_table(ct);
    pretty_table(io, [hdColumn  valueM],  hdRow);
end

"""
	$(SIGNATURES)

Find parameters where scalar deviation changes less than `rtol`.
"""
function find_unchanged_devs(ct :: ChangeTable; rtol = 0.01)
    relDevV = ct.scalarDevV ./ ct.scalarDev0 .- 1;
    return param_names(ct)[relDevV .< rtol]
end

"""
	$(SIGNATURES)

Make table with largest individual deviation changes for each parameter.
Rows are parameters.
Can also show the transposed version: the `n` parameters with the largest effects on each devation.
"""
function make_largest_change_table(ct :: ChangeTable, n :: Integer; transposed :: Bool = false)
    @assert validate_change_table(ct)

    if transposed
        devChangeM = (ct.devM .- ct.dev0V);
        rowNameV = dev_names(ct);
        colNameV = string.(param_names(ct));
    else
        # Deviation changes by [parameter, deviation]
        devChangeM = (ct.devM .- ct.dev0V)';
        colNameV = dev_names(ct);
        rowNameV = string.(param_names(ct));
    end
    nRows = length(rowNameV);
    nCols = length(colNameV);
    @check n <= nCols
    @check isequal((nRows, nCols),  size(devChangeM))

    # Records `n` largest deviation changes for each column of `devChangeM`
    changeM = Matrix{String}(undef, nRows, n);
    for j = 1 : nRows
        # Largest deviation changes
        idxV = sortperm(abs.(devChangeM[j,:]), rev = true);
        for k = 1 : n
            idx = idxV[k];
            devChange = round(devChangeM[j, idx], digits = 3);
            changeM[j,k] = "$(colNameV[idx]): $devChange";
        end
    end

    hdRow = ["$j" for j = 0 : n];
    hdRow[1] = " ";

    return  hdRow,  rowNameV,  changeM
end

"""
	$(SIGNATURES)

Show a table with the largest deviation changes for each parameter.
"""
function show_largest_change_table(ct :: ChangeTable, n :: Integer; io = stdout,
    transposed :: Bool = false)
    hdRow, hdColumn, valueM = make_largest_change_table(ct, n; transposed = transposed);
    pretty_table(io, [hdColumn  valueM],  hdRow);
end

# ---------------


text:
module Hello

using Comm

function main()
    println("Hello, World! This is process $(myproc()) of $(nprocs()).")
end

run_main(main)

end


text:
using POMDPToolbox
using POMDPModels
using Base.Test

############################################################
################# TESTS PARTICLE FILTER ####################
############################################################

pomdp = TigerPOMDP()
s = false
n = 10000 # num particles
bu = SIRParticleUpdater(pomdp, n)

srand(2)
b = initialize_belief(bu, initial_state_distribution(pomdp))

a = 0
o = false

srand(1)
bp = update(bu, b, a, o)
w1 = pdf(bp, s)
srand(1)
bp = update(bu, b, a, o)
bp.keep_dict = false
w2 = pdf(bp, s)

bp.keep_dict = false
os = [true, false, true, false, true, false, true, false, true, false]
b = initialize_belief(bu, initial_state_distribution(pomdp))
for o in os
    b = update(bu, b, a, o)
end
w3 = pdf(b, s)
s

s = rand(MersenneTwister(0), b, s)

# dict and array method should give same pdf
@test isapprox(w1, w2, atol=0.01)
# after one update should have ~0.85 prob
@test isapprox(w1, 0.85, atol=0.01)
# after 10 of each observations should be back to uniform belief
@test isapprox(w3, 0.5, atol=0.05)


text:
"""
$(SIGNATURES)

Stores a convex planar polygon in with at most 8 sides. The polygon should be in 3 dimensions (Cartesian coordinates)
or 4 dimensions (tetrahedral coondinates).
"""
struct poly_eight{N,T}
    n::Int64
    v::NTuple{8,SVector{N,T}}
    function poly_eight{N,T}() where {N,T}
        return new{N,T}(0)
    end
    function poly_eight(n::Int64, v::NTuple{8,SVector{N,T}}) where {N,T}
        return new{N,T}(n, v)
    end
    function poly_eight(v::NTuple{3,SVector{N,T}}) where {N,T}
        return new{N,T}(3, (v[1], v[2], v[3], v[1], v[1], v[1], v[1], v[1]))
    end
    function poly_eight(v::NTuple{4,SVector{N,T}}) where {N,T}
        return new{N,T}(4, (v[1], v[2], v[3], v[4], v[1], v[1], v[1], v[1]))
    end
end

@inline Base.isempty(p::poly_eight) = (p.n == 0)
@inline Base.length(p::poly_eight) = p.n
@inline Base.getindex(p::poly_eight, k::Int64) = p.v[k]

"""
$(SIGNATURES)

Finds the centroid of a 3 dimensional poly_eight by dividing the area into triangles. This funciton requires a reference
normal n̂.
"""
function MathKernel.centroid(p_new::poly_eight{3,T}, n̂::SVector{3,T}) where {T}
	cart_a = p_new.v[1]
	cart_c = p_new.v[2]  # because c becomes b
	cum_sum  = zero(T)
	cum_prod = zeros(SVector{3,T})
	for k = 3:p_new.n
		cart_b = cart_c
		cart_c = p_new.v[k]
		area_ = triangle_area((cart_a, cart_b, cart_c), n̂)
		cum_prod += area_ * centroid(cart_a, cart_b, cart_c)
		cum_sum  += area_
	end
	if cum_sum == 0.0  # just return an arbitrary vertex if the area is zero
		return cum_sum, cart_a
	else
		return cum_sum, cum_prod / cum_sum
	end
end

"""
$(SIGNATURES)

One pads each element of a 3 dimensional polygon and then multiplies it by a 4x4 matrix m. This funciton is used to
convert each element of a Cartesian polygon into tetrahedral coordinates.
"""
function MathKernel.one_pad_then_mul(m::SMatrix{4,4,T1,16}, p::poly_eight{3,T2}) where {T1,T2}
    t1 = one_pad_then_mul(m, p[1])
    t2 = one_pad_then_mul(m, p[2])
    t3 = one_pad_then_mul(m, p[3])
    t4 = one_pad_then_mul(m, p[4])
    length_p = length(p)
    if length_p <= 4
        return poly_eight(length_p, (t1, t2, t3, t4, t1, t1, t1, t1))
    else
        t5 = one_pad_then_mul(m, p[5])
        t6 = one_pad_then_mul(m, p[6])
        t7 = one_pad_then_mul(m, p[7])
        t8 = one_pad_then_mul(m, p[8])
        return poly_eight(length_p, (t1, t2, t3, t4, t5, t6, t7, t8))
    end
end

"""
$(SIGNATURES)

Multiplies a 4 dimensional polygon by a 4x4 matrix and then unpads the 4 dimensional polygon. This function is used to
convert a polygon in tetrahedral coordinates to a polygon in Cartesian coordinates.
"""
function MathKernel.mul_then_un_pad(m::SMatrix{4,4,T1,16}, p::poly_eight{4,T2}) where {T1,T2}
    t1 = mul_then_un_pad(m, p[1])
    t2 = mul_then_un_pad(m, p[2])
    t3 = mul_then_un_pad(m, p[3])
    t4 = mul_then_un_pad(m, p[4])
    length_p = length(p)
    if length_p <= 4
        return poly_eight(length_p, (t1, t2, t3, t4, t1, t1, t1, t1))
    else
        t5 = mul_then_un_pad(m, p[5])
        t6 = mul_then_un_pad(m, p[6])
        t7 = mul_then_un_pad(m, p[7])
        t8 = mul_then_un_pad(m, p[8])
        return poly_eight(length_p, (t1, t2, t3, t4, t5, t6, t7, t8))
    end
end

"""
$(SIGNATURES)

Sets small coondinates of a polygon in tetrahedral coordinates to zero. This function addresses a specific probability
one degeneracy that results when clipping an edge in the same place twice.
"""
function zero_small_coordinates(p::poly_eight{4,T}) where {T}
    function zero_if_small(x::SVector{4,T}) where {T}
        abs_x = 1.0e-14 .< abs.(value.(x))
        return x .* abs_x
    end

    v1 = zero_if_small(p[1])
    v2 = zero_if_small(p[2])
    v3 = zero_if_small(p[3])
    v4 = zero_if_small(p[4])
    length_p = length(p)
    if length_p <= 4
        return poly_eight(length_p, (v1, v2, v3, v4, v1, v1, v1, v1))
    else
        v5 = zero_if_small(p[5])
        v6 = zero_if_small(p[6])
        v7 = zero_if_small(p[7])
        v8 = zero_if_small(p[8])
        return return poly_eight(length_p, (v1, v2, v3, v4, v5, v6, v7, v8))
    end
end


text:
"""
Function to perform a single boostingstep, which is called by the wrapper-function 'boosting'
Arguments:
myscratch::Boostscratch -> Container with data which are relevant to perform the boostinsteps
"""
function boost!(myscratch::Boostscratch)
	myscratch.actualscore = myscratch.actualnom.^2 # ./varx
	myscratch.actualsel = indmax(myscratch.actualscore)
	myscratch.selections[myscratch.actualstepno] = myscratch.pooledunibeta.labels[myscratch.actualsel]
	myscratch.actualupdates[myscratch.actualstepno] = myscratch.nu * myscratch.actualnom[myscratch.actualsel]
	myscratch.actualbeta[myscratch.actualsel] += myscratch.actualupdates[myscratch.actualstepno]
	covar = getcovarfromtriangular(myscratch.pooledcovarmat.covarmat, findfirst(myscratch.pooledcovarmat.labels, myscratch.pooledunibeta.labels[myscratch.actualsel]))
	@simd for i = 1 : length(covar)
		@inbounds myscratch.actualnom[findfirst(myscratch.pooledunibeta.labels, myscratch.pooledcovarmat.labels[i])] -= myscratch.actualupdates[myscratch.actualstepno] * covar[i]
	end
	newselval = myscratch.actualnom[myscratch.actualsel]^2
	return newselval
end
"""
Function to repeat the boostingsteps (which were already performed) with new data resulting from a new data call.
Arguments:
myscratch::Boostscratch -> Container with data which are relevant to perform the boostinsteps
"""
function reboost!(myscratch::Boostscratch)
	for i = 1 : myscratch.actualstepno
		covar = getcovarfromtriangular(myscratch.pooledcovarmat.covarmat, findfirst(myscratch.pooledcovarmat.labels, myscratch.selections[i]))
		@simd for j = 1 : length(myscratch.wantedlabels)
			@inbounds myscratch.actualnom[findfirst(myscratch.pooledunibeta.labels, myscratch.wantedlabels[j])] -= myscratch.actualupdates[i] * covar[findfirst(myscratch.pooledcovarmat.labels, myscratch.wantedlabels[j])]
		end
	end
	myscratch.actualstepno += 1
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
const CONVERGENCE_MEMORY=500

mutable struct GMC_NS_Progress{T<:Real} <: AbstractProgress
    interval::T
    dt::AbstractFloat
    start_it::Integer
    counter::Integer
    triggered::Bool
    tfirst::AbstractFloat
    tlast::AbstractFloat
    tstp::AbstractFloat
    printed::Bool        # true if we have issued at least one status update
    desc::AbstractString # prefix to the percentage, e.g.  "Computing..."
    color::Symbol        # default to green
    output::IO           # output stream into which the progress is written
    numprintedvalues::Integer   # num values printed below progress in last iteration
    offset::Integer             # position offset of progress bar (default is 0)

    e::GMC_NS_Ensemble
    #tuner::GMC_Tuner
    top_m::GMC_NS_Model

    mean_stp_time::AbstractFloat

    upper_displays::Vector{Vector{Function}}
    lower_displays::Vector{Vector{Function}}
    display_rotate_iterates::Int64
    upperidx::Integer
    loweridx::Integer

    convergence_history::Vector{Float64}
    time_history::Vector{Float64}

    function GMC_NS_Progress{T}(e::GMC_NS_Ensemble,
                               top_m::GMC_NS_Model,
                               #tuner::GMC_Tuner,
                               interval::T;
                               dt::Real=0.1,
                               desc::AbstractString="Nested Sampling::",
                               color::Symbol=:green,
                               output::IO=stdout,
                               offset::Int=0,
                               start_it::Int=1,
                               upper_displays=[[evidence_display]],
                               lower_displays=[[ensemble_display]],
                               disp_rot_its=0,
                               ) where T
        tfirst = tlast = time()
        printed = false
        new{T}(interval,
         dt,
         start_it,
         start_it,
         false,
         tfirst,
         tlast,
         0.,
         printed,
         desc,
         color,
         output,
         0,
         offset,
         e,
         #tuner,
         top_m,
         0.,
         upper_displays,
         lower_displays,
         disp_rot_its,
         1,
         1,
         zeros(CONVERGENCE_MEMORY),
         zeros(CONVERGENCE_MEMORY))
    end
end

"""
    GMC_NS_Progress(e, interval; kwargs...)

Return a GMC_NS.jl ProgressMeter for the given ensemble, to be updated at 'interval'. Important keyword arguments:

'upper_displays <: AbstractVector{<:AbstractVector{Function}}': define a vector of above-status-line display functions to rotate through
'lower_displays <: AbstractVector{<:AbstractVector{Function}}': define a vector of below-status-line display functions to rotate through
'disp_rot_its <: Integer': define the number of iterates to rotate to the next upper and lower display function vectors

"""
function GMC_NS_Progress(e::GMC_NS_Ensemble,
    #tuner::GMC_Tuner,
     interval::Real; dt::Real=0.1, desc::AbstractString="GMC-NS::", color::Symbol=:green, output::IO=stderr, offset::Integer=0, start_it::Integer=1, upper_displays::AbstractVector{<:AbstractVector{Function}}=Vector{Vector{Function}}(), lower_displays::AbstractVector{<:AbstractVector{Function}}=Vector{Vector{Function}}(), disp_rot_its::Integer=0)
    top_m = deserialize(e.models[findmax([model.log_Li for model in e.models])[2]].path)

    return GMC_NS_Progress{typeof(interval)}(e, top_m,
    # tuner,
     interval, dt=dt, desc=desc, color=color, output=output, offset=offset, start_it=start_it, upper_displays=upper_displays, lower_displays=lower_displays, disp_rot_its=disp_rot_its)
end

function update!(p::GMC_NS_Progress, val, thresh; options...)
    p.counter += 1

    p.tstp=time()-p.tlast
    popfirst!(p.time_history)
    push!(p.time_history, p.tstp)

    p.interval = val - thresh
    isinf(p.interval) && (p.interval = 0.)
    popfirst!(p.convergence_history)
    push!(p.convergence_history, p.interval)

    string(p.top_m.trajectory,'.',p.top_m.i) != basename(p.e.models[findmax([model.log_Li for model in p.e.models])[2]].path) && (p.top_m=deserialize(p.e.models[findmax([model.log_Li for model in p.e.models])[2]].path))

    updateProgress!(p; options...)
end

function updateProgress!(p::GMC_NS_Progress; offset::Integer = p.offset, keep = (offset == 0))
    p.offset = offset
    t = time()
    p.display_rotate_iterates > 0 && p.counter%p.display_rotate_iterates == 0 && rotate_displays(p)

    if p.interval < 0 && !p.triggered
        p.triggered = true
        if p.printed
            p.triggered = true
            dur = durationstring(t-p.tfirst)
            msg = @sprintf "%s Converged. Time: %s (%d iterations). logZ: %s\n" p.desc dur p.counter p.e.log_Zi[end]

            print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
            move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
            length(p.upper_displays)>0 ? (upper_lines=display_dash(p, p.upper_displays[p.upperidx])) : (upper_lines=0)
            printover(p.output, msg, :magenta)
            length(p.lower_displays)>0 ? (lower_lines=display_dash(p, p.lower_displays[p.loweridx])) : (lower_lines=0)

            p.numprintedvalues=upper_lines + lower_lines + 1

            if keep
                println(p.output)
            else
                print(p.output, "\r\u1b[A" ^ (p.offset + p.numprintedvalues))
            end
        end
        return
    end

    if t > p.tlast+p.dt && !p.triggered
        p.counter < CONVERGENCE_MEMORY ? mean_step_time=mean(p.time_history[end-(p.counter-1):end]) :
                                            mean_step_time=mean(p.time_history)
        msg = @sprintf "%s Iterate: %s Recent step time μ: %s Convergence Interval: %g\n" p.desc p.counter hmss(mean_step_time) p.interval

        print(p.output, "\n" ^ (p.offset + p.numprintedvalues))
        move_cursor_up_while_clearing_lines(p.output, p.numprintedvalues)
        length(p.upper_displays)>0 ? (upper_lines=display_dash(p, p.upper_displays[p.upperidx])) : (upper_lines=0)
        printover(p.output, msg, p.color)
        length(p.lower_displays)>0 ? (lower_lines=display_dash(p, p.lower_displays[p.loweridx])) : (lower_lines=0)

        p.numprintedvalues=upper_lines + lower_lines + 1
        print(p.output, "\r\u1b[A" ^ (p.offset + p.numprintedvalues))

        # Compensate for any overhead of printing. This can be
        # especially important if you're running over a slow network
        # connection.
        p.tlast = t + 2*(time()-t)
        p.printed = true
    end
end
                function rotate_displays(p)
                    p.upperidx < length(p.upper_displays) ? p.upperidx+=1 : p.upperidx=1
                    p.loweridx < length(p.lower_displays) ? p.loweridx+=1 : p.loweridx=1
                end

                function hmss(dt)
                    dt<0 ? (dt=-dt; prfx="-") : (prfx="")
                    isnan(dt) && return "NaN"
                    (h,r) = divrem(dt,60*60)
                    (m,r) = divrem(r, 60)
                    (isnan(h)||isnan(m)||isnan(r)) && return "NaN"
                    string(prfx,Int(h),":",Int(m),":",Int(ceil(r)))
                end

                function display_dash(p::GMC_NS_Progress, funcvec::Vector{Function})
                    lines=0
                    for func in funcvec
                        lines+=func(p)
                    end
                    return lines
                end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   This file contains functions related to key bindings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Dictionary with the default bindings.
const _default_keybindings = Dict{Tuple{Union{Symbol, String}, Bool, Bool, Bool}, Symbol}(
    ("q",       false, false, false) => :quit,
    ("?",       false, false, false) => :help,
    (:up,       false, false, false) => :up,
    ("k",       false, false, false) => :up,
    (:down,     false, false, false) => :down,
    ("j",       false, false, false) => :down,
    (:left,     false, false, false) => :left,
    ("h",       false, false, false) => :left,
    (:right,    false, false, false) => :right,
    ("l",       false, false, false) => :right,
    (:up,       false, false, true ) => :fastup,
    (:down,     false, false, true ) => :fastdown,
    (:left,     false, false, true ) => :fastleft,
    (:right,    false, false, true ) => :fastright,
    (:left,     true,  false, false) => :bol,
    ("0",       false, false, false) => :bol,
    (:right,    true,  false, false) => :eol,
    ("\$",      false, false, false) => :eol,
    (:pageup,   false, false, false) => :pageup,
    ("b",       false, false, false) => :pageup,
    (:pagedown, false, false, false) => :pagedown,
    (" ",       false, false, false) => :pagedown,
    (:home,     false, false, false) => :home,
    (:up,       true,  false, false) => :home,
    ("g",       false, false, false) => :home,
    ("<",       false, false, false) => :home,
    (:end,      false, false, false) => :end,
    (:down,     true,  false, false) => :end,
    ("G",       false, false, false) => :end,
    (">",       false, false, false) => :end,
    ("/",       false, false, false) => :search,
    ("n",       false, false, false) => :next_match,
    ("N",       false, false, false) => :previous_match,
    (:esc,      false, false, false) => :quit_search,
    ("f",       false, false, false) => :change_freeze,
    (:eot,      false, false, false) => :quit_eot,
    ("r",       false, false, false) => :toggle_ruler,
    ("t",       false, false, false) => :change_title_rows,
)

# Dictionary with the current keybindings, it is initialized here with the
# default values to improve startup time.
const _keybindings = copy(_default_keybindings)

"""
    delete_keybinding(key::Union{Char, Symbol}; alt::Bool = false, ctrl::Bool = false, shift::Bool = false)

Delete the keybinding `key`. The modifiers keys can be selected using the
keywords `alt`, `ctrl`, and `shift`.

For more information about how specify `key` see [`set_keybinding`](@ref).

"""
function delete_keybinding(
    key::Union{Char, Symbol};
    alt::Bool = false,
    ctrl::Bool = false,
    shift::Bool = false
)
    dict_key = (key isa Char ? string(key) : key, alt, ctrl, shift)
    delete!(_keybindings, dict_key)
    return nothing
end

"""
    reset_keybindings()

Reset key bindings to the original ones.

"""
function reset_keybindings()
    empty!(_keybindings)
    merge!(_keybindings, _default_keybindings)

    # Key bindings that depends on the mode.
    if get(ENV, "PAGER_MODE", "default") == "vi"
        _keybindings[(:eot,     false, false, false)] = :halfpagedown
        _keybindings[(:shiftin, false, false, false)] = :halfpageup
    else
        _keybindings[(:eot, false, false, false)] = :quit_eot
    end

    return nothing
end

"""
    set_keybinding(key::Union{Char, Symbol}, action::Symbol; alt::Bool = false, ctrl::Bool = false, shift::Bool = false)

Set key binding `key` to the action `action`. The modifiers keys can be selected
using the keywords `alt`, `ctrl`, and `shift`.

`key` can be a `Char` or a `Symbol` indicating one of the following special
keys:

    :up, :down, :right, :left, :home, :end, :F1, :F2, :F3, :F4, :F5, :F6, :F7,
    :F8, :F9, :F10, :F11, :F12, :keypad_dot, :keypad_enter, :keypad_asterisk,
    :keypad_plus, :keypad_minus, :keypad_slash, :keypad_equal, :keypad_0,
    :keypad_1, :keypad_2, :keypad_3, :keypad_4, :keypad_5, :keypad_6, :keypad_7,
    :keypad_8, :keypad_9, :delete, :pageup, :pagedown, :tab,

`action` can be one of the following symbols:

    :quit, :help, :up, :down, :left, :right, :fastup, :fastdown, :fastleft,
    :fastright :bol, :eol, :pageup, :pagedown, :home, :end

"""
function set_keybinding(
        key::Union{Char, Symbol}, action::Symbol;
        alt::Bool = false,
        ctrl::Bool = false,
        shift::Bool = false)
    dict_key = (key isa Char ? string(key) : key, alt, ctrl, shift)
    _keybindings[dict_key] = action
    return nothing
end


text:
using DynamicPolynomials
@testset "DynamicPolynomials basic groebner" begin
    @polyvar x y

    @test Groebner.groebner([x, y]) == Groebner.groebner([y, x]) == [y, x]
    @test Groebner.groebner([5x, y]) == Groebner.groebner([y, 5x]) == [y, x]
    @test Groebner.groebner([x, y, y]) == Groebner.groebner([y, y, x]) == [y, x]
    @test Groebner.groebner([x + 3//2, y]) == [y, x + 3//2]
    @test Groebner.groebner([x + 3, 2y - 6]) == [y - 3, x + 3]

    fs = [x^2 - 5y]
    G = Groebner.groebner(fs)
    @test G == [x^2 - 5y]

    fs = [
        x + y^2,
        x*y - y^2
    ]
    G = Groebner.groebner(fs)
    @test Groebner.isgroebner(G)

end


text:
using Winston

p = FramedPlot(
        title="Title",
        xlabel="X axis",
        ylabel="Y axis")

add(p, Histogram(hist(randn(1000))...))
add(p, PlotLabel(.5, .5, "Histogram", color=0xcc0000))

t1 = Table(1, 2)
t1[1,1] = p
t1[1,2] = p

t2 = Table(2, 1)
t2[1,1] = t1
t2[2,1] = p

file(t2, "example3.png")


text:
# Calcul des distances a vol d'oiseau entre 2 coordonnees GPS
# X. Gandibleux - Octobre 2019

using Printf
using DelimitedFiles

# Donnees dans un fichier au format CSV séparé par des virgules avec :
# - ligne 1 : header
# - colonne 1 : nom du lieu
# - colonne 2 : position x du lieu
# - colonne 3 : position y du lieu
# - colonne 4 : latitude du lieu
# - colonne 5 : longitude du lieu
# donnees extraites de : https://data.education.gouv.fr
m = readdlm("GPSlyceesNantes.csv", ',')

# Fonction de conversion degre vers radian
function deg_rad(v)
    return pi/180*v
end

# Calcul du distancier
lieu = String[]
for i in 2:size(m,1)
    push!(lieu, m[i,1])
end

dist = zeros(size(m,1)-1, size(m,1)-1)
for i in 2:size(m,1)-1
    for j in i+1:size(m,1)

        LatitudeA  = m[i,4]
        LongitudeA = m[i,5]
        LatitudeB  = m[j,4]
        LongitudeB = m[j,5]

        # formule des sinus : http://villemin.gerard.free.fr/aGeograp/Distance.htm
        # https://www.coordonnees-gps.fr/conversion-coordonnees-gps

        d = acos(
             sin(deg_rad(LatitudeA)) * sin(deg_rad(LatitudeB))
             +
             cos(deg_rad(LatitudeA)) * cos(deg_rad(LatitudeB)) * cos(deg_rad(LongitudeB) - deg_rad(LongitudeA))
            ) * 6371

        @printf("%2d - %2d = %6.3f km | %s <-> %s \n",i-1, j-1, d, m[i,1], m[j,1])

        dist[i-1,j-1] = d
        dist[j-1,i-1] = d
    end
end

println("distancier symétrique crée : ", length(lieu), " lieux")

max_stars_repo_name: xgandibleux/distances

max_issues_repo_name: xgandibleux/distances

max_forks_repo_name: xgandibleux/distances


text:
using Test
using LinearAlgebra

using CompScienceMeshes
using BEAST

x = point(1.0, 0.0, 0.0)
y = point(0.0, 1.0, 0.0)
z = point(0.0, 0.0, 1.0)

κ = 1.0
S = SingleLayerTrace(im*κ)

h = 0.25
Γ = meshrectangle(1.0,1.0,h)
γ = meshsegment(1.0,1.0,3)

in_interior = CompScienceMeshes.interior_tpredicate(Γ)
on_junction = CompScienceMeshes.overlap_gpredicate(γ)
edges = skeleton(Γ,1) do edge
    in_interior(edge) ||on_junction(chart(Γ,edge))
end
X = raviartthomas(Γ, edges)

x = divergence(X)
y = ntrace(X,γ)
Z = assemble(S,y,x)

# test for the correct sparsity pattern
# I, J, V = findall(!iszero, Z)
Q = findall(!iszero, Z)
I = getindex.(Q,1)
J = getindex.(Q,2)
@test length(unique(I)) == 4
@test length(unique(J)) == 44

## test the acutal value of the penalty terms
using CompScienceMeshes
using BEAST
using Test

p1 = point(0,0,0)
p2 = point(1,0,0)
p3 = point(0,1,0)

q1 = point(0,0,0)
q2 = point(1,0,0)

m = Mesh([p1,p2,p3],[index(1,2,3)])
n = Mesh([q1,q2], [index(1,2)])
translate!(n, point(0,0,20))

X = lagrangecxd0(m)
Y = lagrangecxd0(n)

x = refspace(X)
y = refspace(Y)

N = BEAST.SingleLayerTrace(0.0)
Nyx = assemble(N,Y,X)

@test size(Nyx) == (1,1)

sx = chart(m, first(cells(m)))
sy = chart(n, first(cells(n)))

cx = neighborhood(sx, [1,1]/3)
cy = neighborhood(sy, [1]/2)

vx = x(cx)
vy = y(cy)

R = norm(cartesian(cx) - cartesian(cy))
estimate = volume(sx) * volume(sy) * vx[1][1] * vy[1][1] / (4π*R)
actual = Nyx[1,1]
@test norm(estimate - actual) / norm(actual) < 5e-4

# test that the trace and divergence work as advetised
X = raviartthomas(m,boundary(m))

D = divergence(X)
Y = ntrace(X,boundary(m))

@test numfunctions(D) == 3
@test isa(refspace(X), BEAST.RTRefSpace)
for _f in D.fns
    @test length(_f) == 1
    @test _f[1].cellid == 1
    @test _f[1].refid == 1
    @test _f[1].coeff ≈ (1 / volume(sx))
end

Σ = geometry(Y)
@test numfunctions(Y) == 3
@test isa(refspace(Y), BEAST.LagrangeRefSpace)
for _f in Y.fns
    @test length(_f) == 1
    @test 0 < _f[1].cellid < 4
    seg = chart(Σ, cells(Σ)[_f[1].cellid])
    @test _f[1].refid == 1
    @test _f[1].coeff ≈ (1 / volume(seg))
end


text:
@enum(Kind,
    ENDMARKER, # EOF
    ERROR,
    NAME, # foo, Σxx
    AT, # @
    PIPE, #|
    AMP, #&

    begin_keywords,
        KEYWORD, # general
    end_keywords,

    begin_literal,
        INT, # 4
        FLOAT, # 3.5, 3.7e+3
        STRING, # "foo"
        TRIPLE_STRING, # """ foo \n """
    end_literal,

    begin_delimiters,
        LSQUARE, # [
        RSQUARE, # ]
        LBRACE, # {
        RBRACE, # }
        LPAREN, # (
        RPAREN,  # )
    end_delimiters,

    begin_ops,
        OP,
        SPREAD, # ...
        EQUALS, # =
        COLON, # :
        DOLLAR, # $
        BANG, # !
    end_ops,
)


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using MelissasModule
using Plots

environment = Environment(5, 100)

function create_shoot_plot(
    gridlength = 10;
    weight_values::Tuple{Float64,Float64} = (100.0, 500.0),
    angle_values = (pi / 10, pi / 3),
    environment = environment,
    kwargs...,
)
    weights = range(weight_values[1], weight_values[2], length = gridlength)
    angles = range(angle_values[1], angle_values[2], length = gridlength)
    distances = [
        shoot_distance(Trebuchet(weight, angle), environment) for
        (angle, weight) in Iterators.product(angles, weights)
    ]
    surface(weights, angles, distances; kwargs...)
end

##
plt = create_shoot_plot(
    30;
    weight_values = (20.0, 1100.0),
    angle_values = (pi / 30, 2pi / 3),
    camera = (15, 45),
    yguidefontrotation = 45,
    xlabel = "Counterweight",
    ylabel = "α",
    zlabel = "Distance",
    dpi = 300
)
savefig(plt, joinpath(@__DIR__, "..", "..", "fig", "shoot_surface.png"))
##

max_stars_repo_name: tkphd/julia-novice

max_issues_repo_name: tkphd/julia-novice

max_forks_repo_name: tkphd/julia-novice


text:
module Metatheory

using RuntimeGeneratedFunctions
using Base.Meta
using Reexport

include("docstrings.jl")

RuntimeGeneratedFunctions.init(@__MODULE__)

include("options.jl")

macro log(args...)
    quote options.verbose && @info($(args...)) end |> esc
end

export options

include("Util/Util.jl")
using .Util
export Util

include("rgf.jl")
export @metatheory_init

include("Patterns/Patterns.jl")
@reexport using .Patterns

include("ematch_compiler.jl")
@reexport using .EMatchCompiler

include("Rules/Rules.jl")
@reexport using .Rules

include("match.jl")

include("EGraphs/EGraphs.jl")
@reexport using .EGraphs

include("Library/Library.jl")
export Library

using Rewriters
function rewrite(expr, theory; order=:outer)
   if order == :inner
      Fixpoint(Prewalk(Fixpoint(Chain(theory))))(expr)
   elseif order == :outer
      Fixpoint(Postwalk(Fixpoint(Chain(theory))))(expr)
   end
end
export rewrite

end # module


text:
export ancestralsampling!

"""
    ancestralsampling!(root::AbstractRegionGraphNode, x::AbstractMatrix{T}; obs = 1:size(x,1))

Sample z_n assignments for each observation using ancestral sampling and update all parameters of the network in-place.
"""
function ancestralsampling!(root::AbstractRegionGraphNode, x::AbstractMatrix{T}; obs = 1:size(x,1)) where {T<:Real}

    (N,D) = size(x)
    nodes = topoligicalorder(root)
    nodes = filter(n -> !(n isa PartitionGraphNode), nodes)
    out = AxisArray(Vector{Matrix}(undef, length(nodes)), map(id, nodes))

    # Bottom-up path.
    for n in nodes
        logpdf!(n, x, out)
    end

    # Top-down path.
    _ancestralsampleassignments!(root, out, [collect(obs)], x)

    return root
end

function _ancestralsampleassignments!(n::RegionGraphNode,
                                      out::AxisArray{<:AbstractArray},
                                      obs::Vector{<:AbstractVector},
                                      x::AbstractMatrix)
    @assert length(obs) == length(n)

    (N, D) = size(x)
    couts = map(c -> _getchildlogpdf(c, out), children(n))
    Ch = size(n.logweights,1)
    T = eltype(out[n.id])
    clusters = Vector{Vector{Int}}(undef, Ch)
    for c in 1:length(clusters)
        clusters[c] = Vector{Int}()
    end

    out_ = _getchildlogpdf.(children(n), Ref(out))

    for k in 1:length(n)
        z = zeros(Int, length(obs[k]))
        Threads.@threads for i in 1:length(obs[k])
            @inbounds begin
                local ni = obs[k][i]
                local p_ = zeros(T, Ch)
                local j = 1
                for ccnt in 1:length(children(n))
                    local out_ = couts[ccnt]
                    for j_ in 1:size(out_,2)
                        p_[j] = hasscope(n) ? exp(out_[ni,j_] + n.logweights[j,k]) : exp(n.logweights[j,k])
                        j += 1
                    end
                end
                z[i] = randomindex(p_)
            end
        end

        for c in 1:length(clusters)
            obs_ = findall(z .== c)
            @assert obs_ isa Vector
            clusters[c] = obs[k][obs_]
        end

        @inbounds for j in 1:size(n.active,1)
            n.active[j,k] = j ∈ z
        end

        # Update weights.
        sstat = DirichletSufficientStats(map(j -> sum(z .== j), 1:size(n.logweights, 1)))
        w = rand(Dirichlet(postparams(n.prior, sstat)))
        w /= sum(w)
        @inbounds n.logweights[:,k] .= log.(w)
    end

    j = 1
    for child in children(n)
        obs_ = clusters[j:(j+length(child)-1)]
        @assert all(g -> all((i) -> (0 < i <= N), g), obs_) "$obs_"
        _ancestralsampleassignments!(child, out, obs_, x)
        j += length(child)
    end
end

function _ancestralsampleassignments!(n::PartitionGraphNode,
                                      out::AxisArray{<:AbstractArray},
                                      obs::Vector{<:AbstractVector},
                                      x::AbstractMatrix)

    (N,D) = size(x)
    clusters = Vector{Vector{Vector{Int}}}(undef, length(children(n)))
    for (c, child) in enumerate(children(n))
        clusters[c] = map(k -> Vector{Int}(), 1:length(child))
    end

    dims = tuple(map(c -> length(c), children(n))...)

    for k in 1:length(n)
        if !isempty(obs[k])
            @inbounds begin
                idxs = Base._ind2sub(dims, k)
                for (c, child) in enumerate(children(n))
                    append!(clusters[c][idxs[c]], obs[k])
                end
            end
        end
    end

    for (c, child) in enumerate(children(n))
        @assert all(g -> all((i) -> (0 < i <= N), g), clusters[c])
        _ancestralsampleassignments!(child, out, clusters[c], x)
    end
end

function _ancestralsampleassignments!(n::FactorizedDistributionGraphNode,
                                      out::AxisArray{<:AbstractArray},
                                      obs::Vector{<:AbstractVector},
                                      x::AbstractMatrix)
    (D,K) = size(n.sstats)
    fill!(n.obsVecs, false)

    #Threads.@threads
    for d in 1:D
        for k in 1:length(n)
            if d == 1
                n.obsVecs[obs[k],k] .= true
            end
            x_ = x[obs[k],d]
            x_ = x_[.!isnan.(x_)]
            update!(n.sstats[d,k], x_)
        end
    end
end

function _ancestralsampleassignments!(n::FactorizedMixtureGraphNode,
                                      out::AxisArray{<:AbstractArray},
                                      obs::Vector{<:AbstractVector},
                                      x::AbstractMatrix)
    T = eltype(out[n.id])
    D = length(n.sstats)
    fill!(n.obsVecs, 0)

    for d in 1:D
        for k in 1:length(n)
            z = zeros(Int, length(obs[k]))
            for i in 1:length(obs[k])
                local ni = obs[k][i]
                @assert ni > 0 "$(obs), k: $k"
                local p_ = zeros(T, size(n.logweights[d],1))
                for m in 1:size(n.logweights[d],1)
                    θ = postparams(n.priors[d][m], n.sstats[d][m,k])
                    p_[m] = d ∈ scope(n) ? exp(_logpostpdf(n.likelihoods[d][m], vec(view(x, ni, d)), θ)[1]) : 1.0
                    p_[m] *= exp(n.logweights[d][m,k])
                end
                z[i] = randomindex(p_)
            end

            @inbounds n.obsVecs[obs[k],k] .= z

            for m in 1:size(n.sstats[d],1)
                obs_ = filter(i -> !isnan(x[i,d]), obs[k][z .== m])
                @inbounds update!(n.sstats[d][m,k], x[obs_,d])
            end

            sstat = DirichletSufficientStats(map(j -> sum(z .== j), 1:size(n.logweights[d], 1)))
            w = postparams(n.weightpriors[d], sstat)
            w /= sum(w)
            @inbounds n.logweights[d][:,k] .= log.(w)
        end
    end
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
"""
    restriction(submesh, supermesh)

Computes the restriction matrix relative to a submesh `submesh` of `supermesh`.

The restriction matrix has size `(m,n)`, where

    m == numcells(submesh)
    n == numcells(supermesh)

It has entries `1` at location `[i,j]` iff cell `i` of submesh equals cell `j` of supermesh.
The remaining entries are zero.

This matrix is referred to as the restriction matrix because if it acts on an array of
samples taken at the cells of `supermesh` is selects out the samples in the cells that
are retained in `submesh`, taking into account any renumbering. Likewise, its transpose
is sometimes referred to as the extension-by-zero operator because it maps arrays of samples
taken in the cells of `submesh` into an array of samples taken in the cells of `supermesh`
by inserting zeros at cells that were not retained in `submesh`.
"""
function restriction(submesh::Mesh, supermesh::Mesh)
    index_in_supermesh = mapper(supermesh)
    R = spzeros(Int, numcells(submesh), numcells(supermesh))
    for (i,c) in enumerate(cells(submesh))
        j = index_in_supermesh[c]
        R[i,j] = 1
    end
    R
end


text:
"""
    InverseFunctions.test_inverse(f, x; compare=isapprox, kwargs...)

Test if [`inverse(f)`](@ref) is implemented correctly.

The function tests (as a `Test.@testset`) if

* `compare(inverse(f)(f(x)), x) == true` and
* `compare(inverse(inverse(f))(x), f(x)) == true`.

`kwargs...` are forwarded to `compare`.
"""
function test_inverse(f, x; compare=isapprox, kwargs...)
    @testset "test_inverse: $f with input $x" begin
        y = f(x)
        inverse_f = inverse(f)
        @test compare(inverse_f(y), x; kwargs...)
        inverse_inverse_f = inverse(inverse_f)
        @test compare(inverse_inverse_f(x), y; kwargs...)
    end
    return nothing
end


text:
using SeisIO, PlotlyJS, JLD2, Dates, Noise, Printf, ORCA

fi = "./dataset/BPnetwork.jld2"

t = jldopen(fi)

DLtimestamplist = t["info/DLtimestamplist"]
stationlist = t["info/stationlist"]

tid  = 1
for tid = 1:length(DLtimestamplist)
    plotspan=1000
    # make beampower timestamp: bptimestamp
    t1 = parse.(Int64, split(DLtimestamplist[tid], ".")[1:2])
    hm = split(DLtimestamplist[1], ".")[3]
    m1, d1 = SeisIO.j2md(t1[1], t1[2])
    st = DateTime(join([t1[1], "-", m1, "-",  d1, hm]))
    tstamp = Noise.timestamp(st)

    layout = Layout(width=800, height=600,
                    xaxis=attr(title="Time [hours]", tickvals =[0, 6, 12, 18, 24]),
                    yaxis=attr(tickvals = collect(1:length(stationlist)),
                    ticktext=stationlist, automargin=true),
                    #font =attr(size=12),
                    showlegend=false,
                    xaxis_range=[0, 24],
                    #yaxis_range=[],
                    title = @sprintf("Time:%s", tstamp)
                    )

    p = PlotlyJS.plot([NaN], layout)

    amp = 1e4
    #plot trace
    for i = 1:length(stationlist)
        station = stationlist[i]
        S = t[joinpath(DLtimestamplist[tid], station)]
        tvec = collect(0:length(S.x)-1)/S.fs / 60 / 60
        trace1 = scatter(;x=tvec[1:plotspan:end], y=amp*S.x[1:plotspan:end].+float(i), mode="lines")
        PlotlyJS.addtraces!(p, trace1)
    end

    PlotlyJS.deletetraces(p,1)

    #display(p)

    mkpath("./fig_dawdata")
    figname = @sprintf("./fig_dawdata/%s.png",tstamp)
    savefig(p, figname)
end

max_stars_repo_name: kura-okubo/SeisBeamforming

max_issues_repo_name: kura-okubo/SeisBeamforming

max_forks_repo_name: kura-okubo/SeisBeamforming


text:
using HorizonSideRobots
include("../basics.jl")

md"""
    ДАНО: Робот - в произвольной клетке поля (без внутренних перегородок и маркеров)

    РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы
"""

"""Двигает робота на север до границы, возвращая кол-во шагов, потребовавшихся ему для этого"""
function move_to_border!(r)
    cnt = 0
    while isborder(r, Nord) == 0
        move!(r, Nord)
        cnt += 1
    end
    return cnt
end

"""
Двигает робота вдоль границы.
Если оказывается в углу: true. Если встретил промаркированную клетку (вернулся в начало, пройдя по периметру): false.
"""
function go_and_mark_border!(r, side)
    while isborder(r, side) == 0
        if ismarker(r)
            return false
        end
        putmarker!(r)
        move!(r, side)
    end
    return true
end

 """Решение задачи 2"""
 function solve_problem_2!(r)
    cnt = move_to_border!(r)
    i = 0
    while go_and_mark_border!(r, HorizonSide(i % 4))
        i+=1
    end
    for i in 1:cnt
        move!(r, Sud)
    end
    show(r)
end


text:
function affinegap_local_align{T}(a, b, submat::AbstractSubstitutionMatrix{T}, gap_open_penalty::T, gap_extend_penalty::T)
    m = length(a)
    n = length(b)
    ge = gap_extend_penalty
    goe = gap_open_penalty + ge
    trace = Matrix{Trace}(m + 1, n + 1)
    H = Vector{T}(m + 1)
    E = Vector{T}(m)
    # run dynamic programming column by column
    @inbounds begin
        H[1] = T(0)
        trace[1,1] = TRACE_NONE
        for i in 1:m
            H[i+1] = T(0)
            trace[i+1,1] = TRACE_NONE
        end
        best_score = typemin(T)
        best_endpos = (0, 0)
        for j in 1:n
            h_diag = H[1]
            H[1] = T(0)
            trace[1,j+1] = TRACE_NONE
            # any value goes well since this will be set in the first iteration
            F = T(0)
            for i in 1:m
                # gap in the sequence A
                e = H[i+1] - goe
                gap_a = TRACE_GAPOPEN_A
                if j > 1
                    e′ = E[i] - ge
                    if e′ == e
                        gap_a |= TRACE_GAPEXTD_A
                    elseif e′ > e
                        gap_a  = TRACE_GAPEXTD_A
                        e = e′
                    end
                end
                # gap in the sequence B
                f = H[i] - goe
                gap_b = TRACE_GAPOPEN_B
                if i > 1
                    f′ = F - ge
                    if f′ == f
                        gap_b |= TRACE_GAPEXTD_B
                    elseif f′ > f
                        gap_b  = TRACE_GAPEXTD_B
                        f = f′
                    end
                end
                # match
                h = h_diag + submat[a[i],b[j]]
                # find the best score and its trace
                best = max(T(0), e, f, h)
                t = TRACE_NONE
                e == best && (t |= gap_a)
                f == best && (t |= gap_b)
                h == best && (t |= TRACE_MATCH)
                # update
                E[i] = e
                F = f
                h_diag = H[i+1]
                H[i+1] = best
                trace[i+1,j+1] = t
                if best > best_score
                    best_score = best
                    best_endpos = (i, j)
                end
            end
        end
    end
    return best_score, trace, best_endpos
end

function affine_local_traceback(a, b, trace, endpos)
    # gap/character counts (reversed order)
    counts_a = [0, 0]
    counts_b = [0, 0]
    # if gap extension is selected in the previous traceback step, either
    # gap extension or gap open in that direction should be selected.
    gapext_a = false
    gapext_b = false
    i, j = endpos
    while trace[i+1,j+1] != TRACE_NONE
        @assert !(gapext_a && gapext_b)
        t = trace[i+1,j+1]
        if gapext_a
            if j ≥ 2 && t & TRACE_GAPEXTD_A > 0
                @gapext a
            elseif t & TRACE_GAPOPEN_A > 0
                @gapopen a
            end
        elseif gapext_b
            if i ≥ 2 && t & TRACE_GAPEXTD_B > 0
                @gapext b
            elseif t & TRACE_GAPOPEN_B > 0
                @gapopen b
            end
        elseif t & TRACE_MATCH > 0
            @match
        elseif j ≥ 2 && t & TRACE_GAPEXTD_A > 0
            @gapext a
        elseif t & TRACE_GAPOPEN_A > 0
            @gapopen a
        elseif i ≥ 2 && t & TRACE_GAPEXTD_B > 0
            @gapext b
        elseif t & TRACE_GAPOPEN_B > 0
            @gapopen b
        end
        # do not come here
        @assert false
    end
    # update counts
    reverse!(counts_a)
    reverse!(counts_b)
    return GappedSequence(a, i + 1, counts_a), GappedSequence(b, j + 1, counts_b)
end


text:
function tcm_predictMMNL(N, K, class_probs, u, assortments)
	M = size(assortments,1)
	v = zeros(size(assortments))
	v_by_class = Array{ Array{Float64,2}, 1}(K)
	for k in 1:K
		v_one_class = zeros(Float64,size(assortments))
	    for m = 1 : M
	        v_one_class[m,N] = 1.0
	        for p = 1 : N-1
	            v_one_class[m,p] = exp(u[k,p]) * assortments[m,p]
	        end
	        v_one_class[m,:] = v_one_class[m,:] / sum(v_one_class[m,:])
	    end
		v += class_probs[k] * v_one_class

		v_by_class[k] = v_one_class
	end

	return v, v_by_class
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using FractionalDiffEq, Plots

function classical!(du, u, p, t)
    A, B = 1, 3
    du[1]=A-(B+1)*u[1]+u[1]^2*u[2]
    du[2]=B*u[1]-u[1]^2*u[2]
end
α = [0.8, 0.7]
u0 = [1.2, 2.8]
h=0.01
prob = FODESystem(classical!, α, u0)
result = solve(prob, h, 10, NonLinearAlg())
tspan = collect(0:h:10)
plot(result[:, 1], result[:, 2])


text:
immutable Rayleigh <: ContinuousUnivariateDistribution
    σ::Float64

    function Rayleigh(σ::Real)
        σ > zero(σ) || error("Rayleigh: σ must be positive")
        @compat new(Float64(σ))
    end

    Rayleigh() = new(1.0)
end

@distr_support Rayleigh 0.0 Inf

#### Parameters

scale(d::Rayleigh) = d.σ
params(d::Rayleigh) = (d.σ,)

#### Statistics

mean(d::Rayleigh) = sqrthalfπ * d.σ
median(d::Rayleigh) = 1.177410022515474691 * d.σ   # sqrt(log(4.0)) = 1.177410022515474691
mode(d::Rayleigh) = d.σ

var(d::Rayleigh) = 0.429203673205103381 * d.σ^2   # (2.0 - π / 2) = 0.429203673205103381
std(d::Rayleigh) = 0.655136377562033553 * d.σ

skewness(d::Rayleigh) = 0.631110657818937138
kurtosis(d::Rayleigh) = 0.245089300687638063

entropy(d::Rayleigh) = 0.942034242170793776 + log(d.σ)

#### Evaluation

function pdf(d::Rayleigh, x::Float64)
	σ2 = d.σ^2
	x > 0.0 ? (x / σ2) * exp(- (x^2) / (2.0 * σ2)) : 0.0
end

function logpdf(d::Rayleigh, x::Float64)
	σ2 = d.σ^2
	x > 0.0 ? log(x / σ2) - (x^2) / (2.0 * σ2) : -Inf
end

logccdf(d::Rayleigh, x::Float64) = - (x^2) / (2.0 * d.σ^2)
ccdf(d::Rayleigh, x::Float64) = exp(logccdf(d, x))

cdf(d::Rayleigh, x::Float64) = 1.0 - ccdf(d, x)
logcdf(d::Rayleigh, x::Float64) = log1mexp(logccdf(d, x))

quantile(d::Rayleigh, p::Float64) = sqrt(-2.0 * d.σ^2 * log1p(-p))

#### Sampling

rand(d::Rayleigh) = d.σ * sqrt(2.0 * randexp())


text:
"""
    compute_all_culverts(path::String)

Parameters:

**`path`**: The file path to an INI file containing run parameters
"""
function compute_all_culverts(path::String)
    ## Parse INI Arguments
    cfg = Omniscape.parse_cfg(path)

    culvert_id_field = Symbol(cfg["culvert_id_field"])
    culvert_shapefile_path = cfg["culvert_shapefile_path"]
    culvert_ids = readdlm(cfg["culvert_ids_path"], ',', String)
    landcover_path = cfg["landcover_path"]
    os_radius_pixels = parse(Int, cfg["os_radius_pixels"])
    summary_radius_pixels = parse(Int, cfg["summary_radius_pixels"])
    os_source_threshold = parse(Float64, cfg["os_source_threshold"])
    reclass_table_path = cfg["reclass_table_path"]
    source_strength_path = cfg["source_strength_path"]
    os_block_size = parse(Int, cfg["os_block_size"])
    mask_values = float.(eval(Meta.parse(cfg["mask_values"])))
    output_folder = cfg["output_folder"]

    ## Side length derivation - assumes 0.5m resolution (hard coded)
    input_resolution = 0.5
    side_length_meters = (os_radius_pixels + summary_radius_pixels) * input_resolution * 2

    ## Setup output file structure and write empty file to store metadata
    # Create directory
    dir_suffix = 2
    while isdir(string(output_folder, "_$(dir_suffix)"))
        dir_suffix+=1
    end
    if isdir(output_folder)
        @warn string("A folder already exists at $(output_folder).",
                     " Creating a new folder",
                     ", $(string(output_folder, "_$(dir_suffix)")), ",
                     "and saving outputs there.")
        output_folder = string(output_folder, "_$(dir_suffix)")
    end
    mkpath(output_folder)

    ## Check that all ids in culvert_ids are in the shapefile
    culvert_shapefile_ids = string.(getproperty(Shapefile.Table(culvert_shapefile_path), culvert_id_field))
    bad_ids = culvert_ids[(!).(map(x -> in(x, culvert_shapefile_ids), culvert_ids))]
    if !isempty(bad_ids)
        @warn("The following user-provided culvert IDs are missing from the culverts shapefile and will be skipped. $(bad_ids)")

        missing_str_file = string(output_folder, "/missing_structures.txt")
        touch(missing_str_file)

        open(missing_str_file, "a") do f
          writedlm(f, bad_ids, ",")
        end

        culvert_ids = culvert_ids[map(x -> in(x, culvert_shapefile_ids), culvert_ids)]
    end

    # Create metadata CSV
    results_file = string(output_folder, "/results_summary.csv")
    touch(results_file)

    # write column headers to CSV
    summary_radius_str = string(Int(round(summary_radius_pixels * input_resolution)), "m")

    column_names = hcat([String(culvert_id_field)],
                        ["x_utm"],
                        ["y_utm"],
                        [string("mean_current_within_", summary_radius_str)],
                        [string("stdev_current_within_", summary_radius_str)],
                        [string("min_current_within_", summary_radius_str)],
                        [string("max_current_within_", summary_radius_str)])

    open(results_file, "w") do io
        writedlm(io, column_names, ",")
    end

    ## Build Dict for Omnsicape ini params
    os_args = Dict{String, String}()
    os_args["radius"] = string(os_radius_pixels)
    os_args["block_size"] = string(os_block_size)
    os_args["source_threshold"] = string(os_source_threshold)
    os_args["source_from_resistance"] = "false"
    os_args["reclassify_resistance"] = "true"
    os_args["solver"] = cfg["solver"]
    reclass_table = Omniscape.read_reclass_table(reclass_table_path, Float64)

    total_solves = length(culvert_ids)
    for i in 1:total_solves
        @info "Solving culvert $(culvert_ids[i]), $(i) of $(total_solves)"
        try
            os_args["project_name"] = string(output_folder, "/", culvert_ids[i])
            local output = compute_connectivity(culvert_shapefile_path,
                                                culvert_id_field,
                                                String(culvert_ids[i]),
                                                landcover_path,
                                                side_length_meters,
                                                source_strength_path,
                                                mask_values,
                                                os_args,
                                                reclass_table)

            cum_current, geotransform, x_coord, y_coord = output
            # Clip current map for summary, should always be an odd number of rows and
            # columns based on internal code
            culvert_row_y = Int(ceil((geotransform[4] - y_coord) / input_resolution)) + 1
            culvert_col_x = Int(floor((x_coord - geotransform[1]) / input_resolution)) + 1
            current_clipped = Omniscape.clip(cum_current, x = culvert_col_x,
                                             y = culvert_row_y,
                                             distance = summary_radius_pixels)

            current_mean = mean(skipmissing(current_clipped))
            current_sd = sqrt(var(skipmissing(current_clipped)))
            current_min = minimum(skipmissing(current_clipped))
            current_max = maximum(skipmissing(current_clipped))

            row_values = hcat([String(culvert_ids[i])],
                              [x_coord],
                              [y_coord],
                              [current_mean],
                              [current_sd],
                              [current_min],
                              [current_max])

            open(results_file, "a") do io
              writedlm(io, row_values, ",")
            end
        catch e
            error_message = sprint(showerror, e, catch_backtrace())
            @warn "Solve failed for culvert with ID $(culvert_ids[i]). Skipping."

            # Keep track of failed structures and their error messages
            if !isfile(string(output_folder, "/failed_structures.csv"))
                touch(string(output_folder, "/failed_structures.csv"))
                open(string(output_folder, "/failed_structures.csv"), "a") do f
                  writedlm(f, ["Culvert ID" "Error Message"], ",")
                end
            end

            open(string(output_folder, "/failed_structures.csv"), "a") do f
              writedlm(f, [culvert_ids[i] String(error_message)], ",")
            end
        end
        println(" ")
    end

    println("Done! Find all the output in $(string(pwd(), "/", output_folder))")

    nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewaySslCertificatePropertiesFormat <: SwaggerModel
    data::Any # spec type: Union{ Nothing, String } # spec name: data
    password::Any # spec type: Union{ Nothing, String } # spec name: password
    publicCertData::Any # spec type: Union{ Nothing, String } # spec name: publicCertData
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ApplicationGatewaySslCertificatePropertiesFormat(;data=nothing, password=nothing, publicCertData=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewaySslCertificatePropertiesFormat, Symbol("data"), data)
        setfield!(o, Symbol("data"), data)
        validate_property(ApplicationGatewaySslCertificatePropertiesFormat, Symbol("password"), password)
        setfield!(o, Symbol("password"), password)
        validate_property(ApplicationGatewaySslCertificatePropertiesFormat, Symbol("publicCertData"), publicCertData)
        setfield!(o, Symbol("publicCertData"), publicCertData)
        validate_property(ApplicationGatewaySslCertificatePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewaySslCertificatePropertiesFormat

const _property_map_ApplicationGatewaySslCertificatePropertiesFormat = Dict{Symbol,Symbol}(Symbol("data")=>Symbol("data"), Symbol("password")=>Symbol("password"), Symbol("publicCertData")=>Symbol("publicCertData"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewaySslCertificatePropertiesFormat = Dict{Symbol,String}(Symbol("data")=>"String", Symbol("password")=>"String", Symbol("publicCertData")=>"String", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ApplicationGatewaySslCertificatePropertiesFormat }) = collect(keys(_property_map_ApplicationGatewaySslCertificatePropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewaySslCertificatePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewaySslCertificatePropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewaySslCertificatePropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewaySslCertificatePropertiesFormat[property_name]

function check_required(o::ApplicationGatewaySslCertificatePropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewaySslCertificatePropertiesFormat }, name::Symbol, val)
end


text:
for (op, fun) in ((:maximum, :max), (:minimum, :min), (:prod, :prod), (:all, :all))
    isdefined(Base, op) || continue
    _op = Symbol("_", op)
    __op = Symbol("__", op)

    @eval begin
        Base.$op(x::AbstractJaxArray; dims=:, keepdims=_autokeep(dims)) = $_op(nothing, x, keepdims, dims)
        Base.$op(T::Type, x::AbstractJaxArray; dims=:, keepdims=_autokeep(dims)) = $_op(T, x, keepdims, dims)

        $_op(T, x, keepdims, dims...) = $__op(T, x, keepdims, _convert_dims(size(x), dims...))
        $__op(T, x, keepdims, dims) = numpy.$fun(x, axis=dims, dtype=T, keepdims=keepdims)
    end
end


text:
function chemical_equilibrium(A, B, C, ::Val{:JuMP})
    model = JuMP.Model(CBLS.Optimizer)

    n = length(C)
    m = length(B)

    # Add the number of moles per compound (continuous interval)
    @variable(model, 0 ≤ X[1:n] ≤ maximum(B))

    # mass_conservation function
    conserve = i -> (x ->
        begin
            δ = abs(sum(A[:, i] .* x) - B[i])
            return δ ≤ 1.e-6 ? 0. : δ
        end
    )

    for i in 1:m
        @constraint(model, X in Error(conserve(i)))
    end

    # computes the total energy freed by the reaction
    free_energy = x -> sum(j -> x[j] * (C[j] + log(x[j] / sum(x))))

    @objective(model, Min, ScalarFunction(free_energy))

    return model, X
end

"""
    chemical_equilibrium(atoms_compounds, elements_weights, standard_free_energy; modeler = :JuMP)

!!! warning

    Even the structure to model problems with continuous domains is available, the default solver is not yet equiped to solve such problems efficiently.

# From Wikipedia

In a chemical reaction, chemical equilibrium is the state in which both the reactants and products are present in concentrations which have no further tendency to change with time, so that there is no observable change in the properties of the system. This state results when the forward reaction proceeds at the same rate as the reverse reaction. The reaction rates of the forward and backward reactions are generally not zero, but they are equal. Thus, there are no net changes in the concentrations of the reactants and products. Such a state is known as dynamic equilibrium.
"""
function chemical_equilibrium(
    atoms_compounds,
    elements_weights,
    standard_free_energy;
    modeler = :JuMP,
)
    return chemical_equilibrium(
        atoms_compounds,
        elements_weights,
        standard_free_energy,
        Val(modeler)
    )
end


text:
JSON.lower(f::Filter) = non_nothing_dict(f)
JSON.lower(sqs::SearchQuerySpec) = non_nothing_dict(sqs)

"""
    Selector(dimension::String, value::String)

dimension = value filter.
"""
struct Selector <: Filter
    type::String
    dimension::String
    value::String
    Selector(dimension, value) = new("selector", dimension, value)
end

"""
    ColumnComparison(dimensions::Tuple{DimensionSpec, DimensionSpec})

    ColumnComparison(dim1::DimensionSpec, dim2::DimensionSpec)

dimension1 = dimension2 filter.
"""
struct ColumnComparison <: Filter
    type::String
    dimensions::Tuple{DimensionSpec, DimensionSpec}
    ColumnComparison(dimensions) = new("columnComparison", dimensions)
end

ColumnComparison(dim1, dim2) = ColumnComparison((dim1, dim2))

"""
    RegexF(dimension::String, pattern::String)

Regex pattern matching filter.
"""
struct RegexF <: Filter
    type::String
    dimension::String
    pattern::String
    RegexF(dimension, pattern) = new("regex", dimension, pattern)
end

"""
    AndF(filters::Vector{<:Filter})

Match all `filters`.
"""
struct AndF <: Filter
    type::String
    fields::Vector{<:Filter}
    AndF(filters) = new("and", filters)
end

"""
    OrF(filters::Vector{<:Filter})

Match at least one of the `filters`.
"""
struct OrF <: Filter
    type::String
    fields::Vector{<:Filter}
    OrF(filters) = new("or", filters)
end

"""
    NotF(filter::Filter)

Do not match `filter`.
"""
struct NotF <: Filter
    type::String
    field::Filter
    NotF(filter) = new("not", filter)
end

"""
    JavaScriptF(dimension::String, jsfunction::String)

Filter using a custom JavaScript function.
"""
struct JavaScriptF <: Filter
    type::String
    dimension::String
    jsfunction::String
    JavaScriptF(dimension, jsfunction) = new("javascript", dimension, jsfunction)
end

JSON.lower(f::JavaScriptF) = Dict("type" => f.type, "dimension" => f.dimension, "function" => f.jsfunction)

"""
    ExtractionFilter(dimension::String, value, extractionFn::ExtractionFunction)

Filter to `value` extracted using `extractionFn`.
"""
struct ExtractionFilter <: Filter
    type::String
    dimension::String
    value
    extractionFn::ExtractionFunction
    ExtractionFilter(dimension, value, extractionFn) = new("extraction", dimension, value, extractionFn)
end

"""
    Contains(value::String; case_sensitive=nothing)

Search query spec for strings containing `value`.

case_sensitive should be a Bool if provided.
"""
struct Contains <: SearchQuerySpec
    type::String
    value::String
    caseSensitive
    function Contains(value; case_sensitive=nothing)
        nothing_or_type(case_sensitive, Bool)
        new("contains", value, case_sensitive)
    end
end

"""
    InsensitiveContains(value::String)

Case insensitive string contains searhc query spec.
"""
struct InsensitiveContains <: SearchQuerySpec
    type::String
    value::String
    InsensitiveContains(value) = new("insensitive_contains", value)
end

"""
    Fragment(value::String; case_sensitive=nothing)

Fragment search query spec.

case_sensitive should be a Bool if provided.
"""
struct Fragment <: SearchQuerySpec
    type::String
    value::String
    caseSensitive
    function Fragment(value; case_sensitive=nothing)
        nothing_or_type(case_sensitive, Bool)
        new("fragment", value, case_sensitive)
    end
end

"""
    SearchF(dimension::String, query::SearchQuerySpec; extractionFn=nothing)

Filter using a query and optionally extraction function.

extractionFn should be an ExtractionFunction if provided.
"""
struct SearchF <: Filter
    type::String
    dimension::String
    query::SearchQuerySpec
    extractionFn
    function SearchF(dimension, query; extractionFn=nothing)
        nothing_or_type(extractionFn, ExtractionFunction)
        new("search", dimension, query, extractionFn)
    end
end

"""
    InF(dimension::String, values::Vector)

Filter if value is in `values`.
"""
struct InF <: Filter
    type::String
    dimension::String
    values::Vector
    InF(dimension, values) = new("in", dimension, values)
end

"""
    Like(dimension::String, pattern::String; escape=nothing, extractionFn=nothing)

Filter strings matching `pattern`.

escape should be a String and extractionFn an ExtractionFunction if provided.
"""
struct Like <: Filter
    type::String
    dimension::String
    pattern::String
    escape
    extractionFn
    function Like(dimension, pattern; escape=nothing, extractionFn=nothing)
        nothing_or_type(escape, AbstractString)
        nothing_or_type(extractionFn, ExtractionFunction)
        new("like", dimension, pattern, escape, extractionFn)
    end
end

"""
    Bound(dimension::String;
        lower=nothing, upper=nothing, lowerStrict=nothing, upperStrict=nothing,
        ordering=nothing, extractionFn=nothing)

Filter values bounded in a range.

lower, upper and ordering should be `String`s if provided. lowerStrict and
upperStrict should be `Bool`s if provided. extractionFn should be an
ExtractionFunction if provided.
"""
struct Bound <: Filter
    type::String
    dimension::String
    lower
    upper
    lowerStrict
    upperStrict
    ordering
    extractionFn
    function Bound(dimension; lower=nothing, upper=nothing, lowerStrict=nothing, upperStrict=nothing, ordering=nothing, extractionFn=nothing)
        nothing_or_type(lower, String)
        nothing_or_type(upper, String)
        lower === upper === nothing && error("At least one of lower and upper must be specified")
        nothing_or_type(lowerStrict, Bool)
        nothing_or_type(upperStrict, Bool)
        nothing_or_type(ordering, String)
        nothing_or_type(extractionFn, ExtractionFunction)
        new("bound", dimension, lower, upper, lowerStrict, upperStrict, ordering, extractionFn)
    end
end

"""
    IntervalF(dimension::String, intervals::Vector{<:Interval}; extractionFn=nothing)

Filter by time falling in one of the `intervals`.

extractionFn should be an ExtractionFunction if provided.
"""
struct IntervalF <: Filter
    type::String
    dimension::String
    intervals::Vector{<:Interval}
    extractionFn
    function IntervalF(dimension, intervals; extractionFn=nothing)
        nothing_or_type(extractionFn, ExtractionFunction)
        new("interval", dimension, intervals, extractionFn)
    end
end

"""
    TrueF()

Filter everything.
"""
struct TrueF <: Filter
    type::String
    TrueF() = new("true")
end


text:
using PyCall

@pyimport keras.models as models
@pyimport yaml as pyyaml

abstract Policy
# Policy must have:
# choose_move(board::Board, policy::Policy)

type RandomPolicy <: Policy end

function choose_move(board::Board, policy::RandomPolicy)
    candidates = legal_moves(board)
    rand(candidates)
end

###
# Keras Policy
# Training code in keras.jl
###

type KerasNetwork <: Policy
    # Ick - is there some way of improving typing here?
    model::Module
    features::Vector{Function}  # Feature extractors to run
    # Make it callable if given a raw object
    KerasNetwork(model::PyCall.PyObject, features) = new(pywrap(model), features)
    KerasNetwork(model::Module, features) = new(model, features)
end

"Load a Keras based policy from `folder` with `name`"
function KerasNetwork(folder::AbstractString, name::AbstractString)
    println(STDERR, "Loading Keras model from $(joinpath(folder, name))")
    open(joinpath(folder, string(name, ".yml"))) do file
        yaml = readall(file)
        model = pywrap(models.model_from_yaml(yaml))
        model.load_weights(joinpath(folder, string(name, ".hf5")))

        # Need to compile in order to predict anything even if we aren't training
        model.compile(loss="categorical_crossentropy", optimizer="adadelta", metrics=["accuracy"])

        # Read in extra data on the feature maps this model was trained on
        # Do this with pyyaml because the yaml has !!python directives
        features = pyyaml.load(yaml)["features"]
        features = map(x -> eval(symbol(x)), features)

        println(STDERR, "Loaded model")
        KerasNetwork(model, features)
    end
end

# A policy takes a board and outputs a probability distribution over moves

function choose_move(board::Board, policy::KerasNetwork)
    if board.last_move == PASS_MOVE
        return PASS_MOVE
    end
    moves = move_probs(board, policy)
    color = current_player(board)
    for move in moves
        point = pointindex(move)
        if is_legal(board, point, color)
            return point
        end
    end
    return PASS_MOVE
end

function move_probs(board::Board, policy::KerasNetwork)
    X = get_features(board, features=policy.features)
    X = to_python_array(X, size(X)..., 1) # Pad it out to create a batch
    probs = policy.model.predict(X)[:]
    sortperm(probs, rev=true)
end


text:
using BlockDiagonalMatrices
using Test

@testset "BlockDiagonalMatrices" begin
    include("blockdiagonal.jl")
    include("base_maths.jl")
    include("linalg.jl")
end  # tests


text:
# # Heat equation
#
# In this example, we numerically solve the 1D [heat
# equation](https://en.wikipedia.org/wiki/Heat_equation)
#
# ```math
# \frac{∂θ}{∂t} = ν \frac{∂^2 θ}{∂x^2},
# ```
#
# in a bounded domain ``x ∈ [-1, 1]`` with homogeneous Neumann boundary
# conditions, ``∂_x θ(±1, t) = 0``.

# ## Defining a B-spline basis

# The general idea is to approximate the unknown solution by a spline of order
# ``k``.
# For this, we first define a B-spline basis ``\{ b_i(x), \, i = 1, …, N \}``, such
# that the solution at a given time $t$ is approximated by
#
# ```math
# θ(x, t) = ∑_{i = 1}^N v_i(t) b_i(x),
# ```
#
# where the $v_i$ are the B-spline coefficients.
#
# A B-spline basis is uniquely defined by its order ``k`` and by a choice of
# *knot* locations within the spatial domain, which form the spatial grid.

# For this example, we take a uniform repartition of knots in ``[-1, 1]``.
knots_in = range(-1, 1; length = 11)

# We then create a B-spline basis of order ``k = 4`` using these knots.
using BSplineKit
B = BSplineBasis(BSplineOrder(4), knots_in)

# Note that the generated basis includes an *augmented* set of knots, in which
# each boundary is repeated ``k`` times:

knots(B)

# In other words, the boundary knots have multiplicity ``k``, while interior
# knots have multiplicity 1.
# This is common practice in bounded domains, and translates the fact that the
# solution does not need to be continuous at the boundaries.
# This provides additional degrees of freedom notably for the boundary
# conditions.
# This behaviour can be disabled via the `augment` argument of
# [`BSplineBasis`](@ref).

# We can now plot the knot locations (crosses) and the generated B-spline basis:

using CairoMakie
using LaTeXStrings
CairoMakie.activate!(type = "svg")

function plot_knots!(ax, ts; knot_offset = 0.03, kws...)
    ys = zero(ts)
    ## Add offset to distinguish knots with multiplicity > 1
    for i in eachindex(ts)[(begin + 1):end]
        if ts[i] == ts[i - 1]
            ys[i] = ys[i - 1] + knot_offset
        end
    end
    scatter!(ax, ts, ys; marker = :x, markersize = 20, color = :gray, kws...)
    ax
end

function plot_basis!(ax, B; eval_args = (), kws...)
    cmap = cgrad(:tab20)
    N = length(B)
    ts = knots(B)
    hlines!(ax, 0; color = :gray)
    for (n, bi) in enumerate(B)
        color = cmap[(n - 1) / (N - 1)]
        i, j = extrema(support(bi))
        lines!(ax, ts[i]..ts[j], x -> bi(x, eval_args...); color, linewidth = 2.5)
    end
    plot_knots!(ax, ts; kws...)
    ax
end

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = L"b_i(x)")
plot_basis!(ax, B)
fig

# ## Imposing boundary conditions

# In BSplineKit, the recommended approach for solving boundary value problems is
# to use the **basis recombination** method.
# That is, to expand the solution onto a new basis consisting on linear
# combinations of B-splines ``b_i(x)``, such that each recombined basis function
# ``ϕ_j(x)`` individually satisfies the required homogeneous boundary conditions
# (BCs).
# Thanks to the local support of B-splines, basis recombination only involves a
# small number of B-splines near the boundaries.
#
# Using the [`RecombinedBSplineBasis`](@ref) type, we can easily define such
# recombined bases for many different BCs.
# In this example we generate a basis satisfying homogeneous Neumann BCs:

R = RecombinedBSplineBasis(Derivative(1), B)

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = L"ϕ_i(x)")
plot_basis!(ax, R)
fig

# We notice that, on each of the two boundaries, the two initial (or final)
# B-splines of the original basis have been combined to produce a single basis
# function that has zero derivative at each respective boundary.
# To verify this, we can plot the basis function derivatives:

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = L"ϕ_i^′(x)")
plot_basis!(ax, R; eval_args = (Derivative(1), ), knot_offset = 0.4)
fig

# Note that the new basis has two less functions than the original one,
# reflecting a loss of two degrees of freedom corresponding to the new
# constraints on each boundary:

length(B), length(R)

# ### Recombination matrix

# As stated above, the basis recombination approach consists in performing
# linear combinations of B-splines ``b_i`` to obtain a derived basis of
# functions ``ϕ_j`` satisfying certain boundary conditions.
# This can be conveniently expressed using a transformation matrix
# ``\mathbf{T}`` relating the two bases:
#
# ```math
# ϕ_j(x) = ∑_{i = 1}^N T_{ij} b_i(x)
# \quad \text{for } j = 1, 2, …, M,
# ```
#
# where ``N`` is the number of B-splines ``b_i``, and ``M < N`` is the number of
# ``ϕ_j`` functions (in this example, ``M = N - 2``).
#
# The recombination matrix associated to the generated basis can be obtained
# using [`recombination_matrix`](@ref):

T = recombination_matrix(R)

# Note that the matrix is almost an identity matrix, since most B-splines are
# kept intact in the new basis.
# This simple structure allows for very efficient computations using this
# matrix.
# The first and last columns indicate that Neumann BCs are imposed by adding the
# two first (and two last) B-splines, i.e.
#
# ```math
# ϕ_1(x) = b_1(x) + b_2(x),
# \qquad
# ϕ_M(x) = b_{N - 1}(x) + b_N(x).
# ```

# ### Representation of the solution
#
# Note that the solution ``θ(x, t)`` can be represented in the original and in
# the recombined B-spline bases as
#
# ```math
# θ(x, t) = ∑_{i = 1}^N v_i(t) b_i(x) = ∑_{j = 1}^M u_j(t) ϕ_j(x),
# ```
#
# where the ``u_j`` are the coefficients in the recombined basis.
#
# The recombination matrix introduced above can be used to transform between the
# coefficients ``u_j`` and ``v_i`` in both bases, via the linear relation
# ``\bm{v} = \mathbf{T} \bm{u}``.

# ## Initial condition

# We come back now to our problem.
# We want to impose the following initial condition:
#
# ```math
# θ(x, 0) = θ_0(x) = 1 + \cos(π x).
# ```

# First, we approximate this initial condition in the recombined
# B-spline basis that we have just constructed.
# This may be easily done using [`approximate`](@ref):

θ₀(x) = 1 + cos(π * x)
θ₀_spline = approximate(θ₀, R, MinimiseL2Error())

# To see that everything went well, we can plot the exact initial condition and
# its spline approximation, which show no important differences.

fig = Figure(resolution = (800, 400))
let ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = L"\theta")
    lines!(ax, -1..1, θ₀; label = L"θ_0(x)", color = :blue)
    lines!(ax, -1..1, x -> θ₀_spline(x); label = "Approximation", color = :orange, linestyle = :dash)
    axislegend(ax; position = :cb)
end
let ax = Axis(fig[1, 2]; xlabel = L"x", ylabel = "Difference")
    lines!(ax, -1..1, x -> θ₀(x) - θ₀_spline(x))
    plot_knots!(ax, knots(R); knot_offset = 0)
end
fig

# Note that we have access to the recombined B-spline coefficients ``u_j``
# associated to the initial condition, which we will use further below:

u_init = coefficients(θ₀_spline)

# ## Solving the heat equation
#
# To solve the governing equation, the strategy is to project the unknown
# solution onto the chosen recombined basis.
# That is, we approximate the solution as
#
# ```math
# θ(x, t) = \sum_{j = 1}^M u_j(t) \, ϕ_j(x).
# ```
#
# Plugging this representation into the heat equation, we find
#
# ```math
# \newcommand{\dd}{\mathrm{d}}
# ∑_j \frac{\dd u_j}{\dd t} \, ϕ_j(x) = ν ∑_j u_j \, ϕ_j''(x),
# ```
#
# where primes denote spatial derivatives.
#
# We can now use the [method of mean weighted
# residuals](https://en.wikipedia.org/wiki/Method_of_mean_weighted_residuals) to
# find the coefficients ``u_j``, by projecting the above equation onto a chosen
# set of *test* functions ``φ_i``:
#
# ```math
# ∑_j \frac{\mathrm{d} u_j}{\mathrm{d} t} \, ⟨φ_i, ϕ_j⟩ = ν ∑_j u_j \, ⟨φ_i, ϕ_j''⟩,
# ```
#
# where ``⟨ f, g ⟩ = ∫_{-1}^1 f(x) \, g(x) \, \mathrm{d} x`` is the inner product
# between functions.
#
# By choosing ``M`` different test functions ``φ_i``, the above problem can be
# written as the linear system
#
# ```math
# \mathbf{A} \frac{\mathrm{d} \bm{u}(t)}{\mathrm{d} t} =
# ν \mathbf{L} \bm{u}(t),
# ```
#
# where the matrices are defined by ``A_{ij} = ⟨ φ_i, ϕ_j ⟩`` and
# ``L_{ij} = ⟨ φ_i, ϕ_j'' ⟩``.
#
# Two of the most common choices of test functions ``φ_i`` are:
#
# - ``φ_i(x) = δ(x - x_i)``, where ``δ`` is Dirac's delta, and ``x_i`` are a set of
#   *collocation* points where the equation will be satisfied.
#   This approach is known as the **collocation** method.
#
# - ``φ_i(x) = ϕ_i(x)``, in which case this is the **Galerkin** method.
#
# We describe the solution using both methods in the following.

# ## Collocation method

# For the collocation method, we need to choose a set of ``M`` grid points
# ``x_j``.
# Since the basis functions implicitly satisfy the boundary conditions, these
# points must be chosen *inside* of the domain.
#
# The collocation points may be automatically generated by calling
# [`collocation_points`](@ref).
# Note that, since we pass the recombined basis `R`, the boundaries are not
# included in the chosen points:

xcol = collocation_points(R)

# We can now construct the matrices ``\mathbf{A}`` and ``\mathbf{L}`` associated
# to the collocation method.
# By definition, these matrices simply contain the evaluations of all basis
# functions ``ϕ_j`` and their derivatives at the collocation points:
# ``A_{ij} = ϕ_j(x_i)`` and ``L_{ij} = ϕ_j''(x_i)``.
# Both these matrices can be constructed in BSplineKit using
# [`collocation_matrix`](@ref).
# Note that both matrices are of type [`CollocationMatrix`](@ref), which wrap
# matrices defined in
# [BandedMatrices.jl](https://github.com/JuliaMatrices/BandedMatrices.jl).

Acol = collocation_matrix(R, xcol)
Lcol = collocation_matrix(R, xcol, Derivative(2))

# For convenience and performance, we can incorporate the heat diffusivity ``ν``
# in the ``\mathbf{L}`` matrix:

ν = 0.01
Lcol *= ν

# Finally, for the time integration, we use OrdinaryDiffEq.jl from the
# [DifferentialEquations.jl suite](https://diffeq.sciml.ai/stable/).

using LinearAlgebra
using OrdinaryDiffEq

function heat_rhs!(du, u, params, t)
    mul!(du, params.L, u)    # du = ν * L * u
    ldiv!(du, params.A, du)  # du = A \ (ν * L * u)
    du
end

## Solver parameters
params_col = (
    A = lu(Acol),  # we pass the factorised matrix A for performance
    L = Lcol,
)

tspan = (0.0, 10.0)
prob = ODEProblem(heat_rhs!, u_init, tspan, params_col)
sol_collocation = solve(prob, Tsit5(); saveat = 0.5)

function plot_heat_solution(sol, R)
    fig = Figure()
    ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = L"θ(x, t)")
    colormap = cgrad(:viridis)
    tspan = sol.prob.tspan
    Δt = tspan[2] - tspan[1]
    for (u, t) in tuples(sol)
        S = Spline(R, u)
        color = colormap[(t - tspan[1]) / Δt]
        lines!(ax, -1..1, x -> S(x); label = string(t), color, linewidth = 2)
    end
    Colorbar(fig[1, 2]; colormap, limits = tspan, label = "Time")
    fig
end

plot_heat_solution(sol_collocation, R)

# ## Galerkin method

# We start by constructing the Galerkin matrices ``\mathbf{A}`` and
# ``\mathbf{L}``.
# The first of these matrices, ``A_{ij} = ⟨ ϕ_i, ϕ_j ⟩``, is usually
# known as the *mass matrix* of the system.
# It is a positive definite symmetric matrix, which enables the use of Cholesky
# factorisation to solve the resulting linear system.
# Moreover, here it is banded thanks to the local support of the B-splines.
# The mass matrix can be constructed by calling [`galerkin_matrix`](@ref):

Agal = galerkin_matrix(R)

# Note that, unlike the collocation method, in the Galerkin method we don't need
# to specify a set of grid points, as functions are not evaluated at collocation
# points (they are instead integrated over the whole domain).
# The integration is performed using Gauss--Legendre quadrature, which can be
# made exact up to numerical precision, taking advantage of the fact that the
# product of two B-splines is a piecewise polynomial.

# As for the matrix ``\mathbf{L}`` representing the second derivative operator,
# we can write it using integration by parts as
#
# ```math
# L_{ij} = ⟨ ϕ_i, ϕ_j'' ⟩
# = -⟨ ϕ_i', ϕ_j' ⟩ + \left[ ϕ_i ϕ_j' \right]_{-1}^1  = -R_{ij},
# ```
#
# where ``R_{ij} = ⟨ ϕ_i', ϕ_j' ⟩`` is a positive definite symmetric matrix.
# Note that the boundary terms all vanish since all basis functions individually
# satisfy homogeneous Neumann boundary conditions, ``ϕ_i'(±1) = 0``.
# (The same result would be obtained with homogeneous Dirichlet boundary
# conditions.)
#
# As can be seen above, one well-known advantage of the Galerkin method is that
# the basis functions can satisfy weaker continuity conditions than in the
# collocation method, as high-order derivatives can be reduced using integration
# by parts.
#
# The matrix ``\mathbf{R}`` can be constructed using [`galerkin_matrix`](@ref):

Rgal = galerkin_matrix(R, (Derivative(1), Derivative(1)))

# Note that, instead, we could have constructed the original matrix
# ``\mathbf{L}``, which, as expected, is equal to ``\mathbf{R}`` up to a sign:

Lgal = galerkin_matrix(R, (Derivative(0), Derivative(2)))

# As in the collocation example, we include the heat diffusivity ``ν`` in the
# ``\mathbf{R}`` matrix:

parent(Rgal) .*= -ν  # we can't directly multiply Rgal, as it's a Hermitian wrapper

# We finally solve using DifferentialEquations.jl.
# Note that not much is changed compared to the collocation example.
# The only difference is that we use a Cholesky factorisation for the mass
# matrix ``\mathbf{A}``.

params_gal = (
    A = cholesky(Agal),
    L = Rgal,
)

prob = ODEProblem(heat_rhs!, u_init, tspan, params_gal)
sol_galerkin = solve(prob, Tsit5(); saveat = 0.5)

plot_heat_solution(sol_galerkin, R)

# ## Result comparison

# The solution of the Galerkin method looks very similar to the one obtained
# with the collocation method.
# However, as seen below, there are non-negligible differences between the two.

fig = Figure(resolution = (800, 400))
let ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = latexstring("θ(x, t = $(tspan[end]))"))
    for pair in (
            "Collocation" => sol_collocation,
            "Galerkin" => sol_galerkin,
        )
        label, sol = pair
        u = last(sol.u)
        S = Spline(R, u)
        lines!(ax, -1..1, x -> S(x); label, linewidth = 2)
    end
    axislegend(ax; position = :cb)
end
let ax = Axis(fig[1, 2]; xlabel = L"x", ylabel = "Difference")
    Sc = Spline(R, last(sol_collocation.u))
    Sg = Spline(R, last(sol_galerkin.u))
    lines!(ax, -1..1, x -> Sc(x) - Sg(x); linewidth = 2)
end
fig

# Compared to the Galerkin method, there seems to be some additional dissipation
# in the domain interior when using the collocation method.
# This hints at the presence of numerical dissipation introduced by this method.

# To finish, we compare the two solutions to a solution at a higher resolution, using
# a higher number of B-spline knots and a higher B-spline order.
# This last solution is obtained using the collocation method to allow for
# better comparisons between both methods.

hi_res = let
    knots_in = range(-1, 1; length = 101)
    B = BSplineBasis(BSplineOrder(6), knots_in)
    R = RecombinedBSplineBasis(Derivative(1), B)
    θ₀_spline = approximate(θ₀, R)
    u_init = coefficients(θ₀_spline)
    xcol = collocation_points(R)
    Acol = collocation_matrix(R, xcol)
    Lcol = ν .* collocation_matrix(R, xcol, Derivative(2))
    params_col = (A = lu(Acol), L = Lcol)
    prob = ODEProblem(heat_rhs!, u_init, tspan, params_col)
    sol = solve(prob, Tsit5(); saveat = 0.5)
    (; R, sol)
end

fig = Figure(resolution = (800, 400))
let ax = Axis(fig[1, 1]; xlabel = L"x", ylabel = latexstring("θ(x, t = $(tspan[end]))"))
    for pair in (
            "Collocation" => sol_collocation,
            "Galerkin" => sol_galerkin,
        )
        label, sol = pair
        u = last(sol.u)
        S = Spline(R, u)
        lines!(ax, -1..1, x -> S(x); label, linewidth = 2)
    end
    let u = last(hi_res.sol.u)
        S = Spline(hi_res.R, u)
        lines!(ax, -1..1, x -> S(x); label = "Hi-res", linewidth = 2, linestyle = :dash, color = :gray)
    end
    axislegend(ax; position = :cb)
end
let ax = Axis(fig[1, 2]; xlabel = L"x", ylabel = "Difference with hi-res solution")
    Sc = Spline(R, last(sol_collocation.u))
    Sg = Spline(R, last(sol_galerkin.u))
    S_hi = Spline(hi_res.R, last(hi_res.sol.u))
    lines!(ax, -1..1, x -> Sc(x) - S_hi(x); label = "Colloc.", linewidth = 2)
    lines!(ax, -1..1, x -> Sg(x) - S_hi(x); label = "Galerkin", linewidth = 2)
    axislegend(ax; position = :rb)
end
fig

# We see that the low-resolution solution with the Galerkin method matches the
# high-resolution solution.
# This confirms that the Galerkin method provides higher accuracy than the
# collocation method when both are used at the same resolution.

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
# TODO: Nice API for all the different forms for scale and shift?
batchnorm_flags() = Lib.dnnl_use_scaleshift
batchnorm_flags(::typeof(identity)) = batchnorm_flags()
batchnorm_flags(::typeof(Flux.relu)) = batchnorm_flags() | Lib.dnnl_fuse_norm_relu

function batchnorm_inference(
    src::Memory,
    scale_shift::Memory;
    epsilon = 1f-10,
    attributes = noattributes(),
    opdesc = Ref{Lib.dnnl_batch_normalization_desc_t}(),
)
    @apicall dnnl_batch_normalization_forward_desc_init(
        opdesc, Inference(), src, epsilon, batchnorm_flags()
    )

    return temp_primitive(opdesc, attributes, global_engine(), noforward()) do p, pd
        dst = similar(src)
        execute!(p, @dnnl_args src scale_shift dst)
        return dst
    end
end

function batchnorm_training(
    src::Memory,
    scale_shift::Memory,
    activation::F = identity;
    epsilon = 1f-10,
    attributes = noattributes(),
    opdesc = Ref{Lib.dnnl_batch_normalization_desc_t}(),
) where {F<:Union{typeof(identity),typeof(Flux.relu)}}
    flags = batchnorm_flags(activation)
    @apicall dnnl_batch_normalization_forward_desc_init(
        opdesc, Training(), src, epsilon, flags
    )

    return temp_primitive(opdesc, attributes, global_engine(), noforward()) do p, pd
        dst = similar(src)
        mean = similar(src, eltype(src), (size(src, ndims(src) - 1),))
        variance = similar(mean)

        # Fuse the Relu activation function is possible.
        if activation === Base.identity
            execute!(p, @dnnl_args src scale_shift dst mean variance)
            return (; dst, mean, variance, workspace = nothing, forward = noforward())
        elseif activation === Flux.relu
            workspace = similar(src)
            execute!(p, @dnnl_args src scale_shift dst mean variance workspace)
            return (; dst, mean, variance, workspace, forward = pd)
        end
    end
end

function batchnorm_backward(
    diff_dst::Memory,
    src::Memory,
    scale_shift::Memory,
    mean::Memory,
    variance::Memory;
    forward = noforward(),
    epsilon = 1f-10,
    attributes = noattributes(),
    workspace = nothing,
    opdesc = Ref{Lib.dnnl_batch_normalization_desc_t}(),
)
    flags = (workspace === nothing) ? batchnorm_flags(identity) : batchnorm_flags(Flux.relu)
    src = maybe_reorder(memorydesc(diff_dst), src)

    @apicall dnnl_batch_normalization_backward_desc_init(
        opdesc, Lib.dnnl_backward, diff_dst, src, epsilon, flags
    )

    return temp_primitive(opdesc, attributes, global_engine(), forward) do p, pd
        diff_src = similar(diff_dst)
        diff_scale_shift = similar(scale_shift)
        args = @dnnl_args diff_dst src mean variance scale_shift diff_src diff_scale_shift
        # Maybe append on the workspace as well.
        if workspace !== nothing
            args = append(args, @dnnl_args workspace)
        end
        execute!(p, args)
        return (; diff_src, diff_scale_shift)
    end
end

#####
##### Higher level API
#####

struct BatchNorm{T<:OneDNN.Memory,F}
    scale_shift::T
    activation::F
    epsilon::Float32
    attributes::Attributes
    opdesc::Base.RefValue{Lib.dnnl_batch_normalization_desc_t}
end

function BatchNorm(scale_shift::AbstractMatrix, activation = identity; epsilon = 1f-5)
    attributes = Attributes()
    postops = PostOps()
    eltwise!(postops, activation)
    append!(attributes, postops)
    opdesc = Ref{Lib.dnnl_batch_normalization_desc_t}()
    return BatchNorm(OneDNN.Memory(scale_shift), activation, epsilon, attributes, opdesc)
end

function (bn::BatchNorm)(x)
    return batchnorm_inference(
        Memory(x), bn.scale_shift; bn.epsilon, bn.attributes, bn.opdesc
    )
end

function set_scale_shift!(x::AbstractArray{T,2}) where {T}
    for i in eachindex(x)
        x[i] = randn(T)
    end
    return nothing
end

_batchnorm_fusable(::Any) = false
_batchnorm_fusable(::typeof(identity)) = true
_batchnorm_fusable(::typeof(Flux.relu)) = true

function ChainRulesCore.rrule(bn::BatchNorm, src::Memory)
    return if _batchnorm_fusable(bn.activation)
        batchnorm_rrule_fused(bn, src)
    else
        batchnorm_rrule_unfused(bn, src)
    end
end

function batchnorm_rrule_fused(bn::BatchNorm, src::Memory)
    nt = batchnorm_training(src, bn.scale_shift, bn.activation; bn.epsilon, bn.opdesc)

    function batchnorm_fused_pullback(diff_data::Memory)
        Δ = batchnorm_backward(
            diff_data,
            src,
            bn.scale_shift,
            nt.mean,
            nt.variance;
            nt.forward,
            bn.epsilon,
            nt.workspace,
            bn.opdesc,
        )
        return ((; scale_shift = Δ.diff_scale_shift), Δ.diff_src)
    end

    return (nt.dst, batchnorm_fused_pullback)
end

function batchnorm_rrule_unfused(bn::BatchNorm, src::Memory)
    nt = batchnorm_training(src, bn.scale_shift, identity; bn.epsilon)
    eltwise_output, eltwise_pullback = ChainRulesCore.rrule(eltwise, bn.activation, nt.dst)

    function batchnorm_unfused_pullback(_diff_data::Memory)
        # Reverse the eltwise function.
        diff_data = eltwise_pullback(_diff_data)[3]
        Δ = batchnorm_backward(
            diff_data, src, bn.scale_shift, nt.mean, nt.variance; bn.epsilon
        )
        return ((; scale_shift = Δ.diff_scale_shift), Δ.diff_src)
    end

    return (eltwise_output, batchnorm_unfused_pullback)
end